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

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'
  3. 'UnAdjustedPrice'
  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')
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])

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="")
[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)).

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="")
[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:

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[1]=res@fit$coef[1]*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)[1]
    alpha = coef(res)[2]
    beta = coef(res)[3]
    nu = coef(res)[4]
    sigma = sqrt(omega + alpha*tail(y,1)^2 + beta*tail(res@fit$var,1))* scale
    res@fit$coef[1]=res@fit$coef[1]*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)