# Profile Likelihood

Consider some simulated data

> set.seed(1)
> x=exp(rnorm(100))

Assume that those data are observed i.id. random variables with distribution$F_{\boldsymbol{\theta}}$, with $\boldsymbol{\theta}=(\alpha,\beta)$. The natural idea is to consider the maximum likelihood estimator

$\widehat{\boldsymbol{\theta}}=\underset{\boldsymbol{\theta}}{\text{argmax}}\left\lbrace\log \mathcal{L}(\boldsymbol{\theta})\right\rbrace$For instance, consider some maximum likelihood estimator,

> library(MASS)
> (F=fitdistr(x,"gamma"))
shape       rate
1.4214497   0.8619969
(0.1822570) (0.1320717)
> F$estimate[1]+c(-1,1)*1.96*F$sd[1]
[1] 1.064226 1.778673

Here, we have an approximated (since the maximum likelihood has an asymptotic Gaussian distribution) confidence interval for $\alpha$. We can use numerical optimization routine to get the maximum of the log-likelihood function

> log_lik=function(theta){
+   a=theta[1]
+   b=theta[2]
+   logL=sum(log(dgamma(x,a,b)))
+   return(-logL)
+ }

> optim(c(1,1),log_lik)
$par [1] 1.4214116 0.8620311$value
[1] 146.5909

And we have the same value.

Now, what if we care only about $\alpha$, and not $\beta$. The we can use profile likelihood. The idea is to solve

$\widehat{\alpha}=\underset{\alpha}{\text{argmax}}\left\lbrace \underset{\beta}{\text{max}}\left\lbrace\log \mathcal{L}(\alpha,\beta)\right\rbrace \right\rbrace$

i.e.

$\widehat{\alpha}=\underset{\alpha}{\text{argmax}}\left\lbrace \log \mathcal{L}_p(\alpha)\right\rbrace$or, equivalently,

$\widehat{\alpha}=\underset{\alpha}{\text{argmax}}\left\lbrace \log \mathcal{L}(\alpha,\widehat{\beta}_{\alpha})\right\rbrace$

> prof_log_lik=function(a){
+   b=(optim(1,function(z) -sum(log(dgamma(x,a,z)))))$par + return(-sum(log(dgamma(x,a,b)))) + } > vx=seq(.5,3,length=101) > vl=-Vectorize(prof_log_lik)(vx) > plot(vx,vl,type="l") > optim(1,prof_log_lik)$par
[1] 1.421094

$value [1] 146.5909 A few weeks ago, we have mentioned the likelihood ratio test, i.e. $2[\log\mathcal{L}(\widehat{\theta})-\log\mathcal{L}(\theta_0)]\sim n[\widehat{\theta}-{\theta}_0]^2I(\theta_0)\sim \chi^2(k)$ The analogous can be obtained here, since $2[\log\mathcal{L}_p(\widehat{\alpha})-\log\mathcal{L}_p(\alpha_0)]\sim\chi^2(1)$ (the 1 comes from the fact that $\alpha$ is a one-dimensional coefficient). The (technical) proof can be found in Suhasini Subba Rao’s notes (see also Section 4.5.2 in Antony Davison’s Statistical Models). From that property, we can easily obtain a confidence interval for $\alpha$ Hence, from our sample, we get the following 95% confidence interval, > abline(v=optim(1,prof_log_lik)$par,lty=2)
> abline(h=-optim(1,prof_log_lik)$value) > abline(h=-optim(1,prof_log_lik)$value-qchisq(.95,1)/2)

