Below we address implementation of volatility models as discussed in Chapter four and five of the book and Chapter 10 in these notes.

Along the way, we will make some functions and put them into `functions.r`.

## 13.1 Libraries

``suppressPackageStartupMessages(library(rugarch))``

## 13.2 Data

``````load('data/data.RData')
sp500=data\$sp500
y=data\$sp500\$y
Return=data\$Return
Ticker=data\$Ticker
names(data)
sd(y)``````
1. 'Return'
2. 'Price'
4. 'sp500'
5. 'sp500tr'
6. 'Ticker'
0.0121610363140039

## 13.3 Specification

Specify the VaR probability, portfolio value, estimation window size and portfolio weight:

``````p = 0.05
value = 1000
WE = 1000
w = rep(1/6,6)``````

In what follows below, we will use several methods for forecast risk. To facilitate the simulation study at end we put each univariate VaR into a function.

## 13.4 Historical Simulation — HS

### 13.4.1 Single asset

``````par(mar=c(2,3,0,0),bty='l')
plot(y,type='l',las=1,ylab="",xlab="")`````` Let’s sort the values using `sort()`:

``````ys = sort(y)
par(mfrow=c(1,2))
plot(ys, type = "l",main="all sorted values",las=1,bty='l') Historical Simulation is a method that assumes that history repeats itself. It uses the empirical distribution of the returns to find the p-th quantile.

``````plot(head(ys,60),main="The lowest 60 sorted values",las=1,bty='l')
segments(WE*0.01,-1,WE*0.01,ys[WE*0.01])
segments(WE*0.01,ys[WE*0.01],0,ys[WE*0.01])
segments(WE*0.05,-1,WE*0.05,ys[WE*0.05])
segments(WE*0.05,ys[WE*0.05],0,ys[WE*0.05])`````` ``````Risk.HS=function(y,WE,p,value,Print=FALSE){
ys=sort(tail(y,WE))
VaR= round(-ys[p*WE] * value,3)
ES= round(- mean(ys[1:(p*WE)]) * value ,3)
if(Print) cat("The ",p*100,"% HS VaR is \$",VaR," while the ES is \$",ES,".\n",sep="")
return(list(risk=list(VaR=VaR,ES=ES),method="HS",p=p,value=value,WE=WE,sigma=sd(y)))
}
res=list()
res\$HS=Risk.HS(y=y,WE=WE,p=p,value=value,Print=TRUE)
res\$HS``````
``The 5% HS VaR is \$21.683 while the ES is \$36.515.``
\$risk
\$VaR
21.683
\$ES
36.515
\$method
'HS'
\$p
0.05
\$value
1000
\$WE
1000
\$sigma
0.0121610363140039

### 13.4.2 Portfolio

``````yw = as.matrix(tail(Return[,Ticker],WE))
dim(yw)
yw=yw %*% w
dim(yw)``````
1. 1000
2. 6
1. 1000
2. 1

Then do the same calculation as for one asset

## 13.5 Univariate parametric VaR with known values

Since we will do the calculations repeatedly, we set them up as functions, and put them into `functions.r`.

``````NormalVaR=function(p,sd=1) return(-qnorm(p,sd=sd))

NormalES=function(p,sd=1) return(sd*dnorm(qnorm(p))/(p))

tVaR=function(p,df,sd=1,Standardized=TRUE){
Scale=1
if(Standardized) Scale=sqrt(df/(df-2))
return(-sd*qt(p=p,df=df)/Scale)
}
tES=function(p,df,sd=1,Standardized=TRUE) {
Scale=1
if(Standardized) Scale=sqrt(df/(df-2))
return((sd*dt(qt(p=p,df=df),df=df)*((df+(qt(p=p,df=df))^2)/(df-1))/p)/Scale)
}``````

## 13.6 EWMA

``````Risk.EWMA=function(y,WE=30,lambda=0.94,p,value,Print=FALSE){

sigma2=var(y)
for (i in 2:length(y)){
sigma2 = lambda*sigma2+(1-lambda)*y[i-1] ^2
}
sigma=sqrt(sigma2)
VaR = round( NormalVaR(p=p,sd=sigma)* value,3)
ES = round( NormalES(p=p,sd=sigma)*value,3)
if(Print){
cat("The forecast volatility is",round(sigma,4),"\n")
cat("The EWMA ",p*100,"% VaR is \$",VaR,", while the ES is \$",ES,"\n",sep="")
}
if(Print) cat("The ",p*100,"% HS VaR is \$",VaR," while the ES is \$",ES,"\n",sep="")
return(list(risk=list(VaR=VaR,ES=ES),method="EWMA",sigma=sigma,p=p,value=value,WE=WE,lambda=lambda,sigma=sigma))
}
res\$EWMA=Risk.EWMA(y=y,p=p,value=value)
res\$EWMA``````
\$risk
\$VaR
22.063
\$ES
27.668
\$method
'EWMA'
\$sigma
0.0134135930446279
\$p
0.05
\$value
1000
\$WE
30
\$lambda
0.94
\$sigma
0.0134135930446279

### 13.6.1 Gaussian

``````sigma=0.01
VaR = round( NormalVaR(p=p,sd=sigma)* value,1)
ES = round( NormalES(p=p,sd=sigma)*value,1)
cat("The ",p*100,"% VaR is \$",VaR," while the ES is \$",ES,"\n",sep="")``````
``The 5% VaR is \$16.4 while the ES is \$20.6``

### 13.6.2 Student-t

We show two cases for the ES, one with a numerical integration of the t density while the other has the exact solution. These should be the same.

Here an important issue is if the student-t distribution is standardised or not, i.e. with variance set to 1. This becomes an issues since many GARCH libraries standardise.

