Risk Measures with Extreme Value Models

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 $F$. The cumulative distribution function of the  Pareto distribution is here

$G_{(\xi,\sigma)}(x) = \begin{cases}\displaystyle{1 - \left(1+ \frac{\xi x}{\sigma}\right)^{-1/\xi}} & \text{ }\xi \neq 0 \\\displaystyle1 - \exp \left(-\frac{x}{\sigma}\right) & \text{ }\xi = 0\end{cases}$

For some threshold $u$, and $x\geq u$, we can write

$F(x)=F(u)+[1-F(u)]\cdot \mathbb{P}(X­\leq x\vert X­>u)$

From Pickands–Balkema–de Haan theorem, if $u$ is large enough, then

$F(x)\sim F(u)+[1-F(u)]\cdot G_{\xi,\sigma}(x-u)$

Given our sample $\{x_1,\cdots,x_n\}$, let $N_u$ denote the number of observations over,  threshold $u$. Then we can write

$\widehat{F}(x)=\frac{N_u}{n}\left(1-\left(1+\frac{\widehat \xi}{\widehat \sigma}[x-u]\right)^{-\frac{1}{\widehat \xi}}\right)+\left(1-\frac{N_u}{n}\right)$

or equivalently

$\widehat{F}(x)=1-\frac{N_u}{n}\left(1+\frac{\widehat \xi}{\widehat \sigma}[x-u]\right)^{-\frac{1}{\widehat \xi}}\right)$

If we invert this function, we get the quantile of level $p\in(0,1)$,

$\widehat{Q}(p)=u+\frac{\widehat \sigma}{\widehat \xi}\left(\left(\frac{n(1-p)}{N_u}\right)^{-\widehat \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

$g_{(\xi,\sigma)}(x) = \frac{1}{\sigma}\left(1 + \frac{\xi x}{\sigma}\right)^{\left(-\frac{1}{\xi} - 1\right)}$

which is here function of two paramters, $\xi$ and $\sigma$.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 $\sigma$ is usually not a parameter of interest, why not considering a reparametrization of our density, as a function of  $\xi$ and $Q(p)$ (for some probability $p$ that will be considered as fixed from now on). We can easily get (assuming that $\xi\neq 0$) that

$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 x\right)^{-\frac{1}{\xi}-1}$

Tis expression is simple, and can be used to derive the likelihood (on the observations exceeding the threshold)

$\log\mathcal{L}(\xi,Q(p);\boldsymbol{x})=\sum_{i=0}^{N_u-1} \log g_{\xi,Q(p)}(x_{n-i:n})$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, $ES(p)=\mathbb{E}(X\vert X>Q(p))$ can be estimated as

$\widehat{ES}(p)=\widehat{Q}(p)+\underbrace{\mathbb{E}\left(X-\widehat{Q}(p)\vert X> \widehat{Q}(p)\right)}_{e(\widehat{Q}(p))}$

where $e(\cdot )$ denotes the mean excess function, which can be writen, with a Generalized Pareto Distribution

$e(x)=\frac{\sigma+\xi x}{1-\xi}$

Thus, a natural estimator for the expected shortfall is

$\widehat{ES}(p)=\frac{\widehat{Q}(p)}{1-\widehat \xi}+\frac{\widehat \sigma-\widehat \xi u}{1-\widehat \xi}$

One more time, it is possible to re-parametrize the density of the Pareto distribution, using $ES(p)$ instead of $\sigma$. Here, we get

$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 x\right)^{-\frac{1}{\xi}-1}$

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