# This file contains the R-code for the empirical analysis in Section 4 # of "Extreme Quantiles and Value-at-Risk" by V. Chernozhukov and S. Du. data<- as.matrix(read.table("Var_data.txt")) oxi.return <- data[,5]; oil.price <- data[,2]; dji.return <- data[,3]; lag.return <- data[,4]; data<- data.frame(oxi.return, lag.return, oil.price, dji.return) formla<- as.formula("oxi.return ~ lag.return + oil.price + dji.return") set.seed(0); est1<- rq(formla, tau=.001, method="pfn") res1.e<- summary.rq.nestextreme(est1,xi=.25, spacing=100, R=300) res1.n<- summary.rq(est1,se="rank") set.seed(0); est2<- rq(formla, tau=.01, method="pfn") res2.e<- summary.rq.nestextreme(est2,xi=.25,spacing=100, R=300) res2.n<- summary.rq(est2,se="rank") set.seed(0); est3<- rq(formla, tau=.05, method="pfn") res3.e<- summary.rq.nestextreme(est3,xi=.25,spacing=100, R=300) res3.n<- summary.rq(est3,se="rank") report.1<- rbind(res1.e$coef[,-2],res2.e$coef[,-2], res3.e$coef[,-2]) report<- report.1[, c(1,4, 2,3)] latex.table(as.matrix(report), file="C:\\aaa\\comput\\Extreme-Var\\Table1", caption="Estimation Results", dec=2) #============================ Tail Analysis set.seed(0); est2<- rq(formla, tau=.01); hill2<- summary.rq.hill(est2, R=300); set.seed(0); est3<- rq(formla, tau=.025); hill3<- summary.rq.hill(est3, R=300); set.seed(0); est4<- rq(formla, tau=.05); hill4<- summary.rq.hill(est4, R=300); set.seed(0); est5<- rq(formla, tau=.1); hill5<- summary.rq.hill(est5, R=300); set.seed(0); est6<- rq(formla, tau=.005); hill6<- summary.rq.hill(est6, R=300); report<- rbind(hill2, hill3, hill4) latex.table(as.matrix(report), file="C:\\aaa\\comput\\Extreme-Var\\Table2", caption="Estimation Results for EV Index") # ============================== # Extrapolated Fits extrap.fit<- extrapolate2.rq(formla, tau=.05, tau.e=.0001, xi=.23) simple.fit<- rq(formla, tau=.0001)$fitted n<- length(extrap.fit); postscript("C:\\aaa\\comput\\Extreme-VaR\\extrapolate.eps", width = 8, height = 5, horizontal = FALSE, onefile = FALSE, family="AvantGarde") plot((n-50):n, extrap.fit[(n-50):n], ylim=c(-.2,-.05), type="l", lty=2, xlab="time", ylab="Conditional .0001-Quantile", main="Extrapolalated vs Ordinary Estimates of Cond Quantiles") lines((n-50):n,simple.fit[(n-50):n], lty=1) legend("bottomright", c("ordinary fit", "extrap fit"), lty=c(1,2)) dev.off() # ================================ extrapolate.rq<- function(formula, tau, tau.e, xi=.2){ # a quick extrapolation estimator est<- rq(formula, tau=tau, method="pfn"); n<- length(est$fitted); m<- 2; est.m<- rq(formula, tau=tau*m, method="pfn"); fit.e<- (( tau.e/tau )^{-xi}-1)/(m^{-xi}-1)*(est.m$fitted -est$fitted) + est$fitted return(fit.e) } extrapolate2.rq<- function(formula, tau, tau.e, xi=.2){ # a quick extrapolation estimator est<- rq(formula, tau=tau, method="pfn"); n<- length(est$fitted); fit.e<- ( tau.e/tau )^{-xi}*est$fitted return(fit.e) } extrap.fit<- extrapolate.rq(formla, tau=.1, tau.e=.001, xi=.2) simple.fit<- rq(formla, tau=.001)$fitted n<- length(extrap.fit); plot(extrap.fit[(n-50):n], ylim=c(-.17,0), type="l", lty=2) lines(simple.fit[(n-50):n]) summary.rq.nestextreme<- function (object, R=100, xi=.3, spacing =20, ...) { mt <- terms(object) m <- model.frame(object) y <- model.response(m) x <- model.matrix(mt, m, contrasts = object$contrasts) wt <- model.weights(m) tau <- object$tau eps <- .Machine$double.eps^(2/3) coef <- coefficients(object) if (is.matrix(coef)) coef <- coef[, 1] vnames <- dimnames(x)[[2]] resid <- object$residuals n <- length(resid) p <- length(coef) rdf <- n - p if (!is.null(wt)) { resid <- resid * wt x <- x * wt y <- y * wt } infer.object <- rq.infer.nestextreme(x, y, tau = tau, R=R, xi=xi, spacing=20,...) coef <- array(coef, c(p, 5)) dimnames(coef) <- list(vnames, c("Value", "Pseudo Std. Er.", "Lower 5% Bound", "Upper 95% Bound", "Bias-Corrected Estimate") ) coef[, 3] <- infer.object$ciL.e coef[, 4] <- infer.object$ciU.e coef[, 2] <- (coef[, 4]- coef[, 3])/(2*1.69) coef[, 5] <- infer.object$BC.betat object <- object[c("call", "terms")] object$R<- infer.object$R object$coefficients <- coef object$tau <- tau class(object) <- "summary.rq" object } # function that does not complain about singular designs and records instead a NaN rq.fit.fn.nostop<- function (x, y, tau = 0.5, beta = 0.99995, eps = 1e-06) { n <- length(y) p <- ncol(x) if (n != nrow(x)) stop("x and y don't match n") if (tau < eps || tau > 1 - eps) stop("tau is outside (0,1)") rhs <- (1 - tau) * apply(x, 2, sum) d <- rep(1, n) u <- rep(1, n) wn <- rep(0, 10 * n) wn[1:n] <- (1 - tau) z <- .Fortran("rqfn", as.integer(n), as.integer(p), a = as.double(t(as.matrix(x))), c = as.double(-y), rhs = as.double(rhs), d = as.double(d), as.double(u), beta = as.double(beta), eps = as.double(eps), wn = as.double(wn), wp = double((p + 3) * p), aa = double(p * p), it.count = integer(2), info = integer(1), PACKAGE = "quantreg") if (z$info != 0) { coefficients <- NA; residuals<- NA; fit<- NA } else { coefficients <- -z$wp[1:p] residuals <- y - x %*% coefficients fitted<- x %*% coefficients } return(list(coefficients = coefficients, tau = tau, residuals = residuals, fit=fitted)) } rq.infer.nestextreme<- function(X,Y, tau=.1, xi=.3, R=R, spacing=20, ...) { if (tau>.5) {Y<- -Y; tau.e<- 1-tau } else tau.e<- tau n<- length(Y); p <- dim(X)[2] m<- (tau.e*n+spacing)/(tau.e*n); muX<- apply(X,2, mean) ; betat<- rq(Y~X[,-1], tau=tau.e, method="fn")$coef ; betatm<- rq(Y~X[,-1], tau=m*tau.e, method="fn")$coef ; An<- (m-1)*tau.e*sqrt(n/(tau.e*(1-tau.e)) )/(muX%*%( betatm - betat ) ) Res<- NULL; set.seed(1); for(i in 1:R) { Xs<- X; Ys<- (rexp(n)^{-xi}-1)/(-xi); betat.e<- rq.fit.fn.nostop(Xs,Ys, tau=tau.e)$coef ; betatm.e<- rq.fit.fn.nostop(Xs, Ys, tau=m*tau.e)$coef ; Ans<- (m-1)*tau.e*sqrt(n/(tau.e* (1-tau.e)) )/(-muX%*%( betat.e - betatm.e)); betat.model<- c( ((-log(1-tau.e))^(-xi) - 1)/(-xi), rep(0,p-1)); Bs.e<- Ans * (betat.e - betat.model); Res<- rbind(Res, c(Bs.e)); } object<- NULL names(betat)<- dimnames(X)[[2]] dimnames(Res)[[2]]<- dimnames(X)[[2]] if(tau <=.5) { object$betat<- betat object$BC.betat<- betat - apply(Res, 2, quantile, .5, na.rm=TRUE)/An ; object$ciL.e<- betat - apply(Res, 2, quantile, .95, na.rm=TRUE)/An ; object$ciU.e<- betat - apply(Res, 2, quantile, .05, na.rm=TRUE)/An ; } else { object$betat<- -betat object$BC.betat<- -(betat - apply(Res, 2, quantile, .5, na.rm=TRUE)/An) ; object$ciL.e<- -(betat - apply(Res, 2, quantile, .05, na.rm=TRUE)/An) ; object$ciU.e<- -(betat - apply(Res, 2, quantile, .95, na.rm=TRUE)/An) ; } object$R<- R-sum(is.na(Res[,1])); return(object) ; } infer.rq.hill.nestextreme<- function(x,y, tau=.1, R=R, ...) { if (tau>.5) {y<- -y; tau.e<- 1-tau } else tau.e<- tau n<- length(y); p <- dim(X)[2] Res<- NULL; set.seed(1); hill<- hill.rq(x,y, tau=tau.e); xi<- hill; for(i in 1:R) { ys<- (rexp(n)^{-xi}-1)/(-xi); hills<- hill.rq(x,ys, tau=tau.e); Res<- rbind(Res, c(hills)); } object<- NULL object$hill<- hill object$BC.hill<- hill- quantile(Res-hill, .5, na.rm=TRUE) ; object$ciL.e<- hill-quantile(Res-hill, .95, na.rm=TRUE) ; object$ciU.e<- hill-quantile(Res-hill, .05, na.rm=TRUE) ; object$se<- sqrt(var(Res)) ; object$R<- R-sum(is.na(Res[,1])); return(object) ; } summary.rq.hill<- function (object, R=300, ...) { mt <- terms(object) m <- model.frame(object) y <- model.response(m) x <- model.matrix(mt, m, contrasts = object$contrasts) wt <- model.weights(m) tau <- object$tau eps <- .Machine$double.eps^(2/3) coef <- coefficients(object) if (is.matrix(coef)) coef <- coef[, 1] vnames <- dimnames(x)[[2]] resid <- object$residuals n <- length(resid) p <- length(coef) rdf <- n - p if (!is.null(wt)) { resid <- resid * wt x <- x * wt y <- y * wt } res <- infer.rq.hill.nestextreme(x,y, tau=tau,R=R, ...) res<- array(res,c(1,4)); dimnames(res) <- list("EV Index", c("Estimate", "Bias-Corrected", "Lower 5%", "Upper 95%") ) return(res); } hill.rq<- function (x,y, tau) { est<- rq.fit.fn.nostop(x,y, tau=tau); fitted<- est$fit; y<- est$fit + est$res; return( sum(log( ifelse(y <=fitted, abs(y/fitted), 1 )))/(length(y)*tau) ) } pickands.rq<- function (x,y, tau) { # a quick pickands-type estimator est<- rq.fit.fn.nostop(x,y, tau=tau); n<- length(est$fit); m<- 2 est.m<- rq.fit.fn.nostop(x,y, tau=tau*m); est.minv<- rq.fit.fn.nostop(x,y,tau=tau/m); fit<- mean(est$fit) fit.m<- mean(est.m$fit) fit.minv<- mean(est.minv$fit) return( - log(( fit.m - fit )/( fit - fit.minv))/ log(m) ) } infer.rq.pickands.nestextreme<- function(x,y, tau=.1, R=R, ...) { if (tau>.5) {y<- -y; tau.e<- 1-tau } else tau.e<- tau n<- length(y); p <- dim(X)[2] Res<- NULL; set.seed(1); pickands<- pickands.rq(x,y, tau=tau.e); xi<- pickands; for(i in 1:R) { ys<- (rexp(n)^{-xi}-1)/(-xi); pickandss<- pickands.rq(x,ys, tau=tau.e); Res<- rbind(Res, c(pickandss)); } object<- NULL object$pickands<- pickands object$BC.pickands<- pickands- quantile(Res-pickands, .5, na.rm=TRUE) ; object$ciL.e<- pickands-quantile(Res-pickands, .95, na.rm=TRUE) ; object$ciU.e<- pickands-quantile(Res-pickands, .05, na.rm=TRUE) ; object$se<- sqrt(var(Res)) ; object$R<- R-sum(is.na(Res[,1])); return(object) ; } summary.rq.pickands<- function (object, R=300, ...) { mt <- terms(object) m <- model.frame(object) y <- model.response(m) x <- model.matrix(mt, m, contrasts = object$contrasts) wt <- model.weights(m) tau <- object$tau eps <- .Machine$double.eps^(2/3) coef <- coefficients(object) if (is.matrix(coef)) coef <- coef[, 1] vnames <- dimnames(x)[[2]] resid <- object$residuals n <- length(resid) p <- length(coef) rdf <- n - p if (!is.null(wt)) { resid <- resid * wt x <- x * wt y <- y * wt } res <- infer.rq.pickands.nestextreme(x,y, tau=tau,R=R, ...) res<- array(res,c(1,4)); dimnames(res) <- list("EV Index", c("Estimate", "Bias-Corrected", "Lower 5%", "Upper 95%") ) return(res); } pickands2<- summary.rq.pickands(est2, R=1);