GARCH Models

The Generalized Auto-Regressive Conditional Heteroscedasticity (GARCH) model of Bollerslev (1986) and the numerous extensions which have followed since, is a framework for modeling the dynamics of the conditional variance. It has proved particularly popular, particularly among financial market practitioners, but also in other areas were the arrival of unexpected new information may lead to non-instantaneous decay (persistence) and/or asymmetric reaction to good and bad news (news impact).

Many programming languages have one or more implementations of GARCH, with R having no less than 3, including the garch function from the tseries package, fGarch and rugarch. A select R package review is provided in Hill and McCullough (2019), and a cross language package review on asymmetric GARCH in Charles and Darné (2019). The suggestions and constructive feedback in these papers were taken into account when developing this new package, as were many of the user comments over the last 10+ years.

The tsgarch package is a partial re-implementation of rugarch, by the same author, with key differences summarized below:

• does not (yet) implement all GARCH models in rugarch. FIGARCH, Multiplicative Component GARCH and Realized GARCH are not currently implemented.

• does not implement joint ARFIMA-GARCH estimation. The conditional mean equation only allows for a constant. With so many options for modelling the conditional mean, many of which are available in the tsmodels framework, it was decided to keep this package simpler and avoid the joint estimation of both conditional mean and variance dynamics. While the 2 step estimation approach, whereby the residuals of the conditional mean are passed to the variance dynamics estimation, may be less efficient for small sized datasets, it is expected to be more flexible in what can be achieved. Additionally, the ARCH-in-mean model is no longer available, as it was found to have very limited value within the tsmodels framework or in general, this author’s experience. A separate tsarma package for ARMA(p,q)-X models is however available.

• makes use of automatic differentiation (autodiff) during estimation, via the TMB package. This is in line with similar approaches in other models written in the tsmodels framework. Using autodiff allows for more confident estimation and more accurate standard errors. Autodiff is also used for the more complex inequality constraints of some models which involve evaluation of an expectation, in order to generate more accurate Jacobians of these inequalities for use during optimization.

• adopts a custom approach in optimization using parameter scaling which avoids problems with convergence and numerical stability.

• fully implements and correctly documents a number of sandwich estimators making use of the sandwich package framework (with methods for bread and estfun and meat/meat_HAC functions).

• makes use of S3 methods and classes, abandoning the S4 approach of rugarch. Additionally, while making use of standard methods from the stats package, some of the methods are based on those exported from tsmethods, consistent with other packages in the tsmodels framework. The Appendix provides a table with a hierarchical overview of the main functions and methods currently implemented in the package.

• provides more accurate and informative output, including standard errors, confidence intervals, and a more detailed summary output including a fancy flextable output for the summary method.

• provides numerous fixes and enhancements to the original rugarch code, including the unconditional variance of the exponential GARCH(1,1) model, and better startup conditions for the initialization of the recursion.

• provides more extensive documentation and high coverage of unit tests, with a focus on the accuracy of the estimation and the correctness of the output as well as catching many corner cases.

This vignette provides a mathematical summary of each type of GARCH model used, general assumptions regarding initialization and other estimation details as well as the forecast equation. Separate vignettes are available with code demonstrations.

The density function, based on the distributions available in the tsdistributions package, is expressed in terms of the location, scale, skew and shape parameters {μ_t, σ_t, ζ, ν}, normalized to give zero mean and unit variance:

where μ_t = E(y_t|x_t) and σ_t² = E((y_t − μ_t)²|x_t), with x_t the information set available at the period immediately before time t and may include both external regressors and lagged values of y. Assuming that the distribution of z_t is independent of the conditioning information set¹, then:

which is related to the distribution of y_t by through the scale term:

The parameter set θ, which includes the conditional mean, conditional variance and distributional skew and shape parameters is estimated by minimizing the negative of the log-likelihood:

