Re-parametrization and Maximum Likelihood

The maximum likelihood estimator is invariant in the sense that for all bijective function  $h(\cdot)$, if  $\widehat{\theta}$ is the maximum likelihood estimator of  $\theta$ then  $h(\theta)$. Let  $\tau=h(\theta)$, then  $f_{\theta}(\cdot)$ is equal to  $f_{h^{-1}(\tau)}$, and the likelihood function in  $\tau$ is  $\ell(\tau)=\mathcal{L}(h^{-1}(\tau))$. And since  $\widehat{\theta}$ is the maximum likelihood estimator of  $\theta$,

$\mathcal{L}(h^{-1}(h(\widehat{\theta})))=\mathcal{L}(\widehat{\theta})\geq\mathcal{L}({\theta})= \mathcal{L}(h^{-1}(h({\theta})))$

hence,  $h(\widehat{\theta})$ is the maximum likelihood estimator of  $h(\theta)$.

For instance, the Bernoulli distribution is  $\mathcal{B}(p)$ with  $p\in[0,1]$ and

$f_p(x)=p^x(1-p)^(1-x),~ x\in\{0,1\}.$

Given sample  $\{x_1,\cdots,x_n\}$, the likelihood is

$\mathcal{L}(p)=\prod_{i=1}^n f_p(x_i)=p^{\sum x_i}[1-p]^{n-\sum x_i}.$

The log-likelihood is then

$\log\mathcal{L}(p)={\sum x_i}\log[p]+\Big({n-\sum x_i}\Big) \log[1-p]$

with ICI

$\frac{\partial}{\partial p}\log\mathcal{L}(p)=\frac{\sum x_i}{p}-\frac{n-\sum x_i}{1-p}.$

Thus, the first order condition

$\frac{\partial}{\partial p}\log\mathcal{L}(p)=0$

is satisfied when  $p=\overline{x}$. In order to illustrate, consider the following data

 > set.seed(1) > X=sample(0:1,size=15,replace=TRUE) > X [1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 

The (negative) log-likelihood is here

 > loglik=function(p){ + -sum(log(dbinom(X,size=1,prob=p))) + } 

that we can visualize below

 > u=seq(0,1,by=.025) > v=-Vectorize(loglik)(u) > plot(u,v,type="l",xlab="",ylab="") 

From calculations above, we know that the maximum likelihood estimator for $p$ is

 > mean(X) [1] 0.5333333 

The numerical version is

 > (opt=optim(.5,loglik)) $par [1] 0.5333008$value
[1] 10.36385

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

$message NULL Somehow, we were lucky here, because we did not say that the optimization was on the interval $[0,1]$. Nevertheless, our estimator for the probability belongs to $[0,1]$. In order to insure that the optimal value is in $[0,1]$, we can consider some constrained optimization routine  > constrOptim(.5, loglik, grad=NULL,ui=matrix(c(1,-1),2,1), ci=c(0,-1))$par [1] 0.5333008

$value [1] 10.36385$counts
20 NA

$convergence [1] 0$message
NULL

$outer.iterations [1] 2$barrier.value
[1] 6.909277e-05

On the previous graph, we did – indeed – reach that maximum of the log-likelihood

 > abline(v=opt$par,col="red")  An alternative is to consider $\theta=\log[p/(1-p)]$ (as in the exponential family). The log-likelihood is then $\log\ell(\theta)={\sum x_i}\theta+n\log(1+e^\theta)$ since $\log[1+e^\theta]=\log\Big[1+\frac{p}{1-p}\Big]=-\log[1-p]$ Here $\frac{\partial}{\partial \theta} \log\ell(p)=\sum x_i -n \frac{e^\theta}{1+e^\theta}$ Thus, the first order condition $\frac{\partial}{\partial p} \log\ell(p)=0$ is satisfied when $\frac{e^\theta}{1+e^\theta}=\overline{x}$ i.e. $\theta=\log\frac{\overline{x}}{1-\overline{x}}$ From a numerical perspective, we have the same optimal value  > loglik=function(theta){ + -sum(log(dbinom(X,size=1,prob=exp(theta)/(1+exp(theta))))) + } > (opt=optim(0,loglik))$par [1] 0.1335938

$value [1] 10.36385$counts
20 NA

$convergence [1] 0$message
NULL
> exp(opt$par)/(1+exp(opt$par))
[1] 0.5333489

Modeling the Marginals and the Dependence separately

When introducing copulas, it is commonly admitted that copulas are interesting because they allow to model the marginals and the dependence structure separately. The motivation is probably Sklar’s theorem, which says that given some marginal cumulative distribution functions (say $F$ and $G$, in dimension 2), and a copula (denoted $C$), then we can generate a multivariate cumulative distribution function with marginals the one specified previously, using

$H(x,y)=C(F(x),G(y))$

But this separability might be misleading. Consider the case of a fully parametric model,

$\left\{\begin{array}{l}F=F_{\alpha_0}\in\{F_\alpha;\alpha\in A\}\\G=G_{\beta_0}\in\{G_\beta;\beta\in B\} \\ C=C_{\gamma_0}\in\{C_\gamma;\gamma\in \Gamma\}\end{array}\right.$

Assume that those distributions are continuous, so that we can write the likelihood using densities,

$\mathcal{L}=\prod_{i=1}^n f_{\alpha}(X_i)\cdot g_{\beta}(Y_i)\cdot c_\gamma(F_{\alpha}(X_i),G_{\beta}(Y_i))$and the log-likelihood is

$\log\mathcal{L}=\underbrace{\sum_{i=1}^n \log f_{\alpha}(X_i)}_{\log\mathcal{L}(\alpha)}+\underbrace{\sum_{i=1}^n \log g_{\beta}(Y_i)}_{\log\mathcal{L}(\beta)}{\color{blue}+\sum_{i=1}^n \log c_\gamma(F_{\alpha}(X_i),G_{\beta}(Y_i))}$

The first part is the log-likelihood if we consider the first marginal (only). The second part is the log-likelihood if we consider the second marginal (only). If the two components are not independent (i.e. the copula density $c$ is not equal to 1 everywhere) the third part cannot be considered as null, and so, in a general context,

$(\alpha^\star,\beta^\star)\neq (\widehat{\alpha},\widehat{\beta})$

where

$(\alpha^\star,\beta^\star,\gamma^\star)=\text{argmax}\{\log \mathcal{L}\}$

while

$\left\{\begin{array}{l}\widehat{\alpha}=\text{argmax}\{\log \mathcal{L}(\alpha)\}\\ \widehat{\beta}=\text{argmax}\{\log \mathcal{L}(\beta)\}\end{array}\right.$

In order to illustrate this point, consider a bivariate lognormal distribution (obtained by taking the exponential of a Gaussian vector)

> mu1=1
> mu2=2
> MU=c(mu1,mu2)
> s1=1
> s2=sqrt(2)
> r=.8
> SIGMA=matrix(c(s1^2,r*s1*s2,r*s1*s2,s2^2),2,2)
> library(mnormt)
> set.seed(1)
> Z=exp(rmnorm(25,MU,SIGMA))

If we believe that marginals and correlations can be treated separately, we can start with marginal distributions.

> library(MASS)
> (p1=fitdistr(Z[,1],"lognormal"))
meanlog      sdlog
1.1686652   0.9309119
(0.1861824) (0.1316508)
> (p2=fitdistr(Z[,2],"lognormal"))
meanlog      sdlog
2.2181721   1.1684049
(0.2336810) (0.1652374)

Based on those marginal distributions, define $\tilde{U}_i=F_{\widehat{\alpha}}(X_i)$ and $\tilde{V}_i=G_{\widehat{\beta}}(Y_i)$, and consider the maximum likelihood estimator $\widehat{\gamma}}$ of the copula parameter, obtained from this pseudo sample,

$\widehat{\gamma}}=\text{argmax}\left\{\sum_{i=1}^n \log c_\gamma(\tilde U_i,\tilde V_i) \right\}$

Numerically, we get (since we consider a Gaussian copula, which is the true copula generated here)

> library(copula)
> Gcop=normalCopula(.3,dim=2)
> U=cbind(plnorm(Z[,1],p1$estimate[1],p1$estimate[2]),
+ plnorm(Z[,2],p2$estimate[1],p2$estimate[2]))
> fitCopula(Gcop,data=U,method="ml")
fitCopula() estimation based on 'maximum likelihood'
and a sample of size 25.
Estimate Std. Error z value Pr(>|z|)
rho.1  0.86530    0.03799   22.77

But clearly, we did not treat the dependence structure separately, since it was a function of marginal distributions,

$\widehat{\gamma}}=\text{argmax}\left\{\sum_{i=1}^n \log c_\gamma(F_{\widehat{\alpha}}(X_i),G_{\widehat{\beta}}(Y_i)) \right\}$

