GARCH(1,1)#

Introduction#

The GARCH(1,1) model is a commonly used model for capturing the time-varying volatility in financial time series data. The model can be defined as follows:

Tip

Check Examples section for code guide and comparison to Stata.

Return equation#

There are many ways to specify return dynamics. Here a constant mean model is used.

(1)#\[r_t = \mu + \epsilon_t\]

where \(r_t\) represents the return at time \(t\), and \(\mu\) is the mean return.

Shock equation#

(2)#\[\epsilon_t = \sigma_t \cdot z_t\]

In this equation, \(\epsilon_t\) is the shock term, \(\sigma_t\) is the conditional volatility, and \(z_t\) is a white noise error term with zero mean and unit variance (\(z_t \sim N(0,1)\)).

Note

We can also assume that the noise term follows a different distribution, such as Student-t, and modify the likelihood function below accordingly.

Volatility equation#

(3)#\[\sigma_t^2 = \omega + \alpha \cdot \epsilon_{t-1}^2 + \beta \cdot \sigma_{t-1}^2\]

Here \(\sigma_t^2\) is the conditional variance at time \(t\), and \(\omega\), \(\alpha\), \(\beta\) are parameters to be estimated. This equation captures how volatility evolves over time.

The unconditional variance and persistence

The unconditional variance, often denoted as \(\text{Var}(\epsilon_t)\) or \(\sigma^2\), refers to the long-run average or steady-state variance of the return series. It is the variance one would expect the series to revert to over the long term, and it doesn’t condition on any past information.

For a GARCH(1,1) model to be stationary, the persistence, sum of \(\alpha\) and \(\beta\), must be less than 1 ( \(\alpha + \beta < 1\) ). Given this condition, the unconditional variance \(\sigma^2\) can be computed as follows:

(4)#\[\sigma^2 = \frac{\omega}{1 - \alpha - \beta}\]

In this formulation, \(\omega\) is the constant or “base” level of volatility, while \(\alpha\) and \(\beta\) determine how shocks to returns and past volatility influence future volatility. The unconditional variance provides a long-run average level around which the conditional variance oscillates.

Log-likelihood function#

The likelihood function for a GARCH(1,1) model is used for the estimation of parameters \(\mu\), \(\omega\), \(\alpha\), and \(\beta\). Given a time series of returns \(\{ r_1, r_2, \ldots, r_T \}\), the likelihood function \(L(\mu, \omega, \alpha, \beta)\) can be written as:

(5)#\[L(\mu, \omega, \alpha, \beta) = \prod_{t=1}^{T} \frac{1}{\sqrt{2\pi \sigma_t^2}} \exp\left(-\frac{(r_t-\mu)^2}{2\sigma_t^2}\right)\]

Taking the natural logarithm of \(L\), we obtain the log-likelihood function \(\ell(\mu, \omega, \alpha, \beta)\):

(6)#\[\ell(\mu, \omega, \alpha, \beta) = -\frac{1}{2} \sum_{t=1}^{T} \left[\ln(2\pi) + \ln(\sigma_t^2) + \frac{(r_t-\mu)^2}{\sigma_t^2} \right]\]

The parameters \(\mu, \omega, \alpha, \beta\) can then be estimated by maximizing this log-likelihood function.

Estimation techniques#

A few tips to improve the estimation and enhance its numerical stability.

Initial value of conditional variance#

Note that the conditional variance in a GARCH(1,1) model is (3):

\[\sigma_t^2 = \omega + \alpha \cdot \epsilon_{t-1}^2 + \beta \cdot \sigma_{t-1}^2\]

We need a good starting value \(\sigma_0^2\) to begin with, which can be estimated via the backcasting technique. Once we have that \(\sigma^2_0\) through backcasting, we can proceed to calculate the entire series of conditional variances using the standard GARCH recursion formula.

To backcast the initial variance, we can use the Exponential Weighted Moving Average (EWMA) method, setting \(\sigma^2_0\) to the EWMA of the sample variance of the first \(n \leq T\) returns:

\[\sigma^2_0 = \sum_{t=1}^{n} w_t \cdot r_t^2\]

where \(w_t\) are the exponentially decaying weights and \(r_t\) are residuals of returns, i.e., returns de-meaned by sample average. This \(\sigma^2_0\) is then used to derive \(\sigma^2_1\) the starting value for the conditional variance series.

Initial value of \(\omega\)#