subject to upper and lower parameters bound constraints as well as the persistence constraint (P < 1) and other non-negativity constraints. Estimation is performed using the nloptr solver with analytic Gradient and Jacobian of the constraints² provided through autodiff in TMB.

When estimating a GARCH model, we need to initialize the variance to some value in order to start the recursion. In the tsgarch package we set all values of σ_t², ∀t ≤ 0 to the sample variance or some other estimator which the user can choose from. Three options are available:

1. unconditional: $\left(\frac{1}{T}\sum_{j=1}^T\varepsilon_j^2\right)^{(\delta/2)}$
2. sample: $\left(\frac{1}{S}\sum_{j=1}^S\varepsilon_j^2\right)^{(\delta/2)}$, where S is the length of the sample to use.
3. backcasting: $\left(\lambda^T\hat\sigma^2 + (1-\lambda)\sum_{j=0}^{T-1}\lambda^j \varepsilon^2{1+j}\right)^{\delta/2},\quad \lambda\in\{0,1\}$

where : $\hat\sigma^2=\frac{1}{T-1}\sum_{j=1}^{T}\varepsilon^2_j$

For the backasting method, setting λ = 1 we obtain the sample variance, else values of λ less than one and greater than zero yields the exponential smoothing backcast estimator. The exponent δ is for the power ARCH type models, otherwise it is set to 2 (variance). One popular commercial econometric software package defaults to backasting using a value of 0.7 for λ.

In addition to the variance, we also need to initialize the value of the ARCH component for t ≤ j, by taking the sample average of the equation or part of the equation. Details of this initialization is provided for each model. Note that for all models, the initialization values of the ARCH(q) components will follow the user choice for the variance initialization.

The persistence (P) of a GARCH model is a measure which quantifies the degree of volatility clustering and rate of decay. It also forms a bounding constraint ( < 1) on the model dynamics in order to ensure stationarity. Another way to think about persistence is by looking at the unconditional variance of a GARCH model which is defined as

where χ̄_k is the sample mean of any external regressors.

Equation illustrates that positivity of the unconditional variance requires P <  = 1, whilst existence of this value requires P < 1, which is not the case for the integrated GARCH model where P = 1 by design. The form that P takes will depend on the type of model, with the formulas provided in Section . Closely related to the persistence is the half-life measure which is defined as the number of periods it takes for a shock to revert half way back to the long run variance, and defined as −log_e(2)/log_e(P).

A special note is warranted for the half-life of the Component GARCH (cgarch) model which is composed of a permanent and transitory component, each of which have a persistence (see ). The permanent component half-life, based on the estimate of ρ, measures the time taken for the long-run influence of a shock in volatility to revert by half towards it’s long run unconditional value, whereas the transitory component half-life accounts for the time taken for a shock’s influence to revert to its long-run rate.

Variance targeting sets the value of the GARCH intercept (ω) to it’s long run estimate as :

where σ̂² is the unconditional variance of ε², consistently estimated by its sample counterpart, P is the model persistence and χ̄_k the mean of the external variance regressors (if present). A review of variance targeting can be found in Francq, Horvath, and Zakoian (2011). In this author’s experience, more than 90% of model estimation problems come from trying to estimate ω as a result of parameter scaling issues. In tsgarch, despite attempts to apply some scaling rules during optimization, failure to converge will sometimes happen, in which case the model will be automatically re-estimated with variance targeting (the output will indicate when this has happened).

Unlike initialization estimators for σ₀² and ε₀², the value of σ̂² and v̂_j is based on the full sample.

Inclusion of additive external regressors has always been a little tricky when it comes to variance modelling due to the non-negativity constraint. This is an issue in all GARCH flavors with the exception of the exponential model. One way to deal with this is to constrain coefficients and regressors to be strictly positive which is not ideal. Another option, which is now offered in tsgarch is to have multiplicative regressors where the intercept is now calculated as follows:

which does not require any bound constraints on either the constant ω or the regressors.