If we consider a global optimization problem, then results are different. The joint density can be derived (see e.g. Mostafa & Mahmoud (1964))

> dbivlognorm=function(x,theta){
+ mu1=theta[1]
+ mu2=theta[2]
+ s1=theta[3]
+ s2=theta[4]
+ r=theta[5]
+ a1=(log(x[,1])-mu1)/s1
+ a2=(log(x[,2])-mu2)/s2
+ d=1/(2*pi*s1*s2*sqrt(1-r^2))*1/(x[,1]*x[,2])*
+ exp(-(a1^2-2*r*a1*a2+a2^2)/(2*(1-r^2)))
+ return(d)
+ }
> LogLik=function(theta){
+ return(-sum(log(dbivlognorm(Z,theta))))}
> optim(par=c(0,0,1,1,0),fn=LogLik)$par [1] 1.1655359 2.2159767 0.9237853 1.1610132 0.8645052 The difference is not huge, but still. The estimators are not identical. From a statistical point of view, we can hardly treat the marginals and the dependence structure separately. Another point we should keep in mind is that the estimation of the copula parameter depends on the margins, not only through the parameters, but more deeply, through the choice of the marginal distributions (that might be misspecified). For instance, if we assume that margins are exponentially distributed, > (p1=fitdistr(Z[,1],"exponential")) rate 0.22288362 (0.04457672) > (p2=fitdistr(Z[,2],"exponential")) rate 0.06543665 (0.01308733) the estimation of the parameter of the Gaussian copula yields > U=cbind(pexp(Z[,1],p1$estimate[1]),
+ pexp(Z[,2],p2$estimate[1])) > fitCopula(Gcop,data=U,method="ml") fitCopula() estimation based on 'maximum likelihood' and a sample of size 25. Estimate Std. Error z value Pr(>|z|) rho.1 0.87421 0.03617 24.17 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 The maximized loglikelihood is 15.4 Optimization converged The problem is that since we misspecify marginal distribution, our pseudo sample is defined on the unit-interval, but there is no chance that we get uniform margins. If we generate a sample of size 500 with the code above, > x <- U[,1]; y <- U[,2] > xhist <- hist(x, plot=FALSE) ; yhist <- hist(y, plot=FALSE) > top <- max(c(xhist$counts, yhist$counts)) > nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE), c(3,1), c(1,3), TRUE) > par(mar=c(3,3,1,1)) > plot(x, y, xlab="", ylab="",col="red",xlim=0:1,ylim=0:1) > par(mar=c(0,3,1,1)) > barplot(xhist$counts, axes=FALSE, ylim=c(0, top),
+ space=0,col="light green")
> par(mar=c(3,0,1,1))
> barplot(yhist$counts, axes=FALSE, xlim=c(0, top), + space=0, horiz=TRUE,col="light blue") If we compare with the previous case, when marginal distribution were well-specified, we can clearly see that the dependence structure depends on marginal distributions, Inference for ARMA(p,q) Time Series As we mentioned in our previous post, as soon as we have a moving average part, inference becomes more complicated. Again, to illustrate, we do not need a two general model. Consider, here, some $ARMA(1,1)$ process, $X_t=\phi X_{t-1}+\varepsilon_t+\theta \varepsilon_{t-1}$ where $(\varepsilon_t)$ is some white noise, and assume further that $\theta+\phi\neq0$. > theta=.7 > phi=.5 > n=1000 > Z=rep(0,n) > set.seed(1) > e=rnorm(n) > for(t in 2:n) Z[t]=phi*Z[t-1]+e[t]+theta*e[t-1] > Z=Z[800:1000] > plot(Z,type="l") • A two step procedure To start with something simple, assume that we did miss the moving average component, $X_t=\phi X_{t-1}+u_t$ The estimator of $\phi$ – by least squares – is not longer consistent. But still. We can still compute it > base=data.frame(Y=Z[2:n],X=Z[1:(n-1)]) > regression=lm(Y~0+X,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X, data = base) Residuals: Min 1Q Median 3Q Max -3.2445 -0.7909 0.0626 0.9707 3.0685 Coefficients: Estimate Std. Error t value Pr(>|t|) X 0.69571 0.05101 13.64 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.225 on 199 degrees of freedom (799 observations deleted due to missingness) Multiple R-squared: 0.4832, Adjusted R-squared: 0.4806 F-statistic: 186 on 1 and 199 DF, p-value: < 2.2e-16 and then, we cancompute the autocorrelation of the noise, > n=200 > cor(residuals(regression)[2:n],residuals(regression)[1:(n-1)]) [1] 0.2663076 or more formally, use Durbin-Watson estimator, to get autocorrelation of the noise (and some significance test) > library(car) > durbinWatsonTest(regression) lag Autocorrelation D-W Statistic p-value 1 0.2656555 1.46323 0 Alternative hypothesis: rho != 0 The point, here, is that we would like to assume that $u_t=\varepsilon_t+\theta\varepsilon_{t-1}$ meaning that $(u_t)$ should be some $MA(1)$ process. And $\rho(1)=\frac{\theta}{1+\theta^2}$ i.e. $\theta$ is a root of this quadratic problem, > polyroot(c(1,-1/cor(residuals(regression)[2:n],residuals(regression)[1:(n-1)]),1)) [1] 0.2884681+0i 3.4665883+0i Here, we do have two positive roots. I would go for the one smaller than one, in order to be able to invert the polynomial, if necessary… • Use of the empirical autocorrelation function An alternative might be to use properties of the autocorrelation function, $\gamma(0)=\frac{1+\theta^2+2\phi\theta}{1-\phi^2}\cdot\sigma^2$ $\gamma(1)=\phi\gamma(0)+\theta\sigma^2$ and $\gamma(2)=\phi\gamma(1)$ Again, we have a set of three equations, with three unknown parameters. Numerically, it is possible to find some roots. If we run the code, we get > v=c(as.numeric(acf(Z)$acf[2:3]),1)*var(Z)
> library(rootSolve)
> seteq=function(x){
+ F1=v[1]-x[3]^2*(x[2]^2+2*x[1]*x[2]+1)/(1-x[1]^2)
+ F2=v[2]-(x[1]*v[1]+x[2]*x[3]^2)
+ F3=v[3]-x[1]*v[2]
+ return(c(F1,F2,F3))}
> multiroot(f=seteq,start=c(.1,.1,1))
$root [1] 3.643734 -3.188145 1.427759$f.root
[1]  1.371170e-11 -3.714573e-11  0.000000e+00