> segments(F$estimate[1]-1.96*F$sd[1],
-170,F$estimate[1]+1.96*F$sd[1],-170,lwd=3,col="blue")
> borne=-optim(1,prof_log_lik)$value-qchisq(.95,1)/2 > (b1=uniroot(function(z) Vectorize(prof_log_lik)(z)+borne,c(.5,1.5))$root)
[1] 1.095726
> (b2=uniroot(function(z) Vectorize(prof_log_lik)(z)+borne,c(1.25,2.5))$root) [1] 1.811809 that can be visualized below, > segments(b1,-168,b2,-168,lwd=3,col="red") In blue the obtained obtained using the asymptotic Gaussian property of the maximum likelihood estimator, and in red, the obtained obtained using the asymptotic chi-square distribution of the log (profile) likelihood ratio. Cite this article as: Arthur Charpentier, "Profile Likelihood," in Freakonometrics, 16/11/2015, http://freakonometrics.hypotheses.org/20573. # Where does that 2 come from in the likelihood ratio test? This afternoon, in class, we’ve seen Wald test, the likelihood-ratio test, and finally the score test. All of them rely on the same idea and then, use that if with , we can write Or – slightly more interesting – if , then Then one can get that Based on that property, we can derive Wald statistics, that can be visualized below The score test is a test on the square of the slope The idea for the likelihood ratio test is to consider Observe that $\log\mathcal{L}(\widehat{\theta})-\log\mathcal{L}({\theta}_0)$ can be written, using Taylor’s expansion $[\widehat{\theta}-{\theta}_0]\frac{\partial \log\mathcal{L}(\widehat{\theta})}{\partial \theta}+\frac{1}{2}[\widehat{\theta}-{\theta}_0]^2\frac{\partial^2 \log\mathcal{L}(\widetilde{\theta})}{\partial \theta^2}$ for some $\widetilde{\theta}\in[\widehat{\theta},\theta_0]$. The first term is null, since the maximum likelihood estimator is precisely at the maximum of the (log) likelihood. So $2[\log\mathcal{L}(\widehat{\theta})-\log\mathcal{L}(\theta_0)]=[\widehat{\theta}-{\theta}_0]^2\frac{\partial^2 \log\mathcal{L}(\widetilde{\theta})}{\partial \theta^2}$ That’s more or less where the 2 comes from. Then observe that $\frac{\partial^2 \log\mathcal{L}(\widetilde{\theta})}{\partial \theta^2}\sim nI(\theta_0)$ and therefore $2[\log\mathcal{L}(\widehat{\theta})-\log\mathcal{L}(\theta_0)]\sim n[\widehat{\theta}-{\theta}_0]^2I(\theta_0)\sim \chi^2(k)$ This test will be discussed further next week (since it is related to Neyman-Pearson’s theorem), but also, that result can be used to derive confidence intervals. With a log-likelihood as follows it is possible to get a confidence interval for the parameter by looking for$\theta$‘s such that $\log\mathcal{L}(\theta)\geq\log\mathcal{L}(\widehat{\theta})-q$ We will discuss that idea later on, in the context of profile likelihood. # 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

# Likelihood Based Methods, for Extremes

This week, in the MAT8595 course, we will start the section on inference for extreme values. To start with something simple, we will use maximum likelihood techniques on a Generalized Pareto Distribution (we’ve seen Monday Pickands-Balkema-de Hann theorem).

• Maximum Likelihood Estimation

In the context of parametric models, the standard technique is to consider the maximum of the likelihood (or the log-likelihod).Let $\boldsymbol{\theta}$ denote the parameter (with $\boldsymbol{\theta}\in\Theta\subset\mathbb{R}^k$). Given some – stnardard – technical assumptions, such as $\boldsymbol{\theta}_0\notin\partial\Theta$, or $f(x;\boldsymbol{\theta})­>0$ on some neighbourhood of $\boldsymbol{\theta}_0$, then

$\sqrt{n}\big(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0\big)\ \xrightarrow{\mathcal{L}}\ \mathcal{N}(0,\,I^{-1})$

where $I$ denotes Fisher information matrix (see any textbook for mathematical statistics courses). Consider here some i.i.d. sample, from a Generalized Pareto Distribution, with parameter $\boldsymbol{\theta}=(\xi,\sigma)$, so that

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

