Chapter 5. Implementing Risk Forecasts (in Python/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 Python
Last updated June 2018

import numpy as np
p = np.loadtxt('stocks.csv',delimiter=',',skiprows=1)
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)
y1 = y1[len(y1)-4100:]
y2 = y2[len(y2)-4100:]
y = np.stack([y1,y2], axis = 1)
T = len(y1)
value = 1000                                          # portfolio value
p = 0.01                                              # probability
		
Listing 5.1/5.2: Download stock prices in Julia
Last updated June 2018

using CSV;
p = CSV.read("stocks.csv",nullable=false);
## 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 Python
Last updated June 2018

ys = np.sort(y1)           # sort returns
op = int(T*p)              # p percent smallest
VaR1 = -ys[op - 1] * 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 Python
Last updated June 2018

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.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 Python
Last updated June 2018

ES1 = -np.mean(ys[:op]) * 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 Python
Last updated June 2018

sigma = np.std(y1, ddof=1)                # estimate volatility
VaR3 = -sigma * stats.norm.ppf(p) * 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 Python
Last updated June 2018

## portfolio volatility
sigma = np.sqrt(np.mat(np.transpose(w))*np.mat(np.cov(y,rowvar=False))*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.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 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.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 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.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 Python
Last updated June 2018

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)[0] / p * 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)[1] / p * value
		

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.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 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.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: Two-asset EWMA VaR in Python
Last updated June 2018

## 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.transpose(w)*s*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.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: GARCH VaR in Python
Last updated June 2018

from arch import arch_model
am = arch_model(y1, mean = 'Zero', vol='Garch', p=1, o=0, q=1, dist='Normal')
res = am.fit(update_freq=5)
omega = res.params[0]
alpha = res.params[1]
beta = res.params[2]
## 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)
## Note: arch_model's GARCH optimization has issues with convergence
		
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