$iter [1] 8$estim.precis
[1] 1.695248e-11

Here, we have a situation…

• Use of least square techniques

We can use, here, the algorithm described in the context of $MA(q)$ processes.

> V=function(p){
+ phi=p[1]
+ theta=p[2]
+ u=rep(0,length(Z))
+ for(t in 2:length(Z)) u[t]=Z[t]-phi*Z[t-1]-theta*u[t-1]
+ return(sum(u^2))
+ }
> p=optim(par=c(.1,.1),V)$par [1] 0.3637783 0.7773845 > coef=c(p,sqrt(V(p)/(length(Z)))) which is not so bad. Actually, if we run that procecure on 1,000 samples, we get the following output • Use of maximum likelihood techniques Last, but not least, one more time, we can use (global) maximum likelihood techniques, since the process is a Gaussian process (all finite dimensional vector will have a joint Gaussian distribution) if we assume that the noise is Gaussian. > library(mnormt) > GlobalLogLik=function(A,TS){ + n=length(TS) + phi=A[1]; theta=A[2] + sigma=A[3] + SIG=matrix(0,n,n) + rho=rep(0,n) + rho[1]=sigma^2*(theta^2+2*phi*theta+1)/(1-phi^2) + rho[2]=phi*rho[1]+theta*sigma^2 + for(h in 3:n) rho[h]=phi*rho[h-1] + for(i in 1:n){for(j in 1:n){ + SIG[i,j]=rho[abs(i-j)+1]}} + return(dmnorm(TS,rep(0,n),SIG,log=TRUE))} > LogL=function(A) -GlobalLogLik(A,TS=Z) > optim(c(.1,.1,1),LogL)$par
[1] 0.3890991 0.7672036 1.0731340

It works fine, one more time. But maybe we got lucky here. We’ve seen in the post on autoregressive time series that the algorithm might fell if the time series is not stationary. In order to avoid such problems, we can consider a constraint optimization problem, where we simply recall that $\phi\in(-1,1)$,

> U=matrix(c(1,-1,0,0,0,0),2,3)
> C=-c(.999,.999)
$par [1] 0.3890991 0.7672036 1.0731340$value
[1] 300.1956

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

$message NULL$outer.iterations
[1] 2

$barrier.value [1] -1.536358e-05 If we run that algorithm 1,000 times, on simulated time series (with the same parameters), we get Inference for AR(p) Time Series Consider a (stationary) autoregressive process, say of order 2, $Y_t =\varphi_1 Y_{t-1}+\varphi_2 Y_{t-2}+\varepsilon_t$ for some white noise with variance . Here is a code to generate such a process, > phi1=.25 > phi2=.7 > n=1000 > set.seed(1) > e=rnorm(n) > Z=rep(0,n) > for(t in 3:n) Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+e[t] > Z=Z[800:1000] > n=length(Z) > plot(Z,type="l") Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . Several techniques can be used to estimate those parameters. • using least square regression A natural idea is to see here a regression model, since (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 -3.0268 -0.7063 0.1065 0.6925 3.2566 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 0.23400 0.05463 4.283 2.88e-05 *** X2 0.62863 0.05476 11.479 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.062 on 197 degrees of freedom Multiple R-squared: 0.6349, Adjusted R-squared: 0.6312 F-statistic: 171.3 on 2 and 197 DF, p-value: < 2.2e-16 so we get the following estimators, for the autocorrelation coefficients, and the volatility of the noise > regression$coefficients
X1        X2
0.2339959 0.6286321
> summary(regression)$sigma [1] 1.061839 • 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 (again, using a matrix expression) 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.2256270 [2,] 0.6315329 Now, we need to extract the estimated innovation process, from this set of parameters > estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2)
> sd(estWN)
[1] 1.058558

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.

An alternative could be to include the variance term in Yule-Walker equations, to get a three dimensional linear equation,