If we solve (numerically) the first order condition of the maximum likelihood, we get an estimator  $\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma }_n)$ which satisfies

$\sqrt{n}\left(\left[\begin{array}{c}\widehat{\xi}_n\\\widehat{\sigma }_n\end{array}\right]-\left[\begin{array}{c}\xi_0\\\sigma_0 \end{array}\right]\right)\rightarrow \mathcal{N}\left(\left[\begin{array}{c}0\\0\end{array}\right],\left[\begin{array}{cc}(1+\xi_0)^2 & \sigma_0[1+\xi_0]\\ \sigma_0 [1+\xi_0] & 2\sigma^2_0(1+\xi_0) \end{array}\right]\right)$

The idea of this asymptotic normality is the following : if the true distribution of the sample is a GPD with parameter $\boldsymbol{\theta}_0$, then, if $n$ is large enough, then $\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma }_n)$ will have a joint normal distribution. So if we generate a lot of sample (sufficently large, say 200 observations), then the scatterplot of the estimator should the same as the scatterplot of a Gaussian distribution,

> library(evir)
> n=200
> param=matrix(NA,1000,2)
> for(s in 1:1000){
+ x=rgpd(n,xi=1/1.5,beta=1)
+ param[s,]=gpd(x,0)$par.ests + } > m=apply(param,2,mean) > S=var(param) > library(mnormt) > x=seq(min(param[,1])-.05,max(param[,1])+.05,length=101) > y=seq(min(param[,2])-.05,max(param[,2])+.05,length=101) > vx=rep(x,each=length(y)) > vy=rep(y,length(x)) > vz=dmnorm(cbind(vx,vy),m,S) > z=matrix(vz,length(y),length(x)) > COL=rev(heat.colors(100)) > image(x,y,z,col=COL) > points(param) and to get a 3d representation > x=seq(min(param[,1])-.05,max(param[,1])+.05,length=31) > y=seq(min(param[,2])-.05,max(param[,2])+.05,length=31) > vx=rep(x,each=length(y)) > vy=rep(y,length(x)) > vz=dmnorm(cbind(vx,vy),m,S) > z=matrix(vz,length(y),length(x)) > persp(x,y,t(z),shade=TRUE,col="green",theta=-30,phi=20,ticktype="detailed", + xlab="xi",ylab="sigma") With 200 observations, if the true underlying distribution is a GPD, then, indeed, the joint distribution of $\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma }_n)$ seems to be normal. That would be interesting to generate some confidence intervals for instance, or define some tests. To go further, see any standard textbook on statistical mathematics, e.g. Casella & Berger (2002). • Delta Method Another important property is the so called delta-method (we’ve seen Monday in class that it was obtained easily using a first order Taylor expansion). The idea is that if $\widehat{\boldsymbol{\theta}}_n$ is asymptotically normal, and if $h:\Theta\subset\mathbb R^k \rightarrow \mathbb R^m$ is sufficently smooth, then $h(\widehat{\boldsymbol{\theta}}_n)$ will also be asymptotically Gaussian. More precicely (see also the header of this blog) $\sqrt{n}\big(h(\widehat{\boldsymbol{\theta}}_n)-h(\boldsymbol{\theta}_0)\big)\ \xrightarrow{\mathcal{L}}\ \mathcal{N}(0,\,\nabla h((\boldsymbol{\theta}_0)'\cdot I^{-1}\cdot \nabla h((\boldsymbol{\theta}_0))$ From this property, we can get the normality of $\widehat{\alpha}_n=\widehat{\xi}_n^{-1}$ (which is another parametrization used in extreme value models), or on any quantile, $\widehat{Q}_u=F^{-1}_{\widehat{\boldsymbol{\theta}}_n}(u)=h_u(\widehat{\xi}_n,\widehat{\sigma}_n)$. Let us run some simulation, one more time to check that we actually have a joint normality. > library(evir) > n=200 > param=riskm=matrix(NA,1000,2) > for(s in 1:1000){ + x=rgpd(n,xi=1/1.5,beta=1) + param[s,]=gpd(x,0)$par.ests
+ xihat=param[s,1]
+ shat=param[s,2]
+ q=shat * (.01^(-xihat) - 1)/xihat
+ tvar=q+(shat + xihat * q)/(1 - xihat)
+ riskm[s,]=c(1/xihat,q)
+ }
> m=apply(riskm,2,mean)
> S=var(riskm)
> library(mnormt)
> x=seq(min(riskm[,1])-.05,max(riskm[,1])+.05,length=101)
> y=seq(min(riskm[,2])-.05,max(riskm[,2])+.05,length=101)
> vx=rep(x,each=length(y))
> vy=rep(y,length(x))
> vz=dmnorm(cbind(vx,vy),m,S)
> z=matrix(vz,length(y),length(x))
> image(x,y,t(z),col=COL)
> points(riskm)

As we can see bellow, with samples of size 200, we cannot use this asymptotical result: it looks like we do not have enought data. But if we run the same code with

> n=5000

We get the joint normality of $\widehat{\alpha}_n$ and $\widehat{Q}_n(u)$. This is what we can get from this result, called delta-method in statistical textbooks. See again Casella & Berger (2002) for more details.

• Profile Likelihood

Another interesting tool is the concept of profile likelihood. This would be interesting here since the main interest is the tail index $\xi$, $\sigma$ being here some kind of auxilary parameter. See Venzon & Moolgavkar (1988) for more details. Here, we will plot

But more generally, it is possible to consider

where  is the set of interesting parameters. Then (under standard suitable conditions) we can prove that

which can be used to derive confidence intervals. In the GPD case, for each $\xi$, we have to find an optimal $\sigma^\star(\xi)$. We compute the (profile) likelihood i.e. $\mathcal{L}(\xi,\sigma^\star(\xi))$. And we can compute the maximum of this profile likelihood. This two-stage optimization is, in general, not equivalent with the (global) maximization of the likelihood, as computed below

>  n=500
>  set.seed(1)
>  x=rgpd(n,xi=1/1.5,beta=1)
>  loglikelihood=function(xi,beta){
+  sum(log(dgpd(x,xi,mu=0,beta))) }
>  XIV=(1:300)/100;L=rep(NA,300)
>  for(i in 1:300){
+  XI=XIV[i]
+  profilelikelihood=function(beta){
+  -loglikelihood(XI,beta) }
+  L[i]=-optim(par=1,fn=profilelikelihood)$value } > plot(XIV,L,type="l") > XIV[which.max(L)] [1] 0.67 > gpd(x,0)$par.ests
xi      beta
0.6730145 0.9725483

We are not far away. Actually, if we want to compute the maximum of the profile likelihood (and not only compute the values of the profile likelihood on a grid, as before), we use

>  PL=function(XI){
+  profilelikelihood=function(beta){
+  -loglikelihood(XI,beta) }
+  return(optim(par=1,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(0,3)))$minimum
[1] 0.6731025

$objective [1] 822.5574 Observe that, indeed, we are not far away from the maximum likelihood estimator of $\xi$ (I believe that it’s mainly a computational issue here, and theat the two are similar, here… actually, I’d be glad to hear about cases where maximum of the profile likelihood is not the same as the maximum of the likelihood). The interesting point is that we can use this technique to compute a confidence interval, and even visualize it on a graph > 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(XIV[I],rep(-up-qchisq(p=.95,df=1),length(I)),
+  lwd=5,col="red")
>  abline(v=range(XIV[I]),lty=2,col="red")

The vertical lines are the lower and the upper bound of a 95% confidence interval for parameter $\xi$.

To go further, see Murphy, S.A & van der Vaart, A.W. (2000). On Profile Likelihood.

# Maximum likelihood estimates for multivariate distributions

Consider our loss-ALAE dataset, and – as in Frees & Valdez (1998) – let us fit a parametric model, in order to price a reinsurance treaty. The dataset is the following,

> library(evd)
> data(lossalae)
> Z=lossalae
> X=Z[,1];Y=Z[,2]

The first step can be to estimate marginal distributions, independently. Here, we consider lognormal distributions for both components,

> Fempx=function(x) mean(X<=x)
> Fx=Vectorize(Fempx)
> u=exp(seq(2,15,by=.05))
> plot(u,Fx(u),log="x",type="l",
+ xlab="loss (log scale)")
> Lx=function(px) -sum(log(Vectorize(dlnorm)(
+ X,px[1],px[2])))
> opx=optim(c(1,5),fn=Lx)
> opx$par [1] 9.373679 1.637499 > lines(u,Vectorize(plnorm)(u,opx$par[1],
+ opx$par[2]),col="red") The fit here is quite good, For the second component, we do the same, > Fempy=function(x) mean(Y<=x) > Fy=Vectorize(Fempy) > u=exp(seq(2,15,by=.05)) > plot(u,Fy(u),log="x",type="l", + xlab="ALAE (log scale)") > Ly=function(px) -sum(log(Vectorize(dlnorm)( + Y,px[1],px[2]))) > opy=optim(c(1.5,10),fn=Ly) > opy$par
[1] 8.522452 1.429645
> lines(u,Vectorize(plnorm)(u,opy$par[1], + opy$par[2]),col="blue")

It is not as good as the fit obtained on losses, but it is not that bad,

Now, consider a multivariate model, with Gumbel copula. We’ve seen before that it worked well. But this time, consider the maximum likelihood estimator globally.

> Cop=function(u,v,a) exp(-((-log(u))^a+
+ (-log(v))^a)^(1/a))
> phi=function(t,a) (-log(t))^a
> cop=function(u,v,a) Cop(u,v,a)*(phi(u,a)+
+ phi(v,a))^(1/a-2)*(
+ a-1+(phi(u,a)+phi(v,a))^(1/a))*(phi(u,a-1)*
+ phi(v,a-1))/(u*v)
> L=function(p) {-sum(log(Vectorize(dlnorm)(
+ X,p[1],p[2])))-
+ sum(log(Vectorize(dlnorm)(Y,p[3],p[4])))-
+ sum(log(Vectorize(cop)(plnorm(X,p[1],p[2]),
+ plnorm(Y,p[3],p[4]),p[5])))}
> opz=optim(c(1.5,10,1.5,10,1.5),fn=L)
> opz$par [1] 9.377219 1.671410 8.524221 1.428552 1.468238 Marginal parameters are (slightly) different from the one obtained independently, > c(opx$par,opy$par) [1] 9.373679 1.637499 8.522452 1.429645 > opz$par[1:4]
[1] 9.377219 1.671410 8.524221 1.428552

And the parameter of Gumbel copula is close to the one obtained with heuristic methods in class.

Now that we have a model, let us play with it, to price a reinsurance treaty. But first, let us see how to generate Gumbel copula… One idea can be to use the frailty approach, based on a stable frailty. And we can use Chambers et al (1976)to generate a stable distribution. So here is the algorithm to generate samples from Gumbel copula

> alpha=opz$par[5] > invphi=function(t,a) exp(-t^(1/a)) > n=500 > x=matrix(rexp(2*n),n,2) > angle=runif(n,0,pi) > E=rexp(n) > beta=1/alpha > stable=sin((1-beta)*angle)^((1-beta)/beta)* + (sin(beta*angle))/(sin(angle))^(1/beta)/ + (E^(alpha-1)) > U=invphi(x/stable,alpha) > plot(U) Here, we consider only 500 simulations, Based on that copula simulation, we can then use marginal transformations to generate a pair, losses and allocated expenses, > Xloss=qlnorm(U[,1],opz$par[1],opz$par[2]) > Xalae=qlnorm(U[,2],opz$par[3],opz$par[4]) In standard reinsurance treaties – see e.g. Clarke (1996) – allocated expenses are splited prorata capita between the insurance company, and the reinsurer. If $X$ denotes losses, and $Y$ the allocated expenses, a standard excess treaty can be has payoff $\begin{cases} 0\cdot \boldsymbol{1}(X\leq R)\\ \displaystyle{\left(X-R+ \frac{X-R}{X}\cdot Y\right) \right )}\cdot\boldsymbol{1}(R where $L$ denotes the (upper) limit, and $R$ the insurer’s retention. Using monte carlo simulation, it is then possible to estimate the pure premium of such a reinsurance treaty. > L=100000 > R=50000 > Z=((Xloss-R)+(Xloss-R)/Xloss*Xalae)* + (R<=Xloss)*(Xloss<L)+ + ((L-R)+(L-R)/R*Xalae)*(L<=Xloss) > mean(Z) [1] 12596.45 Now, play with it… it is possible to find a better fit, I guess… # Inference and autoregressive processes Consider a (stationary) autoregressive process, say of order 2, for some white noise with variance . Here is a code to generate such a process, > phi1=.5 > phi2=-.4 > sigma=1.5 > set.seed(1) > n=240 > WN=rnorm(n,sd=sigma) > Z=rep(NA,n) > Z[1:2]=rnorm(2,0,1) > for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]} Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters. • using least square regression A natural idea is to see here a regression model, and thus, if we consider a matrix formulation, Here we can run (conditional) ordinary least squares estimation, > base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)]) > regression=lm(Y~0+X1+X2,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -4.3491 -0.8890 -0.0762 0.9601 3.6105 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 0.45107 0.05924 7.615 6.34e-13 *** X2 -0.41454 0.05924 -6.998 2.67e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.449 on 236 degrees of freedom Multiple R-squared: 0.2561, Adjusted R-squared: 0.2497 F-statistic: 40.61 on 2 and 236 DF, p-value: 6.949e-16 > regression$coefficients
X1         X2
0.4510703 -0.4145365
> summary(regression)$sigma [1] 1.449276 • using Yule-Walker equations As we’ve seen in class, we can easily get the following equations for the autocovariance functions, which can also be written So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions The code is the following > rho1=cor(Z[1:(n-1)],Z[2:n]) > rho2=cor(Z[1:(n-2)],Z[3:n]) > A=matrix(c(1,rho1,rho1,1),2,2) > b=matrix(c(rho1,rho2),2,1) > (PHI=solve(A,b)) [,1] [1,] 0.4517579 [2,] -0.4155920 Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation) > estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2)
> sd(estWN)
[1] 1.445706

