# 10  Univariate volatility

Univariate volatility modelling is designed to provide an in-sample description of the stochastic processes governing volatility. More important, for our purpose here, it also gives one day ahead volatility forecasts which we use for risk calculations.

The mathematics of univariate volatility models are discussed in Chapter 2 of the financial risk forecasting book

We use maximum likelihood methods for estimating volatility models, and use the R package `rugarch` for the actual implementation. See the package documentation for more details, manual and a more detailed vignette. It is developed by Alexios Galanos and is constantly maintained and updated. The development code can be found on github.

This chapter uses a long data sample and all the estimation works without problems. In Chapter 11, we will see what can go wrong with volatility models and how that can be fixed.

## 10.1 Libraries

We used two libraries, `rugarch` for volatility forecasting and `car` for QQ plots.

``````source('functions.r')
suppressPackageStartupMessages(library(rugarch))
suppressPackageStartupMessages(library(car))``````

## 10.2 Data

The data we use is returns on the S&P500 index, and for convenience, put them into variable `y`.

``````load('data/data.RData')
y=data\$sp500\$y
y=y-mean(y)``````

## 10.3 What to do with the mean?

One should either remove the mean — de-mean — of the returns, or estimate it along with the other parameters. We can even forecast the mean along with volatility, as we show below.

But since the mean is so small, can usually get away with ignoring it if the sample size is large. Chapter 11 shows how that can fail.

## 10.4 EWMA

Start with the exponentially weighted moving average model, EWMA. We estimate the in-sample EWMA volatility and plot the returns along with +/- two standard deviations.

``````EWMA = vector(length=length(y))
lambda = 0.94
EWMA[1] = var(y)
for (i in 2:length(y)){
EWMA[i] = lambda * EWMA[i-1]+
(1-lambda) * y[i-1] ^2
}
EWMA=sqrt(EWMA)

par(mar=c(2,3.5,1,0))
matplot(
cbind(y,2*EWMA,-2*EWMA),
type='l',
lty=1,
col=c("black","red","red"),
bty='l',
ylab="",
main="SP500 returns with += 2 EWMA sd",
las=1
)``````

## 10.5`rugarch`

### 10.5.1 Setup

To estimate a univariate GARCH model with `rugarch` we need to follow two steps:

1. Create an object of the `uGARCHspec` class which is the specification of the model you want to estimate. This includes the type of model (standard, asymmetric, power, etc.), the GARCH order, the distribution, and model for the mean;
2. Fit the specified model to the data.

### 10.5.2`ugarchspec()`

First let’s take a look at `ugarchspec()`, which is the function to create an instance of the `uGARCHspec` class. The complete syntax with its default values is:

``````ugarchspec(
variance.model = list(
model = "sGARCH",
garchOrder = c(1, 1),
submodel = NULL,
external.regressors = NULL,
variance.targeting = FALSE
),
mean.model = list(
armaOrder = c(1, 1),
include.mean = TRUE,
archm = FALSE,
archpow = 1,
arfima = FALSE,
external.regressors = NULL,
archex = FALSE
),
distribution.model = "norm",
start.pars = list(),
fixed.pars = list()
)``````
``````
*---------------------------------*
*       GARCH Model Spec          *
*---------------------------------*

Conditional Variance Dynamics
------------------------------------
GARCH Model     : sGARCH(1,1)
Variance Targeting  : FALSE

Conditional Mean Dynamics
------------------------------------
Mean Model      : ARFIMA(1,0,1)
Include Mean        : TRUE
GARCH-in-Mean       : FALSE

Conditional Distribution
------------------------------------
Distribution    :  norm
Includes Skew   :  FALSE
Includes Shape  :  FALSE
Includes Lambda :  FALSE ``````

You can check the details of this function and what every argument means by running `?ugarchspec`.

For the purposes of our course, we are going to focus on three arguments of the `ugarchspec()` function:

#### 10.5.2.1`variance.model`

