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 written
The 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)
> library(tseries) > X=get.hist.quote("^GSPC") > T=time(X) > D=as.POSIXlt(T) > Y=X$Close > R=diff(log(Y)) > D=D[-1] > X=-R > plot(D,X) > library(evir) > GPD=gpd(X,quantile(X,.975)) > xi=GPD$par.ests[1] > beta=GPD$par.ests[2] > u=GPD$threshold > QpGPD=function(p){ + u+beta/xi*((100/2.5*(1-p))^(-xi)-1) + } > QpGPD(1-1/250) 97.5% 0.04557386 > QpGPD(1-1/2500) 97.5% 0.08925095
This is similar with the following outputs, with the return period of a yearly event (one observation out of 250 trading days)
> 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 CI 0.04172534 0.04557655 0.05086785
or the decennial one
> 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 CI 0.07165395 0.08925558 0.13636620
Note 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 expression
so 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 write
or equivalently
If we substitute estimators to unknown quantities on that expression, we get
The code is here
> EpGPD=function(p){ + u-beta/xi+beta/xi/(1-xi)*(100/2.5*(1-p))^(-xi) + } > EpGPD(1-1/250) 97.5% 0.06426508 > EpGPD(1-1/2500) 97.5% 0.1215077
An 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). Thus
i.e. if we invert that function, we derive an estimator for the quantile function
which can also be written
(which is close to the relation we derived using a GPD model). Here the code is
> k=trunc(length(X)*.025) > Xs=rev(sort(as.numeric(X))) > xiHill=mean(log(Xs[1:k]))-log(Xs[k+1]) > u=Xs[k+1] > QpHill=function(p){ + u+u*((100/2.5*(1-p))^(-xiHill)-1) + }
with the following Hill plot
For yearly and decennial events, we have here
> QpHill(1-1/250) [1] 0.04580548 > QpHill(1-1/2500) [1] 0.1010204
Those quantities seem consistent since they are quite close, but they are different compared with empirical quantiles,
> quantile(X,1-1/250) 99.6% 0.04743929 > quantile(X,1-1/2500) 99.96% 0.09054039
Note that it is also possible to use some functions to derive estimators of those quantities,
> riskmeasures(gpd(X,quantile(X,.975)),1-1/250) p quantile sfall [1,] 0.996 0.04557655 0.06426859 > 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).