p = read.csv('index.csv')
y=diff(log(p$Index)) # get returns
import numpy as np
price = np.loadtxt('index.csv',delimiter=',',skiprows=1)
y = np.diff(np.log(price), n=1, axis=0) # get returns
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
import numpy as np
T = len(y) # number of obs for y
WE = 1000 # estimation window length
p = 0.01 # probability
l1 = int(WE * p) # HS observation
value = 1 # portfolio value
VaR = np.full([T,4], np.nan) # matrix for forecasts
## EWMA setup
lmbda = 0.94
s11 = np.var(y[1:30])
for t in range(1,WE):
s11=lmbda*s11+(1-lmbda)*y[t-1]**2
## 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)
import numpy as np
from scipy import stats
from arch import arch_model
for t in range(WE, T):
t1 = t - WE # start of data window
t2 = t - 1 # end of data window
window = y[t1:t2+1] # data for estimation
s11 = lmbda * s11 + (1-lmbda) * y[t-1]**2
VaR[t,0]=-stats.norm.ppf(p)*np.sqrt(s11)*value # EWMA
VaR[t,1]=-np.std(window,ddof=1)*stats.norm.ppf(p)*value # MA
ys = np.sort(window)
VaR[t,2] = -ys[l1 - 1] * value # HS
am = arch_model(window, mean = 'Zero',vol='Garch',
p=1, o=0, q=1, dist='Normal')
res = am.fit(update_freq=0, disp='off',show_warning=False)
par = [res.params[0], res.params[1], res.params[2]]
s4 = par[0] + par[1] * window[WE - 1]**2 + par[
2] * res.conditional_volatility[-1]**2
VaR[t,3] = -np.sqrt(s4) * stats.norm.ppf(p) * value # GARCH(1,1)
## GARCH optimization in Python has some convergence issues
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")
W1 = WE # Python index starts at 0
m = ["EWMA", "MA", "HS", "GARCH"]
for i in range(4):
VR = sum(y[W1:T] < -VaR[W1:T,i])/(p*(T-WE))
s = np.std(VaR[W1:T, i], ddof=1)
print ([i, m[i], VR, s])
plt.plot(y[W1:T])
plt.plot(VaR[W1:T])
plt.show()
plt.close()
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))
}
def bern_test(p,v):
lv = len(v)
sv = sum(v)
al = np.log(p)*sv + np.log(1-p)*(lv-sv)
bl = np.log(sv/lv)*sv + np.log(1-sv/lv)*(lv-sv)
return (-2*(al-bl))
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))
}
def ind_test(V):
J = np.full([T,4], 0)
for i in range(1,len(V)-1):
J[i,0] = (V[i-1] == 0) & (V[i] == 0)
J[i,1] = (V[i-1] == 0) & (V[i] == 1)
J[i,2] = (V[i-1] == 1) & (V[i] == 0)
J[i,3] = (V[i-1] == 1) & (V[i] == 1)
V_00 = sum(J[:,0])
V_01 = sum(J[:,1])
V_10 = sum(J[:,2])
V_11 = sum(J[:,3])
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 = np.log(1-hat_p)*(V_00+V_10) + np.log(hat_p)*(V_01+V_11)
bl = np.log(p_00)*V_00 + np.log(p_01)*V_01 + np.log(p_10)*V_10 + np.log(p_11)*V_11
return (-2*(al-bl))
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")
}
W1 = WE
ya = y[W1:T]
VaRa = VaR[W1:T,]
m = ['EWMA', 'MA', 'HS', 'GARCH']
for i in range(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])
print ([i, m[i], ber, 1 - stats.chi2.cdf(ber, 1), ind, 1 - stats.chi2.cdf(ind, 1)])
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
}
VaR = np.full([T,2], np.nan) # VaR forecasts
ES = np.full([T,2], np.nan) # ES forecasts
for t in range(WE, T):
t1 = t - WE
t2 = t - 1
window = y[t1:t2+1]
s11 = lmbda * s11 + (1-lmbda) * y[t-1]**2
VaR[t,0] = -stats.norm.ppf(p) * np.sqrt(s11)*value # EWMA
ES[t,0]=np.sqrt(s11)*stats.norm.pdf(stats.norm.ppf(p))/p
ys = np.sort(window)
VaR[t,1] = -ys[l1 - 1] * value # HS
ES[t,1] = -np.mean(ys[0:l1]) * value
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")
}
ESa = ES[W1:T,:]
VaRa = VaR[W1:T,:]
m = ["EWMA", "HS"]
for i in range(2):
q = ya <= -VaRa[:,i]
nES = np.mean(ya[q] / -ESa[q,i])
print ([i, m[i], 'nES', nES])