The starting value of \(\omega\) is relatively straightforward. Notice that earlier we have jotted down the unconditional variance \(\sigma^2 = \frac{\omega}{1-\alpha-\beta}\). Therefore, given a level of persistence (\(\alpha+\beta\)), we can set the initial guess of \(\omega\) to be the sample variance times one minus persistence:

\[\omega = \hat{\sigma}^2 \cdot (1-\alpha-\beta)\]

where we use the known sample variance of residuals \(\hat{\sigma}^2\) as a guess for the unconditional variance \(\sigma^2\). However, we still need to find good starting values for \(\alpha\) and \(\beta\).

Initial value of \(\alpha\) and \(\beta\)#

Unfortunately, there is no better way to find good starting values for \(\alpha\) and \(\beta\) than a grid search. Luckily, this grid search can be relatively small.

First, we don’t know ex ante the persistence level, so we need to vary the persistence level from some low values to some high values, e.g., from 0.1 to 0.98.
Second, generally the \(\alpha\) parameter is not too big, for example, ranging from 0.01 to 0.2.

We can permute combinations of the persistence level and \(\alpha\), which naturally gives the corresponding \(\beta\) and hence \(\omega\). The “optimal” set of initial values of \(\omega, \alpha, \beta\) is the one that gives the highest log-likelihood.

Note

The initial value of \(\mu\) is reasonably set to the sample mean return.

Variance bounds#

Another issue is that we want to ensure that in the estimation, condition variance does not blow up to infinity or becomes zero. Hence, we need to construct bounds for conditional variances during the GARCH(1,1) parameter estimation process.

To do this, we can calculate loose lower and upper bounds for each observation. Specifically, we can use sample variance of the residuals to compute global lower and upper bounds. We then use EWMA to compute the conditional variance for each time point. The EWMA variances are then adjusted to ensure they are within global bounds. Lastly, we scale the adjusted EWMA variances to form the variance bounds at each time.

During the estimation process, whenever we compute the conditional variances based on the prevailing model parameters, we ensure that they are adjusted to be reasonably within the bounds at each time.

References#

Engle, R. F. (1982), “Autoregressive Conditional Heteroskedasticity with Estimates of the Variance of United Kingdom Inflation.” Econometrica, 50(4), 987-1007.
Bollerslev, T. (1986), “Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics, 31(3), 307-327.
arch by Kevin Sheppard, et al.

API#

class frds.algorithms.GARCHModel(returns: ndarray, zero_mean=False)[source]#

GARCH(1,1) model with constant mean and Normal noise

It estimates the model parameters only. No standard errors calculated.

This code is heavily based on the arch. Modifications are made for easier understaing of the code flow.

__init__(returns: ndarray, zero_mean=False) → None[source]#

Parameters:

returns (np.ndarray) – (T,) array of T returns
zero_mean (bool) – whether to use a zero mean returns model. Default to False.

Note

returns is best to be percentage returns for optimization

fit() → Parameters[source]#

Estimates the GARCH(1,1) parameters via MLE

Returns:: frds.algorithms.GARCHModel.Parameters
Return type:: params

preparation() → List[float][source]#

Prepare starting values.

Returns:: list of starting values
Return type:: List[float]

loglikelihood_model(params: ndarray, backcast: float, var_bounds: ndarray, zero_mean=False) → float[source]#

Calculates the negative log-likelihood based on the current params. This function is used in optimization.

Parameters:

params (np.ndarray) – [mu, omega, alpha, (gamma), beta]
backcast (float) – backcast value
var_bounds (np.ndarray) – variance bounds

Returns:

negative log-likelihood

Return type:

float

static loglikelihood(resids: ndarray, sigma2: ndarray) → float[source]#

Computes the log-likelihood assuming residuals are normally distributed conditional on the variance.

Parameters:

resids (np.ndarray) – residuals to use in computing log-likelihood.
sigma2 (np.ndarray) – conditional variance of residuals.

Returns:

log-likelihood

Return type:

float

static compute_variance(params: List[float], resids: ndarray, backcast: float, var_bounds: ndarray) → ndarray[source]#

Computes the variances conditional on given parameters

Parameters:

params (List[float]) – [omega, alpha, beta]
resids (np.ndarray) – residuals from mean equation
backcast (float) – backcast value
var_bounds (np.ndarray) – variance bounds

Returns:

conditional variance

Return type:

np.ndarray

starting_values(resids: ndarray) → List[float][source]#

Finds the optimal initial values for the volatility model via a grid search. For varying target persistence and alpha values, return the combination of alpha and beta that gives the highest loglikelihood.