This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.

• using (conditional) likelihood estimators

Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e.

In that case, the conditional log likelihood (conditional since we set the first two observations here) is

> CondLogLik=function(A,TS){
+ phi1=A[1];  phi2=A[2]
+ sigma=A[3]	; L=0
+ for(t in 3:length(TS)){
+ L=L+dnorm(TS[t],mean=phi1*TS[t-1]+
+ phi2*TS[t-2],sd=sigma,log=TRUE)}
+ return(-L)}

Now, we can run standard optimization procedures,

> LogL=function(A) CondLogLik(A,TS=Z)
> optim(c(0,0,1),LogL)
$par [1] 0.4509685 -0.4144938 1.4430930$value
[1] 425.0164

$counts function gradient 88 NA$convergence
[1] 0

$message NULL Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient for the second one, and finally for the standard deviation of the innovation process All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details). # Tests on tail index for extremes Since several students got the intuition that natural catastrophes might be non-insurable (underlying distributions with infinite mean), I will post some comments on testing procedures for extreme value models. A natural idea is to use a likelihood ratio test (for composite hypotheses). Let denote the parameter (of our parametric model, e.g. the tail index), and we would like to know whether is smaller or larger than (where in the context of finite versus infinite mean ). I.e. either belongs to the set or to its complementary . Consider the maximum likelihood estimator , i.e. Let and denote the constrained maximum likelihood estimators on and respectively, Either and (on the left), or and (on the right) So likelihood ratios are either equal to or If we use the code mentioned in the post on profile likelihood, it is easy to derive that ratio. The following graph is the evolution of that ratio, based on a GPD assumption, for different thresholds, > base1=read.table( + "http://freakonometrics.free.fr/danish-univariate.txt", + header=TRUE) > library(evir) > X=base1$Loss.in.DKM
> U=seq(2,10,by=.2)
> LR=P=ES=SES=rep(NA,length(U))
> for(j in 1:length(U)){
+ u=U[j]
+ Y=X[X>u]-u
+ loglikelihood=function(xi,beta){
+ sum(log(dgpd(Y,xi,mu=0,beta))) }
+ XIV=(1:300)/100;L=rep(NA,300)
+ for(i in 1:300){
+ XI=XIV[i]
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ L[i]=-optim(par=1,fn=profilelikelihood)$value } + plot(XIV,L,type="l") + PL=function(XI){ + profilelikelihood=function(beta){ + -loglikelihood(XI,beta) } + return(optim(par=1,fn=profilelikelihood)$value)}
+ (L0=(OPT=optimize(f=PL,interval=c(0,10)))$objective) + profilelikelihood=function(beta){ + -loglikelihood(1,beta) } + (L1=optim(par=1,fn=profilelikelihood)$value)
+ LR[j]=L1-L0
+ P[j]=1-pchisq(L1-L0,df=1)
+ G=gpd(X,u)
+ ES[j]=G$par.ests[1] + SES[j]=G$par.ses[1]
+ }
>
> plot(U,LR,type="b",ylim=range(c(0,LR)))
> abline(h=qchisq(.95,1),lty=2)

