We’ve seen Monday, in the MAT8595 course how to use the Generalized Pareto Distribution to estimate some downside risk measures, given a sample (assumed to be i.i.d., I will not mention here properties on extremes for stochastic processes) with distribution
. The cumulative distribution function of the Pareto distribution is here
}(x)%20=%20\begin{cases}\displaystyle{1%20-%20\left(1+%20\frac{\xi%20x}{\sigma}\right)^{-1/\xi}}%20&%20\text{%20}\xi%20\neq%200%20\\\displaystyle1%20-%20\exp%20\left(-\frac{x}{\sigma}\right)%20&%20\text{%20}\xi%20=%200\end{cases})
For some threshold
, and
, we can write
=F(u)+[1-F(u)]\cdot%20\mathbb{P}(X%C2%AD\leq%20x\vert%20X%C2%AD%3Eu))
From Pickands–Balkema–de Haan theorem, if
is large enough, then
\sim%20F(u)+[1-F(u)]\cdot%20G_{\xi,\sigma}(x-u))
Given our sample
, let
denote the number of observations over, threshold
. Then we can write
=\frac{N_u}{n}\left(1-\left(1+\frac{\widehat%20\xi}{\widehat%20\sigma}[x-u]\right)^{-\frac{1}{\widehat%20\xi}}\right)+\left(1-\frac{N_u}{n}\right))
or equivalently
=1-\frac{N_u}{n}\left(1+\frac{\widehat%20\xi}{\widehat%20\sigma}[x-u]\right)^{-\frac{1}{\widehat%20\xi}}\right))
If we invert this function, we get the quantile of level
,
=u+\frac{\widehat%20\sigma}{\widehat%20\xi}\left(\left(\frac{n(1-p)}{N_u}\right)^{-\widehat%20\xi}-1\right))
Actually, a threshold and then the implied number of observation exceeding that threshold, it is possible to consider a fixed number of observation, and then the associated threshold will be the associated order statistics.
The density of the Pareto distribution is here

