R and Matlab Chapter 8. Backtesting and Stress Testing

# Chapter 8. Backtesting and Stress Testing

### R and Matlab

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
p = read.csv('index.csv')
y = diff(log(p$Index)) # get returns  ##### Listing 8.1/8.2 ##### % Load data in MATLAB price = csvread('index.csv', 1, 0); y=diff(log(price)); % get returns  ##### Listing 8.3/8.4 ##### Set backtest up in R WE = 1000 # estimation window length p = 0.01 # probability l1 = ceiling(WE*p) # HS quantile portfolio_value = 1 # portfolio value VaR = matrix(nrow=length(y),ncol=4) # matrix for forecasts lambda = 0.94 s11 = var(y) for(t in 2:WE) { s11=lambda*s11+(1-lambda)*y[t-1]^2 }  ##### 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.5/8.6 ##### Running backtest in R library(rugarch) spec = ugarchspec(variance.model = list( garchOrder = c(1, 1)), mean.model = list( armaOrder = c(0,0),include.mean = FALSE)) start_time <- Sys.time() for (t in (WE+1):length(y)){ t1 = t-WE; # start of the data window t2 = t-1; # end of the data window window = y[t1:t2] # data for estimation s11 = lambda*s11 + (1-lambda)*y[t-1]^2 VaR[t,1] = -qnorm(p) * sqrt(s11) * portfolio_value # EWMA VaR[t,2] = - sd(window) * qnorm(p)*portfolio_value # MA ys = sort(window) VaR[t,3] = -ys[l1]*portfolio_value # HS res = ugarchfit(spec = spec, data = window,solver="hybrid") omega = res@fit$coef['omega']
alpha = res@fit$coef['alpha1'] beta = res@fit$coef['beta1']
sigma2 = omega + alpha * tail(window,1)^2 + beta * tail(res@fit\$var,1)
VaR[t,4] = -sqrt(sigma2) * qnorm(p) * portfolio_value       # GARCH(1,1)
}
end_time <- Sys.time()
print(end_time - start_time)

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


##### Backtesting analysis in R
for (i in 1:4){
VR = sum(y[(WE+1):length(y)]< -VaR[(WE+1):length(y),i])/(p*(length(y)-WE))
s = sd(VaR[(WE+1):length(y),i])
cat(i,"VR",VR,"VaR vol",s,"\n")
}
matplot(cbind(y,VaR),type='l',col=1:5,las=1,ylab="",lty=1:5)
legend("topleft",legend=c("Returns","EWMA","MA","HS","GARCH"),lty=1:5,col=1:5,bty="n")

##### % 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,:)])


##### Bernoulli coverage test in R
bern_test=function(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))
}

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


##### Independence test in R
ind_test=function(V){
J=matrix(ncol=4,nrow=length(V))
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
}
V_00=sum(J[,1],na.rm=TRUE)
V_01=sum(J[,2],na.rm=TRUE)
V_10=sum(J[,3],na.rm=TRUE)
V_11=sum(J[,4],na.rm=TRUE)
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))
}

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


##### Backtesting S&P 500 in R
W1=WE+1
ya=y[W1:length(y)]
VaRa=VaR[W1:length(y),]
m=c("EWMA","MA","HS","GARCH")
for (i in 1:4){
q= y[W1:length(y)]< -VaR[W1:length(y),i]
v=VaRa*0
v[q,i]=1
ber=bern_test(p,v[,i])
ind=ind_test(v[,i])
cat(i,m[i], "\n",
"Bernoulli - ","Test statistic:",ber,"  p-value:",1-pchisq(ber,1),"\n",
"Independence - ", "Test statistic:",ind,"  p-value:",1-pchisq(ind,1),"\n \n")
}

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


##### Backtest ES in R
VaR2 = matrix(nrow=length(y), ncol=2)                    # VaR forecasts for 2 models
ES = matrix(nrow=length(y), ncol=2)                      # ES forecasts for 2 models
for (t in (WE+1):length(y)){
t1 = t-WE;
t2 = t-1;
window = y[t1:t2]
s11 = lambda * s11  + (1-lambda) * y[t-1]^2
VaR2[t,1] = -qnorm(p) * sqrt(s11) * portfolio_value    # EWMA
ES[t,1] = sqrt(s11) * dnorm(qnorm(p)) / p
ys = sort(window)
VaR2[t,2] = -ys[l1] * portfolio_value                  # HS
ES[t,2] = -mean(ys[1:l1]) * portfolio_value
}

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


##### Backtest ES in R
ESa = ES[W1:length(y),]
VaRa = VaR2[W1:length(y),]
for (i in 1:2){
q = ya <= -VaRa[,i]
nES = mean(ya[q] / -ESa[q,i])
cat(i,"nES",nES,"\n")
}

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


##### Financial Risk Forecasting
Market risk forecasting with R, Julia, Python and Matlab. Code, lecture slides, implementation notes, seminar assignments and questions.