with on top the values of the ratio (the dotted line is the quantile of a chi-square distribution with one degree of freedom) and below the associated p-value

> plot(U,P,type="b",ylim=range(c(0,P)))
> abline(h=.05,lty=2)

In order to compare, it is also possible to look at confidence interval for the tail index of the GPD fit,

> plot(U,ES,type="b",ylim=c(0,1))
> lines(U,ES+1.96*SES,type="h",col="red")
> abline(h=1,lty=2)

To go further, see Falk (1995), Dietrich, de Haan & Hüsler (2002), Hüsler & Li (2006) with the following table, or Neves & Fraga Alves (2008). See also here or there (for the latex based version) for an old paper I wrote on that topic.

# a short word on profile likelihood

Profile likelihood is an interesting theory to visualize and compute confidence interval for estimators (see e.g. Venzon & Moolgavkar (1988)). As we will use is, we will plot

But more generally, it is possible to consider

where . Then (under standard suitable conditions)

which can be used to derive confidence intervals.

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
> library(evir)
> X=base1$Loss.in.DKM > u=5 The function to draw the profile likelihood for the tail index parameter is then > Y=X[X>u]-u > loglikelihood=function(xi,beta){ + sum(log(dgpd(Y,xi,mu=0,beta))) } > XIV=(1:300)/100;L=rep(NA,300) > for(i in 1:300){ + XI=XIV[i] + profilelikelihood=function(beta){ + -loglikelihood(XI,beta) } + L[i]=-optim(par=1,fn=profilelikelihood)$value }
> plot(XIV,L,type="l")

