Chapter 8. Backtesting and Stress Testing (in R/MATLAB)


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 R
Last updated August 2019

p = read.csv('index.csv')
y=diff(log(p$Index))      # get returns
		
Listing 8.1/8.2: Load data in MATLAB
Last updated August 2016

price = csvread('index.csv', 1, 0);
y=diff(log(price));                 % get returns
		

Listing 8.3/8.4: Set backtest up in R
Last updated August 2019

WE = 1000                                         # estimation window length
p = 0.01                                          # probability
l1 = WE*p                                         # HS observation
portfolio_value = 1;                              # portfolio value
VaR = matrix(nrow=length(y),ncol=4)               # matrix for forecasts
## EWMA setup
lambda = 0.94;
s11 = var(y[1:30]);
for(t in 2:WE) s11=lambda*s11+(1-lambda)*y[t-1]^2
		
Listing 8.3/8.4: Set backtest up in MATLAB
Last updated August 2016

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

Listing 8.5/8.6: Running backtest in R
Last updated August 2019

## this will take several minutes to run
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[1]
  alpha = res@fit$coef[2]
  beta = res@fit$coef[3]
  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)
		
Listing 8.5/8.6: Running backtest in MATLAB
Last updated June 2018

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.7/8.8: Backtesting analysis in R
Last updated August 2019

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")
		
Listing 8.7/8.8: Backtesting analysis in MATLAB
Last updated 2011

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([i VR s])
end
plot([y(WE+1:T) VaR(WE+1:T,:)])
		

Listing 8.9/8.10: Bernoulli coverage test in R
Last updated August 2016

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))
}
		
Listing 8.9/8.10: Bernoulli coverage test in MATLAB
Last updated 2011

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.11/8.12: Independence test in R
Last updated June 2018

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))
}
		
Listing 8.11/8.12: Independence test in MATLAB
Last updated 2011

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.13/8.14: Backtesting S&P 500 in R
Last updated August 2016

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],'Bernoulli',ber,1-pchisq(ber,1),"independence",ind,1-pchisq(ind,1),"\n")
}
		
Listing 8.13/8.14: Backtesting S&P 500 in MATLAB
Last updated June 2018

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([i,ber,1-chi2cdf(ber,1),in,1-chi2cdf(in,1)])
end
		

Listing 8.15/8.16: Backtest ES in R
Last updated August 2019

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) * value   # EWMA
  ES[t,1] = sqrt(s11) * dnorm(qnorm(p)) / p
  ys = sort(window)
  VaR2[t,2] = -ys[l1]*value                   # HS
  ES[t,2] = -mean(ys[1:l1]) * value
}
		
Listing 8.15/8.16: Backtest ES in MATLAB
Last updated 2011

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.17/8.18: Backtest ES in R
Last updated August 2019

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")
}
		
Listing 8.17/8.18: ES in MATLAB
Last updated 2011

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([i,nES])
end