suppressPackageStartupMessages(library(rugarch))
13 Implementing risk forecasts
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
13.2 Data
load('data/data.RData')
=data$sp500
sp500=data$sp500$y
y=data$Return
Return=data$Ticker
Tickernames(data)
sd(y)
- 'Return'
- 'Price'
- 'UnAdjustedPrice'
- 'sp500'
- 'sp500tr'
- 'Ticker'
13.3 Specification
Specify the VaR probability, portfolio value, estimation window size and portfolio weight:
= 0.05
p = 1000
value = 1000
WE = rep(1/6,6) w
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()
:
= sort(y)
ys par(mfrow=c(1,2))
plot(ys, type = "l",main="all sorted values",las=1,bty='l')
plot(head(ys,50),main="the lowest 50 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])
=function(y,WE,p,value,Print=FALSE){
Risk.HS=sort(tail(y,WE))
ys= round(-ys[p*WE] * value,3)
VaR= round(- mean(ys[1:(p*WE)]) * value ,3)
ESif(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)))
}=list()
res$HS=Risk.HS(y=y,WE=WE,p=p,value=value,Print=TRUE)
res$HS res
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
= as.matrix(tail(Return[,Ticker],WE))
yw dim(yw)
=yw %*% w
ywdim(yw)
- 1000
- 6
- 1000
- 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
.
=function(p,sd=1) return(-qnorm(p,sd=sd))
NormalVaR
=function(p,sd=1) return(sd*dnorm(qnorm(p))/(p))
NormalES
=function(p,df,sd=1,Standardized=TRUE){
tVaR=1
Scaleif(Standardized) Scale=sqrt(df/(df-2))
return(-sd*qt(p=p,df=df)/Scale)
}=function(p,df,sd=1,Standardized=TRUE) {
tES=1
Scaleif(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
=function(y,WE=30,lambda=0.94,p,value,Print=FALSE){
Risk.EWMA
=var(y)
sigma2for (i in 2:length(y)){
= lambda*sigma2+(1-lambda)*y[i-1] ^2
sigma2
}=sqrt(sigma2)
sigma= round( NormalVaR(p=p,sd=sigma)* value,3)
VaR = round( NormalES(p=p,sd=sigma)*value,3)
ES 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))
}$EWMA=Risk.EWMA(y=y,p=p,value=value)
res$EWMA res
- $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
=0.01
sigma= round( NormalVaR(p=p,sd=sigma)* value,1)
VaR = round( NormalES(p=p,sd=sigma)*value,1)
ES 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
=4
nu= round( tVaR(p=p,sd=sigma,df=nu,Standardized=FALSE)* value,1)
VaR = round( tES(p=p,sd=sigma,df=nu,Standardized=FALSE)* value,5)
ES
print(ES)
= function(x) qt(p=x,df=nu)
f = -value*sigma*integrate(f, 0, p)$value/p
ES.integrate
print(ES.integrate)
cat("The ",p*100,"% VaR is $",VaR,", while the ES is $",ES,"\n",sep="")
[1] 32.0287
[1] 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))
.
=4
nu= round( tVaR(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,1)
VaR = round( tES(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,5)
ES
print(ES)
= function(x) qt(p=x,df=nu)
f = -value*sigma*integrate(f, 0, p)$value/p / sqrt(nu/(nu-2))
ES.integrate print(ES.integrate)
cat("The ",p*100,"% VaR is $",VaR,", while the ES is $",ES,"\n",sep="")
[1] 22.64771
[1] 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:
=function(y,p,value,WE,Print=FALSE){
Risk.GARCH=tail(y,WE)
ys=y-mean(y)
y=sd(y)
scale=y/scale
y= ugarchspec(
spec variance.model = list(
garchOrder= c(1,1)),
mean.model= list(
armaOrder = c(0,0),
include.mean=FALSE)
)= ugarchfit(spec = spec, data = y, solver = "hybrid")
res = coef(res)['omega']
omega = coef(res)['alpha1']
alpha = coef(res)['beta1']
beta @fit$coef[1]=res@fit$coef[1]*scale^2
res= sqrt(omega + alpha*tail(y,1)^2 + beta*tail(res@fit$var,1))* scale
sigma names(sigma)=NULL
= round( NormalVaR(p=p,sd=sigma)* value,3)
VaR = round( NormalES(p=p,sd=sigma)*value,3)
ES 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.
$GARCH=Risk.GARCH(y=y,p=p,WE=WE,value=value,Print=TRUE)
res$GARCH res
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
=function(y,p,WE,value,Print=FALSE){
Risk.tGARCH=tail(y,WE)
ys=y-mean(y)
y=sd(y)
scale=y/scale
y= ugarchspec(
spec variance.model = list(
garchOrder= c(1,1)),
mean.model= list(
armaOrder = c(0,0),
include.mean=FALSE),
distribution.model = "std"
)= ugarchfit(spec = spec, data = y, solver = "hybrid")
res = coef(res)[1]
omega = coef(res)[2]
alpha = coef(res)[3]
beta = coef(res)[4]
nu = sqrt(omega + alpha*tail(y,1)^2 + beta*tail(res@fit$var,1))* scale
sigma @fit$coef[1]=res@fit$coef[1]*scale^2
resnames(sigma)=NULL
= round( tVaR(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,3)
VaR = round( tES(p=p,sd=sigma,df=nu,Standardized=TRUE)* value,3)
ES 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.
$tGARCH=Risk.tGARCH(y=y,p=p,WE=WE,value=value,Print=TRUE)
res$tGARCH res
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)