GJR-GARCH(1,1)#

Introduction#

The GJR-GARCH model extends the basic GARCH(1,1) by accounting for leverage effects, where bad news (negative returns) has a greater impact on volatility than good news.

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\]

Here, \(r_t\) is the asset return at time \(t\), \(\mu\) is the mean of the asset return, and \(\epsilon_t\) is the shock term.

Shock equation#

The shock term \(\epsilon_t\) can be decomposed as:

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

where \(\sigma_t\) is the conditional volatility and \(z_t \sim N(0, 1)\) a standard normal noise term.

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#

The conditional volatility in GJR-GARCH(1,1) is given by:

(3)#\[\sigma^2_t = \omega + \alpha \epsilon_{t-1}^2 + \gamma \epsilon_{t-1}^2 I_{t-1} + \beta \sigma_{t-1}^2\]

Here, \(I_{t-1} = 1\) if \(\epsilon_{t-1} < 0\) and 0 otherwise. \(\omega, \alpha, \gamma, \beta\) are parameters to be estimated. The term \(\gamma \epsilon_{t-1}^2 I_{t-1}\) captures the leverage effect.

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 GJR-GARCH(1,1) model to be stationary, the persistence must be less than 1 ( \(\alpha + \beta + \gamma/2 < 1\) ). Given this condition, the unconditional variance \(\sigma^2\) can be computed as follows:

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

Note that the division by 2 for \(\gamma\) assumes that the leverage effect occurs about half the time, given that \(I_{t-1}\) takes the value 1 if the return is negative and 0 otherwise.

Log-likelihood function#

The log-likelihood function for the GJR-GARCH(1,1) \(\ell(\mu, \omega, \alpha, \beta)\) is:

(5)#\[\ell(\mu, \omega, \alpha, \gamma, \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, \gamma, \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. These are basically the same as in estimating GARCH(1,1).

Initial value of conditional variance#

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

\[\sigma^2_t = \omega + \alpha \epsilon_{t-1}^2 + \gamma \epsilon_{t-1}^2 I_{t-1} + \beta \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. Note that the unconditional variance is given by (4): \(\sigma^2 = \frac{\omega}{1-\alpha-\gamma/2-\beta}\). Therefore, given a level of persistence (\(\alpha+\gamma/2+\beta\)), we can set the initial guess of \(\omega\) to be the sample variance times one minus persistence:

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

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\), \(\gamma\) and \(\beta\).

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

Finding good starting values for \(\alpha\) \(\gamma\) and \(\beta\) can be done by a small grid search.

  • 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.

  • Third, the leverage effect \(\gamma\) is generally not big too. We can set it to be in the same range as \(\alpha\).

We can permute combinations of the persistence level, \(\alpha\) and \(\gamma\), which naturally gives the corresponding \(\beta\) and hence \(\omega\). The “optimal” set of initial values of \(\omega, \alpha, \gamma, \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 GJR-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#

API#

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

GJR-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.

fit() → Parameters[source]#

Estimates the GJR-GARCH(1,1) parameters via MLE

Returns:

[mu, omega, alpha, gamma, beta, loglikelihood]

Return type:

List[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, gamma, 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

static forecast_variance(params: Parameters, resids: ndarray, initial_variance: float) → ndarray[source]#

Forecast the variances conditional on given parameters and residuals.

Parameters:
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, gamma, beta]

Return type:

List[float]

__init__(returns: ndarray, zero_mean=False) → None#
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

static backcast(resids: ndarray) → float#

Computes the starting value for estimating conditional variance.

Parameters:

resids (np.ndarray) – residuals

Returns:

initial value from backcasting

Return type:

float

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

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 ewma(resids: ndarray, initial_value: float, lam=0.94) → ndarray#

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

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

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

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

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

preparation() → List[float]#

Prepare starting values.

Returns:

list of starting values

Return type:

List[float]

static variance_bounds(resids: ndarray) → ndarray#

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

class frds.algorithms.GJRGARCHModel.Parameters(mu: float = nan, omega: float = nan, alpha: float = nan, gamma: 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.GJRGARCHModel to estimate a GJR-GARCH(1,1).

>>> from frds.algorithms import GJRGARCHModel
>>> from pprint import pprint
>>> model = GJRGARCHModel(returns)
>>> res = model.fit()
>>> pprint(res)
Parameters(mu=0.010528449295629098,
           omega=0.05512898468355955,
           alpha=0.07700974411970742,
           gamma=0.021814015760057957,
           beta=0.9013499076166999,
           loglikelihood=-4085.741514140086)

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, o=1, q=1)
>>> model.fit(disp=False)
                   Constant Mean - GJR-GARCH Model Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.000
Mean Model:             Constant Mean   Adj. R-squared:                  0.000
Vol Model:                  GJR-GARCH   Log-Likelihood:               -4085.74
Distribution:                  Normal   AIC:                           8181.48
Method:            Maximum Likelihood   BIC:                           8209.52
                                        No. Observations:                 2015
Date:                Fri, Sep 29 2023   Df Residuals:                     2014
Time:                        09:59:37   Df Model:                            1
                                  Mean Model
=============================================================================
                 coef    std err          t      P>|t|       95.0% Conf. Int.
-----------------------------------------------------------------------------
mu             0.0105  3.632e-02      0.290      0.772 [-6.066e-02,8.171e-02]
                               Volatility Model
=============================================================================
                 coef    std err          t      P>|t|       95.0% Conf. Int.
-----------------------------------------------------------------------------
omega          0.0551  2.901e-02      1.900  5.743e-02   [-1.740e-03,  0.112]
alpha[1]       0.0770  3.428e-02      2.247  2.467e-02    [9.821e-03,  0.144]
gamma[1]       0.0218  2.214e-02      0.985      0.324 [-2.158e-02,6.522e-02]
beta[1]        0.9014  3.159e-02     28.532 4.682e-179      [  0.839,  0.963]
=============================================================================

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

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

Important! \(\gamma\) in Stata

In Stata, the estimate of \(\gamma\) is the negative of the ones obtained in arch and frds.algorithms.GJRGARCHModel!

This is because in Stata, GJR-GARCH is specified via the tarch option, or threshold ARCH. It defines the indicator \(I_{t-1}\) in equation (3) the opposite way, \(I_{t-1} = 1\) if \(\epsilon_{t-1} > 0\) and 0 otherwise.

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.GJRGARCHModel are marginally better.

Notes#

I did many tests, and in 99% cases frds.algorithms.GJRGARCHModel 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.