GARCH(1,1) Model in Python

In my previous article GARCH(p,q) Model and Exit Strategy for Intraday Algorithmic Traders we described the essentials of GARCH(p,q) model and provided an exemplary implementation in Matlab. In general, we apply GARCH model in order to estimate the volatility one time-step forward, where:
$$
\sigma_t^2 = \omega + \alpha r_{t-1}^2 + \beta \sigma_{t-1}^2
$$ based on the most recent update of $r$ and $\sigma$, where $r_{t-1} = \ln({P_{t-1}}/{P_{t-2}})$ and $P$ corresponds to an asset price. For any financial time-series, $\{r_j\}$, the estimation of $(\omega,\alpha,\beta)$ parameters can be conducted utilising the maximum likelihood method. The latter is an iterative process by looking for the maximum value of the sum among all sums defined as:
$$
\sum_{i=3}^{N} \left[ -\ln(\sigma_i^2) – \frac{r_i^2}{\sigma_i^2} \right] $$ where $N$ denotes the length of the return series $\{r_j\}$ ($j=2,…,N$) under study.

GARCH(1,1) Model

The GARCH(1,1) model is widely regarded as a peculiar and distinctive model in the realm of financial time series analysis due to its ability to capture two significant characteristics commonly observed in financial markets: volatility clustering and conditional heteroscedasticity.

Volatility clustering refers to the tendency of financial data to exhibit periods of high volatility followed by periods of low volatility. This phenomenon is pervasive in various financial assets and markets. The GARCH(1,1) model’s uniqueness lies in its capacity to effectively account for volatility clustering by incorporating a lagged term of the conditional variance itself. By doing so, the model acknowledges that large shocks today increase the likelihood of subsequent large shocks in the future. This feature contributes to capturing the observed volatility clustering pattern that is frequently encountered in financial markets.

Moreover, the GARCH(1,1) model addresses another crucial characteristic called conditional heteroscedasticity. In financial markets, the variance of asset returns tends to change over time, indicating that the assumption of constant variance often does not hold. To accommodate this feature, the GARCH(1,1) model represents the conditional variance as a combination of the conditional mean (typically modeled by an autoregressive process) and a term that accounts for the conditional variance itself. This modeling approach allows the variance to vary dynamically based on past shocks, thereby capturing the time-varying nature of volatility observed in financial data.

The peculiarity of the GARCH(1,1) model lies not only in its ability to effectively capture volatility clustering and conditional heteroscedasticity but also in its simplicity. Despite its relatively straightforward formulation, the GARCH(1,1) model has proven to be remarkably useful in finance and econometrics. It has become a standard tool for modeling and forecasting volatility, finding applications in areas such as risk management, option pricing, portfolio optimization, and market microstructure analysis.

However, it is worth noting that the GARCH(1,1) model is just one variant within the broader GARCH family. Over time, researchers have developed more complex variants, such as GARCH-M, EGARCH, and GJR-GARCH, to address certain limitations or incorporate additional features observed in financial data. These advanced models introduce various refinements and extensions, offering a more nuanced representation of volatility dynamics in financial markets.

Back to action. By putting $p=1, q=1$, we are able to fit the most fundamental GARCH model:

r = np.array([0.945532630498276,
              0.614772790142383,
              0.834417758890680,
              0.862344782601800,
              0.555858715401929,
              0.641058419842652,
              0.720118656981704,
              0.643948007732270,
              0.138790608092353,
              0.279264178231250,
              0.993836948076485,
              0.531967023876420,
              0.964455754192395,
              0.873171802181126,
              0.937828816793698])

from arch import arch_model
garch11 = arch_model(r, p=1, q=1)
res = garch11.fit(update_freq=10)
print(res.summary())

Above we have used the functionality of the ARCH: a Python library containing, inter alia, coroutines for the analysis of univariate volatility models. The result of the GARCH(1,1) model to our data are summarised as follows:

Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.118198462057
            Iterations: 9
            Function evaluations: 60
            Gradient evaluations: 9
                     Constant Mean - G ARCH Model Results                      
==============================================================================
Dep. Variable:                      y   R-squared:                      -0.004
Mean Model:             Constant Mean   Adj. R-squared:                 -0.004
Vol Model:                     G ARCH   Log-Likelihood:               0.118198
Distribution:                  Normal   AIC:                           7.76360
Method:            Maximum Likelihood   BIC:                           10.5958
                                        No. Observations:                   15
Date:                Wed, Dec 14 2016   Df Residuals:                       11
Time:                        10:25:30   Df Model:                            4
                                  Mean Model                                  
==============================================================================
                 coef    std err          t      P>|t|        95.0% Conf. Int.
------------------------------------------------------------------------------
mu             0.7187      0.116      6.188  6.080e-10       [  0.491,  0.946]
                               Volatility Model                               
==============================================================================
                 coef    std err          t      P>|t|        95.0% Conf. Int.
------------------------------------------------------------------------------
omega          0.0485  6.877e-02      0.706      0.480    [-8.626e-02,  0.183]
alpha[1]       0.1733      0.688      0.252      0.801       [ -1.176,  1.522]
beta[1]    5.0529e-16      2.172  2.326e-16      1.000       [ -4.257,  4.257]
==============================================================================

Covariance estimator: robust

Python vs. Matlab Solution

Programming requires caution. It is always a good practice to test the outcome of one algorithm against alternative solutions. Let’s run the G ARCH(1,1) model estimation for the same input array and compare Python and Matlab results:

clear all; close all; clc;

