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

Copyright 2011 - 2020 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 July 2020

## Calculate returns. Note that first column is dates
y = apply(log(p),2,diff)
## Specify portfolio value and VaR probability
portfolio_value = 1000
p = 0.01

##### Listing 5.1/5.2: Download stock prices in MATLAB Last updated July 2020

p1 = stocks(:,1); % consider first two stocks
p2 = stocks(:,2);
y1=diff(log(p1)); % convert prices to returns
y2=diff(log(p2));
y=[y1 y2];
T=length(y1);
value = 1000;     % portfolio value
p = 0.01;  % probability


##### Listing 5.3/5.4: Univariate HS in R Last updated July 2020

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)

##### Listing 5.3/5.4: Univariate HS VaR in MATLAB Last updated July 2020

ys = sort(y1);  % sort returns
op = ceil(T*p); % p percent smallest, rounded up to meet VaR probability requirement
VaR1 = -ys(op)*value


##### Listing 5.5/5.6: Multivariate HS in R Last updated July 2020

w = matrix(c(0.3,0.2,0.5)) # vector of portfolio weights
## The number of columns of the left matrix must be the same as the number of rows of the right matrix
yp = y %*% w         # obtain portfolio returns
yps = sort(yp)
VaR2 = -yps[op]*portfolio_value
print(VaR2)

##### Listing 5.5/5.6: Multivariate HS VaR in MATLAB Last updated 2011

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 MATLAB Last updated 2011

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 MATLAB Last updated 2011

sigma = std(y1); % estimate volatility
VaR3 = -sigma * norminv(p) * value


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

sigma = sqrt(t(w) %*% cov(y) %*% w)[1] # portfolio volatility
## Note: the trailing [1] 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 MATLAB Last updated 2011

sigma = sqrt(w' * cov(y) * w); % portfolio volatility
VaR4 = - sigma * norminv(p) *  value


##### Listing 5.13/5.14: Student-t VaR in R Last updated July 2020

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'])
## Note: We are removing the names of the fit.st output using unname()
VaR5 = - sigma1 * qt(df=nu,p=p) *  portfolio_value
print(VaR5)

##### Listing 5.13/5.14: Student-t VaR in MATLAB Last updated 2011

scy1=y1*100;         % scale the returns
res=mle(scy1,'distribution','tlocationscale');
sigma1 = res(2)/100; % rescale the volatility
nu = res(3);
VaR5 = - sigma1 * tinv(p,nu) * value


##### 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 MATLAB Last updated June 2018

sigma = std(y1);
ES2=sigma*normpdf(norminv(p))/p * value


##### Listing 5.17/5.18: Direct integration ES in R Last updated July 2020

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 MATLAB Last updated 2011  VaR = -norminv(p); ES = -sigma*quad(@(q) q.*normpdf(q),-6,-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 MATLAB Last updated June 2018  WE=20; for t=T-5:T t1=t-WE+1; window=y1(t1:t); % estimation window sigma=std(window); VaR6 = -sigma * norminv(p) * value end  ##### Listing 5.21/5.22: EWMA VaR in R Last updated July 2020  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 MATLAB Last updated 2011  lambda = 0.94; s11 = var(y1(1:30)); % initial variance for t = 2:T s11 = lambda * s11 + (1-lambda) * y1(t-1)^2; end VaR7 = -norminv(p) * sqrt(s11) * 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)[1]){ s = lambda*s + (1-lambda)*y[t-1,] %*% t(y[t-1,]) } sigma = sqrt(t(w) %*% s %*% w)[1] # portfolio vol ## Note: trailing[1] 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 MATLAB Last updated 2011  s = cov(y); % initial covariance for t = 2:T s = lambda * s + (1-lambda) * y(t-1,:)' * y(t-1,:); end sigma = sqrt(w' * s * w); % portfolio vol VaR8 = - sigma * norminv(p) * value  ##### Listing 5.25/5.26: Univariate GARCH in R Last updated July 2020  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['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)

##### Listing 5.25/5.26: GARCH in MATLAB Last updated August 2016

[parameters,ll,ht]=tarch(y1,1,0,1);
omega = parameters(1)
alpha = parameters(2)
beta = parameters(3)
sigma2 = omega + alpha*y1(end)^2 + beta*ht(end) % calc sigma2 for t+1
VaR9 = -sqrt(sigma2) * norminv(p) * value