{"word_count": 436, "character_count": 3650}
This week, we conclude the part on extremes with an application of extreme value theory to risk measures. We have seen last week that, if we assume that above a threshold , a Generalized Pareto Distribution will fit nicely, then we can use it to derive an estimator of the quantile function (for percentages such that the quantile is larger than the threshold)It the threshold is , i.e. we keep the  largest observations to fit a GPD, then this estimator can be writtenThe code we wrote last week was the following (here based on log-returns of the SP500 index, and we focus on large losses, i.e. large values of the opposite of log returns, plotted below)&#62; library(tseries)&#62; X=get.hist.quote("^GSPC")&#62; T=time(X)&#62; D=as.POSIXlt(T)&#62; Y=X$Close&#62; R=diff(log(Y))&#62; D=D[-1]&#62; X=-R&#62; plot(D,X)&#62; library(evir)&#62; GPD=gpd(X,quantile(X,.975))&#62; xi=GPD$par.ests[1]&#62; beta=GPD$par.ests[2]&#62; u=GPD$threshold&#62; QpGPD=function(p){+ u+beta/xi*((100/2.5*(1-p))^(-xi)-1)+ }&#62; QpGPD(1-1/250)97.5%0.04557386&#62; QpGPD(1-1/2500)97.5%0.08925095This is similar with the following outputs, with the return period of a yearly event (one observation out of 250 trading days)&#62; gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/250, ci.type =+ "likelihood", ci.p = 0.95,like.num = 50)Lower CI   Estimate   Upper CI0.04172534 0.04557655 0.05086785or the decennial one&#62; gpd.q(tailplot(gpd(X,quantile(X,.975))), 1-1/2500, ci.type =+ "likelihood", ci.p = 0.95,like.num = 50)Lower CI   Estimate   Upper CI0.07165395 0.08925558 0.13636620Note that it is also possible to derive an estimator for another population risk measure (the quantile is simply the so-called Value-at-Risk), the expected shortfall (or Tail Value-at-Risk), i.e.The idea is to write that expressionso that we recognize the mean excess function (discussed earlier). Thus, assuming again that above  (and therefore above that high quantile) a GPD will fit, we can writeor equivalentlyIf we substitute estimators to unknown quantities on that expression, we getThe code is here&#62; EpGPD=function(p){+ u-beta/xi+beta/xi/(1-xi)*(100/2.5*(1-p))^(-xi)+ }&#62; EpGPD(1-1/250)97.5%0.06426508&#62; EpGPD(1-1/2500)97.5%0.1215077An alternative is to use Hill's approach (used to derive Hill's estimator). Assume here that , where  is a slowly varying function. Then, for all ,Since  is a slowly varying function, it seem natural to assume that this ratio is almost 1 (which is true asymptotically). Thusi.e. if we invert that function, we derive an estimator for the quantile functionwhich can also be written(which is close to the relation we derived using a GPD model). Here the code is&#62; k=trunc(length(X)*.025)&#62; Xs=rev(sort(as.numeric(X)))&#62; xiHill=mean(log(Xs[1:k]))-log(Xs[k+1])&#62; u=Xs[k+1]&#62; QpHill=function(p){+ u+u*((100/2.5*(1-p))^(-xiHill)-1)+ }with the following Hill plotFor yearly and decennial events, we have here&#62; QpHill(1-1/250)[1] 0.04580548&#62; QpHill(1-1/2500)[1] 0.1010204Those quantities seem consistent since they are quite close, but they are different compared with empirical quantiles,&#62; quantile(X,1-1/250)99.6%0.04743929&#62; quantile(X,1-1/2500)99.96%0.09054039Note that it is also possible to use some functions to derive estimators of those quantities,&#62; riskmeasures(gpd(X,quantile(X,.975)),1-1/250)p   quantile      sfall[1,] 0.996 0.04557655 0.06426859&#62; riskmeasures(gpd(X,quantile(X,.975)),1-1/2500)p   quantile     sfall[1,] 0.9996 0.08925558 0.1215137(in this application, we have assumed that log-returns were independent and identically distributed... which might be a rather strong assumption).