$\left\{\begin{array}{l} \gamma_0 = \varphi_1 \gamma_1+\varphi_2 \gamma_2+\sigma^2\\ \gamma_1=\varphi_1 \gamma_0+\varphi_2 \gamma_1 \\ \gamma_2=\varphi_1 \gamma_1+\varphi_2 \gamma_0\end{array}\right.$

It is not much more complicated to solve, actually,

> gamma0=var(Z[1:n])
> gamma1=var(Z[1:(n-1)],Z[2:n])
> gamma2=var(Z[1:(n-2)],Z[3:n])
> A=matrix(c(gamma1,gamma0,gamma1,gamma2,gamma1,gamma0,1,0,0),3,3)
> b=matrix(c(gamma0,gamma1,gamma2),3,1)
> (PHISIGMA=solve(A,b))
[,1]
[1,] 0.2283151
[2,] 0.6283431
[3,] 1.1335501
• using (conditional) likelihood estimators

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

$Y_t\vert Y_{t-1}=y_{t-1},Y_{t-2}=y_{t-2}$

has a Gaussian distribution

$\mathcal{N}(\varphi_1y_{t-1}+\varphi_2y_{t-2},\sigma^2)$

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.2339589 0.6285002 1.0565613$value
[1] 293.3042

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

$message NULL It is also possible to consider a global maximum likelihood optimisation problem, since the variance matrix of vector $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ has a know form. • using (unconditional) likelihood estimators The variance matrix of $\boldsymbol{Y}=(Y_1,\cdots,Y_t)$ is $\boldsymbol{\Gamma}=[\gamma(\vert i-j\vert)]$, where autocovariances are not not know, be can easily be computed using a recursive relationship. > library(mnormt) > GlobalLogLik=function(A,TS){ + n=length(TS) + phi1=A[1]; phi2=A[2] + sigma=A[3] + SIG=matrix(0,n,n) + rho=rep(0,n) + rho[1]=1 + rho[2]=phi1/(1-phi2) + for(h in 3:n) rho[h]=phi1*rho[h-1]+phi2*rho[h-2] + for(i in 1:n){for(j in 1:n){ + SIG[i,j]=rho[abs(i-j)+1]}} + gamma0=(1-phi2)*sigma^2/((1+phi2)*((1-phi2)^2-phi1^2)) + SIG=gamma0*SIG + return(dmnorm(TS,rep(0,n),SIG,log=TRUE))} > LogL=function(A) -GlobalLogLik(A,TS=Z) > optim(c(.1,.1,1),LogL) Error in chol.default(x, pivot = FALSE) : Error in pd.solve(varcov, log.det = TRUE) : x appears to be not positive definite The problem is that there is a strong constraint on the pair $(\varphi_1,\varphi_2)$ to get a stationary process (we are not far away, here, from the border of the triangle, where the process become non stationary). To be more specific (this was mentioned in a previous post), we should have $\left\{\begin{array}{l} \phi_2-\phi_1<1 \\\phi_2+\phi_1<1\\ \vert\phi_2\vert<1\end{array}\right.$ i.e. in a standard matrix form $\left[\begin{array}{cc} +1 & -1 \\ -1 & -1 \\ 0 & +1\end{array}\right]\left[\begin{array}{c} \varphi_1 \\ \varphi_2\end{array}\right] > \left[\begin{array}{c} -1 \\ -1 \\ -1\end{array}\right]$ (we can add an additional constraint on the variance parameter, to insure that it will be positive). To run a contrained optimization routine, consider > U=matrix(c(1,0,0,-1,0,1,0,-1,0,0,1,0),4,3) > C=c(0,0,0,-.99999) > constrOptim(c(.1,.1,1),LogL,grad=NULL,ui=U,ci=C)$par
[1] 0.2238892 0.6342850 1.0613388

$value [1] 297.9202$counts
108       NA

$convergence [1] 0$message
NULL

$outer.iterations [1] 2$barrier.value
[1] 0.000189892

(here, to faster, we restrain the parameters so that they will be positive).

• comparing those estimates

Here, our five estimators are rather close. Let us run more samples to see more precisely how they behave. For the first parameter $\widehat{\varphi_1}$, we get

and for the second one, $\widehat{\varphi_2}$, we have

The bias we observe is probably coming from the fact that, with this numerical example, we are not far away from the non-stationary case (the sum of the true parameters should be less than 1, and it is 0.95). When we estimate the parameters, we force them to be inside the triangle, since those parameters can be estimated only if the process is stationary.

Observe that the standard-deviation of the innovation process $\widehat{\sigma}$ is here, well estimated,

(with clearly some estimators that perform better than others).

Bias and MLE

Before leaving the office, this evening, JP decided to knock at my door to ask me a “quick and very basic question” (as he put it). This is JP’s stategy, and he knows it works. His question was – more or less – what do we know about the bias in maximum likelihood estimation when we have a small sample, from a Gamma distribution. He was surprised by some results he got. If I wanted to be naughty, too, I would say that he was suprised to see how long his student spent to code that in SAS. So he wanted to challenge me, and see how fast I could give him a valuable answer. Given the fact that I had to leave early because my elder son had a fencing competition, I tried to write a simple code to “visualize” the bias of the parameter (the first one) of a Gamma distribution, with MLE.

Before showing the graph, I wanted to add that I hate one thing about mathematical statistics courses: we learn nothing interesting there. I mean, we can see nice mathematical concepts, but after this class, you can hardly say anything when you see your first dataset. Like with real data. For instance, this course usually emphasize asymptotical results, using limiting theorem. When you take this course, you learn a lot of thing about maximumum likelihood for instance. You can compute the asymptotic variance and derive asymptotic confidence intervals. But are those results relevant when you have 50 observations? Is it possible, with 50 observations, to have a bias which has the same size as the parameter?

As usual, one possible answer is “if you don’t have a lot of observations, be Bayesian!“. Maybe. Someday. What I tried, here, is to run simulations to see how MLE estimators behave. Given a $n$-i.i.d. sample, from a $\mathcal{G}(\alpha,\beta)$ distribution, let $\widehat{\alpha}_n$ and $\widehat{\beta}_n$ denote the maximum likelihood estimators of the two parameters.

library(fitdistrplus)
maxl=function(x) fitdist(x,"gamma",method="mle")$estimate VK=floor(exp(seq(log(5),log(200),length=25))) V=NULL for(k in 1:length(VK)){ n=VK[k] N=5000 m=matrix(rgamma(n*N,1,2),n,N) ss=apply(m,2,maxl) V=rbind(V,ss)} y=as.vector(V[seq(1,length(VK)*2,by=2),]) x=rep(c(VK),ncol(V)) boxplot(y~x, xlab="Nb. observations (log scale)",ylim=c(0,6)) abline(h=1,lty=2,col="blue") Here, in our simulations, the shape parameter was 1. On the graph, we have boxplots of $\widehat{\alpha}_n$ obtained on several scenarios. We clearly see the positive bias of the MLE. And the bias reduces with $n$ (as expected, since the MLE is asymptotically unbiased). We can also visualize the distribution of $\widehat{\alpha}_n$ (instead of boxplots) It is also possible to derive analytical results. David Cox and Joyce Snell did the maths in 1968 and actually did obtain analytical expressions for the biases. More recently, David Giles (a.k.a. @deagiles on Twitter) and Hui Feng did look at the behavior of bias-adjusted estimators, a few years ago. For instance, one can get that $\text{bias}({\alpha}_n)=\frac{\kappa}{n}+O\left(\frac{1}{n^2}\right)$ where $\kappa=\frac{\alpha \Psi'(\alpha)-\alpha^2 \Psi''(\alpha)-2}{2[\alpha \Psi'(\alpha)-1]^2}$ $\Psi$ being the so-called digamma function, $\Psi(z)=\frac{d \log \Gamma(z)}{dz}$ and where $\Psi'$ and $\Psi''$ are the first and second order derivatives, see e.g. Bowman and Shenton (1982) – yes, there is an book on the topic of estimating parameters of the Gamma distribution… Observe that the bias of $\widehat{\alpha}_n$ does not depend on $\beta$, while the bias of $\widehat{\beta}_n$ will depend on $\alpha$. d1digamma=function(x,h=1e-7) return((digamma(x+h)-digamma(x-h))/(2*h)) d2digamma=function(x,h=1e-7) return((d1digamma(x+h)-d1digamma(x-h))/(2*h)) biasalpha=function(a,n){ return((a*d1digamma(a)-a^2*d2digamma(a) -2)/(2*n*(a*d1digamma(a)-1)^2)) } The way I compute it is probably not optimal, so if you want to improve it, please, go ahead ! If we compare the average bias obtained on our simulation, and the one obtained this first order approximation, we get m=apply(V,1,mean) plot(VK,m[seq(1,length(VK)*2,by=2)],type="b",col="red",xlab="Nb. observations (log scale)",log="x") abline(h=1,lty=2,col="blue") B=Vectorize(function(n) biasalpha(a=1,n))(1:200) lines(1:200,B+1,col="orange") Observe here that neglecting the $1/n^2$ factor yield an underestimation of the real biais… Fun, isn’t it? Maximum Likelihood versus Goodness of Fit Thursday, I got an interesting question from a colleague of mine (JP). I mean, the way I understood the question turned out to be a nice puzzle (but I have to confess I might have misunderstood). The question is the following : consider a i.i.d. sample $\{X_1,\cdots,X_n\}$ of continuous variables. We would like to choose between two (parametric) families for the distribution, $\mathcal{F}=\{F_{\boldsymbol \alpha};\boldsymbol \alpha\in\mathcal{A}\}$ and $\mathcal{G}=\{G_{\boldsymbol \beta};\boldsymbol \beta\in\mathcal{B}\}$. If we use maximum likelihood techniques, we get two estimators, one for each family, $\widehat{\boldsymbol \alpha}$ and $\widehat{\boldsymbol \beta}$. Clearly, $F_{\widehat{\boldsymbol \alpha}}(\cdot)$ is a much better than $G_{\widehat{\boldsymbol \beta}}(\cdot)$, in the sense of a standard goodness of fit test (e.g. Kolmogorov-Smirnov since the sample is assumed to be obtained from a continuous variable). Does that mean that family is $\mathcal{F}$ (somehow) better than family $\mathcal{G}$ ? This is my interpretation of the question, and I found it amusing. So I will try to show (using simulated samples) that some odd situations can easily be obtained. Consider a sample from a mixture of log-normal distributions, > set.seed(228) > X=exp(c(rnorm(50,1,1),rnorm(50,2,1.2))) Consider two standard families for positive random variables: a Gamma distribution and a lognormal distribution. > library(MASS) > ab=fitdistr(X,"gamma") > ms=fitdistr(X,"lognormal") If we want to visualized those two distributions, let us use > vab=pgamma(u,ab$estimate[1],ab$estimate[2]) > vms=plnorm(u,ms$estimate[1],ms$estimate[2]) > plot(ecdf(X)) > lines(u,vab,col="red") > lines(u,vms,col="blue") Here, we get What else can we say ? actually, we can also compute Kolmogorov-Smirnov statistic, $D_n=\sup_x |\widehat F_n(x)-F_\star(x)|$where $\widehat F_n(x)={1 \over n}\sum_{i=1}^n \boldsymbol{1}_{X_i\leq x}$ This can be done using > ks.test(X,"plnorm",ms$estimate[1],ms$estimate[2]) One-sample Kolmogorov-Smirnov test data: X D = 0.0693, p-value = 0.7231 alternative hypothesis: two-sided > ks.test(X,"pgamma",ab$estimate[1],ab$estimate[2]) One-sample Kolmogorov-Smirnov test data: X D = 0.148, p-value = 0.02507 alternative hypothesis: two-sided From a theoretical point of view, we should not look at the p-values, since the null-distribution is based on a fixed distribution, not a fitted one (see the Lilliefors tests for normal samples). But still. The Gamma distribution seems to be very far away from the true distribution. The statistics is twice the one we have with our lognormal distribution. And one p-value is 72%, while the other one is 2.5%. Here, we should prefer this lognormal distribution to that Gamma one. But here, we did consider only one distribution in each family. Does that mean that we cannot find one Gamma distribution that will be better than all possible lognormal distributions ? Better, for instance, according to Kolmogorov-Smirnov statistics… Well, it is possible to use another strategy to find appropriate parameters. We can minimize this statistic actually. Define > ks1=function(ms) {m=ms[1];s=ms[2];ks.test(X,"plnorm",m,s)$statistic}
> ks2=function(ab) {a=ab[1];b=ab[2];ks.test(X,"pgamma",a,b)$statistic} and compute > n1=nlm(ks1,c(ms$estimate[1],ms$estimate[2])) > n1$minimum
[1] 0.05252692

$estimate [1] 1.547437 1.121864 > n2=nlm(ks2,c(ab$estimate[1],ab$estimate[2])) > n2$minimum
[1] 0.04737725

$estimate [1] 1.1449041 0.167041 So here, it is possible to find a distribution much closer to the empirical sample, within the Gamma family actually. > vab=pgamma(u,n2$estimate[1],n2$estimate[2]) > vms=plnorm(u,n1$estimate[1],n1$estimate[2]) > lines(u,vab,col="red",lwd=2) > lines(u,vms,col="blue",lwd=2) What would be the point here? Maybe just the idea that the maximum likelihood estimator is only one estimator among a lot of them. And if it has interesting asymptotic properties, on small samples, it might not be the best estimator to consider… And to be completely honest, I’ve been cheating here… I mean, not really cheating (not more than any researcher using a statistical test to validate the findings). But here, I did fix the seed of the random number generator. Actually, such example does not occur that frequently. Here, out of 1000 samples, I got this odd conclusion almost 15 times. And the smaller the sample, the more likely we can observe that story, where the maximum likelihood estimator can be far away from the best fit. Here is the proportion of opposite conclusions, as a function of the sample size, > SIM=function(ns=1000,n=100){ + t=0 + for(s in 1:ns){ + set.seed(s) + X=exp(c(rnorm(n/2,1,1),rnorm(n/2,2,1.2))) + ks1=function(ms) {m=ms[1];s=ms[2];ks.test(X,"plnorm",m,s)$statistic}
+  ks2=function(ab) {a=ab[1];b=ab[2];ks.test(X,"pgamma",a,b)$statistic} + library(MASS) + ab=fitdistr(X,"gamma") + ms=fitdistr(X,"lognormal") + n1=nlm(ks1,c(ms$estimate[1],ms$estimate[2])) + n2=nlm(ks2,c(ab$estimate[1],ab$estimate[2])) + if( (ks.test(X,"plnorm",ms$estimate[1],ms$estimate[2])$statistic-
+  ks.test(X,"pgamma",ab$estimate[1],ab$estimate[2])$statistic) + *(n1$minimum-n2$minimum)<=0 ) t=t+1 + } + return(t/ns)} > VM=rep(NA,20) > VS=seq(10,200,by=10) > for(i in 1:20){VM[i]=SIM(n=VS[i],ns=1000)} > plot(VS,VM,type="p") So to provide a more complete answer to JP’s question, with a very large sample, I guess that your intuition should be valid…. but clearly not on a small sample. 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…

Tail index estimation

These data were collected at Copenhagen Reinsurance and comprise 2167 fire losses over the period 1980 to 1990, They have been adjusted for inflation to reflect 1985 values and are expressed in millions of Danish Kron. Note that it is possible to work with the same data as above but the total claim has been divided into a building loss, a loss of contents and a loss of profits.

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ "http://freakonometrics.free.fr/danish-multivariate.txt",
+ header=TRUE)

