R and Julia Chapter 5. Implementing Risk Forecasts

# Chapter 5. Implementing Risk Forecasts

### R and Julia

Copyright 2011 - 2023 Jon Danielsson. This code is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. The GNU General Public License is available at: www.gnu.org/licenses.

##### Listing 5.1/5.2
p = read.csv('stocks.csv')
y = apply(log(p),2,diff)
portfolio_value = 1000
p = 0.01

##### Listing 5.1/5.2
using CSV, DataFrames;
y1 = diff(log.(p[:,1]));
y2 = diff(log.(p[:,2]));
y = hcat(y1,y2);
T = size(y,1);
value = 1000; # portfolio value
p = 0.01;     # probability


##### Univariate HS in R
y1 = y[,1]                     # select one asset
ys = sort(y1)                  # sort returns
op = ceiling(length(y1)*p)     # p percent smallest, rounded up
VaR1 = -ys[op]*portfolio_value
print(VaR1)

##### Univariate HS in Julia
ys = sort(y1)            # sort returns
op = ceil(Int, T*p)      # p percent smallest, rounding up
VaR1 = -ys[op] * value
println("Univariate HS VaR ", Int(p*100), "%: ", round(VaR1, digits = 3), " USD")


##### Multivariate HS in R
w = matrix(c(0.3,0.2,0.5))   # vector of portfolio weights
yp = y %*% w                 # obtain portfolio returns
yps = sort(yp)
VaR2 = -yps[op]*portfolio_value
print(VaR2)

##### Multivariate HS in Julia
w = [0.3; 0.7]          # vector of portfolio weights
yp = y * w              # portfolio returns
yps = sort(yp)
VaR2 = -yps[op] * value
println("Multivariate HS VaR ", Int(p*100), "%: ", round(VaR2, digits = 3), " USD")


##### Univariate ES in R
ES1 = -mean(ys[1:op])*portfolio_value
print(ES1)

##### Univariate ES in Julia
using Statistics;
ES1 = -mean(ys[1:op]) * value
println("ES: ", round(ES1, digits = 3), " USD")


##### Normal VaR in R
sigma = sd(y1) # estimate volatility
VaR3 = -sigma * qnorm(p) * portfolio_value
print(VaR3)

##### Normal VaR in Julia
using Distributions;
sigma = std(y1);      # estimate volatility
VaR3 = -sigma * quantile(Normal(0,1),p) * value
println("Normal VaR", Int(p*100), "%: ", round(VaR3, digits = 3), " USD")


##### Portfolio normal VaR in R
sigma = sqrt(t(w) %*% cov(y) %*% w) # portfolio volatility
VaR4 = -sigma * qnorm(p)*portfolio_value
print(VaR4)