This argument takes a `list` with the specifications of the GARCH model. Its most important components are:

• `model`: Specify the type of model, currently implemented models are “sGARCH”, “fGARCH”, “eGARCH”, “gjrGARCH”, “apARCH” and “iGARCH” and “csGARCH”
• `garchOrder`: Specify the ARCH(q) and GARCH(p) orders

Other components include options to do `variance.targeting`or `external regressors`.

#### 10.5.2.2`mean.model`

This argument takes a `list` with the specifications of the mean model, if assumed to have one. A traditional assumption is that the model has zero mean, in which case we specify
`armaOrder = c(0,0), include.mean = FALSE`
However, it is also common to assume that the mean follows a certain ARMA process.

#### 10.5.2.3`distribution.model`

Here we can specify the conditional density to use for innovations. The default is a normal distribution, but we can specify some others, perhaps the Student-t distribution `std`, or the asymmetric Student-t distribution `sstd`.

### 10.5.3`ugarchfit()`

Once the specification has been created, you can fit this to the data using
`ugarchfit( spec = ..., data = ...)`
The result will be an object of the class `uGARCHfit`, which is a list that contains useful information which will be shown below.

### 10.5.4 Optimiser failure — Hybrid

`ugarchfit()` performs maximum likelihood to fit the specified model to our data. This optimisation problem required a `solver`, which is the numerical method used to maximise the likelihood function. See discussion in Section 11.6. In some cases, the default solver in `ugarchfit()` can fail to converge. Then we need the option `solver = "hybrid"` when fitting the model, since this will try different optimisation methods in case the default one does not converge.

## 10.6 The models

We will put the results into a list called `results()` so we can automatically process them after.

``Results=list()``

### 10.6.1 Gaussian ARCH(1) and no mean

The `garchOrder` in `variance.model` will be set as `c(1,0)`. Note that even if we only have one component, we need to put it inside a `list`.

``````spec.1 = ugarchspec(
variance.model = list(
garchOrder = c(1,0)
),
mean.model = list(
armaOrder = c(0,0),
include.mean = FALSE
)
)``````

We can call `spec.1` to see what is inside:

``spec.1``
``````
*---------------------------------*
*       GARCH Model Spec          *
*---------------------------------*

Conditional Variance Dynamics
------------------------------------
GARCH Model     : sGARCH(1,0)
Variance Targeting  : FALSE

Conditional Mean Dynamics
------------------------------------
Mean Model      : ARFIMA(0,0,0)
Include Mean        : FALSE
GARCH-in-Mean       : FALSE

Conditional Distribution
------------------------------------
Distribution    :  norm
Includes Skew   :  FALSE
Includes Shape  :  FALSE
Includes Lambda :  FALSE ``````

Now we can fit the specified model to our data using the `ugarchfit()` function:

