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.
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:
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:
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:
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:
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):
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:
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:
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#
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.
Glosten, L. R., Jagannathan, R., & Runkle, D. E. (1993), “On the Relation Between the Expected Value and the Volatility of the Nominal Excess Return on Stocks.” The Journal of Finance, 48(5), 1779-1801.
arch
by Kevin Sheppard, et al.
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:
params (Parameters) –
frds.algorithms.GJRGARCHModel.Parameters
resids (np.ndarray) – residuals to use
initial_variance (float) – starting value of variance forecasts
- 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 ofT
returnszero_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.