r=[0.945532630498276, ...
   0.614772790142383, ...
   0.834417758890680, ...
   0.862344782601800, ...
   0.555858715401929, ...
   0.641058419842652, ...
   0.720118656981704, ...
   0.643948007732270, ...
   0.138790608092353, ...
   0.279264178231250, ...
   0.993836948076485, ...
   0.531967023876420, ...
   0.964455754192395, ...
   0.873171802181126, ...
   0.937828816793698]';

% (p,q) parameter estimation
model = garch(1,1) % define model
[fit,VarCov,LogL,Par] = estimate(model,r);
% extract model parameters
parC=Par.X(1); % omega
parG=Par.X(2); % beta (G ARCH)
parA=Par.X(3); % alpha (ARCH)
% estimate unconditional volatility
gamma=1-parA-parG;
VL=parC/gamma;
volL=sqrt(VL);
% redefine model with estimatated parameters
model=garch('Constant',parC,'GARCH',parG,'ARCH',parA)

what returns:

model = 

    GARCH(1,1) Conditional Variance Model:
    --------------------------------------  
    Distribution: Name = 'Gaussian'
               P: 1
               Q: 1
        Constant: NaN
           GARCH: {NaN} at Lags [1]
            ARCH: {NaN} at Lags [1]

____________________________________________________________
   Diagnostic Information

Number of variables: 3

Functions 
Objective:       @(X)OBJ.nLogLikeGaussian(X,V,E,Lags,1,maxPQ,T,nan,trapValue)
Gradient:        finite-differencing
Hessian:         finite-differencing (or Quasi-Newton)

Constraints
Nonlinear constraints:                do not exist
 
Number of linear inequality constraints:    1
Number of linear equality constraints:      0
Number of lower bound constraints:          3
Number of upper bound constraints:          3

Algorithm selected
   sequential quadratic programming


____________________________________________________________
   End diagnostic information
                                                          Norm of First-order
 Iter F-count            f(x) Feasibility  Steplength        step  optimality
    0       4    1.748188e+01   0.000e+00                           5.758e+01
    1      27    1.723863e+01   0.000e+00   1.140e-03   6.565e-02   1.477e+01
    2      31    1.688626e+01   0.000e+00   1.000e+00   9.996e-01   1.510e+00
    3      35    1.688234e+01   0.000e+00   1.000e+00   4.099e-02   1.402e+00
    4      39    1.686305e+01   0.000e+00   1.000e+00   1.440e-01   8.889e-01
    5      44    1.685246e+01   0.000e+00   7.000e-01   2.379e-01   5.088e-01
    6      48    1.684889e+01   0.000e+00   1.000e+00   9.620e-02   1.379e-01
    7      52    1.684835e+01   0.000e+00   1.000e+00   2.651e-02   2.257e-02
    8      56    1.684832e+01   0.000e+00   1.000e+00   8.389e-03   7.046e-02
    9      60    1.684831e+01   0.000e+00   1.000e+00   1.953e-03   7.457e-02
   10      64    1.684825e+01   0.000e+00   1.000e+00   7.888e-03   7.738e-02
   11      68    1.684794e+01   0.000e+00   1.000e+00   3.692e-02   7.324e-02
   12      72    1.684765e+01   0.000e+00   1.000e+00   1.615e-01   5.862e-02
   13      76    1.684745e+01   0.000e+00   1.000e+00   7.609e-02   8.429e-03
   14      80    1.684740e+01   0.000e+00   1.000e+00   2.368e-02   4.072e-03
   15      84    1.684739e+01   0.000e+00   1.000e+00   1.103e-02   3.142e-03
   16      88    1.684739e+01   0.000e+00   1.000e+00   1.183e-03   2.716e-04
   17      92    1.684739e+01   0.000e+00   1.000e+00   9.913e-05   1.378e-04
   18      96    1.684739e+01   0.000e+00   1.000e+00   6.270e-05   2.146e-06
   19      97    1.684739e+01   0.000e+00   7.000e-01   4.327e-07   2.146e-06

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are 
satisfied to within the selected value of the constraint tolerance.

<stopping criteria details>

 
    GARCH(1,1) Conditional Variance Model:
    ----------------------------------------
    Conditional Probability Distribution: Gaussian

                                  Standard          t     
     Parameter       Value          Error       Statistic 
    -----------   -----------   ------------   -----------
     Constant       0.278061       26.3774      0.0105417
     GARCH{1}       0.457286       49.4915      0.0092397
      ARCH{1}      0.0328433       1.65576      0.0198358

model = 

    GARCH(1,1) Conditional Variance Model:
    --------------------------------------  
    Distribution: Name = 'Gaussian'
               P: 1
               Q: 1
        Constant: 0.278061
           GARCH: {0.457286} at Lags [1]
            ARCH: {0.0328433} at Lags [1]

id est
$$
(\omega,\beta,\alpha)_{\rm Matlab} = (0.278061,0.457286,0.0328433) \ .
$$This slightly differs itself from the Python solution which was
$$
(\omega,\beta,\alpha)_{\rm Py} =(0.0485,5.0529e-16,0.1733) \ .
$$At this stage it is difficult to assess which solution is “better”. Both algorithms and applied methodologies simply may be different what is usually the case. Having that in mind, further extensive tests are required, for example, a dependance of Python solution on trial input $(\omega,\alpha,\beta)$ values.

References

John Hull, Options, Futures, and Other Derivatives, Global Edition, 11th Ed (2021; updated link)

Dive Deeper

GARCH(p,q) Model and Exit Strategy for Intraday Algorithmic Traders