``Results\$ARCH1 = ugarchfit(spec = spec.1, data = y)``
``````Warning message in .sgarchfit(spec = spec, data = data, out.sample = out.sample, :
“
ugarchfit-->warning: solver failer to converge.”``````

That fails, need the hybdid solver.

``Results\$ARCH1 = ugarchfit(spec = spec.1, data = y,solver = "hybrid")``

We can check the class of this new object:

``class(Results\$ARCH1)``
'uGARCHfit'
##### 10.6.1.0.1 Printing and plotting

To get all the results, simply do `print(Results\$ARCH1)`
and `plot(Results\$ARCH1)`

Objects from the `uGARCHfit` class have two `slots` (an R term for an object of what is known as the S4). You can think of every slot as a list of components. The two slots are:

1. `@model`;
2. `@fit`.

To access a slot you need to use the syntax: `object@slot`

##### 10.6.1.0.2`@model`

The `@model` slot includes all the information that was needed to estimate the GARCH model. This includes both the model specification and the data. To see every that is included in this slot, you can use the `names()` function:

``names(Results\$ARCH1@model)``
1. 'modelinc'
2. 'modeldesc'
3. 'modeldata'
4. 'pars'
5. 'start.pars'
6. 'fixed.pars'
7. 'maxOrder'
8. 'pos.matrix'
9. 'fmodel'
10. 'pidx'
11. 'n.start'

And you can access each element with the `\$` sign. For example, let’s see what is inside `pars`: (only rows 7-10).

``Results\$ARCH1@model\$pars[7:10,]``
A matrix: 4 × 6 of type dbl
LevelFixedIncludeEstimateLBUB
omega8.889786e-050112.220446e-160.1478908
alpha14.407729e-010110.000000e+001.0000000
beta0.000000e+00000 NA NA
gamma0.000000e+00000 NA NA

We have a matrix with all the possible parameters to be used in fitting a GARCH model. We see there is a value of 1 for the parameters that were included in our GARCH specification (`omega`, `alpha1` and `beta1`). The matrix also includes what are the lower and upper bounds for these parameters, which has nothing to do with our particular data, only with what is expected from the model.

We can also see the model description:

``Results\$ARCH1@model\$modeldesc``
\$distribution
'norm'
\$distno
1
\$vmodel
'sGARCH'

And inside `Results\$ARCH1@model@modeldata` we will find the data vector we have used when fitting the model. You can think of the `@model` slot as everything that R needs to know in order to be able to estimate the model.

##### 10.6.1.0.3`@fit`

Inside the `fit` slot we will find the estimated model, including the coefficients, likelihood, fitted conditional variance, and more. Let’s check everything that is included:

``names(Results\$ARCH1@fit)``
1. 'hessian'
2. 'cvar'
3. 'var'
4. 'sigma'
5. 'condH'
6. 'z'
7. 'LLH'
8. 'log.likelihoods'
9. 'residuals'
10. 'coef'
11. 'robust.cvar'
12. 'A'
13. 'B'
14. 'scores'
15. 'se.coef'
16. 'tval'
17. 'matcoef'
18. 'robust.se.coef'
19. 'robust.tval'
20. 'robust.matcoef'
21. 'fitted.values'
22. 'convergence'
23. 'message'
24. 'kappa'
25. 'persistence'
26. 'timer'
27. 'ipars'
28. 'solver'

Let’s see the coefficients, along with their standard errors, in `matcoef`. There are 2 ways to do that:

``````Results\$ARCH1@fit\$matcoef
coef(Results\$ARCH1)``````
A matrix: 2 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega8.889786e-052.506444e-0635.467730
alpha14.407729e-013.359505e-0213.120170
omega
8.8897858096543e-05
alpha1
0.440772863748013

We can also see the log-likelihood. Again two ways

``````Results\$ARCH1@fit\$LLH
likelihood(Results\$ARCH1)``````
15519.6114054872
15519.6114054872

This makes it easy to extract the estimated conditional volatility:

``plot(Results\$ARCH1@fit\$sigma, type = "l")``

### 10.6.2 Gaussian GARCH(1,1)

``````spec.2 = ugarchspec(
variance.model = list(garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE)
)
Results\$GARCH11 = ugarchfit(spec = spec.2, data = y)
Results\$GARCH11@fit\$matcoef``````
A matrix: 3 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega2.459018e-067.645007e-07 3.2165020.001297638
alpha11.238915e-019.913236e-0312.4975850.000000000
beta18.557904e-011.073671e-0279.7069520.000000000

### 10.6.3 tGARCH(1,1)

We can see that the coefficients now include `shape`, which is the estimated degrees of freedom:

``````spec.3 = ugarchspec(
variance.model = list(garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "std"
)
Results\$tGARCH11 = ugarchfit(spec = spec.3, data = y)
Results\$tGARCH11@fit\$matcoef``````
A matrix: 4 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega1.505701e-061.959878e-06 0.76826274.423311e-01
alpha11.240716e-013.451712e-02 3.59449393.250227e-04
beta18.706782e-013.112403e-0227.97446870.000000e+00
shape6.061799e+008.255397e-01 7.34283142.091660e-13

### 10.6.4 Asymmetric tGARCH(1,1)

We can have a skew Student-t conditional distribution.

``````spec.4 = ugarchspec(
variance.model = list(garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "sstd"
)
Results\$atGARCH11 = ugarchfit(spec = spec.4, data = y)
Results\$atGARCH11@fit\$matcoef``````
A matrix: 5 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega1.494685e-061.844169e-06 0.81049234.176573e-01
alpha11.227309e-013.209377e-02 3.82413551.312318e-04
beta18.710876e-012.926484e-0229.76567200.000000e+00
skew8.770661e-011.669387e-0252.53819910.000000e+00
shape6.481471e+009.600333e-01 6.75129831.465272e-11

### 10.6.5 Gaussian apARCH model

The parameter list will now include `gamma1` and `delta`, which are part of the apARCH model specification:

``````spec.5 = ugarchspec(
variance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
)

Results\$APGARCH11 = ugarchfit(spec = spec.5, data = y)
Results\$APGARCH11@fit\$matcoef``````
A matrix: 5 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega1.453484e-074.735252e-07 0.30694960.75888173
alpha16.084138e-023.590656e-02 1.69443630.09018241
beta18.728381e-014.150876e-0221.02780640.00000000
gamma14.460231e-011.840133e-01 2.42386270.01535641
delta2.554894e+005.051892e-0250.57300290.00000000

### 10.6.6 t-apARCH model

``````spec.6 = ugarchspec(
variance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "std"
)
Results\$tAPGARCH11 = ugarchfit(spec = spec.6, data = y)
Results\$tAPGARCH11@fit\$matcoef``````
A matrix: 6 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega0.00025254320.0001215195 2.0782110.03768989
alpha10.10275240460.0092548707 11.1025220.00000000
beta10.89697157000.0097044547 92.4288480.00000000
gamma10.99999999000.00035352162828.6814360.00000000
delta1.00658913500.0907703454 11.0894050.00000000
shape6.51335738660.5743629275 11.3401420.00000000

### 10.6.7 Asymmetric t-apARCH model

``````spec.7 = ugarchspec(
variance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "sstd"
)
Results\$atAPGARCH11 = ugarchfit(spec = spec.7, data = y)
Results\$atAPGARCH11@fit\$matcoef``````
A matrix: 7 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
omega0.00029023210.0001322079 2.195270.02814421
alpha10.10350778950.0095115201 10.882360.00000000
beta10.89752094540.0101043456 88.825240.00000000
gamma10.99999999000.00038335872608.523020.00000000
delta0.97806184290.0842650532 11.606970.00000000
skew0.84873256260.0161120763 52.676800.00000000
shape7.03822139200.6802275855 10.346860.00000000

### 10.6.8 Asymmetric t-apARCH model with ARMA(1,1) mean

``````spec.8 = ugarchspec(
variance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
distribution.model = "sstd"
)
Results\$Mean.atAPGARCH11 = ugarchfit(spec = spec.8, data = y,solver = "hybrid")
Results\$Mean.atAPGARCH11@fit\$matcoef``````
A matrix: 10 × 4 of type dbl
Estimate Std. Error t valuePr(>|t|)
mu-8.298863e-179.366160e-05-8.860476e-131.000000e+00
ar1 4.146306e-015.708996e-02 7.262759e+003.792522e-13
ma1-4.871041e-015.568657e-02-8.747245e+000.000000e+00
omega 3.511539e-082.485181e-06 1.412991e-029.887263e-01
alpha1 1.034993e-011.069228e-02 9.679819e+000.000000e+00
beta1 9.220933e-017.358088e-03 1.253170e+020.000000e+00
gamma1 5.290305e-017.168919e-02 7.379502e+001.589839e-13
delta 1.294503e+001.373357e-01 9.425829e+000.000000e+00
skew 8.448373e-011.631212e-02 5.179198e+010.000000e+00
shape 5.772091e+004.773574e-01 1.209176e+010.000000e+00

## 10.7 Summary

We can put all the results into a dataframe

``````r=Results[[length(Results)]]

df=data.frame(c(likelihood(r),r@fit\$timer,coef(r)))
rownames(df)[1:2]=c("log.likelihood","time")
rownames(df)[3:dim(df)[1]]=names(coef(r))
names(df)=names(Results)[length(Results)]

for(i in 1:length(Results)){
r=Results[[i]]
x=c(likelihood(r),r@fit\$timer,coef(r))
x=data.frame(c(likelihood(r),r@fit\$timer,coef(r)))
rownames(x)[1:2]=c("log.likelihood","time")
rownames(x)[3:dim(x)[1]]=names(coef(r))

df[[names(Results)[i]]]=NA
df[rownames(x),names(Results)[i]]=x

}
df=round(df,3)
df=df[,names(Results)]

df``````
A data.frame: 12 × 8
ARCH1GARCH11tGARCH11atGARCH11APGARCH11tAPGARCH11atAPGARCH11Mean.atAPGARCH11
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
log.likelihood15519.61116455.26316575.75916603.89216518.89816690.21216729.52016670.514
time 0.121 0.065 0.150 0.395 0.137 0.807 1.459 0.985
mu NA NA NA NA NA NA NA 0.000
ar1 NA NA NA NA NA NA NA 0.415
ma1 NA NA NA NA NA NA NA -0.487
omega 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
alpha1 0.441 0.124 0.124 0.123 0.061 0.103 0.104 0.103
beta1 NA 0.856 0.871 0.871 0.873 0.897 0.898 0.922
gamma1 NA NA NA NA 0.446 1.000 1.000 0.529
delta NA NA NA NA 2.555 1.007 0.978 1.295
skew NA NA NA 0.877 NA NA 0.849 0.845
shape NA NA 6.062 6.481 NA 6.513 7.038 5.772

### 10.7.1 Issues

#### 10.7.1.1 Sample size

It is important to consider that the more parameter a model has, the more likely it runs into estimation problems, and the more data we need. For example, as a rule of thumb, for a GARCH(1,1) we need at least 500 observations, while for a Student-t GARCH that minimum number can be around 3,000 observations. In general, the more parameters are estimated, the more data is needed to be confident in the estimation. We see an example of this in Chapter 11.

## 10.8 Diagnostics

### 10.8.1 Likelihood ratio test

The likelihood ratio (LR) test compares two models based on the ratio of their likelihoods,

``````# Perform the LR test
LR_statistic = 2*(likelihood(Results\$GARCH11)-likelihood(Results\$ARCH1))
p_value  = 1 - pchisq(LR_statistic, df = 1)

cat(" Likelihood of ARCH: ", round(likelihood(Results\$ARCH1),1), "\n",
"Likelihood of GARCH:", round(likelihood(Results\$GARCH11),1), "\n",
"2 * (Lu - Lr): ", round(LR_statistic,2), "\n",
"p-value:", p_value)``````
`````` Likelihood of ARCH:  15519.6
Likelihood of GARCH: 16455.3
2 * (Lu - Lr):  1871.3
p-value: 0``````

We find a p-value of 0, meaning that we have enough evidence to reject the null hypothesis that the two models are the same. The GARCH model has a significantly larger likelihood than the ARCH model.

### 10.8.2 Residual analysis

#### 10.8.2.1 GARCH

``````residuals=y/Results\$GARCH11@fit\$sigma
par(mar=c(2,4,2,0))
plot(residuals,type='l',bty='l',las=1)
myACF(residuals,main="ACF of SP-500 returns")
myACF(residuals^2,main="ACF of squared SP-500 returns")
x=qqPlot(residuals, distribution = "norm", envelope = FALSE,main="Normal QQ of SP-500 returns")``````

## 10.9 Monte Carlo analysis of GARCH models

### 10.9.1 Using `ugarchsim()`

#### 10.9.1.1 Gaussian

``````s=ugarchsim(Results\$GARCH11, n.sim = 5000)
par(mar=c(2,4,2,0))
plot(y,type='l',bty='l',las=1,main="SP500")
plot(s@simulation\$seriesSim,type='l',bty='l',las=1,main="Simulated GARCH SP500")``````

``````s=ugarchsim(Results\$tGARCH11, n.sim = 5000)
par(mar=c(2,3,2,0))
plot(y,type='l',bty='l',las=1,main="SP500")
plot(s@simulation\$seriesSim,type='l',bty='l',las=1,main="Simulated tGARCH SP500")``````

#### 10.9.1.2 Simulation and estimation

``````N=5
for(n in 1:N){
s=ugarchsim(Results\$GARCH11, n.sim = 5000)
Results\$GARCH11 = ugarchfit(spec = spec.2, data = s@simulation\$seriesSim)
print(Results\$GARCH11@fit\$matcoef)
}``````
``````           Estimate   Std. Error   t value    Pr(>|t|)
omega  2.710645e-06 1.040379e-06  2.605439 0.009175673
alpha1 1.228770e-01 1.304027e-02  9.422883 0.000000000
beta1  8.541387e-01 1.494085e-02 57.168002 0.000000000
Estimate   Std. Error  t value    Pr(>|t|)
omega  3.130311e-06 1.050228e-06  2.98060 0.002876842
alpha1 1.204353e-01 1.032526e-02 11.66414 0.000000000
beta1  8.515699e-01 1.219078e-02 69.85361 0.000000000
Estimate   Std. Error   t value   Pr(>|t|)
omega  3.285156e-06 1.212014e-06  2.710494 0.00671831
alpha1 9.802658e-02 9.062237e-03 10.817040 0.00000000
beta1  8.697615e-01 1.148797e-02 75.710660 0.00000000
Estimate   Std. Error  t value    Pr(>|t|)
omega  3.018835e-06 1.026916e-06  2.93971 0.003285201
alpha1 1.037525e-01 9.367526e-03 11.07576 0.000000000
beta1  8.666474e-01 1.163049e-02 74.51512 0.000000000
Estimate   Std. Error   t value    Pr(>|t|)
omega  2.962311e-06 1.039697e-06  2.849205 0.004382858
alpha1 1.112850e-01 1.116922e-02  9.963539 0.000000000
beta1  8.603203e-01 1.352397e-02 63.614493 0.000000000``````

### 10.9.2 Manual simulation

#### 10.9.2.1 GARCH(1,1)

``````S=1e3
omega=1e-6
alpha=0.15
beta=0.8

N=10000
set.seed(888)
Res=NULL
for(n in 1:N){
eps=rnorm(S)
y=rep(NA,S)
sigma2 = omega/(1-alpha-beta)

y[1]=eps[1]*sqrt(sigma2)
for(i in 2:S){
sigma2= omega+alpha*y[i-1]^2 + beta * sigma2
y[i]=eps[i]*sqrt(sigma2)
}

}
par(mar=c(2,3.5,0,0))
plot(y,type='l',bty='l',las=1)``````

#### 10.9.2.2 tGARCH

``````df=6
omega=0.000001
alpha=0.15
beta=0.75
N=1
S=1e5
set.seed(888)
for(n in 1:N){
eps=rt(S,df=df)
y=rep(NA,S)
sigma2 = omega/(1-alpha-beta)

y[1]=eps[1]*sqrt(sigma2)
for(i in 2:S){
sigma2 = omega + alpha*y[i-1]^2 + beta * sigma2
y[i]=eps[i]*sqrt(sigma2)
}
}
par(mar=c(2,3.5,0,0))
plot(y,type='l',bty='l',las=1)  ``````