which is here function of two paramters,
and
.As discussed in the course, it is possible to use the Delta method to derive the asymptotic distribution of any quantile, and get then an approximated (asymptotic) confidence interval.
But since
is usually not a parameter of interest, why not considering a reparametrization of our density, as a function of
and
(for some probability
that will be considered as fixed from now on). We can easily get (assuming that
) that
![https://latex.codecogs.com/gif.latex?g_{\xi,Q(p)}(x)=\frac{\displaystyle{\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{\xi[Q(p)-u]}\left(1+\frac{\displaystyle{\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{[Q(p)-u]}\cdot%20x\right)^{-\frac{1}{\xi}-1}](https://latex.codecogs.com/gif.latex?g_%7B%5Cxi,Q%28p%29%7D%28x%29=%5Cfrac%7B%5Cdisplaystyle%7B%5Cleft%28%5Cfrac%7Bn%7D%7BN_u%7D%281-p%29%5Cright%29%5E%7B-%5Cxi%7D-1%7D%7D%7B%5Cxi[Q%28p%29-u]%7D%5Cleft%281+%5Cfrac%7B%5Cdisplaystyle%7B%5Cleft%28%5Cfrac%7Bn%7D%7BN_u%7D%281-p%29%5Cright%29%5E%7B-%5Cxi%7D-1%7D%7D%7B[Q%28p%29-u]%7D%5Ccdot%20x%5Cright%29%5E%7B-%5Cfrac%7B1%7D%7B%5Cxi%7D-1%7D)
Tis expression is simple, and can be used to derive the likelihood (on the observations exceeding the threshold)
Numerically, let us write (and plot) that function. Consider some real data here
> X=as.numeric(danish)
> Xs=sort(X,decreasing=TRUE)
> n=length(X)
> u=10
> nu=sum(X>u)
Consider, say, the 99.9% quantile,
> p=.999
The empirical quantile is here
> quantile(X,p)
99.9%
131.5519
The density and the loglikelihood functions are here
> gq=function(x,xi,q){
+ ( (n/nu*(1-p) ) ^ (-xi)-1)/(xi*(q-u))*
+ (1+((n/nu*(1-p))^(-xi)-1)/(q-u)*x)^(-1/xi-1)}
> loglik=function(param){
+ xi=param[2];q=param[1]
+ lg=function(i) log(gq(Xs[i],xi,q))
+ return(-sum(Vectorize(lg)(1:nu)))
+ }
We can try to plot this likelihood using
> h=201
> Q=seq(50,300,length=h)
> XI=seq(.1,1,length=h)
> XIQ=as.matrix(expand.grid(Q,XI))
> M=mapply(loglik,XIQ)
Unfortunately, it was not working, so I used the old style
> M=matrix(NA,h,h)
> for(i in 1:h){for(j in 1:h){M[i,j]=loglik(c(Q[i],XI[j]))}}
The level curves of the log-likelihood are here
> hc=heat.colors(100)
> image(Q,XI,-M,col=hc)
> contour(Q,XI,-M,add=TRUE)
Again, since our interest is in the quantile, we can draw the profile likelihood and get the maximum of that function
> PL=function(Q){
+ profilelikelihood=function(xi){
+ loglik(c(Q,xi))}
+ return(optim(par=.8,fn=profilelikelihood)$value)}
> (OPT=optimize(f=PL,interval=c(100,500)))
$minimum
[1] 111.1055
and the graph is
> XQ=seq(50,300,length=101)
> L=Vectorize(PL)(XQ)
> plot(XQ,-L,type="l")
> up=OPT$objective
> abline(h=-up)
> abline(h=-up-qchisq(p=.95,df=1),col="red")
> I=which(-L>=-up-qchisq(p=.95,df=1))
> lines(XQ[I],rep(-up-qchisq(p=.95,df=1),length(I)),
+ lwd=5,col="red")
> abline(v=range(XQ[I]),lty=2,col="red")

which can be seen as an alternative to
> gpd.q(tailplot(gpd(X,u)),.999)
Lower CI Estimate Upper CI
64.66184 94.28956 188.91752
$objective
[1] 454.6481

If we want to focus on another downside risk measure, that shouldn’t be too difficult. For instance, the expected shortfall,
can be estimated as
=\widehat{Q}(p)+\underbrace{\mathbb{E}\left(X-\widehat{Q}(p)\vert%20X%3E%20\widehat{Q}(p)\right)}_{e(\widehat{Q}(p))})
where
denotes the mean excess function, which can be writen, with a Generalized Pareto Distribution
=\frac{\sigma+\xi%20x}{1-\xi})
Thus, a natural estimator for the expected shortfall is
=\frac{\widehat{Q}(p)}{1-\widehat%20\xi}+\frac{\widehat%20\sigma-\widehat%20\xi%20u}{1-\widehat%20\xi})
One more time, it is possible to re-parametrize the density of the Pareto distribution, using
instead of
. Here, we get
![https://latex.codecogs.com/gif.latex?g_{\xi,ES(p)}(x)=\frac{\displaystyle{\xi+\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{\xi(1-\xi)[ES(p)-u]}\left(1+\frac{\displaystyle{\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{(1-\xi)[ES(p)-u]}\cdot%20x\right)^{-\frac{1}{\xi}-1}](https://latex.codecogs.com/gif.latex?g_%7B%5Cxi,ES%28p%29%7D%28x%29=%5Cfrac%7B%5Cdisplaystyle%7B%5Cxi+%5Cleft%28%5Cfrac%7Bn%7D%7BN_u%7D%281-p%29%5Cright%29%5E%7B-%5Cxi%7D-1%7D%7D%7B%5Cxi%281-%5Cxi%29[ES%28p%29-u]%7D%5Cleft%281+%5Cfrac%7B%5Cdisplaystyle%7B%5Cleft%28%5Cfrac%7Bn%7D%7BN_u%7D%281-p%29%5Cright%29%5E%7B-%5Cxi%7D-1%7D%7D%7B%281-%5Cxi%29[ES%28p%29-u]%7D%5Ccdot%20x%5Cright%29%5E%7B-%5Cfrac%7B1%7D%7B%5Cxi%7D-1%7D)
The code to get the associated log-likelihood is here
> ge=function(x,xi,es){
+ (xi+(n/nu*(1-p))^(-xi)-1)/(xi*(1-xi)*(es-u))*(1+(xi+(n/nu*(1-p))^(-xi)
+ -1)/((es-u)*(1-xi))*x)^(-1/xi-1)
+ }
> loglik=function(param){
+ xi=param[2];es=param[1]
+ lg=function(i) log(ge(Xs[i],xi,es))
+ return(-sum(Vectorize(lg)(1:nu)))
+ }
and again, we can plot it

and the profile (log) likelihood is here (for the 99.9% expected shortfall)
> PL=function(ES){
+ profilelikelihood=function(xi){
+ loglik(c(ES,xi))}
+ return(optim(par=.8,fn=profilelikelihood)$value)}
> (OPT=optimize(f=PL,interval=c(100,500)))
$minimum
[1] 143.66
$objective
[1] 454.6481

which could be compared with
> gpd.sfall(tailplot(gpd(X,u)),.999)
Lower CI Estimate Upper CI
96.64625 191.36972 394.87555
