 # Chapter 8. Backtesting and Stress Testing (in Python/Julia)

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 8.1/8.2: Load data in Python Last updated June 2018

import numpy as np
y = np.diff(np.log(price), n=1, axis=0) # get returns

##### Listing 8.1/8.2: Load data in Julia Last updated June 2022

using CSV, DataFrames;
y = diff(log.(price[:,1])); # get returns


##### Listing 8.3/8.4: Set backtest up in Python Last updated July 2020

from math import ceil
T = len(y)           # number of obs for y
WE = 1000  # estimation window length
p = 0.01   # probability
l1 = ceil(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)
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 July 2020

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


##### Listing 8.5/8.6: Running backtest in Python Last updated July 2020

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', rescale = False)
res = am.fit(update_freq=0, disp = 'off', show_warning=False)
par = [res.params, res.params, res.params]
s4 = par + par * 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)

##### Listing 8.5/8.6: Running backtest in Julia Last updated July 2020

using Distributions, ARCHModels;
for t in WE+1:T
t1 = t - WE      # start of data window
t2 = t - 1       # end of data window
window = y[t1:t2]          # data for estimation
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
garch1_1 = fit(GARCH{1,1}, window; meanspec = NoIntercept);
cond_vol = predict(garch1_1, :volatility)
VaR[t,4]= -cond_vol * quantile(Normal(0,1),p) * value    # GARCH(1,1)
end


##### Listing 8.7/8.8: Backtesting analysis in Python Last updated July 2020

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 (m[i], "\n",
"Violation ratio:", round(VR, 3), "\n",
"Volatility:", round(s,3), "\n")
plt.plot(y[W1:T])
plt.plot(VaR[W1:T])
plt.title("Backtesting")
plt.show()
plt.close()

##### Listing 8.7/8.8: Backtesting analysis in Julia Last updated July 2020

using Plots;
W1 = WE + 1
names = ["Returns" "EWMA" "MA" "HS" "GARCH"]
for i in 1:4
VR = sum(y[W1:T] .< -VaR[W1:T, i]) / (p * (T - WE))
s = std(VaR[W1:T,i])
println(names[i+1], "\n", "VR: ", round(VR, digits = 3), " VaR volatility: ", round(s, digits = 4))
end
plot([y, VaR[:,1], VaR[:,2], VaR[:,3], VaR[:,4]], label = names, title = "Backtesting")


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

function ind_test(V)
J = fill(0, (T,4))
for i in 2:length(V)
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 July 2020

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 (m[i], "\n",
"Bernoulli:", "Test statistic =", round(ber,3), "p-value =", round(1 - stats.chi2.cdf(ber, 1),3), "\n",
"Independence:", "Test statistic =", round(ind,3), "p-value =", round(1 - stats.chi2.cdf(ind, 1),3), "\n")

##### Listing 8.13/8.14: Backtesting S&P 500 in Julia Last updated July 2020

W1 = WE+1
ya = y[W1:T]
VaRa = VaR[W1:T,:]
names = ["EWMA", "MA", "HS", "GARCH"]
for i in 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(names[i], "\n",
"Bernoulli statistic: ", round(ber, digits = 2), " p-value: ", round(1-cdf(Chisq(1), ber), digits = 4), "\n",
"Independence statistic: ", round(ind, digits = 2), " p-value: ", round(1-cdf(Chisq(1), ind), digits = 4), "\n")
end


##### Listing 8.15/8.16: Backtest ES in Python Last updated July 2020

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

VaR = fill(NaN, (T,2))         # VaR forecasts
ES = fill(NaN, (T,2))          # ES forecasts
for t in WE+1:T
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 July 2020

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 (m[i], 'nES = ', round(nES,3))

##### Listing 8.17/8.18: ES in Julia Last updated July 2020

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