Consider here the first dataset (we deal – so far – with univariate extremes),

> X=base1$Loss.in.DKM > D=as.Date(as.character(base1$Date),"%m/%d/%Y")
> plot(D,X,type="h")

The graph is the following,

A natural idea is then to plot

i.e.

> Xs=sort(X)
> logXs=rev(log(Xs))
> n=length(X)
> plot(log(Xs),log((n:1)/(n+1)))

Points are on a straight line here. The slope can be obtained using a linear regression,

> B=data.frame(X=log(Xs),Y=log((n:1)/(n+1)))
> reg=lm(Y~X,data=B)
> summary(reg)

Call:
lm(formula = Y ~ X, data = B)

Residuals:
Min       1Q   Median       3Q      Max
-0.59999 -0.00777  0.00878  0.02461  0.20309

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.089442   0.001572   56.88   <2e-16 ***
X           -1.382181   0.001477 -935.55   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04928 on 2165 degrees of freedom
Multiple R-squared: 0.9975,	Adjusted R-squared: 0.9975
F-statistic: 8.753e+05 on 1 and 2165 DF,  p-value: < 2.2e-16

> reg=lm(Y~X,data=B[(n-500):n,])
> summary(reg)

Call:
lm(formula = Y ~ X, data = B[(n - 500):n, ])

Residuals:
Min       1Q   Median       3Q      Max
-0.48502 -0.02148 -0.00900  0.01626  0.35798

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.186188   0.010033   18.56   <2e-16 ***
X           -1.432767   0.005105 -280.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07751 on 499 degrees of freedom
Multiple R-squared: 0.9937,	Adjusted R-squared: 0.9937
F-statistic: 7.878e+04 on 1 and 499 DF,  p-value: < 2.2e-16

> reg=lm(Y~X,data=B[(n-100):n,])
> summary(reg)

Call:
lm(formula = Y ~ X, data = B[(n - 100):n, ])

Residuals:
Min       1Q   Median       3Q      Max
-0.33396 -0.03743  0.02279  0.04754  0.62946

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.67377    0.06777   9.942   <2e-16 ***
X           -1.58536    0.02240 -70.772   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1299 on 99 degrees of freedom
Multiple R-squared: 0.9806,	Adjusted R-squared: 0.9804
F-statistic:  5009 on 1 and 99 DF,  p-value: < 2.2e-16