#### 13.6.2.1 Non-standardised

``````nu=4
VaR = round( tVaR(p=p,sd=sigma,df=nu,Standardized=FALSE)* value,1)
ES = round( tES(p=p,sd=sigma,df=nu,Standardized=FALSE)* value,5)

print(ES)

f = function(x) qt(p=x,df=nu)
ES.integrate =  -value*sigma*integrate(f, 0, p)\$value/p

print(ES.integrate)

cat("The ",p*100,"% VaR is \$",VaR,", while the ES is \$",ES,"\n",sep="")``````
`` 32.0287``
`` 32.0287``
``The 5% VaR is \$21.3, while the ES is \$32.0287``

#### 13.6.2.2 Standardised

Note in the direct integration how we divide by the standard error of the Student-t
`sqrt(nu/(nu-2))`.

``````nu=4
VaR = round( tVaR(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,1)
ES = round( tES(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,5)

print(ES)

f = function(x) qt(p=x,df=nu)
ES.integrate =  -value*sigma*integrate(f, 0, p)\$value/p / sqrt(nu/(nu-2))
print(ES.integrate)
cat("The ",p*100,"% VaR is \$",VaR,", while the ES is \$",ES,"\n",sep="")``````
`` 22.64771``
`` 22.64771``
``The 5% VaR is \$15.1, while the ES is \$22.64771``

## 13.7 Univariate parametric VaR with estimated parameters

### 13.7.1 Gaussian GARCH

We set the GARCH up as a function as well:

``````Risk.GARCH=function(y,p,value,WE,Print=FALSE){
ys=tail(y,WE)
y=y-mean(y)
scale=sd(y)
y=y/scale
spec = ugarchspec(
variance.model = list(
garchOrder= c(1,1)),
mean.model= list(
armaOrder = c(0,0),
include.mean=FALSE)
)
res = ugarchfit(spec = spec, data = y, solver = "hybrid")
omega = coef(res)['omega']
alpha = coef(res)['alpha1']
beta = coef(res)['beta1']
res@fit\$coef=res@fit\$coef*scale^2
sigma = sqrt(omega + alpha*tail(y,1)^2 + beta*tail(res@fit\$var,1))* scale
names(sigma)=NULL
VaR = round( NormalVaR(p=p,sd=sigma)* value,3)
ES = round( NormalES(p=p,sd=sigma)*value,3)
if(Print){
cat("The forecast volatility is",round(sigma,4),"\n")
cat("The normal GARCH ",p*100,"% VaR is \$",VaR,", while the ES is \$",ES,"\n",sep="")
}
return(list(risk=list(VaR=VaR,ES=ES),method="GARCH",par=coef(res),p=p,value=value,WE=WE,sigma=sigma))
}``````

We can then run it.

``````res\$GARCH=Risk.GARCH(y=y,p=p,WE=WE,value=value,Print=TRUE)
res\$GARCH``````
``````The forecast volatility is 0.0115
The normal GARCH 5% VaR is \$18.944, while the ES is \$23.757``````
\$risk
\$VaR
18.944
\$ES
23.757
\$method
'GARCH'
\$par
omega
2.45918650421497e-06
alpha1
0.123865389470534
beta1
0.855799127959958
\$p
0.05
\$value
1000
\$WE
1000
\$sigma
0.0115173263130952

### 13.7.2 tGARCH

``````Risk.tGARCH=function(y,p,WE,value,Print=FALSE){
ys=tail(y,WE)
y=y-mean(y)
scale=sd(y)
y=y/scale
spec = ugarchspec(
variance.model = list(
garchOrder= c(1,1)),
mean.model= list(
armaOrder = c(0,0),
include.mean=FALSE),
distribution.model = "std"
)
res = ugarchfit(spec = spec, data = y, solver = "hybrid")
omega = coef(res)
alpha = coef(res)
beta = coef(res)
nu = coef(res)
sigma = sqrt(omega + alpha*tail(y,1)^2 + beta*tail(res@fit\$var,1))* scale
res@fit\$coef=res@fit\$coef*scale^2
names(sigma)=NULL
VaR = round( tVaR(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,3)
ES = round( tES(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,3)
if(Print) {
cat("The forecast volatility is",round(sigma,4),"while nu is",nu,"\n")
cat("The t-GARCH ",p*100,"% VaR is \$",VaR,", while the ES is \$",ES,"\n",sep="")
}
return(list(risk=list(VaR=VaR,ES=ES),method="tGARCH",par=coef(res),p=p,value=value,WE=WE,sigma=sigma))
}``````

We can then run it.

``````res\$tGARCH=Risk.tGARCH(y=y,p=p,WE=WE,value=value,Print=TRUE)
res\$tGARCH``````
``````The forecast volatility is 0.012 while nu is 6.058236
The t-GARCH 5% VaR is \$19.018, while the ES is \$26.497``````
\$risk
\$VaR
shape: 19.018
\$ES
shape: 26.497
\$method
'tGARCH'
\$par
omega
1.51909326692045e-06
alpha1
0.124319142338583
beta1
0.870313461282666
shape
6.05823629039183
\$p
0.05
\$value
1000
\$WE
1000
\$sigma
0.0119787135614369

### 13.7.3 Summary

``````for(i in res){
cat(i\$method,i\$risk\$VaR,i\$sigma,"\n")
}``````
``````HS 21.683 0.01216104
EWMA 22.063 0.01341359
GARCH 18.944 0.01151733
tGARCH 19.018 0.01197871 ``````

## 13.8 How this can be improved

Ideally we would put all the risk functions, like `Risk.tGARCH()` into a single function that would look like
`Risk = function(y,p,value,Model="GARCH11",Print=FALSE,Plot=FALSE)`