It is possible to use it that profile likelihood function to derive a confidenceinterval,

> PL=function(XI){
+ profilelikelihood=function(beta){
+ -loglikelihood(XI,beta) }
+ return(optim(par=1,fn=profilelikelihood)$value)} > (OPT=optimize(f=PL,interval=c(0,3)))$minimum
[1] 0.6315989

$objective [1] 754.1115 > up=OPT$objective
> abline(h=-up)
> abline(h=-up-qchisq(p=.95,df=1)/2,col="red")
> I=which(L>=-up-qchisq(p=.95,df=1)/2)
> lines(XIV[I],rep(-up-qchisq(p=.95,df=1)/2,length(I)),
+ lwd=5,col="red")
> abline(v=range(XIV[I]),lty=2,col="red")

This is done with the following code

> library(ismev)
> gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)

# Wald, score et rapport de vraisemblance

Vendredi, nous devrions aborder un peu la problématique des tests asymptotiques. Avant, rappelons que l’estimateur du maximum de vraisemblance de  vérifie la propriété asymptotique suivante.

Or on a une propriété intéressante sur les convergences de vecteurs Gaussiens (pour rester assez général). Soit  une suite de vecteurs aléatoires telle que  où . Alors

où  (pour un paramètre univarié, ça sera alors ). Mais ce résultat n’est intéressant que si  est connue. En fait, si on a la suite des variances vérifie , alors