The slope here is somehow related to the tail index of the distribution. Consider some heavy tailed distribution, i.e. , so that , where  is some slowly varying function. Equivalently, the exists a slowly varying function  such that . Then

i.e. since a natural estimator for  is the order statistic , the slope of the straight line is the opposite of tail index . The estimator of the slope is (considering only the  largest observations)

Hill‘s estimator is based on the assumption that the denominator above is almost 1 (which means that  , as ), i.e.

Note that, if , but not two fast, i.e.  as , then  (one can even get   with stronger convergence assumptions). Further

Based on that (asymptotic) distribution, it is possible to get a (asymptotic) confidence interval for

> xi=1/(1:n)*cumsum(logXs)-logXs
> xise=1.96/sqrt(1:n)*xi
> plot(1:n,xi,type="l",ylim=range(c(xi+xise,xi-xise)),
+ xlab="",ylab="",)
> polygon(c(1:n,n:1),c(xi+xise,rev(xi-xise)),
+ border=NA,col="lightblue")
> lines(1:n,xi+xise,col="red",lwd=1.5)
> lines(1:n,xi-xise,col="red",lwd=1.5)
> lines(1:n,xi,lwd=1.5)
> abline(h=0,col="grey")

It is also possible to work with , then . And similarly  as  (and again  with additional assumptions on the rate of convergence), and

(obtained using the delta-method). Again, we can use that result to derive (asymptotic) confidence intervals

> alpha=1/xi
> alphase=1.96/sqrt(1:n)/xi
> YL=c(0,3)
> plot(1:n,alpha,type="l",ylim=YL,xlab="",ylab="",)
> polygon(c(1:n,n:1),c(alpha+alphase,rev(alpha-alphase)),
+ border=NA,col="lightblue")
> lines(1:n,alpha+alphase,col="red",lwd=1.5)
> lines(1:n,alpha-alphase,col="red",lwd=1.5)
> lines(1:n,alpha,lwd=1.5)
> abline(h=0,col="grey")

The Deckers-Einmahl-de Haan estimator is

where for

Then (given again conditions on the speed of convergence i.e. , with  as ),

Finally, Pickands‘ estimator

it is possible to prove that, as ,

Here the code is

> Xs=rev(sort(X))
> xi=1/log(2)*log( (Xs[seq(1,length=trunc(n/4),by=1)]-
+ Xs[seq(2,length=trunc(n/4),by=2)])/
+ (Xs[seq(2,length=trunc(n/4),by=2)]-Xs[seq(4,
+ length=trunc(n/4),by=4)]) )
> xise=1.96/sqrt(seq(1,length=trunc(n/4),by=1))*
+sqrt( xi^2*(2^(xi+1)+1)/((2*(2^xi-1)*log(2))^2))
> plot(seq(1,length=trunc(n/4),by=1),xi,type="l",
+ ylim=c(0,3),xlab="",ylab="",)
> polygon(c(seq(1,length=trunc(n/4),by=1),rev(seq(1,
+ length=trunc(n/4),by=1))),c(xi+xise,rev(xi-xise)),
+ border=NA,col="lightblue")
> lines(seq(1,length=trunc(n/4),by=1),
+ xi+xise,col="red",lwd=1.5)
> lines(seq(1,length=trunc(n/4),by=1),
+ xi-xise,col="red",lwd=1.5)
> lines(seq(1,length=trunc(n/4),by=1),xi,lwd=1.5)
> abline(h=0,col="grey")

It is also possible to use maximum likelihood techniques to fit a GPD distribution over a high threshold.

> library(evd)
> library(evir)
> gpd(X,5)
$n [1] 2167$threshold
[1] 5

$p.less.thresh [1] 0.8827873$n.exceed
[1] 254

$method [1] "ml"$par.ests
xi      beta
0.6320499 3.8074817

$par.ses xi beta 0.1117143 0.4637270$varcov
[,1]        [,2]
[1,]  0.01248007 -0.03203283
[2,] -0.03203283  0.21504269

$information [1] "observed"$converged
[1] 0

$nllh.final [1] 754.1115 attr(,"class") [1] "gpd" or equivalently (or almost) > gpd.fit(X,5)$threshold
[1] 5

$nexc [1] 254$conv
[1] 0

$nllh [1] 754.1115$mle
[1] 3.8078632 0.6315749

$rate [1] 0.1172127$se
[1] 0.4636270 0.1116136

The interest of the latest function is that it is possible to visualize the profile likelihood of the tail index,

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

or

> gpd.profxi(gpd.fit(X,20),xlow=0,xup=3)

Hence, it is possible to plot the maximum likelihood estimator of the tail index, as a function of the threshold (including a confidence interval),

> GPDE=Vectorize(function(u){gpd(X,u)$par.ests[1]}) > GPDS=Vectorize(function(u){ + gpd(X,u)$par.ses[1]})
> u=c(seq(2,10,by=.5),seq(11,25))
> XI=GPDE(u)
> XIS=GPDS(u)
> plot(u,XI,ylim=c(0,2))
> segments(u,XI-1.96*XIS,u,XI+
+ 1.96*XIS,lwd=2,col="red")

Finally, it is possible to use block-maxima techniques.

> gev.fit(X)
$conv [1] 0$nllh
[1] 3392.418

$mle [1] 1.4833484 0.5930190 0.9168128$se
[1] 0.01507776 0.01866719 0.03035380

The estimator of the tail index is here the last coefficient, on the right.
Since it is rather difficult to install a package in class rooms, here is the source of rcodes used here (to fit a GPD for exceedances)

> source("http://freakonometrics.blog.free.fr/public/code/gpd.R")

Next time, we will discuss how to use those estimators.

Multivariate probit regression using (direct) maximum likelihood estimators

Consider a random pair of binary responses, i.e. with taking values 1 or 2. Assume that probability can be function of some covariates.

• The Gaussian vector latent structure

A standard model is based a latent Gaussian structure, i.e. there exists some random vector such that if is lower than a given threshold, and 1 otherwise.
As in standard probit models, assume that

where we can assume that is a Gaussian random vector. This assumption can be used to derive the likelihood of a sample .