Robert F. Engle and Ng (1993) defined the news impact curve as a way to analyze the effect of news on the conditional variance by keeping constant the information dated t − 2 and earlier. Therefore, the long run variance of the model is used in place of σ_t − 1 and a range of values for ε_t are chosen to show how news of different sign and size impact the current conditional variance. A most interesting example of this found in the Family GARCH model of Hentschel (1995) which accommodates both shifts and rotations in the news impact curve, with the shift factor being the the main source of asymmetry for small shocks, and rotation driving larger shocks. A more detailed exposition, and compact representation of the news impact can be found in Caporin and Costola (2019) who define asymmetry of a GARCH model, for all shocks θ, as the case where:

The tsgarch package makes use of the methods available in the sandwich package to deliver a number of different estimators for the parameter covariance matrix (S).

Define the objective function (Ψ) as the log-likelihood of the GARCH model with respect to the data and parameters:

where $\frac{1}{\sigma_t}g\left(z_t|\zeta,\nu\right)$ is defined in Section . The estimating (or score) function of the objective function is then:

Inference about the parameter set θ relies on a central limit theorem (CLT) with $\sqrt{n}$ consistency:

where $\overset{d}{\to}$ indicates convergence in distribution. The sandwich package defines the sandwich estimator S(θ) as:

where the meat (M) of the sandwich is the variance of the estimating function:

and the bread (B) is the inverse of the expectation of its first derivative (ψ^′):

In tsgarch, the following 4 estimators for the covariance matrix are defined:

This makes use of the analytic hessian (H) at the optimal solution.

The estimating function is essentially the Jacobian of the likelihood at each time step with respect to the parameters. Currently, this is based on the jacobian function from numDeriv, but will be replaced in a future version by the analytic solution from TMB once I figure out how to do this.

In the presence of residual autocorrelation, the HAC estimator is based on the weighted empirical autocorrelations of the empirical estimating functions:

where the weights w can either be fixed or estimated using an appropriate choice of strategies. In the tsgarch package, the bandwidth of Newey and West (1994) is used to automatically calculate the lag and then the weights, based on the bwNeweyWest function from the sandwich package.

The estimation strategy involves 2 passes through the optimizer. During the first pass, the parameters are first estimated using no-scaling. In the second pass, the problem is reformulated based on re-scaling the n parameters and their bounds by the vector $s = \sqrt{\{s^{-1}_{1,1},s^{-1}_{2,2},\ldots,s^{-1}_{n,n}\}}$, where s_i, i are the diagonals of the hessian at the previous solution.

The rescaled parameters and their bounds are passed to the optimizer, whilst the underlying C++ TMB code scales them back in order to perform the likelihood calculations and correctly adjust the analytic derivatives in the presence of this scaling vector.

The reason for performing this extra step, is to avoid some bad solutions as a result of large differences in scaling (particularly with respect to the constant ω), and to avoid issues with the derivatives. The optimal hessian and scores (the estimating function) in the output object are based on the scaled versions from the second pass which are then re-scaled back. This approach was found to offer substantial stability at a relatively low cost. Some justification for this method can be found in Yang and Lee (2010) and Rao (2019).

The tsgarch package implements a number of different flavors of GARCH, including the option of 10 different distributions. The next subsections provide details of each model’s formulation. For each model, the dynamics of the conditional mean take the following form:

where D is one of the available distributions from tsdistributions, ζ the skew parameter, ν the shape parameter and λ an additional shape parameter used in the Generalized Hyperbolic distribution. These additional distributional parameters may be present in some combination or not at all depending on the distribution. It is also possible to set μ to zero and pass a series of pre-filtered residual values from some other model as discussed in Section .

Note that the parameter symbols for the model equations presented in the following subsections will be exactly the same as those output by the package, with the summary method on estimated objects having the option to replace the names of parameters with symbols when transforming to a flextable object (as_flextable method on summary object).

