Python and Julia Chapter 5. Implementing Risk Forecasts

# Chapter 5. Implementing Risk Forecasts

### Python 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
import numpy as np
from scipy import stats
p = p[:,[0,1]]      # consider two stocks
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.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 Python
from math import ceil
ys = np.sort(y1) # sort returns
op = ceil(T*p)    # p percent smallest
VaR1 = -ys[op - 1] * 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 Python
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)

##### 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 Python
ES1 = -np.mean(ys[:op]) * 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 Python
sigma = np.std(y1, ddof=1) # estimate volatility
VaR3 = -sigma * stats.norm.ppf(p) * 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 Python
sigma = np.sqrt(np.mat(w)*np.mat(np.cov(y,rowvar=False))*np.transpose(np.mat(w)))[0,0]
VaR4 = -sigma * stats.norm.ppf(p) * 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 Python
scy1 = y1 * 100         # scale the returns
res = stats.t.fit(scy1)
sigma = res/100      # rescale volatility
nu = res
VaR5 = -sigma*stats.t.ppf(p,nu)*value
print(VaR5)

##### Student-t VaR in Julia



##### Normal ES in Python
sigma = np.std(y1, ddof=1)
ES2 = sigma * stats.norm.pdf(stats.norm.ppf(p)) / p * 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 Python
from scipy.integrate import quad
VaR = -stats.norm.ppf(p)
integrand = lambda q: q * stats.norm.pdf(q)
ES = -sigma * quad(integrand, -np.inf, -VaR) / p * value
print(ES)

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


##### MA normal VaR in Python
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)

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


##### EWMA VaR in Python
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)

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


##### Two-asset EWMA VaR in Python
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]
VaR8 = -sigma * stats.norm.ppf(p) * value
print(VaR8)

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


##### GARCH VaR in Python
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']
beta = res.params.loc['beta']
sigma2 = omega + alpha*y1[T-1]**2 + beta * res.conditional_volatility[-1]**2
VaR9 = -np.sqrt(sigma2) * stats.norm.ppf(p) * value
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.