> logV=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+ sigma=matrix(c(1,CORRELATION,CORRELATION,1),2,2)
+     a11=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
+ for(i in 2:nrow(z)){a11=c(a11,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
+     a10=pnorm(z[1,1],sd=sqrt(sigma[1,1]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a10=c(a10,pnorm(z[i,1],sd=sqrt(sigma[1,1]))-pmnorm(z[i,],varcov=sigma))}
+     a01=pnorm(z[1,2],sd=sqrt(sigma[2,2]))-pmnorm(z[1,],varcov=sigma)
+ for(i in
+ 2:nrow(z)){a01=c(a01,pnorm(z[i,2],sd=sqrt(sigma[2,2]))-pmnorm(z[i,],varcov=sigma))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     1*log(a01) +
+     2*log(a10) +
+     3*log(a00) )
+ }
> OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par (the code is a bit long since I had trouble working properly with matrices – or more precisely to vectorize my functions – so I used loops… I am sure it is possible to write a better code). It is possible to generate samples (based on that specific model) to check that we can actually derive proper maximum likelihood estimators, > library(mnormt) > set.seed(1) > n=1000 > r=0.5 > X1=runif(n) > X2=rnorm(n) > Y1S=1+5*X1 > Y2S=8-5*X1 > RES=rmnorm(n,mean=c(0,0),varcov=matrix(c(1,r,r,1),2,2)) > YS=cbind(Y1S,Y2S)+RES > Y1=(YS[,1]>quantile(YS[,1],.5))*1 > Y2=(YS[,2]>quantile(YS[,2],.5))*1 > base=data.frame(i,Y1,Y2,X1,X2,YS) > head(base) i Y1 Y2 X1 X2 Y1S Y2S 1 1 0 0 0.2655087 0.07730312 3.177587 5.533884 2 2 0 0 0.3721239 -0.29686864 1.935307 5.089524 3 3 1 0 0.5728534 -1.18324224 4.757848 5.172584 4 4 1 0 0.9082078 0.01129269 4.600029 3.878225 5 5 0 1 0.2016819 0.99160104 2.547362 6.743714 6 6 1 0 0.8983897 1.59396745 5.309974 4.421523 (the two columns on the right are latent observations, that cannot be used since theoretically they are unobservable). Note that it is a simple regression, one of the component is here only to bring some noise. First of all, let us look at marginal probit regression > reg1=glm(Y1~X1+X2,data=base,family=binomial) > reg2=glm(Y2~X1+X2,data=base,family=binomial) > summary(reg1) Call: glm(formula = Y1 ~ X1 + X2, family = binomial, data = base) Deviance Residuals: Min 1Q Median 3Q Max -2.90570 -0.50126 -0.00266 0.49162 2.78256 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -4.291725 0.267149 -16.065 <2e-16 X1 8.656836 0.510153 16.969 <2e-16 *** X2 0.007375 0.090530 0.081 0.935 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1386.29 on 999 degrees of freedom Residual deviance: 726.48 on 997 degrees of freedom AIC: 732.48 Number of Fisher Scoring iterations: 5 > summary(reg2) Call: glm(formula = Y2 ~ X1 + X2, family = binomial, data = base) Deviance Residuals: Min 1Q Median 3Q Max -2.74682 -0.51814 -0.00001 0.57969 2.58565 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.91709 0.24399 16.054 <2e-16 *** X1 -7.89703 0.46277 -17.065 <2e-16 *** X2 0.18360 0.08758 2.096 0.036 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 1386.29 on 999 degrees of freedom Residual deviance: 777.61 on 997 degrees of freedom AIC: 783.61 Number of Fisher Scoring iterations: 5 Here, the optimization yields, > OPT=optim(fn=logV,par=c(0,1,1,1,1,1,1),method="BFGS")$par
> OPT[1]
[1] 0.5261382
> matrix(OPT[2:7],2,3)
[,1]      [,2]       [,3]
[1,] -2.451721  4.908633 0.01600769
[2,]  2.241962 -4.539946 0.10614807

Note that the coefficients we have obtained are almost identical to the ones obtained with R standard procedure,

>  library(Zelig)
>  REG= zelig(list(mu1=Y1~X1+X2,
+             mu2=Y2~X1+X2,
+     rho=~1),
+     model="bprobit",data=base)
>  summary(REG)

Call:
zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2,
rho = ~1), model = "bprobit", data = base)

Pearson Residuals:
Min        1Q     Median      3Q     Max
probit(mu1) -10.5442 -0.377243  0.0041803 0.36709 8.60398
probit(mu2)  -7.8547 -0.376888  0.0083715 0.42923 5.88264
rhobit(rho) -13.8322 -0.091502 -0.0080544 0.37218 0.85101

Coefficients:
Value Std. Error   t value
(Intercept):1 -2.451699   0.135369 -18.11116
(Intercept):2  2.241964   0.125072  17.92536
(Intercept):3  1.169461   0.189771   6.16249
X1:1           4.908617   0.252683  19.42602
X1:2          -4.539951   0.233632 -19.43203
X2:1           0.015992   0.050443   0.31703
X2:2           0.106154   0.049092   2.16235

Number of linear predictors:  3

Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho)
&n
bsp;
Dispersion Parameter for binom2.rho family:   1

Residual Deviance: 1460.355 on 2993 degrees of freedom

Log-likelihood: -730.1774 on 2993 degrees of freedom

Number of Iterations: 3

> matrix(coefficients(REG)[c(1:2,4:7)],2,3)
[,1]      [,2]       [,3]
[1,] -2.451699  4.908617 0.01599183
[2,]  2.241964 -4.539951 0.10615443

The correlation here is also the same

> (exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1)
[1] 0.5260951

