rm(list=ls()) source("MRFzR_FunkcjeBlok2.R") require(zoo) require(ggplot2) require(rugarch) require(forecast) require(binom) require(knitr) require(car) source("LoadFundData.R") # source("LoadWigData.R") par(mfrow=c(2,1), cex = 0.7, bty="l") plot(P, main="Level", xlab="", ylab="") plot(r, main="returns", xlab="", ylab="") ##################################### # VaR model settings # ##################################### T <- length(r) # full sample, nobs N <- 1250 # we shorten the sample to the last 5 years r <- tail(r,N) # log-returns R <- coredata(r) p <- 0.05 # tolerance level backN <- 250 # backtest sample w_length <- 1000 # rolling window length ## Historical simulation varHS <- rollapply(r, width=w_length, function(w){quantile(w,p)}, by=1, align="right") varHS <- lag(varHS, -1) # lags so that VaR(t+1) = quantile(t) rr <- r[index(varHS)] # rr - realized returns etaHS <- tail(rr<=varHS, backN) # VaR violations) nHS <- sum(etaHS); nHS # number of exceedances piHS <- mean(etaHS); piHS # share of exceedances dev.off() VaRplot(alpha=p, actual=tail(rr,backN), VaR=tail(varHS,backN)) title(main="Historical simulation: VaR violations") # substitute for selected method (from Topic7a.R) var = varHS eta = etaHS pi = piHS n1 = nHS ############################### # 1. Regression test # ############################### c(pi,p) Acf(as.numeric(coredata(eta))) # 1a. Unconditional coverage. For eta[t] = pi + eps[t] H0: pi = p reg1a <- lm(eta~1) summary(reg1a) linearHypothesis(reg1a, hypothesis.matrix=1, rhs=p, test="Chisq") # H0: pi = p # 1b. Independence. For eta[t] = pi + rho*eta[t-1] + eps[t] H0: rho = 0 reg1b <- lm(eta[-1]~lag(eta,-1)) summary(reg1b) linearHypothesis(reg1b, hypothesis.matrix=c(0,1), rhs=0, test="Chisq") # 1c. Conditional coverage. For eta[t] = pi + rho*eta[t-1] + eps[t] H0: pi = p, rho = 0 reg1c <- lm(eta[-1]~lag(eta,-1)) summary(reg1c) linearHypothesis(reg1c, hypothesis.matrix=diag(1,2), rhs=c(p,0), test="Chisq") ###################################################### # 2. Unconditional coverage, Binomial distribution # # H0: n1 ~ Binom(n,p), 95% confidence interval # ###################################################### 95% confidence interval kable(cbind(seq(0,20,2),pbinom(seq(0,20,2),backN,p)),col.names = c("n1","cdf"),digits=3) LB = qbinom(0.025,backN,p) UB = qbinom(0.975,backN,p) print(paste0("LB=",LB,"; UB=",UB,", n. violations =",n1)) # Kupiec test n <- length(eta); n n1 <- sum(eta); n1 n0 <- n - n1; n0 pi <- n1/n; pi # test stat. LR_uc = (p/pi)^n1 *((1-p)/(1-pi))^n0 stat_uc = -2*log(LR_uc) prob_uc = 1 - pchisq(q=stat_uc,df=1, lower.tail=TRUE) paste0("Kupiec test stat.: ", round(stat_uc,3),"; p-value:", round(prob_uc,3)) # with rugarch package temp <- VaRTest(alpha=p, actual=rr, VaR=var) paste0("Kupiec test stat.:", round(temp$uc.LRstat,3),"; p-value:", round(temp$uc.LRp,3)) ########################################################### # 3. Christofersen independence test # # H0: p(eta[t]=1|eta[t-1]=1) = p(eta[t]=1|eta[t-1]=0) # ########################################################### Acf(as.numeric(coredata(tail(eta,backN)))) ## Ljung-Box test, H0: independence Box.test(as.numeric(coredata(tail(eta, backN))), 1, type="Ljung-Box") # Christofersen independence test eta1 <- coredata(eta[-length(eta)]) eta0 <- coredata(eta[-1]) n00 = sum(!eta1 & !eta0) # no exceedance after no exceedance n01 = sum(!eta1 & eta0) # exceedance after no exceedance n10 = sum( eta1 & !eta0) # no exceedance after exceedance n11 = sum( eta1 & eta0) # exceedance after exceedance #n0 = n00 + n10 #n1 = n01 + n11 pi0 = n01 / (n00+n01) # prob. of exceedance after no exceedance pi1 = n11 / (n10+n11) # prob. of exceedance after exceedance pi = (n01+n11) / (n00+n01+n10+n11) c(pi0,pi1,pi) #LR_ind = pi^n1*(1-pi)^n0 / (pi0^n01 * pi1^n11 * (1-pi0)^n00 * (1-pi1)^n10) LR_ind = (pi/pi0)^n01 * (pi/pi1)^n11 * ((1-pi)/(1-pi0))^n00 * ((1-pi)/(1-pi1))^n10 stat_ind = -2*log(LR_ind) prob_ind = 1 - pchisq(stat_ind,1) paste0("Ch1 test stat.: ", round(stat_ind,3),"; p-value:", round(prob_ind,3)) ############################################################## # 4. Christofersen conditional coverage test # # H0: p(eta[t]=1|eta[t-1]=1) = p(eta[t]=1|eta[t-1]=0) = p # ############################################################## # Christofersen conditional coverage test LR_cc = p^n1*(1-p)^n0 / (pi0^n01 * pi1^n11 * (1-pi0)^n00 * (1-pi1)^n10) stat_cc = -2*log(LR_cc) prob_cc = 1 - pchisq(stat_cc,1) prob_cc paste0("Ch2 test stat.: ", round(stat_cc,3),"; p-value:", round(prob_cc,3)) # stat_cc = stat_uc + stat_ind c(stat_uc+stat_ind,stat_cc) # with rugarch package temp <- VaRTest(alpha=p, actual=coredata(rr), VaR=coredata(var)) paste0("Ch2 test stat.: ", round(temp$cc.LRstat,3),"; p-value:", round(temp$cc.LRp,3)) ############################## # Parametric models # ############################## MAmean <- rollapply(r, width=w_length, mean, by=1, align="right") MAstd <- rollapply(r, width=w_length, sd, by=1, align="right") var <- MAmean + MAstd*qdist("norm", p); # var <- MAmean + MAstd*qdist("std", p, shape=5); var <- lag(var, -1) dev.off() VaRplot(alpha=p, actual=rr, VaR=var) title(main="Parametric model: VaR violations") temp <- VaRTest(alpha=p, actual=coredata(rr), VaR=coredata(var)) paste0("Kupiec test stat.: ", temp$uc.LRstat,"; p-value:", temp$uc.LRp) paste0("Ch2 test stat.: ", temp$cc.LRstat,"; p-value:", temp$cc.LRp) ################# # EWMA model # ################# lambda <- 0.94 EWMAspec <- ugarchspec(mean.model=list(armaOrder=c(0,0), include.mean=FALSE), variance.model=list(model="iGARCH", garchOrder=c(1,1)), fixed.pars=list(alpha1=1-lambda, omega=0), distribution.model = "norm") #fixed.pars=list(alpha1=1-lambda, omega=0,shape=5), distribution.model = "std") ## rolling VaR varEWMA <- rollapply(r, width=w_length, function(w) { # here we don't estimate but calibrate # forecast on the basis of specification frc <- ugarchforecast(EWMAspec,data=w, n.ahead=1) quantile(frc, p) }, by=1, align="right") varEWMA <- lag(varEWMA, -1) VaRplot(alpha=p, actual=rr, VaR=varEWMA) title(main="EWMA: VaR violations") temp <- VaRTest(alpha=p, actual=coredata(rr), VaR=coredata(varEWMA)) paste0("Kupiec test stat.: ", temp$uc.LRstat,"; p-value:", temp$uc.LRp) paste0("Ch2 test stat.: ", temp$cc.LRstat,"; p-value:", temp$cc.LRp) ################## # GARCH model # ################## GARCHspec <- ugarchspec(mean.model=list(armaOrder=c(0,0), include.mean=TRUE), variance.model=list(model="sGARCH", garchOrder=c(1,1)), distribution.model = "norm") # fixed.pars=list(shape=5), distribution.model = "std") varGARCH <- rollapply(r, width=w_length, function(w) { fit <- ugarchfit(data=w, spec=GARCHspec, solver="hybrid") frc <- ugarchforecast(fit, n.ahead=1) quantile(frc, p) }, by=1, align="right") varGARCH <- lag(varGARCH, -1) VaRplot(alpha=p, actual=rr, VaR=varGARCH) title(main="GARCH: VaR violation") temp <- VaRTest(alpha=p, actual=coredata(rr), VaR=coredata(varGARCH)) paste0("Kupiec test stat: ", temp$uc.LRstat,"; p-value:", temp$uc.LRp) paste0("Ch2 test stat: ", temp$cc.LRstat,"; p-value:", temp$cc.LRp) ############################################# # Backtesing with ugarchroll from rugarch # ############################################# step = 10; # how often we reestimate the GARCH model varGARCHRoll <- ugarchroll(spec=GARCHspec, data=r, refit.every=step, forecast.length=backN, refit.window="moving", window.size=w_length, calculate.VaR=TRUE, VaR.alpha=c(0.01, 0.025, 0.05)) report(varGARCHRoll, VaR.alpha = 0.05) report(varGARCHRoll, VaR.alpha = 0.01) ######################### # Appendix: tests power # ######################### backN <- 250 alpha <- 0.05 # significance level (donot confuse with tolerance level) ############################### # Traffic lights method # ############################### nexceed <- list(0:4,5:9,10:backN) pval <- 0.01*c(1,2,3,5) tab1 <- matrix(0, nrow=4, ncol=3) colnames(tab1) <- c("n1<=4", "4