Traduit sur notre estimateur du maximum de vraisemblance, cela signifie que

ou encore

Trois tests peut alors être mis en œuvre (et tant qu’à faire, on peut essayer de l’illustrer sur un exemple: encore et toujours du pile/face).

> set.seed(1)
> n=20
> X=sample(0:1,size=n,replace=TRUE)
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
> LL=sapply(p,logL)
> plot(p,LL,type="l",col="red",lwd=2)
> p0=.5
> points(p0,logL(p0),pch=3,cex=1.5,lwd=2)
> abline(v=p0,lty=2)

On a alors 20 tirages de pile ou face, on a obtenu 11 piles, et on se demande si la pièce est “juste“.Mais avant toute chose, commençons par calculer l’estimateur du maximum de vraisemblance, ainsi que le score et l’information de Fisher. Dans ce modèle de variables de Bernoulli, on connaît des formes explicites de ces quantités. Mais ici, on va utiliser la méthode la plus générale qui soit, i.e. on va maximiser la log-vraisemblance, puis on va calculer le score et l’information de Fisher à la valeur de  que l’on souhaite tester.

> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=p0,method="BFGS")\$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> Df = D(f, "p")
> Df2 = D(Df, "p")
> p=p0
> score=eval(Df)
> (IF=-eval(Df2))
[1] 80
> 1/(p0*(1-p0)/n)
[1] 80