##### Portfolio normal VaR in Julia
sigma = sqrt(w'*cov(y)*w)   # portfolio volatility
VaR4 = -sigma * quantile(Normal(0,1), p) * value
println("Portfolio normal VaR", Int(p*100), "%: ", round(VaR4, digits = 3), " USD")


##### Student-t VaR in R
library(QRM)
scy1 = (y1)*100                     # scale the returns
res = fit.st(scy1)
sigma1 = unname(res$par.ests['sigma']/100) # rescale the volatility nu = unname(res$par.ests['nu'])
VaR5 = - sigma1 * qt(df=nu,p=p) *  portfolio_value
print(VaR5)

##### Student-t VaR in Julia



##### Normal ES in R
sigma = sd(y1)
ES2 = sigma*dnorm(qnorm(p))/p * portfolio_value
print(ES2)

##### Normal ES in Julia
sigma = std(y1)
ES2 = sigma * pdf(Normal(0,1), (quantile(Normal(0,1), p))) / p * value
println("Normal ES: ", round(ES2, digits = 3), " USD")


##### Direct integration ES in R
VaR = -qnorm(p)
integrand = function(q){q*dnorm(q)}
ES = -sigma*integrate(integrand,-Inf,-VaR)$value/p*portfolio_value print(ES)  ##### Listing 5.17/5.18 ##### Direct integration ES in Julia using QuadGK; VaR = -quantile(Normal(0,1), p) integrand(x) = x * pdf(Normal(0,1), x) ES3 = -sigma * quadgk(integrand, -Inf, -VaR) / p * value println("Normal integrated ES: ", round(ES3, digits = 3), " USD")  ##### Listing 5.19/5.20 ##### MA normal VaR in R WE=20 for (t in seq(length(y1)-5,length(y1))){ t1=t-WE+1 window= y1[t1:t] # estimation window sigma=sd(window) VaR6 = -sigma * qnorm(p) * portfolio_value print(VaR6) }  ##### Listing 5.19/5.20 ##### MA normal VaR in Julia WE = 20 for t in T-5:T t1 = t-WE window = y1[t1+1:t] # estimation window sigma = std(window) VaR6 = -sigma*quantile(Normal(0,1),p)*value println("MA Normal VaR", Int(p*100), "% using observations ", t1, " to ", t, ": ", round(VaR6, digits = 3), " USD") end  ##### Listing 5.21/5.22 ##### EWMA VaR in R lambda = 0.94 s11 = var(y1) # initial variance, using unconditional for (t in 2:length(y1)){ s11 = lambda * s11 + (1-lambda) * y1[t-1]^2 } VaR7 = -qnorm(p) * sqrt(s11) * portfolio_value print(VaR7)  ##### Listing 5.21/5.22 ##### EWMA VaR in Julia lambda = 0.94 s11 = var(y1) # initial variance for t in 2:T s11 = lambda * s11 + (1-lambda) * y1[t-1]^2 end VaR7 = -sqrt(s11) * quantile(Normal(0,1), p) * value println("EWMA VaR ", Int(p*100), "%: ", round(VaR7, digits = 3), " USD")  ##### Listing 5.23/5.24 ##### Three-asset EWMA VaR in R s = cov(y) # initial covariance for (t in 2:dim(y)){ s = lambda*s + (1-lambda)*y[t-1,] %*% t(y[t-1,]) } sigma = sqrt(t(w) %*% s %*% w) # portfolio vol VaR8 = -sigma * qnorm(p) * portfolio_value print(VaR8)  ##### Listing 5.23/5.24 ##### Two-asset EWMA VaR in Julia s = cov(y) # initial covariance for t in 2:T s = lambda * s + (1-lambda) * y[t-1,:] * (y[t-1,:])' end sigma = sqrt(w'*s*w) # portfolio vol VaR8 = -sigma * quantile(Normal(0,1), p) * value println("Two-asset EWMA VaR ", Int(p*100), "%: ", round(VaR8, digits = 3), " USD")  ##### Listing 5.25/5.26 ##### Univariate GARCH in R library(rugarch) spec = ugarchspec(variance.model = list( garchOrder = c(1, 1)), mean.model = list( armaOrder = c(0,0),include.mean = FALSE)) res = ugarchfit(spec = spec, data = y1, solver = "hybrid") omega = res@fit$coef['omega']
alpha = res@fit$coef['alpha1'] beta = res@fit$coef['beta1']
sigma2 = omega + alpha * tail(y1,1)^2 + beta * tail(res@fit\$var,1)
VaR9 = -sqrt(sigma2) * qnorm(p) * portfolio_value
names(VaR9)="VaR"
print(VaR9)

##### GARCH VaR in Julia
using ARCHModels;
garch1_1 = fit(GARCH{1,1}, y1; meanspec = NoIntercept);
garch_VaR_in = VaRs(garch1_1, :0.01)
cond_vol = predict(garch1_1, :volatility)     # 1-day-ahead conditional volatility
garch_VaR_out = -cond_vol * quantile(garch1_1.dist, p) * value
println("GARCH VaR ", Int(p*100), "%: ", round(garch_VaR_out, digits = 3), " USD")


##### Financial Risk Forecasting
Market risk forecasting with R, Julia, Python and Matlab. Code, lecture slides, implementation notes, seminar assignments and questions.