Chapter 8. Backtesting and Stress Testing (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 8.1/8.2: Load data in Python
Last updated June 2018

import numpy as np
price = np.loadtxt('index.csv',delimiter=',',skiprows=1)
y = np.diff(np.log(price), n=1, axis=0)                  # get returns
		
Listing 8.1/8.2: Load data in Julia
Last updated June 2018

using CSV;
price = CSV.read("index.csv", nullable = false);
y = diff(log.(price[:,1]));                      # get returns
		

Listing 8.3/8.4: Set backtest up in Python
Last updated June 2018

import numpy as np
T = len(y)                            # number of obs for y
WE = 1000                             # estimation window length
p = 0.01                              # probability
l1 = int(WE * p)                      # HS observation
value = 1                             # portfolio value
VaR = np.full([T,4], np.nan)          # matrix for forecasts
## EWMA setup
lmbda = 0.94
s11 = np.var(y[1:30])
for t in range(1,WE):
    s11=lmbda*s11+(1-lmbda)*y[t-1]**2
		
Listing 8.3/8.4: Set backtest up in Julia
Last updated June 2018

T = length(y)                          # number of obs for return y
WE = 1000                              # estimation window length
p = 0.01                               # probability
l1 = convert(Int, WE*p)                # HS observation
value = 1                              # portfolio value
VaR = fill!(Array{Float64}(T,4), NaN)  # matrix for forecasts
## EWMA setup
lambda = 0.94
s11 = var(y[1:30])
for t in range(2,WE-1)
    s11=lambda*s11+(1-lambda)*y[t-1]^2
end
		

Listing 8.5/8.6: Running backtest in Python
Last updated June 2018

import numpy as np
from scipy import stats
from arch import arch_model
for t in range(WE, T):
    t1 = t - WE                                                # start of data window
    t2 = t - 1                                                 # end of data window
    window = y[t1:t2+1]                                        # data for estimation
    s11 = lmbda * s11 + (1-lmbda) * y[t-1]**2
    VaR[t,0]=-stats.norm.ppf(p)*np.sqrt(s11)*value             # EWMA
    VaR[t,1]=-np.std(window,ddof=1)*stats.norm.ppf(p)*value    # MA
    ys = np.sort(window)
    VaR[t,2] = -ys[l1 - 1] * value                             # HS
    am = arch_model(window, mean = 'Zero',vol='Garch',
                    p=1, o=0, q=1, dist='Normal')
    res = am.fit(update_freq=0, disp='off',show_warning=False)
    par = [res.params[0], res.params[1], res.params[2]]
    s4 = par[0] + par[1] * window[WE - 1]**2 + par[
        2] * res.conditional_volatility[-1]**2
    VaR[t,3] = -np.sqrt(s4) * stats.norm.ppf(p) * value        # GARCH(1,1)
## GARCH optimization in Python has some convergence issues
		
Listing 8.5/8.6: Running backtest in Julia
Last updated June 2018

using Distributions;
using FRFGarch;
for t in range(WE+1, T-WE)
    t1 = t - WE                                         # start of data window
    t2 = t - 1                                          # end of data window
    window = y[t1:t2]                                   # data for estimation
                                                        ## EWMA
    s11 = lambda * s11 + (1-lambda) * y[t-1]^2
    VaR[t,1]=-quantile(Normal(0,1),p)*sqrt(s11)*value   # EWMA
    VaR[t,2]=-std(window)*quantile(Normal(0,1),p)*value # MA
    ys = sort(window)
    VaR[t,3] = -ys[l1] * value                          # HS
    res = GARCHfit(window)
    s4 = res.seForecast
    VaR[t,4]=-s4*quantile(Normal(0,1),p)*value          # GARCH(1,1)
end
## GARCH VaR estimation will be slightly different from other languages
## this is due to GARCHfit choosing initial conditional vol = sample vol
		

Listing 8.7/8.8: Backtesting analysis in Python
Last updated June 2018

W1 = WE                                         # Python index starts at 0
m = ["EWMA", "MA", "HS", "GARCH"]
for i in range(4):
    VR = sum(y[W1:T] < -VaR[W1:T,i])/(p*(T-WE))
    s = np.std(VaR[W1:T, i], ddof=1)
    print ([i, m[i], VR, s])
plt.plot(y[W1:T])
plt.plot(VaR[W1:T])
plt.show()
plt.close()
		
Listing 8.7/8.8: Backtesting analysis in Julia
Last updated June 2018

W1 = WE + 1
for i in range(1,4)
    VR = sum(y[W1:T] .< -VaR[W1:T, i]) / (p * (T - WE))
    s = std(VaR[W1:T,i])
    println([i, "VR", VR, "VaR vol", s])
end
using Plots;
return plot([y, VaR[:,1], VaR[:,2], VaR[:,3], VaR[:,4]])
		

Listing 8.9/8.10: Bernoulli coverage test in Python
Last updated June 2018

def bern_test(p,v):
    lv = len(v)
    sv = sum(v)
    al = np.log(p)*sv + np.log(1-p)*(lv-sv)
    bl = np.log(sv/lv)*sv + np.log(1-sv/lv)*(lv-sv)
    return (-2*(al-bl))
		
Listing 8.9/8.10: Bernoulli coverage test in Julia
Last updated June 2018

function bern_test(p,v)
    lv = length(v)
    sv = sum(v)
    al = log(p)*sv + log(1-p)*(lv-sv)
    bl = log(sv/lv)*sv + log(1-sv/lv)*(lv-sv)
    return (-2*(al-bl))
end
		

Listing 8.11/8.12: Independence test in Python
Last updated June 2018

def ind_test(V):
    J = np.full([T,4], 0)
    for i in range(1,len(V)-1):
        J[i,0] = (V[i-1] == 0) & (V[i] == 0)
        J[i,1] = (V[i-1] == 0) & (V[i] == 1)
        J[i,2] = (V[i-1] == 1) & (V[i] == 0)
        J[i,3] = (V[i-1] == 1) & (V[i] == 1)
    V_00 = sum(J[:,0])
    V_01 = sum(J[:,1])
    V_10 = sum(J[:,2])
    V_11 = sum(J[:,3])
    p_00=V_00/(V_00+V_01)
    p_01=V_01/(V_00+V_01)
    p_10=V_10/(V_10+V_11)
    p_11=V_11/(V_10+V_11)
    hat_p = (V_01+V_11)/(V_00+V_01+V_10+V_11)
    al = np.log(1-hat_p)*(V_00+V_10) + np.log(hat_p)*(V_01+V_11)
    bl = np.log(p_00)*V_00 + np.log(p_01)*V_01 + np.log(p_10)*V_10 + np.log(p_11)*V_11
    return (-2*(al-bl))
		
Listing 8.11/8.12: Independence test in Julia
Last updated June 2018

function ind_test(V)
    J = fill!(Array{Float64}(T,4), 0)
    for i in range(2,length(V)-1)
        J[i,1] = (V[i-1] == 0) & (V[i] == 0)
        J[i,2] = (V[i-1] == 0) & (V[i] == 1)
        J[i,3] = (V[i-1] == 1) & (V[i] == 0)
        J[i,4] = (V[i-1] == 1) & (V[i] == 1)
    end
    V_00 = sum(J[:,1])
    V_01 = sum(J[:,2])
    V_10 = sum(J[:,3])
    V_11 = sum(J[:,4])
    p_00=V_00/(V_00+V_01)
    p_01=V_01/(V_00+V_01)
    p_10=V_10/(V_10+V_11)
    p_11=V_11/(V_10+V_11)
    hat_p = (V_01+V_11)/(V_00+V_01+V_10+V_11)
    al = log(1-hat_p)*(V_00+V_10) + log(hat_p)*(V_01+V_11)
    bl = log(p_00)*V_00 + log(p_01)*V_01 + log(p_10)*V_10 + log(p_11)*V_11
    return (-2*(al-bl))
end
		

Listing 8.13/8.14: Backtesting S&P 500 in Python
Last updated June 2018

W1 = WE
ya = y[W1:T]
VaRa = VaR[W1:T,]
m = ['EWMA', 'MA', 'HS', 'GARCH']
for i in range(4):
    q = y[W1:T] < -VaR[W1:T,i]
    v = VaRa*0
    v[q,i] = 1
    ber = bern_test(p, v[:,i])
    ind = ind_test(v[:,i])
    print ([i, m[i], ber, 1 - stats.chi2.cdf(ber, 1), ind, 1 - stats.chi2.cdf(ind, 1)])
		
Listing 8.13/8.14: Backtesting S&P 500 in Julia
Last updated June 2018

using Distributions;
W1 = WE+1
ya = y[W1:T]
VaRa = VaR[W1:T,:]
m = ["EWMA", "MA", "HS", "GARCH"]
for i in range(1,4)
    q = y[W1:T] .< -VaR[W1:T,i]
    v = VaRa .* 0
    v[q,i] = 1
    ber = bern_test(p, v[:,i])
    ind = ind_test(v[:,i])
    println([i, m[i], ber, 1-cdf(Chisq(1), ber), ind, 1-cdf(Chisq(1), ind)])
end
		

Listing 8.15/8.16: Backtest ES in Python
Last updated June 2018

VaR = np.full([T,2], np.nan)                                 # VaR forecasts
ES = np.full([T,2], np.nan)                                  # ES forecasts
for t in range(WE, T):
    t1 = t - WE
    t2 = t - 1
    window = y[t1:t2+1]
    s11 = lmbda * s11 + (1-lmbda) * y[t-1]**2
    VaR[t,0] = -stats.norm.ppf(p) * np.sqrt(s11)*value       # EWMA
    ES[t,0]=np.sqrt(s11)*stats.norm.pdf(stats.norm.ppf(p))/p
    ys = np.sort(window)
    VaR[t,1] = -ys[l1 - 1] * value                           # HS
    ES[t,1] = -np.mean(ys[0:l1]) * value
		
Listing 8.15/8.16: Backtest ES in Julia
Last updated June 2018

using Distributions;
VaR = fill!(Array{Float64}(T,2), NaN)                            # VaR forecasts
ES = fill!(Array{Float64}(T,2), NaN)                             # ES forecasts
for t in range(WE+1, T-WE)
    t1 = t - WE
    t2 = t - 1
    window = y[t1:t2]
    s11 = lambda*s11+(1-lambda)*y[t-1]^2
    VaR[t,1]=-quantile(Normal(0,1),p)*sqrt(s11)*value            # EWMA
    ES[t,1]=sqrt(s11)*pdf(Normal(0,1),quantile(Normal(0,1),p))/p
    ys = sort(window)
    VaR[t,2] = -ys[l1] * value                                   # HS
    ES[t,2] = -mean(ys[1:l1]) * value
end
		

Listing 8.17/8.18: ES in Python
Last updated June 2018

ESa = ES[W1:T,:]
VaRa = VaR[W1:T,:]
m = ["EWMA", "HS"]
for i in range(2):
    q = ya <= -VaRa[:,i]
    nES = np.mean(ya[q] / -ESa[q,i])
    print ([i, m[i], 'nES', nES])
		
Listing 8.17/8.18: ES in Julia
Last updated June 2018

ESa = ES[W1:T,:]
VaRa = VaR[W1:T,:]
m = ["EWMA", "HS"]
for i in range(1,2)
    q = ya .<= -VaRa[:,i]
    nES = mean(ya[q] ./ -ESa[q,i])
    println([i, m[i], "nES", nES])
end