Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Consider some ARCH($p$) process, say ARCH($1$),

$\varepsilon_t=\sigma_t\cdot \eta_t$

where

$\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2$

with a Gaussian (strong) white noise $(\eta_t)$.

```> n=500
> a1=0.8
> a2=0.0
> w= 0.2
> set.seed(1)
> eta=rnorm(n)
> epsilon=rnorm(n)
> sigma2=rep(w,n)
> for(t in 3:n){
+ sigma2[t]=w+a1*epsilon[t-1]^2+a2*epsilon[t-2]^2
+ epsilon[t]=eta[t]*sqrt(sigma2[t])
+ }
> par(mfrow=c(1,1))
> plot(epsilon,type="l",ylim=c(min(epsilon)-.5,max(epsilon)))
> lines(min(epsilon)-1+sqrt(sigma2),col="red")```

(the red line is the conditional variance process).

```> par(mfrow=c(1,2))
> acf(epsilon,lag=50,lwd=2)
> acf(epsilon^2,lag=50,lwd=2)```

We did mention in class that if $(\varepsilon_t)$ a ARCH($1$), then $(\varepsilon_t^2)$ is an AR($1$) process. So a first idea is to consider a regression, as we did for Gaussian AR($1$)

```> db=data.frame(Y=epsilon[2:n]^2,X1=epsilon[1:(n-1)]^2)
> summary(lm(Y~X1,data=db))

Call:
lm(formula = Y ~ X1, data = db)

Residuals:
Min      1Q  Median      3Q     Max
-2.4538 -0.3618 -0.2626  0.0935  9.3667

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.34963    0.04342   8.052 6.08e-15 ***
X1           0.31123    0.04262   7.303 1.13e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8413 on 497 degrees of freedom
Multiple R-squared:  0.0969,	Adjusted R-squared:  0.09508
F-statistic: 53.33 on 1 and 497 DF,  p-value: 1.129e-12```

There is some significant autocorrelation here. But since our vectors cannot be considered as Gaussian, using least squares is perhaps not the best strategy. Actually, if our series is not Gaussian, it is still conditionally Gaussian, since we assumed that $(\eta_t)$ is a Gaussian (strong) white noise,

$\varepsilon_t\vert \underline{\boldsymbol{\epsilon}}_{t-1} \sim \mathcal{N}(0,\sigma^2_t)$

The likelihood is then

$\mathcal{L}=\prod_{t=1}^T \frac{1}{\sqrt{2\pi \sigma^2_t}}\exp\left(-\frac{\varepsilon_t^2}{2\sigma_t}\right)$

and the log-likelihood is

$\log \mathcal{L}=-\frac{1}{2} \sum_{t=1}^T \log (\sigma_t)-\frac{1}{2}\sum_{t=1}^T \frac{\varepsilon_t^2}{\sigma_t}$

And a natural idea is to define

$(\widehat{\omega},\widehat{\alpha})\in\text{argmax}\{\log \mathcal{L}(\omega,\alpha)\}$

The code is simply

```> X=epsilon
> loglik=function(param){
+ w=exp(param[1])
+ a1=exp(param[2])
+ s2=rep(w,n)
+ for(t in 2:length(X)){s2[t]=w+a1*X[t-1]^2}
+ logL=-.5*sum(log(s2))-.5*sum(X^2/s2)
+ return(-logL)
+ }
> OPT=optim(par=
+ coefficients(lm(Y~X1,data=db)),fn=loglik)
> exp(OPT\$par)
(Intercept)          X1
0.2482241   0.5858578```

(since the parameters have to be positive, we assume here that they can be written as the exponential of some real values). Observe that those values are closer to the one used to generate our time series.

If we use R functions to estimate those parameters, we get

```> library(tseries)
> summary(garch(epsilon,c(0,1)))
...

Call:
garch(x = epsilon, order = c(0, 1))

Model:
GARCH(0,1)

Residuals:
Min       1Q   Median       3Q      Max
-2.87023 -0.60836 -0.03426  0.66648  3.48443

Coefficient(s):
Estimate  Std. Error  t value Pr(>|t|)
a0   0.24959     0.02470   10.104  < 2e-16 ***
a1   0.58306     0.09737    5.988 2.13e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1```

so that the confidence interval for $\alpha$ is

```> summary(garch(epsilon,c(0,1)))\$coef[2,1]+
+ c(-1.96,1.96)*summary(garch(epsilon,c(0,1)))\$coef[2,2]
[1] 0.3922088 0.7739088```