The vanilla GARCH model of Bollerslev (1986) extended the ARCH model of R. F. Engle (1982) to include a moving average term making it more closely aligned with the literature on ARMA processes, and allowing for a wider range of behavior and more persistent volatility.

In the presence of variance regressors, the choice to use multiplicative type intercept means that:

in which case the bounds on ω and ξ_k are mostly free and we don’t have to worry about the positivity of the variance.

where the first equation is used to initialize the GARCH recursion and the second one (identical to the first) is used for the ARCH recursion initialization.

This is the default choice, but other choices such as less than full sample and backcasting are also available and described in Section .

where $\bar\sigma^2 = \frac{\omega + \sum_{k=1}^n\xi_k\bar\chi_k}{1 - P}$³ and represents the long run unconditional variance in the (optional) presence of m variance regressors χ with coefficients ξ (see Section ).

• 1 > α_j > 0
• 1 > β_j > 0
• ω > 0
• P < 1

For higher order models, it is suggested that the constraints on α_j be relaxed to allow for non positive values This can be achieved by changing the lower value in the parmatrix slot of the specification.

The integrated GARCH model of R. F. Engle and Bollerslev (1986) assumes that the persistence P = 1, hence shocks are permanent and the unconditional variance infinite. The motivation behind this model was to capture the long memory behavior observed in some financial time series.⁴ However, Nelson (1990) showed that the IGARCH process with no drift (ω = 0) would converge to zero with probability one. In the presence of a drift term (ω > 0) the process is neither covariance stationary nor does it have well-defined unconditional variance, though it still remains strictly stationary and ergodic. For truly long memory processes, other GARCH models should be considered such as the Fractionally Integrated GARCH (FIGARCH) or Hyperbolic GARCH (HYGARCH) which may be included at a later time.

The persistence in the igarch model is set to 1 and forms a binding constraint on the parameters.

Not defined.

• 1 > α_j > 0
• 1 > β_j > 0
• ω > 0
• P = 1

Same as the the vanilla GARCH forecast in Section .

The Exponentially Weighted Moving Average (ewma) model is a restricted igarch model where the drift term (ω) is set to zero. It has proven popular among some practitioners for it’s simplicity and speed, with coefficients most often hard coded rather than estimated, based on prior knowledge. However, as mentioned in Section the variance will converge to zero in a finite number of steps so it is unlikely to be a good model for anything but very short term forecasting.

The ewma equation is usually written as σ_t² = (1 − λ)ε_t − 1² + λσ_t − 1², but we present below the more general model which shows that it is a restricted igarch model with no drift (although it is always possible to include regressors).

Since this is simply a restricted igarch model, the persistence is set to 1 and forms a binding constraint on the parameters as in .

Not defined.

Same as the the vanilla GARCH forecast in Section with ω set to zero.

The exponential GARCH model of Nelson (1991) allows for asymmetric effects between positive and negative returns, and does not require specific parameter restrictions to ensure positivity of the variance since the modelling is performed on the log variance.

where $z_t = \frac{\varepsilon_t}{\sigma_t}$, with expectation of the absolute moment given by:

For symmetric distributions, the absolute moment is usually available in closed form (such as in the Gaussian case where it is $\sqrt{\frac{2}{\pi}}$). For other distributions this is calculated using Gauss-Kronrod quadrature in the C++ TMB code so that it is forms part of the autodiff tape.

where the first equation is used in the initialization of the GARCH recursion, whilst the second and third equations are for the ARCH initialization.

This is the default choice, but other choices such as less than full sample and backcasting are also available and described in Section .

The persistence has a rather simple form based on the sum of the moving average coefficients.

The unconditional variance of an EGARCH(1,1) model, is given by the following equation:

where g(z_t) = αz_t + γ(|z_t| − E|z_t|) (see He, Terasvirta, and Malmsten (2002)), and E[exp (β^i − 1g(z_t))] is calculated by numerical quadrature:

