# Algorithm to fit AR(1)/GARCH(1,1) model of log-returns

I am fitting numerically an AR(1)/GARCH(1,1) process to index and stock log-returns, $r_t=\log(P_t/P_{t-1})$, where $P_t$ is the price at time $t$, and thus far am not clear on where the observed log returns would be used in an algorithm. Several author groups have described (some in part) the components of the AR(1)/GARCH(1,1) approach, for example:

E. Zivot: $$r_t=\mu + \phi (r_{t-1} - \mu ) + \epsilon_t$$

Rachev et al: $$\begin{split} r_t&=\mu+\phi r_{t-1}\\ \epsilon_t&=\sigma_t \delta_t \quad \quad (\delta_t \textrm{ is an innovation})\\ \sigma_t &= \sqrt{\alpha_0 + \alpha_1 \epsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2}\\ \end{split}$$

Brummelhuis & Kaufman: $$\begin{split} X_t&=\mu_t+ \sigma_t \epsilon_t\\ \mu_t&=\lambda X_{t-1}\\ \sigma_t &= \sqrt{\alpha_0 + \alpha_1 (X_{t-1}-\mu_{t-1})^2 + \beta_1 \sigma_{t-1}^2}\\ \end{split}$$

Jalal & Rockinger: $$\begin{split} \mu_t&=\phi X_{t-1}\\ \epsilon_t &=X_t - \mu_t\\ \sigma_t &= \sqrt{\alpha_0 + \alpha_1 \epsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2}\\ \end{split}$$

My approach:

$$\begin{split} \mu_t&=\mu+\phi r_{t-1}\\ \epsilon_t&=r_t - \mu_t\\ \sigma_t &= \sqrt{\alpha_0 + \alpha_1 \epsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2}\\ \end{split}$$

Given the multiple descriptions above, my interpretation for an algorithm would be:

Algorithm for AR(1)/GARCH(1,1):

1. Initialize $\mu=0$, $\sigma_1 = 1$, $\epsilon_1=0$, $\mu=\phi=\alpha_0=\alpha_1=\beta_1=U(0,1)* 0.01$
2. For $t$ = 2 to $T$:
3. $\quad \mu_t = \mu + \phi r_{t - 1}\quad \quad (\textrm{Log-returns lag-1 input here})$
4. $\quad \epsilon_t = r_t - \mu_t \quad \quad (\textrm{Log-returns lag-0 input here})$
5. $\quad \sigma_t= \sqrt{\alpha_0 + \alpha_1 \epsilon_{t-1}^ 2 + \beta_1 \sigma_{t - 1} ^ 2}$

6. $\quad \hat{r}_t = \phi \mu_{t - 1} + \sigma_t \epsilon_t \quad \quad \textrm{OR}: \hat{r}_t = \mu + \phi \mu_{t - 1} + \sigma_t \epsilon_t \quad \quad ???$

7. Next $t$
8. Calculate residual, $e_t=r_t-\hat{r}_t$
9. Determine $MSE=\frac{1}{T}\sum_t e_t^2$

The algorithm proposed above is essentially the recursive part to calculate the predicted log-returns $\hat{r}_t$ from the input observed returns $r_t$. Innovations or random quantiles from a probability distribution [such as N(0,1) or $t(\nu)$] would not be employed here since we are fitting a model, not simulating. During each iteration, the goodness-of-fit based on proposed parameters for the objective function would be based on $MSE=\frac{1}{T}\sum_t e_t$, which would be minimized via an optimization technique using non-linear regression, finite differencing, or MLE. Metaheuristics could be used as well where initialization of chromosome (particle) values for $(\mu,\phi,\alpha_0,\alpha_1,\beta_1)$ would occur at the first generation.

In terms of comparing results based R, MATLAB, SAS, etc. the parameterization would be:

mu=$\mu$

ar1=$\phi$

garch0=$\alpha_0$

garch1=$\alpha_1$

garch2=$\beta_1$

I am not sure whether the unconditional mean $\mu$ would be needed in line 6 of the algorithm, however. Please comment on correctness of the algorithm, and possibly suggest coding changes. Again, the goal is to algorithmically use numerical methods to solve for the parameters, not R, MATLAB, or SAS, etc.

• Nice post! Very useful. Dec 29, 2013 at 16:35
• Hi PEL, could you split this post into a question and an answer? Now the question seems unanswered while it actually has one. Aug 1, 2016 at 20:35
• The first line of Rachev's formula seems to be missing an error term. Nov 8, 2018 at 12:22
• Would anyone happen to know how I could code this loglike in MATLAB, is thre a function that I can use etc? Jan 21, 2019 at 9:47
• You'd have to implement some type of maximization method in MATLAB, perhaps a derivative free method. Either that, or take the partial derivatives of $loglik$ w.r.t. $\sigma_t$ and $\epsilon_t$, set to zero, then iterate. Lookup derivative-free function maximization in MATLAB. If you can access the genetic algorithm (GA) methods for MATLAB, then treat $loglik$ as the fitness, and ensure it increases with increasing iterations.
– user6430
Jan 22, 2019 at 16:04

SOLUTION: Let $r_t$ be the log-return at time $t$, and $\hat{r}_t$ be the predicted log-return from the regression model.

1. Initialize $loglik(0:T)=0$,$\epsilon_1=0$, $\sigma_1 = 0$, $\mu=U(0,1)*0.0001,\phi=U(0,1)*0.01$, $\alpha_0=U(0,1)*0.00002,\alpha_1=U(0,1)*0.01,\beta_1=0.9 + U(0,1)* 0.01$, $B=10,000$
2. For $b$ = 1 to $B$
3. $\quad$ For $t$ = 2 to $T$:
4. $\quad\quad \hat{r}_t = \mu + \phi r_{t - 1}$
5. $\quad\quad \epsilon_t = r_t - \hat{r}_t$
6. $\quad\quad \sigma^2_t = \alpha_0 + \alpha_1 \epsilon_{t-1} ^ 2 + \beta_1 \sigma_{t-1} ^ 2$
7. $\quad$ Next $t$
8. $\quad loglik_b = \sum_t^T \left\{ - \frac{1}{2} \log(\sigma^2_t) - \frac{1}{2}\epsilon_t^2 / \sigma^2_t \right\}$
9. Next $b$

At each iteration $b$, the log-likelihood function value ($loglik_b$) is determined, which needs to be maximized over the $B$ iterations, that is, ensure it increases during iterations.

Parameter results after $B=10,000$ iterations:

$\mu=0.000745264$

$\phi=-0.110093524$

$\alpha_0=0.00000401$

$\alpha_1=0.117396192$

$\beta_1=0.849342676$

Below is a plot of the fitted SP500 index (Yahoo ^GSPC) volatility, $\sigma_t$, for dates 1/2/2001 to 10/31/2013. The subprime mortgage crisis in 2008 and debt ceiling in 2011 are very visible. In addition, below is a plot of the observed (input) log-returns (blue), $r_t$, and predicted GARCH-filtered returns (red), $\hat{r}_t$, revealing that volatility clustering was removed. 