17  Implementing risk forecasts

Below, we address the implementation risk forecast models as discussed in Chapters four and fGARCHive of Financial Risk Forecasting.

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

17.1 Libraries

library(reshape2)
source("common/functions.r",chdir=TRUE)
library(rugarch)

17.2 Data

data=ProcessRawData()
sp500=data$sp500

17.3 Specification

When backtesting, we have a number of parameters we have to keep track of, such as probability, portfolio value, sample size, testing window size and estimation winter size. In order to minimise the chance of mistakes and make the coding as fast as possible, we use the techniques discussed in Section Section 12.6 and keep track of them in a list that we call par for parameters.

par         = list()
par$p       = 0.05
par$value   = 1000
par$WE      = 2000
par$T       = length(sp500$y)
par$WT      = par$T-par$WE
par$lambda  = 0.94

cat("The backtesting parameters are:",
  "\n\tp =     ",par$p,
  "\n\tvalue = ",par$value,
  "\n\tpar$WE =",par$WE,
  "\n\tT =     ",par$T,
  "\n\tWT =    ",par$WT,
  "\n\tlambda =",par$lambda,
  "\n"
)

The backtesting parameters are: 
    p =      0.05 
    value =  1000 
    par$WE = 2000 
    T =      5034 
    WT =     3034 
    lambda = 0.94 

print(par)

$p
[1] 0.05

$value
[1] 1000

$WE
[1] 2000

$T
[1] 5034

$WT
[1] 3034

$lambda
[1] 0.94

17.4 Univariate parametric ES and VaR

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

17.4.1 Normal

It is quite straightforward to get ES and VaR for the normal distribution:

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

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

17.4.2 Student-t

An important issue here is whether the student-t distribution is standardised, i.e., with variance set to 1. This becomes an issue since many GARCH libraries standardise.

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))
 ES=(sd*dt(qt(p=p,df=df),df=df)*
    ((df+(qt(p=p,df=df))^2)/
    (df-1))/p)/Scale
 return(ES)
}

17.5 Functions

While we could program the functionality for backtesting without functions, because we may want to repeat the calculations multiple times and with different parameters, it is much better to use functions.

We label their functions as Risk.* where * refers to any of the methods we are using, like HS or EWMA. We also put these functions into functions.rso we can use them later.

Each function takes three arguments: the returns y, the parameter list par and Print, which determines whether the function prints out the parameters and results it achieves. This can be very useful for debugging.

17.6 Historical simulation — Risk.HS

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.

par(mar=c(2,3,0,0),bty='l')
plot(sp500$y,type='l',las=1,ylab="",xlab="")

Let’s sort the values using sort():

ys = sort(sp500$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'
)

plot(head(ys,60),main="The lowest 60 sorted values",las=1,bty='l')
segments(par$WE*0.01,-1,par$WE*0.01,ys[par$WE*0.01])
segments(par$WE*0.01,ys[par$WE*0.01],0,ys[par$WE*0.01])
segments(par$WE*0.05,-1,par$WE*0.05,ys[par$WE*0.05])
segments(par$WE*0.05,ys[par$WE*0.05],0,ys[par$WE*0.05])

17.6.1 Risk.HS

The function of implementing HS, Risk.HS, is quite straightforward. It takes a vector of returns and the parameter list, along with an optional argument, Print, and returns a list with both the ES and VaR, along with the parameter list.

Risk.HS=function(y,par,Print=FALSE){
    ys=sort(tail(y,par$WE))
    VaR= -ys[par$p*par$WE] * par$value
    ES= - mean(ys[1:(par$p*par$WE)]) * par$value
    if(Print) cat("The HS ",par$p*100,"% VaR is $",round(VaR,1),", while the ES is $",round(ES,1),"\n",sep="") 

    return(list(VaR=VaR,ES=ES,method="HS",par=par))
}

17.7 EWMA — Risk.EWMA