We approximate the infinite product by truncating i to 1000 which we found is more than sufficient for convergence.

For higher order GARCH model, i.e. p, q > 1, the unconditional variance is approximated via simulation, averaging the variance over 100 simulations of length 25000.

where σ̄², represents the unconditional variance as given in Equation .

• P < 1

For the max(p, q) <  = 1, the forecast for t + h, h > 1 is:

where g(x) = (α₁x + γ₁(|x| − κ)), and κ as in . We approximate the infinite product in the equation by truncating it to 1000 terms. For higher order models, the forecast is approximated via simulation.

The GJR GARCH model of Glosten, Jagannathan, and Runkle (1993) models positive and negative shocks on the conditional variance asymmetrically using a leverage term for past squared, negative innovations via the indicator function I.

where γ_j now represents the leverage term. The indicator function I takes on value 1 for ε ≤ 0 and 0 otherwise. Because of the presence of the indicator function, the persistence of the model now crucially depends on the asymmetry of the conditional distribution.

where κ is the expected value of ε_t ≤ 0. Since this represents the probability of being less than zero, we can work directly with the standardized innovations when using the location scale distributions⁵, z_t:

For symmetric distributions this value is always 0.5, but for skewed distributions this is calculated using Gauss-Kronrod quadrature in the C++ TMB code so that it is forms part of the autodiff tape allowing to also extract the Jacobian of this function for use with the inequality constraint in the nloptr solver.

where $\bar\sigma^2 = \frac{\omega + \sum\limits_{k=1}^{m}{\xi_{k}\bar\chi_k}}{1 - P}$, represents the long run unconditional variance.

• 1 > α_j > 0
• 1 > β_j > 0
• α_j + γ_j > 0
• ω > 0
• P < 1

Note that we also constrain γ_j > −1. The Jacobian of the inequality constraints is calculated either analytically or using autodiff. For higher order models, it is suggested that the constraints on α_j and γ_j be relaxed to allow for non positive values. This can be achieved by changing the lower value in the parmatrix slot of the specification.

where κ is the expected value of ε_t ≤ 0, see .

The asymmetric power ARCH model of Ding, Granger, and Engle (1993) allows for both leverage and the Taylor effect, named after Taylor (1986) who observed that the sample autocorrelation of absolute returns was usually larger than that of squared returns.

where δ ∈ ℝ⁺, being a Box-Cox transformation of σ_t, and γ_j the coefficient in the leverage term.

See Laurent (2004) page 52:

where κ_j is the expected value of the standardized residuals z_t under the Box-Cox transformation of the term which includes the leverage coefficient:

For symmetric distributions there are closed form solutions for this expression, but for skewed distributions this is calculated using Gauss-Kronrod quadrature in the C++ TMB code so that it is forms part of the autodiff tape allowing to also extract the Jacobian of this function for use with the inequality constraint in the nloptr solver.

where $\bar\sigma^{\delta} = \frac{\omega + \sum\limits_{k=1}^{m}{\xi_{k}\bar\chi_k}}{1 - P}$, represents the long run unconditional volatility raised to the power of δ.

• 1 > α_j > 0
• 1 > β_j > 0
• |γ_j| < 1
• δ > 0
• ω > 0
• P < 1

Non-negativity constraints on α_j may be relaxed for higher order models by changing the lower parameter in the parmatrix of the specification object prior to estimation.

where κ_j = E(|z_t − j| − γ_jz_t − j)^δ, see .

The family GARCH model of Hentschel (1995) is a large omnibus model which subsumes some of the most popular GARCH models. It allows for both shifts and rotations in the news impact curve, where the shift is the main source of asymmetry for small shocks while rotation drives larger shocks. The following restrictions in the parameters leads to specific models:

