# Chapter 8. Backtesting and Stress Testing

##### 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');
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