Chapter 8. Backtesting and Stress Testing


Copyright 2016 Jon Danielsson. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at. http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Listing 8.1: Load data in R
Last edited: August 2016

library("tseries")    
library(zoo)
p = get.hist.quote(instrument = "^gspc", start = "1994-02-11",end="2016-08-30",quote="AdjClose",quiet=T) 
y=diff(log(p))      
y=coredata(y)      
		
Listing 8.2: Load data in Matlab
Last edited: August 2016

stocks = hist_stock_data('11021994','31082016','^gspc');
price=stocks.AdjClose(end:-1:1);  
y=diff(log(price));    
		

Listing 8.3: Set backtest up in R
Last edited: August 2016

TT = length(y)      
WE = 1000          
p = 0.01          
l1 = WE*p         
value = 1;         
VaR = matrix(nrow=TT,ncol=4)  
lambda = 0.94;
s11 = var(y[1:30]);    
for(t in 2:WE) s11 = lambda * s11 + (1-lambda) * y[t-1]^2
library(fGarch)
		
Listing 8.4: Set backtest up in Matlab
Last edited: August 2016

T = length(y);      
WE = 1000;          
p = 0.01;            
value = 1;          
l1 = WE*p ;          
VaR = NaN(T,4);      
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: Running backtest in R
Last edited: August 2016

for (t in (WE+1):TT){
  t1 = t-WE;                  
  t2 = t-1;                  
  window = y[t1:t2]         
  s11   = lambda * s11  + (1-lambda) * y[t-1]^2
  VaR[t,1] = -qnorm(p) * sqrt(s11) * value 
  VaR[t,2] = - sd(window) * qnorm(p)*value
  ys = sort(window)          
  VaR[t,3] = -ys[l1]*value  
  g=garchFit(formula = ~ garch (1,1), window ,trace=FALSE, include.mean=FALSE)
  par=g@fit$matcoef        
  s4=par[1]+par[2]*window[WE]^2+par[3]*g@h.t[WE]
  VaR[t,4] = -qnorm(p) * sqrt(s4) * value
}
		
Listing 8.6: Running backtest in Matlab
Last edited: August 2016

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 ;
  VaR(t,2) = -std(window)*norminv(p)*value;
  ys = sort(window);  
  VaR(t,3) = -ys(l1)*value;
    [par,ll,ht]=tarch(y1,1,0,1)
    h = par(1) + par(2) * y1(end)^2 + par(3) * ht(end)^2  
  VaR(t,4) = -norminv(p) * sqrt(h)  *value;
end
		

Listing 8.7: Backtesting analysis in R
Last edited: August 2016

W1=WE+1
for (i in 1:4){
  VR = sum(y[W1:TT]< -VaR[W1:TT,i])/(p*(T-WE))
  s = sd(VaR[W1:TT,i])                
  cat(i,"VR",VR,"VaR vol",s,"\n")    
}
matplot(cbind(y[W1:TT],VaR[W1:T,]),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.8: Backtesting analysis in Matlab
Last edited: 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: Bernoulli coverage test in R
Last edited: 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.10: Bernoulli coverage test in Matlab
Last edited: 2011

function res=bern_test(p,v)
  a=p^(sum(v))*(1-p)^(length(v)-sum(v));
  b=(sum(v)/length(v))^(sum(v))*(1-(sum(v)/length(v)))^(length(v)-sum(v));
  res=-2*log(a/b);
end  
		

Listing 8.11: Independence coverage test in R
Last edited: 2011

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)
  a=(1-hat_p)^(V_00+V_10)*(hat_p)^(V_01+V_11)
  b=(p_00)^(V_00)*(p_01)^(V_01)*(p_10)^(V_10)*p_11^(V_11)
  return(-2*log(a/b))
}
		
Listing 8.12: Independence coverage test in Matlab
Last edited: 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);
  a=(1-hat_p)^(V_00+V_10)*(hat_p)^(V_01+V_11);
  b=(p_00)^(V_00)*(p_01)^(V_01)*(p_10)^(V_10)*p_11^(V_11);
  res= -2*log(a/b);
end
		

Listing 8.13: Backtesting S&P-500 in R
Last edited: August 2016

W1=WE+1
ya=y[W1:TT]
VaRa=VaR[W1:TT,]
m=c("EWMA","MA","HS","GARCH")
for (i in 1:4){
  q= y[W1:TT]< -VaR[W1:TT,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.14: Backtesting S&P-500 in Matlab
Last edited: 2011

function test()
  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
end
		

Listing 8.15: Backtest ES in R
Last edited: Augst 2016

ES = matrix(nrow=TT,ncol=2)  
VaR = matrix(nrow=TT,ncol=2)  
for (t in (WE+1):TT){
  t1 = t-WE;    
  t2 = t-1;    
  window = y[t1:t2]     
  s11   = lambda * s11  + (1-lambda) * y[t-1]^2
  VaR[t,1] = -qnorm(p) * sqrt(s11) * value  
  ES[t,1] = sqrt(s11) * dnorm(qnorm(p)) / p
  ys = sort(window)      
  VaR[t,2] = -ys[l1]*value    
  ES[t,2] = -mean(ys[1:l1]) * value    
}
		
Listing 8.16: ES in Matlab
Last edited: 2011

VaR = NaN(T,2);    
ES = NaN(T,2);    
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 ;
  ES(t,1) = sqrt(s11) * normpdf(norminv(p)) / p;
  ys = sort(window);  
  VaR(t,2) = -ys(l1) * value;
  ES(t,2) = -mean(ys(1:l1)) * value;    
end
		

Listing 8.17: Backtest ES in R
Last edited: August 2016

ESa = ES[W1:TT,]
VaRa = VaR[W1:TT,]
for (i in 1:2){
  q = ya <= -VaRa[,i]
  nES = mean(ya[q] / -ESa[q,i])
  cat(i,"nES",nES,"\n")
}
		
Listing 8.18: ES in Matlab
Last edited: 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