14  Implementing risk forecasts

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

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

14.1 Libraries

suppressPackageStartupMessages(library(rugarch))

14.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

14.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.

14.4 Historical Simulation — HS

14.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

14.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

14.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)
}

14.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

14.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

14.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.

14.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

14.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

14.7 Univariate parametric VaR with estimated parameters

14.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

14.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

14.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 

14.8 Exercise

Combine all risk functions discussed above into a single function Risk = function(y,p,value,Model="GARCH",Print=FALSE) that computes VaR based on a set of inputs. This function should support HS, EWMA, GARCH, tGARCH models, and should be robust to handle improper inputs.

# Answer
Risk = function(y, p=0.05, value=NULL, WE=NULL, lambda=NULL, Model="HS", Print=FALSE){
  if (Model=="HS"){
    if(is.null(value)+is.null(WE)>0){
      print("Error. Please check the inputs.")
    }else{
    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)))}
  }
  else if (Model=="EWMA"){
    if(is.null(value)+is.null(WE)+is.null(lambda)>0){
      print("Error. Please check the inputs.")
    }else{
    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="")
    }
    return(list(risk=list(VaR=VaR,ES=ES),method="EWMA",sigma=sigma,p=p,value=value,WE=WE,
                lambda=lambda,sigma=sigma))}
  }
  else if (Model=="GARCH"){
    if(is.null(value)+is.null(WE)>0){
      print("Error. Please check the inputs.")
    }else{
    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))}
  }
  else if(Model=="tGARCH"){
    if(is.null(value)+is.null(WE)>0){
      print("Error. Please check the inputs.")
    }else{
    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))}
  }else{
    print("Error. Please check the inputs.")
  }
}
Risk(y=y,p=0.05,value=1000,WE=30,lambda=0.94, Model="EWMA", Print=TRUE)

The forecast volatility is 0.0134 
The EWMA 5% VaR is $22.063, while the ES is $27.668

$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