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

Copyright 2011 - 2022 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 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.1/5.2: Download stock prices in Python Last updated July 2020

import numpy as np
from scipy import stats
p = p[:,[0,1]] # consider two stocks
## convert prices to returns, and adjust length
y1 = np.diff(np.log(p[:,0]), n=1, axis=0)
y2 = np.diff(np.log(p[:,1]), n=1, axis=0)
y = np.stack([y1,y2], axis = 1)
T = len(y1)
value = 1000   # portfolio value
p = 0.01   # probability


##### 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.3/5.4: Univariate HS in Python Last updated July 2020

from math import ceil
ys = np.sort(y1) # sort returns
op = ceil(T*p)   # p percent smallest
VaR1 = -ys[op - 1] * value
print(VaR1)


##### 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.5/5.6: Multivariate HS in Python Last updated July 2020

w = [0.3, 0.7]       # vector of portfolio weights
yp = np.squeeze(np.matmul(y, w)) # portfolio returns
yps = np.sort(yp)
VaR2= -yps[op - 1] * value
print(VaR2)


##### Listing 5.7/5.8: Univariate ES in MATLAB Last updated 2011

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

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

ES1 = -np.mean(ys[:op]) * value
print(ES1)


##### Listing 5.9/5.10: Normal VaR in MATLAB Last updated 2011

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

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

sigma = np.std(y1, ddof=1) # estimate volatility
VaR3 = -sigma * stats.norm.ppf(p) * value
print(VaR3)


##### 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.11/5.12: Portfolio normal VaR in Python Last updated June 2022

## portfolio volatility
sigma = np.sqrt(np.mat(w)*np.mat(np.cov(y,rowvar=False))*np.transpose(np.mat(w)))[0,0]
## Note: [0,0] is to pull the first element of the matrix out as a float
VaR4 = -sigma * stats.norm.ppf(p) * value
print(VaR4)


##### 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.13/5.14: Student-t VaR in Python Last updated June 2018

scy1 = y1 * 100    # scale the returns
res = stats.t.fit(scy1)
sigma = res[2]/100 # rescale volatility
nu = res[0]
VaR5 = -sigma*stats.t.ppf(p,nu)*value
print(VaR5)


##### 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.15/5.16: Normal ES in Python Last updated June 2018

sigma = np.std(y1, ddof=1)
ES2 = sigma * stats.norm.pdf(stats.norm.ppf(p)) / p * value
print(ES2)


##### Listing 5.17/5.18: Direct integration ES in MATLAB Last updated 2011

VaR = -norminv(p);

##### Listing 5.17/5.18: Direct integration ES in Python Last updated June 2018

VaR = -stats.norm.ppf(p)
integrand = lambda q: q * stats.norm.pdf(q)
ES = -sigma * quad(integrand, -np.inf, -VaR)[0] / p * value
print(ES)


##### 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.19/5.20: MA normal VaR in Python Last updated June 2018

WE = 20
for t in range(T-5,T+1):
t1 = t-WE
window = y1[t1:t] # estimation window
sigma = np.std(window, ddof=1)
VaR6 = -sigma*stats.norm.ppf(p)*value
print (VaR6)


##### 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.21/5.22: EWMA VaR in Python Last updated June 2018

lmbda = 0.94
s11 = np.var(y1[0:30], ddof = 1) # initial variance
for t in range(1, T):
s11 = lmbda*s11 + (1-lmbda)*y1[t-1]**2
VaR7 = -np.sqrt(s11)*stats.norm.ppf(p)*value
print(VaR7)


##### 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.23/5.24: Two-asset EWMA VaR in Python Last updated June 2022

## s is the initial covariance
s = np.cov(y, rowvar = False)
for t in range(1,T):
s = lmbda*s+(1-lmbda)*np.transpose(np.asmatrix(y[t-1,:]))*np.asmatrix(y[t-1,:])
sigma = np.sqrt(np.mat(w)*s*np.transpose(np.mat(w)))[0,0]
## Note: [0,0] is to pull the first element of the matrix out as a float
VaR8 = -sigma * stats.norm.ppf(p) * value
print(VaR8)


##### 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

##### Listing 5.25/5.26: GARCH VaR in Python Last updated July 2020

from arch import arch_model
am = arch_model(y1, mean = 'Zero', vol='Garch', p=1, o=0, q=1, dist='Normal', rescale = False)
res = am.fit(update_freq=5, disp = "off")
omega = res.params.loc['omega']
alpha = res.params.loc['alpha[1]']
beta = res.params.loc['beta[1]']
## computing sigma2 for t+1
sigma2 = omega + alpha*y1[T-1]**2 + beta * res.conditional_volatility[-1]**2
VaR9 = -np.sqrt(sigma2) * stats.norm.ppf(p) * value
print(VaR9)