Matlab and Julia Chapter 8. Backtesting and Stress Testing

Chapter 8. Backtesting and Stress Testing

Matlab 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 8.1/8.2
% Load data in MATLAB
price = csvread('index.csv', 1, 0);
y=diff(log(price)); % get returns
Listing 8.1/8.2
Load data in Julia
using CSV, DataFrames;
price = CSV.read("index.csv", DataFrame);
y = diff(log.(price[:,1])); # get returns
Listing 8.3/8.4
% Set backtest up in MATLAB
T = length(y);  % number of obs for return y
WE = 1000;      % estimation window length 
p = 0.01;       % probability
l1 = ceil(WE*p) ;     % HS observation
value = 1;      % portfolio value
VaR = NaN(T,4); % matrix for forecasts
lambda = 0.94;
s11 = var(y);
for t = 2:WE
    s11=lambda*s11+(1-lambda)*y(t-1)^2;
end
Listing 8.3/8.4
Set backtest up in Julia
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
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 MATLAB
for t = WE+1:T
    t1 = t-WE;          % start of the 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) = -norminv(p) * sqrt(s11)  *value; % EWMA
    VaR(t,2) = -std(window)*norminv(p)*value; % MA
    ys = sort(window);
    VaR(t,3) = -ys(l1)*value; % HS
    [par,ll,ht]=tarch(window,1,0,1);
    h=par(1)+par(2)*window(end)^2+par(3)*ht(end);
    VaR(t,4) = -norminv(p)*sqrt(h)*value; % GARCH(1,1)
end
Listing 8.5/8.6
Running backtest in Julia
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 MATLAB
names = ["EWMA", "MA", "HS", "GARCH"];
for i=1:4
    VR = length(find(y(WE+1:T)<-VaR(WE+1:T,i)))/(p*(T-WE)); 
    s = std(VaR(WE+1:T,i));         
    disp([names(i), "Violation Ratio:", VR, "Volatility:", s])          
end
plot([y(WE+1:T) VaR(WE+1:T,:)])
Listing 8.7/8.8
Backtesting analysis in Julia
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 MATLAB
function res=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)
	res=-2*(al-bl);
end
Listing 8.9/8.10
Bernoulli coverage test in Julia
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 MATLAB
function res=ind_test(V)
	T=length(V);
	J=zeros(T,4);
	for i = 2:T
		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;
	res= -2*(al-bl);
end
Listing 8.11/8.12
Independence test in Julia
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 MATLAB
names = ["EWMA", "MA", "HS", "GARCH"];
ya=y(WE+1:T);
VaRa=VaR(WE+1:T,:);
for i=1:4
	q=find(y(WE+1:T)<-VaR(WE+1:T,i));
	v=VaRa*0;
	v(q,i)=1;
	ber=bern_test(p,v(:,i));
	in=ind_test(v(:,i));
	disp([names(i), "Bernoulli Statistic:", ber, "P-value:", 1-chi2cdf(ber,1),...
    "Independence Statistic:", in, "P-value:", 1-chi2cdf(in,1)])
end
Listing 8.13/8.14
Backtesting S&P 500 in Julia
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 MATLAB
VaR = NaN(T,2);  % VaR forecasts for 2 models 
ES = NaN(T,2);   % ES forecasts for 2 models 
for t = 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) = -norminv(p) * sqrt(s11)  *value; % EWMA
    ES(t,1) = sqrt(s11) * normpdf(norminv(p)) / p;
    ys = sort(window);
    VaR(t,2) = -ys(l1) * value;          % HS
    ES(t,2) = -mean(ys(1:l1)) * value;  
end
Listing 8.15/8.16
Backtest ES in Julia
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 MATLAB
names = ["EWMA", "HS"];
VaRa = VaR(WE+1:T,:);
ESa = ES(WE+1:T,:);
for i = 1:2
	q = find(ya <= -VaRa(:,i));
	nES = mean(ya(q) ./ -ESa(q,i));
	disp([names(i), nES])
end
Listing 8.17/8.18
ES in Julia
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
Financial Risk Forecasting
Market risk forecasting with R, Julia, Python and Matlab. Code, lecture slides, implementation notes, seminar assignments and questions.
© All rights reserved, Jon Danielsson,