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:


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)\)).


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.


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.


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


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]#
  • returns (np.ndarray) – (T,) array of T returns

  • zero_mean (bool) – whether to use a zero mean returns model. Default to False.


returns is best to be percentage returns for optimization

fit() Parameters[source]#

Estimates the GARCH(1,1) parameters via MLE



Return type:


preparation() List[float][source]#

Prepare starting values.


list of starting values

Return type:


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.

  • params (np.ndarray) – [mu, omega, alpha, (gamma), beta]

  • backcast (float) – backcast value

  • var_bounds (np.ndarray) – variance bounds


negative log-likelihood

Return type:


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

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

  • resids (np.ndarray) – residuals to use in computing log-likelihood.

  • sigma2 (np.ndarray) – conditional variance of residuals.



Return type:


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

Computes the variances conditional on given parameters

  • params (List[float]) – [omega, alpha, beta]

  • resids (np.ndarray) – residuals from mean equation

  • backcast (float) – backcast value

  • var_bounds (np.ndarray) – variance bounds


conditional variance

Return type:


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.


resids (np.ndarray) – residuals from the mean model


[omega, alpha, beta]

Return type:


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.


resids (np.ndarray) – residuals from the mean model.


an array where each row contains the lower and upper bounds for the conditional variance at each time point.

Return type:


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

Adjust the conditional variance at time t based on its bounds

  • sigma2 (float) – conditional variance

  • var_bounds (np.ndarray) – lower and upper bounds


adjusted conditional variance

Return type:


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

Computes the starting value for estimating conditional variance.


resids (np.ndarray) – residuals


initial value from backcasting

Return type:


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

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

  • resids (np.ndarray) – Residuals from the model.

  • initial_value (float) – Initial value for the conditional variance.

  • lam (float) – Decay factor for the EWMA.


Array containing the conditional variance estimates.

Return type:


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


Let’s import the dataset.

>>> import pandas as pd
>>> data_url = ""
>>> 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 =
>>> pprint(res)

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


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.