 # Chapter 5. Implementing Risk Forecasts (in R/Julia)

Copyright 2011 - 2019 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: https://www.gnu.org/licenses/.

##### Listing 5.1/5.2: Download stock prices in R Last updated August 2019

y=apply(log(p),2,diff)     # calculate returns. note first column is dates
portfolio_value = 1000
p = 0.01                   # probability

##### Listing 5.1/5.2: Download stock prices in Julia Last updated June 2018

using CSV;
## convert prices of first two stocks to returns, and adjust length
y1 = diff(log.(p[:,1]));
y2 = diff(log.(p[:,2]));
y1 = y1[length(y1)-4100+1:length(y1)];
y2 = y2[length(y2)-4100+1:length(y2)];
y = hcat(y1,y2);
T = size(y,1)
value = 1000;                              # portfolio value
p = 0.01;                                  # probability


##### Listing 5.3/5.4: Univariate HS in R Last updated August 2016

ys = sort(y1)                  # sort returns
op = length(y1)*p              # p percent smallest
VaR1 = -ys[op]*portfolio_value
print(VaR1)

##### Listing 5.3/5.4: Univariate HS in Julia Last updated June 2018

ys = sort(y1)            # sort returns
op = convert(Int64, T*p) # p percent smallest
VaR1 = -ys[op] * value


##### Listing 5.5/5.6: Multivariate HS in R Last updated August 2019

w = matrix(c(0.3,0.2,0.5)) # vector of portfolio weights

##### Listing 5.5/5.6: Multivariate HS in Julia Last updated June 2018

w = [0.3; 0.7]          # vector of portfolio weights
yp = y * w              # portfolio returns
yps = sort(yp)
VaR2 = -yps[op] * value


##### Listing 5.7/5.8: Univariate ES in R Last updated August 2019

ES1 = -mean(ys[1:op])*portfolio_value
print(ES1)

##### Listing 5.7/5.8: Univariate ES in Julia Last updated June 2018

ES1 = -mean(ys[1:op]) * value


##### Listing 5.9/5.10: Normal VaR in R Last updated August 2019

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

##### Listing 5.9/5.10: Normal VaR in Julia Last updated June 2018

sigma = std(y1);                                # estimate volatility
using Distributions;
VaR3 = -sigma * quantile(Normal(0,1),p) * value


##### Listing 5.11/5.12: Portfolio normal VaR in R Last updated August 2019

sigma = sqrt(t(w) %*% cov(y) %*% w)   # portfolio volatility
## Note: the trailing  is to convert a single element matrix to float
VaR4 = -sigma * qnorm(p)*portfolio_value
print(VaR4)

##### Listing 5.11/5.12: Portfolio normal VaR in Julia Last updated June 2018

sigma = sqrt(w'*cov(y)*w)                        # portfolio volatility
VaR4 = -sigma * quantile(Normal(0,1), p) * value


##### Listing 5.13/5.14: Student-t VaR in R Last updated August August 2019

library(QRM)
scy1=(y1)*100                                      # scale the returns
res=fit.st(scy1)
sigma1=res$par.ests/100 # rescale the volatility nu=res$par.ests
VaR5 = - sigma1 * qt(df=nu,p=p) *  portfolio_value
print(VaR5)

##### Listing 5.13/5.14: Student-t VaR in Julia Last updated June 2018

## using Distributions;
## res = fit_mle(TDist, y1)
## nu = res.ν (this is the Greek letter nu, not Latin v)
## sigma = sqrt(nu/(nu-2))
## VaR5 = -sigma * quantile(TDist(nu), p) * value
## Julia does not have a function for fitting Student-t data yet
## Currently: there exists Distributions.jl with fit_mle
## usage: Distributions.fit_mle(Dist_name, data[, weights])


##### Listing 5.15/5.16: Normal ES in R Last updated June August 2019

sigma = sd(y1)
ES2 = sigma*dnorm(qnorm(p))/p * portfolio_value
print(ES2)

##### Listing 5.15/5.16: Normal ES in Julia Last updated June 2018

sigma = std(y1)
ES2 = sigma * pdf(Normal(0,1), (quantile(Normal(0,1), p))) / p * value


##### Listing 5.17/5.18: Direct integration ES in R Last updated August 2019

VaR = -qnorm(p)
integrand = function(q){q*dnorm(q)}
ES = -sigma*integrate(integrand,-Inf,-VaR)$portfolio_value/p*portfolio_value print(ES)  ##### Listing 5.17/5.18: Direct integration ES in Julia Last updated June 2018  using QuadGK; VaR = -quantile(Normal(0,1), p) integrand(x) = x*pdf(Normal(0,1), x) ES = -sigma * quadgk(integrand, -Inf, -VaR) / p * value  ##### Listing 5.19/5.20: MA normal VaR in R Last updated June August 2019  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 Last updated June 2018  WE = 20 for t in range(T-5, 6) t1 = t-WE window = y1[t1+1:t] # estimation window sigma = std(window) VaR6 = -sigma*quantile(Normal(0,1),p)*value println(VaR6) end  ##### Listing 5.21/5.22: EWMA VaR in R Last updated August 2019  lambda = 0.94; s11 = var(y1[1:30]); # initial variance 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 Last updated June 2018  lambda = 0.94 s11 = var(y1[1:30]) # initial variance for t in range(2, T-1) s11 = lambda * s11 + (1-lambda) * y1[t-1]^2 end VaR7 = -sqrt(s11) * quantile(Normal(0,1), p) * value  ##### Listing 5.23/5.24: Three-asset EWMA VaR in R Last updated August 2019  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 ## Note:  is to convert single element matrix to float VaR8 = -sigma * qnorm(p) * portfolio_value print(VaR8)  ##### Listing 5.23/5.24: Two-asset EWMA VaR in Julia Last updated June 2018  s = cov(y) # initial covariance for t in range(2, T-1) 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  ##### Listing 5.25/5.26: Univariate GARCH in R Last updated August 2019  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) omega = res@fit$coef
alpha = res@fit$coef beta = res@fit$coef
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)

##### Listing 5.25/5.26: GARCH VaR in Julia Last updated June 2018

## We use the FRFGarch mini-package again (refer to Chapter 2 above)
using FRFGarch;
res = GARCHfit(y1)
sigma_fc = res.seForecast
## seForecast contains the next-day conditional volatility forecast
VaR9 = - sigma_fc * quantile(Normal(0,1), p) * value
## GARCH estimation will be slightly different from other languages
## this is due to GARCHfit choosing initial conditional vol = sample vol