Parameters:: resids (np.ndarray) – residuals from the mean model
Returns:: [omega, alpha, beta]
Return type:: List[float]

static variance_bounds(resids: ndarray) → ndarray[source]#

Compute bounds for conditional variances using EWMA.

This function calculates the lower and upper bounds for conditional variances based on the residuals provided. The bounds are computed to ensure numerical stability during the parameter estimation process of GARCH models. The function uses Exponentially Weighted Moving Average (EWMA) to estimate the initial variance and then adjusts these estimates to lie within global bounds.

Parameters:: resids (np.ndarray) – residuals from the mean model.
Returns:: an array where each row contains the lower and upper bounds for the conditional variance at each time point.
Return type:: np.ndarray

static bounds_check(sigma2: float, var_bounds: ndarray) → float[source]#

Adjust the conditional variance at time t based on its bounds

Parameters:

sigma2 (float) – conditional variance
var_bounds (np.ndarray) – lower and upper bounds

Returns:

adjusted conditional variance

Return type:

float

static backcast(resids: ndarray) → float[source]#

Computes the starting value for estimating conditional variance.

Parameters:: resids (np.ndarray) – residuals
Returns:: initial value from backcasting
Return type:: float

static ewma(resids: ndarray, initial_value: float, lam=0.94) → ndarray[source]#

Compute the conditional variance estimates using Exponentially Weighted Moving Average (EWMA).

Parameters:

resids (np.ndarray) – Residuals from the model.
initial_value (float) – Initial value for the conditional variance.
lam (float) – Decay factor for the EWMA.

Returns:

Array containing the conditional variance estimates.

Return type:

np.ndarray

class frds.algorithms.GARCHModel.Parameters(mu: float = nan, omega: float = nan, alpha: float = nan, beta: float = nan, loglikelihood: float = nan)#

Examples#

Let’s import the dataset.

>>> import pandas as pd
>>> data_url = "https://www.stata-press.com/data/r18/stocks.dta"
>>> df = pd.read_stata(data_url, convert_dates=["date"])

Scale returns to percentage returns for better optimization results

>>> returns = df["nissan"].to_numpy() * 100

Use frds.algorithms.GARCHModel to estimate a GARCH(1,1).

>>> from frds.algorithms import GARCHModel
>>> from pprint import pprint
>>> model = GARCHModel(returns)
>>> res = model.fit()
>>> pprint(res)
Parameters(mu=0.019315543596552513,
           omega=0.05701047522984261,
           alpha=0.0904653253307871,
           beta=0.8983752570013462,
           loglikelihood=-4086.487358003049)

These estimates are identical to the ones produced by arch.

>>> from arch import arch_model
>>> model = arch_model(returns, mean='Constant', vol='GARCH', p=1, q=1)
>>> model.fit(disp=False)
                     Constant Mean - GARCH Model Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.000
Mean Model:             Constant Mean   Adj. R-squared:                  0.000
Vol Model:                      GARCH   Log-Likelihood:               -4086.49
Distribution:                  Normal   AIC:                           8180.97
Method:            Maximum Likelihood   BIC:                           8203.41
                                        No. Observations:                 2015
Date:                Thu, Sep 28 2023   Df Residuals:                     2014
Time:                        13:04:30   Df Model:                            1
                                  Mean Model
=============================================================================
                 coef    std err          t      P>|t|       95.0% Conf. Int.
-----------------------------------------------------------------------------
mu             0.0193  3.599e-02      0.536      0.592 [-5.124e-02,8.985e-02]
                             Volatility Model
==========================================================================
                 coef    std err          t      P>|t|    95.0% Conf. Int.
--------------------------------------------------------------------------
omega          0.0570  2.810e-02      2.029  4.245e-02 [1.943e-03,  0.112]
alpha[1]       0.0905  2.718e-02      3.328  8.744e-04 [3.719e-02,  0.144]
beta[1]        0.8984  2.929e-02     30.670 1.426e-206   [  0.841,  0.956]
==========================================================================

Additionaly, in Stata, we can estimate the same model as below:

webuse stocks, clear
replace nissan = nissan * 100
arch nissan, arch(1) garch(1) vce(robust)

It would produce very similar estimates. The discrepencies are likly due to the different optimization algorithms used. Based on loglikelihood, the estimates from arch and frds.algorithms.GARCHModel are marginally better.

Notes#

I did many tests, and in 99% cases frds.algorithms.GARCHModel performs equally well with arch, simply because it’s adapted from arch. In some rare cases when the return series does not behave well, though, the two would produce very different estimates despite having almost identical log-likelihood.