# This file contains R-functions used to implement the empirical analysis in Section 4 # of "Extreme Quantiles and Value-at-Risk" by V. Chernozhukov and S. Du. 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) } 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); }