• GARCH: δ = 2, η_j = γ_j = 0
• Absolute Value GARCH: δ = 1, |γ_j| ≤ 1
• GJR GARCH: δ = 2,η_j = 0
• Threshold GARCH: δ = 1,η_j = 0,|γ_j| ≤ 1
• Nonlinear GARCH: γ_j = η_j = 0
• Nonlinear Asymmetric GARCH: δ = 2,γ_j = 0
• APARCH: η_j = 0, |γ_j| ≤ 1

The original formulation, takes a slightly different form:

which accommodates the exponential GARCH model as well as a more flexible transformation, albeit more computationally demanding. The tsgarch package does not adopt this more general formulation, instead setting λ = δ.

where κ_j is the expected value of the standardized residuals z_t under the Box-Cox transformation of the absolute value asymmetry term

There are no simple closed form solutions for this expression so it is calculated using Gauss-Kronrod quadrature in the C++ TMB code so that it is forms part of the autodiff tape allowing to also extract the Jacobian of this function for use with the inequality constraint in the nloptr solver.

where $\bar\sigma^{\delta} = \frac{\omega + \sum\limits_{k=1}^{m}{\xi_{k}\bar\chi_k}}{1 - P}$, represents the long run unconditional volatility raised to the power of δ.

• 1 > α_j > 0
• 1 > β_j > 0
• |γ_j| < 1
• −10 < |η_j| < 10
• δ > 0
• ω > 0
• P < 1

Non-negativity constraints on α_j may be relaxed for higher order models by changing the lower parameter in the parmatrix of the specification object prior to estimation.

where κ_j is the expected value of the standardized residuals z_t under the Box-Cox transformation of the absolute value asymmetry term, see .

The model of Lee and Engle (1999) decomposes the conditional variance into a permanent and transitory component so as to investigate the long- and short-run movements of volatility affecting securities.

The process can be re-written in an alternative form to better highlight the decomposition of the permanent and transitory components, shown below for the Component GARCH(1,1) model:

The parameters α and ϕ correspond to immediate impacts of volatility shocks (ε_t − j² − σ_t − j²) on the transitory and permanent components, whereas (α + β) and ρ measure the persistence of the transitory and permanent components, respectively. The model as currently implemented allows for higher order in the transitory component but no higher orders in the permanent component. The regressors only enter through the permanent component equation.⁶

The transitory component is initialized to 0, whilst the squared residuals, variance and permanent component are initialized to the sample variance. As a result, the initial squared residuals and variance cancel out during initialization.

Since the component GARCH model can be represented as a restricted GARCH(2,2) model, we derive the news impact curve using this representation to arrive at the following equation:

where the unconditional variance $\bar\sigma^2 = \frac{\omega + \sum\limits_{k=1}^{m}{\xi_{k}\bar\chi_k}}{1-\rho}$.

• 1 > α_j > 0
• 1 > β_j > 0
• 1 > ϕ > 0
• 1 > ρ > 0
• ω > 0
• 1 > P^P > P^T > 0
• β > ϕ

where P_T is the transitory persistence as defined in . As h grows larger, the forecast converges to $\frac{\omega+\sum\limits_{k=1}^{m}{\xi_{k}\chi_{k,t+h}}}{1-\rho}$.

The tsgarch package includes a number of methods which can be used with either the tsgarch.spec or tsgarch.estimate objects. These include standard methods from the stats package such residuals, fitted, sigma, logLik, AIC, BIC, confint, vcov, coef, predict and simulate. Additionally, the following methods from the tsmethods package are included:

• tsfilter: online filtering
• tsbacktest: backtesting
• tsprofile: parameter profiling through simulation/estimation
• unconditional: unconditional variance
• tsequation: model equations (used mainly by the as_flextable method)
• pit: probability integral transform
• persistence: model persistence

Also included from the sandwich package are methods for estfun and bread.

The package is still in active development. Additional models may be added, porting some of the remaining rugarch models, but the focus is likely to be on models for high frequency data and porting the multivariate models from rmgarch.

A separate package called tstests includes a number of tests which were previously part of rugarch and will be released to CRAN soon.