Risk.EWMA=function(y,par,Print=FALSE){
    y=tail(y,par$WE)
    sigma2=var(y)
    for (i in 2:length(y)){
        sigma2 = par$lambda*sigma2+(1-par$lambda)*y[i-1] ^2
    }
    sigma=sqrt(sigma2)
    VaR = NormalVaR(p=par$p,sd=sigma)* par$value
    ES =  NormalES(p=par$p,sd=sigma)*par$value
    if(Print){
        cat("The EWMA forecast volatility is",round(sigma,4),"\n")
        cat("The EWMA ",par$p*100,"% VaR is $",round(VaR,1),", while the ES is $",round(ES,1),"\n",sep="") 
    }
    return(list(VaR=VaR,ES=ES,method="EWMA",sigma=sigma))
}

17.8 Gaussian GARCH

We set the GARCH up as a function as well:

Risk.GARCH=function(y,par,Print=FALSE){
    ys=tail(y,par$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 =  NormalVaR(p=par$p,sd=sigma)* par$value
    ES =  NormalES(p=par$p,sd=sigma)*par$value

    if(Print){
        cat("The GARCH forecast volatility is",round(sigma,4),"\n")
        cat("The GARCH ",par$p*100,"% VaR is $",round(VaR,1),", while the ES is $",round(ES,1),"\n",sep="") 
    }
    return(list(VaR=VaR,ES=ES,method="GARCH",parameters=coef(res),sigma=sigma))
}

17.9 t-GARCH

Risk.tGARCH=function(y,par,Print=FALSE){
    ys=tail(y,par$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 =  tVaR(p=par$p,sd=sigma,df=nu,Standardized=TRUE)* par$value
    ES =  tES(p=par$p,sd=sigma,df=nu,Standardized=TRUE)* par$value
    if(Print){
        cat("The tGARCH forecast volatility is",round(sigma,4),"\n")
        cat("The tGARCH ",par$p*100,"% VaR is $",round(VaR,1),", while the ES is $",round(ES,1),"\n",sep="") 
    }
    return(list(VaR=VaR,ES=ES,method="tGARCH",parameters=coef(res),sigma=sigma))
}

17.10 Summary

res=list()
res$HS=Risk.HS(y=sp500$y,par=par,Print=TRUE)

The HS 1% VaR is $35.8, while the ES is $51.6

res$EWMA=Risk.EWMA(y=sp500$y,par=par,Print=TRUE)

The EWMA forecast volatility is 0.0134 
The EWMA 1% VaR is $31.2, while the ES is $35.8

res$GARCH=Risk.GARCH(y=sp500$y,par=par,Print=TRUE)

The GARCH forecast volatility is 0.0115 
The GARCH 1% VaR is $26.8, while the ES is $30.7

res$tGARCH=Risk.tGARCH(y=sp500$y,par=par,Print=TRUE)

The tGARCH forecast volatility is 0.012 
The tGARCH 1% VaR is $30.7, while the ES is $39.4

for(l in res){
    cat(l$method,l$VaR,l$ES,"\n")
}

HS 35.75759 51.60686 
EWMA 31.20468 35.7501 
GARCH 26.793 30.696 
tGARCH 30.712 39.353 

for(l in res){
    if(l$method==res[[names(res)[1]]]$method){
        df=t(data.frame(c(l$method,l$VaR,l$ES)))
    }else{
        df=rbind(df,c(l$method,l$VaR,l$ES))
    }
}
df=as.data.frame(df)
names(df)=c("method","VaR","ES")
rownames(df)=NULL
df$VaR=as.numeric(df$VaR)
df$ES=as.numeric(df$ES)

df

method	VaR	ES
HS	35.75759	51.60686
EWMA	31.20468	35.75010
GARCH	26.79300	30.69600
tGARCH	30.71200	39.35300

17.11 Exercise

1. Create a generic function for obtaining VaR and ES from GARCH models, where the particular model specification, normal, t, skew t, lags, and apARCH, are arguments to a single function.
2. Combine all risk functions discussed above into a single function, Risk = function(y,par,Model,Print), that computes ES and VaR based on a set of inputs. This function should support HS, EWMA, GARCH, and tGARCH models and be robust enough to handle improper inputs.