On note que l’information de Fisher calculée par double différenciation de la log-vraisemblance est identique à la formule analytique. Pour rappels, on a

Trois tests (équivalent asymptotiquement comme on le verra en cours) sont alors possibles.

Tout d’abord le test de Wald propose de travailler sur la différence entre l’estimateur du maximum de vraisemblance, et la valeur que l’on cherche à tester (comme le montre le dessin ci-dessous). Cette différence doit être “petite” si  est la vraie valeur.On peut alors utiliser la statistique suivante

Asymptotiquement, en utilisant le résultat évoqué au début, cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> alpha=0.05
> (T=(pml-p0)^2*IF)
[1] 0.1999970
> T<qchisq(1-alpha,df=1)
[1] TRUE

Ici, on accepte l’hypothèse que la pièce n’est pas pipée.

Ensuite le test du rapport de vraisemblance propose de travailler sur les valeurs de la log-vraisemblance. Là encore, cette différence doit être “petite” si  est la vraie valeur. Graphiquement, on mesure une distance normalisée

La statistique est celle d’un test bilatéral,

On pose alors

(le 2 sera détaillé un peu en cours, c’est mis là de manière à avoir trois tests équivalents). Asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=2*(logL(pml)-logL(p0)))
[1] 0.2003347
> T<qchisq(1-alpha,df=1)
[1] TRUE

Là encore, on accepte l’hypothèse que la pièce n’est pas pipée.

Enfin, le test du score propose de travailler sur la pente de la log-vraisemblance en  . Cette pente doit être “petite” si  est la vraie valeur.

La pente correspond au score,

La statistique de test est alors ici

Et là encore, asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=slope^2/IF)
[1] 0.2
> T<qchisq(1-alpha,df=1)
[1] TRUE

Et le test accepte ici encore l’hypothèse que la pièce n’est pas pipée.