Actually, since our main interest is this  $\alpha$ parameter, it is possible to use profile likelihood techniques,

```> proflik=function(a){
+ loglik=function(w){
+ s2=rep(w,n)
+ for(t in 2:length(X)){s2[t]=w+a*X[t-1]^2}
+ logL=-.5*sum(log(s2))-.5*sum(X^2/s2)
+ return(-logL)}
+ return(-optim(par=.3,fn=loglik)\$value)}

> A=seq(0,2,by=.05)
> P=Vectorize(proflik)(A)
> par(mfrow=c(1,1))
> plot(A,P,type="l")
> OPT=optimize(function(x) -proflik(x), interval=c(0,2))
> t=-OPT\$objective-qchisq(.95,df=1)
> abline(h=t,col="red")
> ainf=uniroot(function(x) proflik(x)-t,c(0,OPT\$minimum))\$root
> asup=uniroot(function(x) proflik(x)-t,c(OPT\$minimum,2))\$root
>  abline(v=ainf,lty=2)
>  abline(v=asup,lty=2)```

Of course, all those techniques can be extended to higher order ARCH processes. For instance, if we assume that we have a ARCH($2$) time series

$\varepsilon_t=\sigma_t\cdot \eta_t$

where now

$\sigma^2_t=\omega+\alpha_1 \varepsilon_{t-1}^2+\alpha_2 \varepsilon_{t-2}^2$

with a Gaussian (strong) white noise $(\eta_t)$. The log-likelihood is still

$\log \mathcal{L}=-\frac{1}{2} \sum_{t=1}^T \log (\sigma_t)-\frac{1}{2}\sum_{t=1}^T \frac{\varepsilon_t^2}{\sigma_t}$

and we can define

$(\widehat{\omega},\widehat{\alpha}_1,\widehat{\alpha}_2)\in\text{argmax}\{\log \mathcal{L}(\omega,\alpha_1,\alpha_2)\}$

The code above can be changed, to take into account this additional component,

```> db=data.frame(Y=epsilon[3:n]^2,
+ X1=epsilon[2:(n-1)]^2,
+ X2=epsilon[1:(n-2)]^2)
> X=epsilon
> loglik=function(param){
+ w=exp(param[1])
+ a1=exp(param[2])
+ a2=exp(param[3])
+ s2=rep(w,n)
+ for(t in 3:length(X)){s2[t]=w+a1*X[t-1]^2+a2*X[t-2]^2}
+ logL=-.5*sum(log(s2))-.5*sum(X^2/s2)
+ return(-logL)
+ }
> OPT=optim(par=
+ coefficients(lm(Y~X1+X2,data=db)),fn=loglik)
> exp(OPT\$par)
(Intercept)          X1          X2
0.22710526  0.59475474  0.04741294```

We can also consider some Generalized ARCH process, e.g. a GARCH($1$,$1$),

$\varepsilon_t=\sigma_t\cdot \eta_t$

where now

$\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2+\beta \sigma_{t-1}$

Again, maximum likelihood techniques can be used. Actually, we can also code Fisher-Scoring algorithm, since (in a very general context)

$\frac{\partial \log \mathcal{L}}{\partial \boldsymbol{\theta}}=-\frac{1}{2}\sum_{t=1}^T \frac{\partial \sigma_t^2}{\partial \boldsymbol{\theta}}\cdot \frac{1}{\sigma_t^2}\left(1-\frac{\varepsilon^2_t}{\sigma_t^2}\right)$

with here $\boldsymbol{\theta}=(\omega,\alpha,\beta)$. Using a standard gradient descent algorithm, we get the following estimate for our GARCH process,

```> X=epsilon
> theta=c(.2,.2,.2)
> G=rep(1,3)
> n=length(X)
> j=1
> while(sum(G^2)>1e-12){
+ s2=rep(theta[1],n)
+ for (i in 2:n){s2[i]=theta[1]+theta[2]*X[(i-1)]^2+theta[3]*s2[(i-1)]}
+ z=(X^2-s2)/s2^2
+ V=cbind(z[2:n],z[2:n]*X[1:(n-1)]^2,z[2:n]*s2[1:(n-1)])
+ H=(t(V)%*%V)
+ G=apply(V,2,sum)
+ theta=theta+solve(H)%*%G
+ j=j+1}
> as.numeric(theta)
[1] 0.20372918 0.59183911 0.08936159```

The interesting point, here, is that we also derive the (asymptotic) variance

```> (stdev=sqrt(diag(solve(H))))
[1] 0.01849067 0.04950477 0.02937233```