library(reshape2)
source("common/functions.r",chdir=TRUE)
library(rugarch)
library(car)
13 Univariate volatility
Univariate volatility modelling is designed to provide an in-sample description of the stochastic processes governing volatility. More importantly, 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 Financial Risk Forecasting.
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 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 14, we will see what can go wrong with volatility models and how to fix it.
13.1 Libraries
13.2 Data
The data we use is returns on the S&P500 index, and for convenience, we put them into the variable y
.
=ProcessRawData()
data=data$sp500$y
y=y-mean(y) y
13.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.
However, since the mean is so small, we can usually ignore it if the sample size is large. Chapter 14 shows how that can fail.
13.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.
= vector(length=length(y))
EWMA = 0.94
lambda 1] = var(y)
EWMA[for (i in 2:length(y)){
=
EWMA[i] * EWMA[i-1]+
lambda 1-lambda) * y[i-1] ^2
(
}=sqrt(EWMA)
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
)
13.5 rugarch
13.5.1 Setup
To estimate a univariate GARCH model with rugarch
, we need to follow two steps:
- 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 the model for the mean; - Fit the specified model to the data.
13.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:
13.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”, “iGARCH”, and “csGARCH”garchOrder
: Specify the ARCH(q) and GARCH(p) orders
Other components include options to do variance.targeting
or external regressors
.
13.5.2.2 mean.model
This argument takes a list
with the specifications of the mean model and is 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.
13.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
.
13.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, as shown below.
13.5.4 Optimiser failure — Hybrid
ugarchfit()
performs maximum likelihood to fit the specified model to our data. This optimisation problem requires a solver
, which is the numerical method used to maximise the likelihood function. See discussion in Section 15.3. 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.
13.6 The models
We will put the results into a list called results()
so we can automatically process them afterwards.
=list() Results
13.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
.
.1 = ugarchspec(
specvariance.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:
.1 spec
*---------------------------------*
* 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:
$ARCH1 = ugarchfit(spec = spec.1, data = y) Results
Warning in .sgarchfit(spec = spec, data = data, out.sample = out.sample, :
ugarchfit-->warning: solver failer to converge.
That fails; we need the hybrid solver.
$ARCH1 = ugarchfit(spec = spec.1, data = y,solver = "hybrid") Results
We can check the class of this new object:
class(Results$ARCH1)
[1] "uGARCHfit"
attr(,"package")
[1] "rugarch"
13.6.1.0.1 Printing and plotting
To get all the results, 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:
@model
;@fit
.
To access a slot, you need to use the syntax: object@slot
13.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 everything that is included in this slot, you can use the names()
function:
names(Results$ARCH1@model)
[1] "modelinc" "modeldesc" "modeldata" "pars" "start.pars"
[6] "fixed.pars" "maxOrder" "pos.matrix" "fmodel" "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).
$ARCH1@model$pars[7:10,] Results
Level Fixed Include Estimate LB UB
omega 8.889786e-05 0 1 1 2.220446e-16 0.1478908
alpha1 4.407729e-01 0 1 1 0.000000e+00 1.0000000
beta 0.000000e+00 0 0 0 NA NA
gamma 0.000000e+00 0 0 0 NA NA
We have a matrix with all the possible parameters to fit a GARCH model. The parameters included in our GARCH specification (omega
, alpha1
, and beta1
) have a value of 1. The matrix also includes the lower and upper bounds for these parameters, which have nothing to do with our particular data but with what is expected from the model.
We can also see the model description:
$ARCH1@model$modeldesc Results
$distribution
[1] "norm"
$distno
[1] 1
$vmodel
[1] "sGARCH"
Inside Results$ARCH1@model@modeldata
, we will find the data vector we used when fitting the model. The @model
slot contains everything that R needs to know to estimate the model.
13.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" "cvar" "var" "sigma"
[5] "condH" "z" "LLH" "log.likelihoods"
[9] "residuals" "coef" "robust.cvar" "A"
[13] "B" "scores" "se.coef" "tval"
[17] "matcoef" "robust.se.coef" "robust.tval" "robust.matcoef"
[21] "fitted.values" "convergence" "message" "kappa"
[25] "persistence" "timer" "ipars" "solver"
Let’s see the coefficients, along with their standard errors, in matcoef
. There are 2 ways to do that:
$ARCH1@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 8.889786e-05 2.506444e-06 35.46773 0
alpha1 4.407729e-01 3.359505e-02 13.12017 0
coef(Results$ARCH1)
omega alpha1
8.889786e-05 4.407729e-01
We can also see the log-likelihood. Again two ways
$ARCH1@fit$LLH Results
[1] 15519.61
likelihood(Results$ARCH1)
[1] 15519.61
This makes it easy to extract the estimated conditional volatility:
plot(Results$ARCH1@fit$sigma, type = "l")
13.6.2 Gaussian GARCH(1,1)
.2 = ugarchspec(
specvariance.model = list(garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE)
)$GARCH11 = ugarchfit(spec = spec.2, data = y)
Results$GARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 2.458870e-06 7.647997e-07 3.215051 0.001304212
alpha1 1.238025e-01 9.901141e-03 12.503865 0.000000000
beta1 8.558381e-01 1.073389e-02 79.732345 0.000000000
13.6.3 tGARCH(1,1)
We can see that the coefficients now include shape
, which is the estimated degrees of freedom:
.3 = ugarchspec(
specvariance.model = list(garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "std"
)$tGARCH11 = ugarchfit(spec = spec.3, data = y)
Results$tGARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 1.507668e-06 1.949716e-06 0.7732756 4.393593e-01
alpha1 1.240843e-01 3.430002e-02 3.6176161 2.973289e-04
beta1 8.706493e-01 3.094681e-02 28.1337309 0.000000e+00
shape 6.054042e+00 8.182193e-01 7.3990465 1.372236e-13
13.6.4 Asymmetric tGARCH(1,1)
We can have a skew Student-t conditional distribution.
.4 = ugarchspec(
specvariance.model = list(garchOrder = c(1,1)),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "sstd"
)$atGARCH11 = ugarchfit(spec = spec.4, data = y)
Results$atGARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 1.495222e-06 1.851522e-06 0.807564 4.193416e-01
alpha1 1.228609e-01 3.224600e-02 3.810114 1.389026e-04
beta1 8.710221e-01 2.938361e-02 29.643124 0.000000e+00
skew 8.770726e-01 1.669690e-02 52.529069 0.000000e+00
shape 6.467815e+00 9.572251e-01 6.756838 1.410361e-11
13.6.5 Gaussian apARCH model
The parameter list will now include gamma1
and delta
, which are part of the apARCH model specification:
.5 = ugarchspec(
specvariance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
)
$APGARCH11 = ugarchfit(spec = spec.5, data = y)
Results$APGARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 1.453484e-07 4.735252e-07 0.3069496 0.75888173
alpha1 6.084138e-02 3.590656e-02 1.6944363 0.09018241
beta1 8.728381e-01 4.150876e-02 21.0278064 0.00000000
gamma1 4.460231e-01 1.840133e-01 2.4238627 0.01535641
delta 2.554894e+00 5.051892e-02 50.5730029 0.00000000
13.6.6 t-apARCH model
.6 = ugarchspec(
specvariance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "std"
)$tAPGARCH11 = ugarchfit(spec = spec.6, data = y)
Results$tAPGARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 0.0002526539 0.0001222254 2.067115 0.03872335
alpha1 0.1027769088 0.0092727659 11.083738 0.00000000
beta1 0.8969539656 0.0097020705 92.449748 0.00000000
gamma1 0.9999999900 0.0003536942 2827.301056 0.00000000
delta 1.0065472477 0.0911995498 11.036757 0.00000000
shape 6.5128558586 0.5744228117 11.338087 0.00000000
13.6.7 Asymmetric t-apARCH model
.7 = ugarchspec(
specvariance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(0,0), include.mean = FALSE),
distribution.model = "sstd"
)$atAPGARCH11 = ugarchfit(spec = spec.7, data = y)
Results$atAPGARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
omega 0.0002902321 0.0001322079 2.19527 0.02814421
alpha1 0.1035077895 0.0095115201 10.88236 0.00000000
beta1 0.8975209454 0.0101043456 88.82524 0.00000000
gamma1 0.9999999900 0.0003833587 2608.52302 0.00000000
delta 0.9780618429 0.0842650532 11.60697 0.00000000
skew 0.8487325626 0.0161120763 52.67680 0.00000000
shape 7.0382213920 0.6802275855 10.34686 0.00000000
13.6.8 Asymmetric t-apARCH model with ARMA(1,1) mean
.8 = ugarchspec(
specvariance.model = list(model = "apARCH"),
mean.model = list(armaOrder = c(1,1), include.mean = TRUE),
distribution.model = "sstd"
)$Mean.atAPGARCH11 = ugarchfit(spec = spec.8, data = y,solver = "hybrid")
Results$Mean.atAPGARCH11@fit$matcoef Results
Estimate Std. Error t value Pr(>|t|)
mu -8.298863e-17 9.366143e-05 -8.860492e-13 1.000000e+00
ar1 4.146313e-01 5.708975e-02 7.262799e+00 3.792522e-13
ma1 -4.871048e-01 5.568636e-02 -8.747291e+00 0.000000e+00
omega 3.511556e-08 2.485179e-06 1.412999e-02 9.887263e-01
alpha1 1.034993e-01 1.069227e-02 9.679820e+00 0.000000e+00
beta1 9.220934e-01 7.358085e-03 1.253170e+02 0.000000e+00
gamma1 5.290307e-01 7.168920e-02 7.379503e+00 1.589839e-13
delta 1.294503e+00 1.373358e-01 9.425821e+00 0.000000e+00
skew 8.448373e-01 1.631212e-02 5.179199e+01 0.000000e+00
shape 5.772092e+00 4.773572e-01 1.209177e+01 0.000000e+00
13.7 Summary
We can put all the results into a dataframe
=Results[[length(Results)]]
r
=data.frame(c(likelihood(r),r@fit$timer,coef(r)))
dfrownames(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)){
=Results[[i]]
r=c(likelihood(r),r@fit$timer,coef(r))
x=data.frame(c(likelihood(r),r@fit$timer,coef(r)))
xrownames(x)[1:2]=c("log.likelihood","time")
rownames(x)[3:dim(x)[1]]=names(coef(r))
names(Results)[i]]]=NA
df[[rownames(x),names(Results)[i]]=x
df[
}=round(df,3)
df=df[,names(Results)]
df
=df df.summary
13.7.1 Issues
13.7.1.1 Sample size
It is important to consider that the more parameters a model has, the more likely it is to run 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 14.
13.8 Diagnostics
13.8.1 Likelihood ratio test
The likelihood ratio (LR) test compares two models based on the ratio of their likelihoods,
# Perform the LR test
= 2*(likelihood(Results$GARCH11)-likelihood(Results$ARCH1))
LR_statistic = 1 - pchisq(LR_statistic, df = 1)
p_value
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.
13.8.2 Residual analysis
13.8.2.1 GARCH
=y/Results$GARCH11@fit$sigma
residualspar(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")
=qqPlot(residuals, distribution = "norm", envelope = FALSE,main="Normal QQ of SP-500 returns") x
13.9 Monte Carlo analysis of GARCH models
13.9.1 Using ugarchsim()
13.9.1.1 Gaussian
=ugarchsim(Results$GARCH11, n.sim = 5000)
spar(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")
=ugarchsim(Results$tGARCH11, n.sim = 5000)
spar(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")
13.9.1.2 Simulation and estimation
=5
Nfor(n in 1:N){
=ugarchsim(Results$GARCH11, n.sim = 5000)
s$GARCH11 = ugarchfit(spec = spec.2, data = s@simulation$seriesSim)
Resultsprint(Results$GARCH11@fit$matcoef)
}
Estimate Std. Error t value Pr(>|t|)
omega 2.183211e-06 0.0000011277 1.935986 5.286943e-02
alpha1 9.983926e-02 0.0147237274 6.780841 1.194778e-11
beta1 8.812737e-01 0.0165986574 53.093072 0.000000e+00
Estimate Std. Error t value Pr(>|t|)
omega 2.595323e-06 1.051599e-06 2.467978 1.358787e-02
alpha1 9.688202e-02 1.171525e-02 8.269737 2.220446e-16
beta1 8.767996e-01 1.430066e-02 61.311843 0.000000e+00
Estimate Std. Error t value Pr(>|t|)
omega 2.733132e-06 9.760227e-07 2.800275 0.005105917
alpha1 1.007865e-01 1.056384e-02 9.540711 0.000000000
beta1 8.699549e-01 1.283992e-02 67.753930 0.000000000
Estimate Std. Error t value Pr(>|t|)
omega 2.249089e-06 9.081018e-07 2.476693 0.0132606
alpha1 9.579799e-02 1.134245e-02 8.445970 0.0000000
beta1 8.769719e-01 1.365328e-02 64.231579 0.0000000
Estimate Std. Error t value Pr(>|t|)
omega 2.906452e-06 1.178785e-06 2.465633 0.01367715
alpha1 9.754828e-02 9.270425e-03 10.522525 0.00000000
beta1 8.631767e-01 1.233798e-02 69.960923 0.00000000
13.9.2 Manual simulation
13.9.2.1 GARCH(1,1)
=1e3
S=1e-6
omega=0.15
alpha=0.8
beta
=10000
Nset.seed(888)
=NULL
Resfor(n in 1:N){
=rnorm(S)
eps=rep(NA,S)
y= omega/(1-alpha-beta)
sigma2
1]=eps[1]*sqrt(sigma2)
y[for(i in 2:S){
= omega+alpha*y[i-1]^2 + beta * sigma2
sigma2=eps[i]*sqrt(sigma2)
y[i]
}
}par(mar=c(2,3.5,0,0))
plot(y,type='l',bty='l',las=1)
13.9.2.2 tGARCH
=6
df=0.000001
omega=0.15
alpha=0.75
beta=1
N=1e5
Sset.seed(888)
for(n in 1:N){
=rt(S,df=df)
eps=rep(NA,S)
y= omega/(1-alpha-beta)
sigma2
1]=eps[1]*sqrt(sigma2)
y[for(i in 2:S){
= omega + alpha*y[i-1]^2 + beta * sigma2
sigma2 =eps[i]*sqrt(sigma2)
y[i]
}
}par(mar=c(2,3.5,0,0))
plot(y,type='l',bty='l',las=1)
13.10 Exercise
In this chapter, we estimated univariate models using eight different specifications. Conduct pairwise likelihood ratio tests for all eight model specifications. Comment on the results.
# Answer
# df.summary refers to the data frame generated in the summary section
library(combinat)
=combn(1:8, 2)
cb
=c()
rstfor (i in 1:ncol(cb)){
=cb[,i]
idx=2*abs(df.summary[1,idx[1]]-df.summary[1,idx[2]])
lr.tmp=1-pchisq(lr.tmp, df=1)
pval.tmp=c(lr.tmp, pval.tmp)
rst.tmp=rbind(rst, rst.tmp)
rst
}colnames(rst)=c("LR Stat", "p-value")
rownames(rst)=1:nrow(rst)
rstwhich(rst[,2]>0.05)
# We reject the null hypothesis that the two models are the same for all pairs.