{"word_count": 492, "character_count": 4476}
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 hereFor some threshold , and , we can writeFrom Pickands–Balkema–de Haan theorem, if  is large enough, thenGiven our sample , let  denote the number of observations over,  threshold . Then we can writeor equivalentlyIf we invert this function, we get the quantile of level ,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 herewhich 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 ) thatTis 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&#62; X=as.numeric(danish)&#62; Xs=sort(X,decreasing=TRUE)&#62; n=length(X)&#62; u=10&#62; nu=sum(X&#62;u)Consider, say, the 99.9% quantile,&#62; p=.999The empirical quantile is here&#62; quantile(X,p)   99.9% 131.5519The density and the loglikelihood functions are here&#62; 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)}&#62; 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&#62; h=201&#62; Q=seq(50,300,length=h)&#62; XI=seq(.1,1,length=h)&#62; XIQ=as.matrix(expand.grid(Q,XI))&#62; M=mapply(loglik,XIQ)Unfortunately, it was not working, so I used the old style&#62; M=matrix(NA,h,h)&#62; 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&#62; hc=heat.colors(100)&#62; image(Q,XI,-M,col=hc)&#62; 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&#62; PL=function(Q){+ profilelikelihood=function(xi){+ loglik(c(Q,xi))}+ return(optim(par=.8,fn=profilelikelihood)$value)}&#62; (OPT=optimize(f=PL,interval=c(100,500)))$minimum[1] 111.1055and the graph is&#62; XQ=seq(50,300,length=101)&#62; L=Vectorize(PL)(XQ)&#62; plot(XQ,-L,type="l")&#62; up=OPT$objective&#62; abline(h=-up)&#62; abline(h=-up-qchisq(p=.95,df=1),col="red")&#62; I=which(-L&#62;=-up-qchisq(p=.95,df=1))&#62; lines(XQ[I],rep(-up-qchisq(p=.95,df=1),length(I)),+ lwd=5,col="red")&#62; abline(v=range(XQ[I]),lty=2,col="red")which can be seen as an alternative to&#62; gpd.q(tailplot(gpd(X,u)),.999) Lower CI  Estimate  Upper CI  64.66184  94.28956 188.91752 $objective[1] 454.6481If we want to focus on another downside risk measure, that shouldn't be too difficult. For instance, the expected shortfall,  can be estimated aswhere  denotes the mean excess function, which can be writen, with a Generalized Pareto DistributionThus, a natural estimator for the expected shortfall isOne more time, it is possible to re-parametrize the density of the Pareto distribution, using  instead of . Here, we getThe code to get the associated log-likelihood is here&#62; 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)+ }&#62; 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 itand the profile (log) likelihood is here (for the 99.9% expected shortfall)&#62; PL=function(ES){+ profilelikelihood=function(xi){+ loglik(c(ES,xi))}+ return(optim(par=.8,fn=profilelikelihood)$value)}&#62; (OPT=optimize(f=PL,interval=c(100,500)))$minimum[1] 143.66$objective[1] 454.6481which could be compared with&#62; gpd.sfall(tailplot(gpd(X,u)),.999) Lower CI  Estimate  Upper CI  96.64625 191.36972 394.87555