That procedure works well an can be extended to ordinal responses (not only binary ones, or to three dimensional problems,

logV=function(beta){
BETA=matrix(beta[4:(3+ncol(Y)*ncol(X))],ncol(Y),ncol(X))
z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]),X%*%(BETA[3,]))
r12=beta[1]
r23=beta[2]
r31=beta[3]
s1=s2=s3=1
sigma=matrix(c(s1^2,r12*s1*s2,r31*s1*s3,
r12*s1*s2,s2^2,r23*s2*s3,
r31*s1*s3,r23*s2*s3,s3^2),3,3)
sigma1=matrix(c(s2^2,r23*s2*s3,
r23*s2*s3,s3^2),2,2)
sigma2=matrix(c(s1^2,r31*s1*s3,
r31*s1*s3,s3^2),2,2)
sigma3=matrix(c(s1^2,r12*s1*s2,
r12*s1*s2,s2^2),2,2)
a111=pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a111=c(a111,pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a011=pmnorm(z[1,2:3],varcov=sigma1)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a011=c(a011,pmnorm(z[i,2:3],varcov=sigma1)-pmnorm(z[i,],varcov=sigma))}
a101=pmnorm(z[1,c(1,3)],varcov=sigma2)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a101=c(a101,pmnorm(z[i,c(1,3)],varcov=sigma2)-pmnorm(z[i,],varcov=sigma))}
a110=pmnorm(z[1,1:2],varcov=sigma3)-pmnorm(z[1,],varcov=sigma)
for(i in 2:nrow(z)){a110=c(a110,pmnorm(z[i,1:2],varcov=sigma3)-pmnorm(z[i,],varcov=sigma))}
a100=pnorm(z[1,1],sd=s1)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a100=c(a100,pnorm(z[i,1],sd=s1)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a010=pnorm(z[1,2],sd=s2)-pmnorm(z[1,c(1,2)],varcov=sigma3)-pmnorm(z[1,c(2,3)],varcov=sigma1)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a010=c(a010,pnorm(z[i,2],sd=s2)-pmnorm(z[i,c(1,2)],varcov=sigma3)-pmnorm(z[i,c(2,3)],varcov=sigma1)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a001=pnorm(z[1,3],sd=s3)-pmnorm(z[1,c(2,3)],varcov=sigma1)-pmnorm(z[1,c(1,3)],varcov=sigma2)+pmnorm(z[1,],rep(0,ncol(Y)),varcov=sigma)
for(i in 2:nrow(z)){a001=c(a001,pnorm(z[i,3],sd=s3)-pmnorm(z[i,c(2,3)],varcov=sigma1)-pmnorm(z[i,c(1,3)],varcov=sigma2)+pmnorm(z[i,],rep(0,ncol(Y)),varcov=sigma))}
a000=1-a111-a011-a101-a110-a001-a010-a100

a111[a111<=0]=1e-50
a110[a110<=0]=1e-50
a101[a101<=0]=1e-50
a011[a011<=0]=1e-50
a100[a100<=0]=1e-50
a010[a010<=0]=1e-50
a001[a001<=0]=1e-50
a000[a000<=0]=1e-50

-sum(((Y[,1]==0)&(Y[,2]==0)&(Y[,3]==0))*log(a111) +
4*log(a011) +
5*log(a101) +
6*log(a110) +
7*log(a001) +
8*log(a010) +
9*log(a100) +
10*log(a000) )
}

A strong assumption in that bivariate model is that residuals have a Gaussian structure. It is possible to change that assumption

• marginally: for instance if we use a logistic cumulative distribution function, then we will have a bivariate logit regression
• in terms of dependence structure: it is possible to consider another copula than the gaussian one, e.g. Gumbel’s copula (also called the bivariate logistic copula), or Clayton’s

Here, the following code can be used to extend the model to non Gaussian structures,

> F=function(x,r){pmnorm(x,rep(0,length(x)),
+                 varcov=matrix(c(1,r,r,1),2,2))}
> Fx=function(x1){F(c(x1,1e40),0)}
> Fy=function(x2){Fx(x2)}
>
> logVgen=function(parameter){
+ CORRELATION=parameter[1]
+ BETA=matrix(parameter[2:length(parameter)],ncol(Y),ncol(X))
+ z=cbind(X%*%(BETA[1,]),X%*%(BETA[2,]))
+     a11=F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a11=c(a11,F(z[i,],r=CORRELATION))}
+     a10=Fx(z[1,1])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a10=c(a10,Fx(z[i,1])-F(z[i,],r=CORRELATION))}
+     a01=Fy(z[1,2])-F(z[1,],r=CORRELATION)
+ for(i in 2:nrow(z)){a01=c(a01,Fy(z[i,2])-F(z[i,],r=CORRELATION))}
+     a00=1-a10-a01-a11
+ -sum(((Y[,1]==1)&(Y[,2]==1))*log(a11) +
+     11*log(a01) +
+     12*log(a10) +
+     13*log(a00) )
+ }
>
> beta0=c(0,1,1,1,1,1,1)
> (OPT=optim(fn=logVgen,par=beta0,method="BFGS")$par) [1] 0.52613820 -2.45172059 2.24196154 4.90863292 -4.53994592 0.01600769 [7] 0.10614807 There were 23 warnings (use warnings() to see them) E.g. > library(copula) > F=function(x,r){pcopula(pnorm(x), claytonCopula(2, r))} > Fx=function(x1){F(c(x1,1e40),0) } > Fy=function(x2){Fx(x2)} • An application to school tests Consider the following dataset, hsb2=read.table("http://freakonometrics.free.fr/hsb2.csv", header=TRUE, sep=",") math_male=hsb2$math[female==0]
write_male=hsb2$write[female==0] math_female=hsb2$math[female==1]
write_female=hsb2$write[female==1] plot(math_female, write_female, type="p", pch=19,col="red",xlab="maths",ylab="writing",cex=.8) points(math_male, write_male, cex=1.2, col="blue") with here maths versus writing, with girls in red and boys in blue, where variables here are  female : 0: male 1: female race : 1: hispanic 2: asian 3: african-amer 4: white ses : 1: low 2: middle 3: high schtyp : type of school 1: public 2: private prog : type of program 1: general 2: academic 3: vocation read : reading score write : writing score math : math score science : science score socst : social studies score  We can try to understand correlation between math and writing skills. Covariates can be the sex of the child, and his reading skills. The question will then be: are good students in maths and writing simply students that can read well ? Here the code is simply > W=hsb2$write>=50
> M=hsb2$math>=50 > base=data.frame(Y1=W,Y2=M, + X1=hsb2$female,X2=hsb2$read) > > library(Zelig) > REG= zelig(list(mu1=Y1~X1+X2, + mu2=Y2~X1+X2, + rho=~1), + model="bprobit",data=base) > summary(REG) Call: zelig(formula = list(mu1 = Y1 ~ X1 + X2, mu2 = Y2 ~ X1 + X2, rho = ~1), model = "bprobit", data = base) Pearson Residuals: Min 1Q Median 3Q Max probit(mu1) -4.7518 -0.502594 0.15038 0.53038 1.8592 probit(mu2) -3.4243 -0.653537 0.23673 0.67011 2.6072 rhobit(rho) -4.9821 0.010481 0.13500 0.40776 2.9171 Coefficients: Value Std. Error t value (Intercept):1 -5.484711 0.787101 -6.96825 (Intercept):2 -4.061384 0.633781 -6.40818 (Intercept):3 1.332187 0.322175 4.13497 X1:1 1.125924 0.233550 4.82092 X1:2 0.167258 0.202498 0.82598 X2:1 0.103997 0.014662 7.09286 X2:2 0.082739 0.012026 6.88017 Number of linear predictors: 3 Names of linear predictors: probit(mu1), probit(mu2), rhobit(rho) Dispersion Parameter for binom2.rho family: 1 Residual Deviance: 364.51 on 593 degrees of freedom Log-likelihood: -182.255 on 593 degrees of freedom Number of Iterations: 3 > (exp(summary(REG)@coef3[3])-1)/(exp( summary(REG)@coef3[3])+1) [1] 0.5824045 with a remaining correlation among residuals of 0.58. So with only the sex of the student, and his or her reading skill, we cannot explain the correlation between maths and writing skills. With our previous code, we have here > beta0=c((exp(summary(REG)@coef3[3])-1)/(exp(summary(REG)@coef3[3])+1), + summary(REG)@coef3[c(1:2,4:7),1]) > beta0 (Intercept):1 (Intercept):2 X1:1 X1:2 0.58240446 -5.48471133 -4.06138412 1.12592427 0.16725842 X2:1 X2:2 0.10399668 0.08273879 > (OPT=optim(fn=logV,par=beta0,method="BFGS")$par)
(Intercept):1 (Intercept):2          X1:1          X1:2
0.5824045    -5.4847113    -4.0613841     1.1259243     0.1672584
X2:1          X2:2
0.1039967     0.0827388

i.e. we obtain (almost) exactly the same estimators. But here I have used as starting values for the optimization procedure the estimators given by R. If we change them, hopefully we have a robust maximum likelihood estimator,

> (OPT=optim(fn=logV,par=beta0/2,method="BFGS")\$par)
(Intercept):1 (Intercept):2          X1:1          X1:2
0.58233360   -5.49428984   -4.06839571    1.12696594    0.16760347
X2:1          X2:2
0.10417767    0.08287409
There were 12 warnings (use warnings() to see them)

So once again, it is possible to optimize numerically a likelihood function, and it works.

1. Y[,1]==0)&(Y[,2]==1 []
2. Y[,1]==1)&(Y[,2]==0 []
3. Y[,1]==0)&(Y[,2]==0 []
4. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==0 []
5. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==0 []
6. Y[,1]==0)&(Y[,2]==0)&(Y[,3]==1 []
7. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==0 []
8. Y[,1]==1)&(Y[,2]==0)&(Y[,3]==1 []
9. Y[,1]==0)&(Y[,2]==1)&(Y[,3]==1 []
10. Y[,1]==1)&(Y[,2]==1)&(Y[,3]==1 []
11. Y[,1]==0)&(Y[,2]==1 []
12. Y[,1]==1)&(Y[,2]==0 []
13. Y[,1]==0)&(Y[,2]==0 []