Category Archives: MAT8595

Conditional Distributions from some Elliptical Vectors

18/06/2014 Arthur Charpentier 6 Comments

This winter, in my ACT8595 course, I asked my students (that was some homework) to prove that it was possible to derive the conditional distribution when we have a Student-t random vector (and to get the analytical expression of the later). But before, let us recall a standard result about the Gaussian vector. If $\boldsymbol{X}=(\boldsymbol{X}_1,\boldsymbol{X}_2)$ is a Gaussian random vector, i.e.

$\begin{pmatrix}\boldsymbol{X}_1\\ \boldsymbol{X}_2\end{pmatrix}\sim \mathcal{N}\left(\begin{pmatrix}\boldsymbol{\mu}_1\\ \boldsymbol{\mu}_2\end{pmatrix},\left(\begin{array}{cc}\boldsymbol{\Sigma}_{11}& \boldsymbol{\Sigma}_{12}\\\boldsymbol{\Sigma}_{21}& \boldsymbol{\Sigma}_{22}\end{array}\right)\right)$

then $\boldsymbol{X}_1\vert \boldsymbol{X}_2=\boldsymbol{x}_2$ has a Gaussian distribution. More precisely, it is a $\mathcal{N}({\boldsymbol\mu}_{1\vert 2},{\boldsymbol\Sigma}_{1\vert 2})$ distribution, with

${\boldsymbol\mu}_{1\vert 2}=\boldsymbol\mu_1 + \boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1}\left(\boldsymbol{x}_2 - \boldsymbol\mu_2\right)$

and ${\boldsymbol\Sigma}_{1\vert 2}$ is the Schur complement of the block ${\boldsymbol\Sigma}_{11}$ of the matrix ${\boldsymbol\Sigma}$ ,

${\boldsymbol\Sigma}_{1\vert 2}=\boldsymbol\Sigma_{11} - \boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1} \boldsymbol\Sigma_{21}$

Observe that $\boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1}$ is also related to well known quantity: in the bivariate case, where $\boldsymbol{X}_{1}$ and $\boldsymbol{X}_{2}$ are univariate Gaussian variables,

$\boldsymbol\Sigma_{12} \boldsymbol\Sigma_{22}^{-1}=\frac{\text{cov}(X_1,X_2)}{\text{var}(X_2)}=\text{corr}(X_1,X_2)\sqrt{\frac{\text{var}(X_1)}{\text{var}(X_2)}}$

which is the slope in the linear regression of $X_1$ on $X_2$ .

In the case of the Student-t distribution, the conditional distrubution will not be a Student-t distribution anymore, but it will still be an elliptical distribution, and some interpretations of various quantities can actually be obtained.

The density of the multivariate centred Student-t distribution, with unit variance, and parameters $\nu$ and $\boldsymbol{R}$ is

$https://latex.codecogs.com/gif.latex?f(\boldsymbol{x})=%20\frac{\Gamma([d+\nu]/2)}{(\nu\pi)^{d/2}%20\Gamma(\nu/2)\vert\boldsymbol{R}\vert^{1/2}}%20\left(%201+\frac{1}{\nu}\boldsymbol{x}%27\boldsymbol{R}^{-1}\boldsymbol{x}%20\right)^{-(d+\nu)/2}$

If we consider the following blocks,

$https://latex.codecogs.com/gif.latex?\boldsymbol{R}=%20\left(%20\begin{array}{cc}%20\boldsymbol{R}_{11}&%20\boldsymbol{R}_{12}\\%20\boldsymbol{R}_{21}&%20\boldsymbol{R}_{22}%20\end{array}%20\right)$

then we can get that marginal distributions have a centred Student-t distribution, with unit variance, and parameters $\nu$ and $\boldsymbol{R}_{22}$ ,

$https://latex.codecogs.com/gif.latex?f_2(\boldsymbol{x}_2)=%20\frac{\Gamma([d_2+\nu]/2)}{(\nu\pi)^{d_2/2}%20\Gamma(\nu/2)\vert\boldsymbol{R}_{22}\vert^{1/2}}%20\left(%201+\frac{1}{\nu}\boldsymbol{x}_2%27\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_2%20\right)^{-(d_2+\nu)/2}$

Then, to derive the conditional density, we can use Bayes formula,

$https://latex.codecogs.com/gif.latex?f_{1\vert%202}(\boldsymbol{x}_1\vert%20\boldsymbol{x}_2)=%20\frac{f(\boldsymbol{x}_1,\boldsymbol{x}_2)}{f_2(\boldsymbol{x}_2)}$

One can write (as in Section 9.1 in Tong, 1990, The Multivariate Normal Distribution, but other expressions can be found in Section 2.5 in Fang, Ng and Kotz, 1989, Symmetric multivariate and related distributions, or in Section 1.11 in Kotz and Nadarajah, 2004, Multivariate t distributions and their applications) this conditional density as

$https://latex.codecogs.com/gif.latex?f_{1\vert%202}(\boldsymbol{x}_1\vert%20\boldsymbol{x}_2)=\kappa%20\left(1+\frac{1}{\nu}\boldsymbol{x}_2%27\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_2\right)^{(d_2+\nu)/2}%20\left(1+\frac{1}{\nu}\left[\boldsymbol{x}_2%27\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_2+\alpha(\boldsymbol{x}_1,\boldsymbol{x}_2)\right]\right)^{-(d_1+\nu)/2}$

with

$https://latex.codecogs.com/gif.latex?\kappa=\frac{\Gamma([d+\nu]/2)}{(\nu\pi)^{d_1/2}%20\Gamma([d_2+\nu]/2)}\frac{1}{\vert\boldsymbol{R}_{11}-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{R}_{21}\vert^{1/2}}$

and

$https://latex.codecogs.com/gif.latex?\alpha(\boldsymbol{x}_1,\boldsymbol{x}_2)=(\boldsymbol{x}_1-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_{2})%27%20[\boldsymbol{R}_{11}-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{R}_{21}]^{-1}(\boldsymbol{x}_1-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_{2})$

This conditional distribution is elliptical, ~~but it is not a Student-t distribution, except in the case where $\boldsymbol{x}_2\in\{(\pm1,\cdots,\pm 1)\}$ , or when the correlation matrix $\boldsymbol{R}$~~ ~~is the identity~~.

[June 2016] actually, as shown in Ding (2016), this is a Student-t distribution. “Kotz & Nadarajah (2004) and Nadarajah & Kotz (2005) failed to recognize that the conditional distribution of the MVT distribution is also a MVT distribution due to the complexity of the conditional density function […] Conditional distributions of elliptically contoured distributions are also elliptically contoured distributions. But this does not immediately guarantee that conditional distributions of the MVT distributions are also MVT distributions without some further algebra.“

Now, if we look at the components of this density, we can observe that we have

$https://latex.codecogs.com/gif.latex?(\boldsymbol{x}_1-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{x}_{2})$

which was mentioned previously, in the Gaussian case: the term on the right is the conditional mean,

$\mathbb{E}(\boldsymbol{X}_1\vert \boldsymbol{X}_2=\boldsymbol{x}_2)=\boldsymbol{\mu}_{1\vert 2}$

and the bloc that appears at several places is the conditional variance,

$\text{Var}(\boldsymbol{X}_1\vert \boldsymbol{X}_2=\boldsymbol{x}_2)=\boldsymbol{R}_{1\vert 2}=\boldsymbol{R}_{11}-\boldsymbol{R}_{12}\boldsymbol{R}_{22}^{-1}\boldsymbol{R}_{21}$

Now, if we want to visualize that conditional density, let us plot it. The code below is based on Bayes formula

> library(mnormt)
> r=.6
> R=matrix(c(1,r,r,1),2,2)
> nu=4
> f2=function(x2) dt(x2,df=nu)
> f =function(x) dmt(x,S=R,df=nu)
> f1.2=function(x1,x2) f(c(x1,x2))/f2(x2)

In order to compare that conditional density with a Student-t one, let us define the density of a non-centred Student-t random variable,

> dstd=function(x,mu,s,nu) gamma((nu+1)/2)/
+ (gamma(nu/2)*s*sqrt(pi*nu))*
+ (1+1/nu*(x-mu)^2/(s^2))^(-(nu+1)/2)

Here is the function we can use to plot those two densities,

> graphdensity=function(x2=-1.5){
+ vectx1=seq(-3,3,length=251)
+ y=Vectorize(function(x) f1.2(x,x2))(vectx1)
+ plot(vectx1,y,type="l",col="red",ylim=c(0,.5),
+ xlab="",ylab="")
+ abline(v=r*x2,lty=2)
+ lines(vectx1,dstd(vectx1,x2*r,sqrt(1-r^2),nu),col="blue",lty=2)}
> graphdensity(-1.5)

In the case where $x_2=\pm1$ , the two lines are rather close (the difference migth come from computational issues)

> graphdensity(-1)

and just to conclude, a last one

> graphdensity(0)

Copulas, MAT8595, Maths, Probability

On Hoeffding’s identity

17/06/2014 Arthur Charpentier 1 Comment

In 1940, Wassily Hoeffding published Masstabinvariante Korrelationstheorie, which was an impressive paper. For those (like me) who unfortunately barely speak German, an English translation could be found in The Collected Works of Wassily Hoeffding, published a few years ago. As I keep saying in my courses about copulas, almost everything was in that paper, by Wassily Hoeffding. For instance, we can see the following graph, of a cumulative distribution function,

What is the difference with a copula? A copula (in dimension 2) is the cumulative distribution function of a random pair with uniform on $[0,1]\times[0,1]$ , as defined by Abe Sklar

But Wassily Hoeffding considered a random pair with uniform on $[-1/2,+1/2]\times[-1/2,+1/2]$ . But everything else is the same. He can even derive the level curves of the density of the Gaussian copula,

$c(u,v)=\frac{\varphi_r(\Phi^{-1}(u),\Phi^{-1}(v))}{\varphi(\Phi^{-1}(u))\cdot \varphi(\Phi^{-1}(v))}$

> library(mnormt)
> r=.6
> dc=function(u,v) return(
+ as.numeric(dmnorm(cbind(qnorm(u),qnorm(v)),varcov=
+ matrix(c(1,r,r,1),2,2))/dnorm(qnorm(u))/dnorm(qnorm(v))))
> n=500
> vectu=seq(1/n,1-1/n,length=n-1)
> matdc=outer(vectu,vectu,dc)
> contour(vectu,vectu,matdc,levels=
+ c(.325,.944,1.212,1.250,1.290,1.656,3.85),lwd=2)

But another interesting point is that there is the so-called Hoeffding’s equality

$\text{cov}(X,Y)=\int_{\mathbb{R}\times\mathbb{R}} [F_{XY}(x,y) - F_X(x)F_Y(y)]dxdy$

which is interesting, and quite important, actually, to understand that the covariance (or the correlation) can be seen as some ‘distance‘ to the independence. More precisely, observe that

$\text{cov}(X,Y)=\int_{\mathbb{R}\times\mathbb{R}} [F_{XY}(x,y) - F_{X,Y}^\perp(x,y)]dxdy$

where $F_{X,Y}^\perp$ would be the joint cumulative distribution function of some independent variables, with the same marginal distributions.

Of course, it is not exactly a distance, since it can be negative. But still. Now, the thing is that the proof is not trivial. But it is using interesting identities. For instance, in 1885, Franklin wrote a nice paper, Proof of a Theorem of Tchebycheff’s on Definite Integrals, in the American Journal of Mathematics. To get some heuristics about the identity, consider some (finite) sequences $(\boldsymbol{u})=\{u_1,\cdots,u_n\}$ and $(\boldsymbol{v})=\{v_1,\cdots,v_n\}$ , then one can prove that

$\sum_{i,j=1}^n [(u_i-u_j)(v_i-v_j)]=2\left[n\sum_{i=1}^n u_iv_i -\sum_{i=1}^n u_i\cdot \sum_{i=1}^n v_i\right]$

And there is a continuous version of that identity. Consider two bounded functions $u(\cdot)$ and $v(\cdot)$ , on some interval, $[a,b]$ then

$\int_{[a.b]\times[a,b]}[(u(x)-u(y))\cdot(v(x)-v(y))]dxdy$

is equal to

$2\left[(b-a)\int_{[a,b]} u(x)v(x)dx -\int_{[a,b]} u(x)dx\cdot \int_{[a,b]} v(x)dx\right]$

In 1979, in Monotone Regression and Covariance Structure, Gerald Shea gave a more probabilistic interpretation of that results, using the fact that

$(b-a)=\int_{[a,b]}1dx$

and using a different measure. More precisely, assume now that functions $u(\cdot)$ and $v(\cdot)$ are integrable, with respect to some measure $\mu$ , on some set $\mathcal{S}\subset\mathbb{R}$ . Then

$\int_{\mathcal{S}\times\mathcal{S}}[(u(x)-u(y))\cdot(v(x)-v(y))]\mu(dx)\mu(dy)$

is equal to

$2\left[\mu(\mathcal{S})\int_{\mathcal{S}} u(x)\cdotv(x)\mu(dx) -\int_{\mathcal{S}} u(x)\mu(dx)\cdot \int_{\mathcal{S}} v(x)\mu(dx)\right]$

In the case where $\mu$ is a probability measure of $\mathcal{S}\subset\mathbb{R}$ , i.e. $\mu(\mathcal{S})=1$ , this equality is the one used by Wassily Hoeffding, in 1940. The interpretation in terms of random variable is simple that

$\mathbb{E}[(X_1-X_2)\cdot(Y_1-Y_2)]=2\left[\mathbb{E}(X_1Y_1)-\mathbb{E}(X_1)\mathbb{E}(Y_1)\right]$

(with standard assuptions of existence of those quantitites) where $(X_1,Y_1)$ and $(X_2,Y_2)$ are two independent vectors, with identical distribution, $F_{XY}$ . Actually, this relationship can also be found in Some Concepts of Dependence, by E. L. Lehmann, published in 1966. Oh, and by the way, the connection with Chebyshev inequality (claimed in the title of seminal paper by Franklin) come from the fact that if $u(\cdot)$ and $v(\cdot)$ are monotonic, then the left part of the identity is positive, and thus,

$(b-a)\int_{[a,b]} u(x)\cdot v(x)dx \geq\int_{[a,b]} u(x)dx\cdot \int_{[a,b]} v(x)dx$

But let’s get back to Hoeffiding’s result. How do we get it from that lemma. The idea is to write

$\mathbb{E}[(X_1-X_2)(Y_1-Y_2)]$

$\mathbb{E}\left(\int_{\mathbb{R}\times\mathbb{R}} [\boldsymbol{1}_{u\leq X_1}-\boldsymbol{1}_{u\leq X_2}]\cdot[\boldsymbol{1}_{v\leq Y_1}-\boldsymbol{1}_{v\leq Y_2}] dudv\right)$

i.e.

$\mathbb{E}\left(\int_{\mathbb{R}\times\mathbb{R}}[\boldsymbol{1}_{u\leq X_1}\boldsymbol{1}_{u\leq Y_1}-\boldsymbol{1}_{u\leq X_1}\boldsymbol{1}_{u\leq Y_2}-\boldsymbol{1}_{u\leq X_2}\boldsymbol{1}_{u\leq Y_1}+\boldsymbol{1}_{u\leq X_2}\boldsymbol{1}_{u\leq Y_2}]dudv\right)$

We can then intervert the integral and the expectation, use the fact that

$\mathbb{E}(\boldsymbol{1}_{u\leq X_1})=\mathbb{P}(X_1\geq u)$

and then, and some integral calculus can be used to rewrite that expression as

$2\left[\int_{\mathbb{R}\times\mathbb{R}}[F_{X_1Y_1}(u,v)-F_{X_1}(u)F_{Y_1}(v)]dudv\right]$

So we get here Hoeffding’s identity. Actually, as mentioned by Ben Derrett about the equality above, it can be observed (see http://math.stackexchange.com/105713) that $2\text{cov}(X,Y)=2\big(\mathbb{E}[XY]-\mathbb{E}[X]\mathbb{E}[Y]\big)$ can also be written

$\mathbb E[X_1Y_1]-\mathbb E[X_2]\mathbb E[Y_1]+\mathbb E[X_2Y_2]-\mathbb E[X_1]\mathbb E[Y_2]$

where again, $(X_1,Y_1)$ and $(X_2,Y_2)$ are two independent vectors, with identical distribution, $F_{XY}$ . The later can be writen

$\mathbb E((X_1-X_2)(Y_1-Y_2))$

Copulas, MAT8595

Examen, Copules et Valeurs Extrêmes

26/04/2014 Arthur Charpentier Leave a comment

Après les exposés des dernières séances, l’examen du cours MAT8595, Copules et Valeurs Extrêmes avait lieu hier matin. L’énoncé est en ligne, et j’ai aussi écrit quelques éléments de correction, là. En cas de désaccord, merci de me le faire savoir rapidement ! Sinon, bonne fin de session..

Copulas, Inference, MAT8595

Modeling the Marginals and the Dependence separately

01/04/2014 Arthur Charpentier 3 Comments

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,

MAT8181, MAT8595

Independence and correlation

29/03/2014 Arthur Charpentier Leave a comment

A short post to get back on a property I gave briefly in the MAT8595 class in January, and again in the MAT8181 class this week (to illustrate the distinction between weak and strong white noises). Recall that (real-valued) random variables $X$ and $Y$ are independent if for all $A, B\subset\mathbb{R}$ , $\mathbb{P}((X,Y)\in A\times B)=\mathbb{P}(X\in A)\cdot \mathbb{P}(Y\in B)$ Another characterization, for integrable variable is that $\mathbb{E}(f(X)\cdot g(Y)) = \mathbb{E}(f(X))\cdot \mathbb{E}( g(Y))$ for all $f,g:\mathbb{R}\mapsto \mathbb{R}$ , which can be written, if variables are square integrable $\text{corr}(f(X), g(Y))=0$ The idea to prove this characterization is to observe that if $X$ and $Y$ are independent $\mathbb{P}((X,Y)\in A\times B)=\mathbb{P}(X\in A)\cdot \mathbb{P}(Y\in B)$ can be written $\mathbb{E}(\boldsymbol{1}_A(X)\cdot \boldsymbol{1}_B(Y)) = \mathbb{E}(\boldsymbol{1}_A(X))\cdot \mathbb{E}( \boldsymbol{1}_B(Y))$ Using a standard argument in integration theory, equality $\mathbb{E}(f(X)\cdot g(Y)) = \mathbb{E}(f(X))\cdot \mathbb{E}( g(Y))$ is valid for step functions (not only indicators), and then to positive measurable functions, and finally to integrable functions. Proving this result is not that difficult. Observe that Rényi (1959) – inspired by Gebelein (1947) – followed by Sarmanov (1958) introduced the concept of maximal correlation, that can be related to this result, $\text{corr}^{\star}(X,Y)=\max\{\text{corr}(f(X),g(Y))\}$ where the maximum is taken over all functions $f$ and $g$ such that the correlation exist. Actually, it is possible to consider only transformations such that $\mathbb{E}(f(X))=0$ and $\text{Var}(f(X))=1$ (and similarly for $g$ , the idea is that we simple center and scale, which does not impact the correlation.Thus, $X$ and $Y$ are independent if and only if $\text{corr}^{\star}(X,Y)=0$ Algorithm to estimate that coefficient are interesting. The problem can be written, equivalently $(f^\star,g^\star)=\text{argmin}\{\mathbb{E}[f(X)-g(Y)]^2\}$ And if the minimization is considered over $g$ , assuming that $f$ is fixed, then the optimal transformation is $g(y)=\mathbb{E}[f(X)\vert Y=y]$ And similarly for $f$ . So using an iterative algorithm, it is possible to get $f^\star$ and $g^\star$ (see Breiman and Friedman (1985) for more details). Actually, those functions appear in nonlinear canonical analysis. As mentioned in Lancaster (1957), for a Gaussian random vector $(X,Y)$ $\text{corr}^{\star}(X,Y)=\vert \text{corr}(X,Y)\vert$ and in that case $f^\star$ and $g^\star$ are affine functions. This can be related to Hermite’s polynomial and to the expansion of the bivariate Gaussian density. I still hope that someone will go further for the project in the MAT8181 course.

Copulas, MAT8595, Multivariate Analysis, Probability

Bivariate Densities with N(0,1) Margins

18/02/2014 Arthur Charpentier 1 Comment

This Monday, in the ACT8595 course, we came back on elliptical distributions and conditional independence (here is an old post on de Finetti’s theorem, and the extension to Hewitt-Savage’s). I have shown simulations, to illustrate those two concepts of dependent variables, but I wanted to spend some time to visualize densities. More specifically what could be the joint density is we assume that margins are $\mathcal{N}(0,1)$ distributions.

The Bivariate Gaussian distribution

Here, we consider a Gaussian random vector, with margins $\mathcal{N}(0,1)$ , and with correlation $r\in[-1,+1]$ . This is the standard graph, with elliptical isodensity curves

r=.5
library(mnormt)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmnorm(cbind(x,y),varcov=S)
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
X=rmnorm(1500,varcov=S)
xhist <- hist(X[,1], plot=FALSE)
yhist <- hist(X[,2], plot=FALSE)
top <- max(c(xhist$density, yhist$density,dnorm(0)))
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))
image(vx,vy,z,col=rev(heat.colors(101)))
contour(vx,vy,z,col="blue",add=TRUE)
points(X,cex=.2)
par(mar=c(0,3,1,1))
barplot(xhist$density, axes=FALSE, ylim=c(0, top), space=0,col="light green")
lines((density(X[,1])$x-xhist$breaks[1])/diff(xhist$breaks)[1],
dnorm(density(X[,1])$x),col="red")
par(mar=c(3,0,1,1))
barplot(yhist$density, axes=FALSE, xlim=c(0, top), space=0, 
horiz=TRUE,col="light green")
lines(dnorm(density(X[,2])$x),(density(X[,2])$x-yhist$breaks[1])/
diff(yhist$breaks)[1],col="red")

That was the simple part.

The Bivariate Student-t distribution

Consider now another elliptical distribution. But we want here to normalize the margins. Thus, instead of a pair $(X,Y)$ , we would like to consider the pair $(\Phi^{-1}(T_\nu(X)),\Phi^{-1}(T_\nu(Y)))$ , so that the marginal distributions are $\mathcal{N}(0,1)$ . The new density is obtained simply since the transformation is a one-to-one increasing transformation. Here, we have

k=3
r=.5
G=function(x) qnorm(pt(x,df=k))
dg=function(x) dt(x,df=k)/dnorm(qnorm(pt(x,df=k)))
Ginv=function(x) qt(pnorm(x),df=k)
S=matrix(c(1,r,r,1),2,2)
f=function(x,y) dmt(cbind(Ginv(x),Ginv(y)),S=S,df=k)/(dg(x)*dg(y))
vx=seq(-3,3,length=201)
vy=seq(-3,3,length=201)
z=outer(vx,vy,f)
set.seed(1)
Z=rmt(1500,S=S,df=k)
X=G(Z)

Because we considered a nonlinear transformation of the margins, the level curves are no longer elliptical. But there is still some kind of symmetry.

The Exchangeable Case with Conditionally Independent Random Variables

We did consider the case where $X$ and $Y$ with independent random variables, given $\Theta$ , and that both variables are exponentially distributed, with parameter $\Theta$ . As we’ve seen in class, it might be difficult to visualize that sample, unless we have log scales on both axis. But instead of a log transformation, why not consider a transformation so that margins will be $\mathcal{N}(0,1)$ . The only technical problem is that we do not have the (nonconditional) distributions of the margins. Well, we have them, but they are integral based. From a computational point of view, that’s not a bit deal… Computations might take a while, but we can visualize the density using the following code (here, we assume that is Gamma distributed)

a=.6
b=1
h=.0001
G=function(x) qnorm(ifelse(x<0,0,integrate(function(z) pexp(x,z)*
dgamma(z,a,b),lower=0,upper=Inf)$value))
Ginv=function(x) uniroot(function(z) G(z)-x,lower=-40,upper=1e5)$root
dg=function(x) (Ginv(x+h)-Ginv(x-h))/2/h
H=function(xy) integrate(function(z) dexp(xy[2],z)*dexp(xy[1],z)*
dgamma(z,a,b),lower=0,upper=Inf)$value
f=function(x,y) H(c(Ginv(x),Ginv(y)))*(dg(x)*dg(y))
vx=seq(-3,3,length=151)
vy=seq(-3,3,length=151)
z=matrix(NA,length(vx),length(vy))
for(i in 1:length(vx)){
for(j in 1:length(vy)){
z[i,j]=f(vx[i],vy[j])}}
set.seed(1)
Theta=rgamma(1500,a,b)
Z=cbind(rexp(1500,Theta),rexp(1500,Theta))
X=cbind(Vectorize(G)(Z[,1]),Vectorize(G)(Z[,2]))

There is a small technical problem, but no big deal.

Here, the joint distribution is quite different. Margins are – one more time – standard Gaussian, but the shape of the joint distribution is quite different, with an asymmetry from the lower (left) tail to the upper (right) tail. More details when we’ll introduce copulas. The only difference will be that the margins will be uniform on the unit interval, and not standard Gaussian.

Extreme Value, Inference, MAT8595

Risk Measures with Extreme Value Models

05/02/2014 Arthur Charpentier Leave a comment

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 $https://latex.codecogs.com/gif.latex?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 $https://latex.codecogs.com/gif.latex?x\geq%20u$ , 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 $https://latex.codecogs.com/gif.latex?\{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

$https://latex.codecogs.com/gif.latex?%20%20%20%20%20g_{(\xi,\sigma)}(x)%20=%20\frac{1}{\sigma}\left(1%20+%20\frac{\xi%20x}{\sigma}\right)^{\left(-\frac{1}{\xi}%20-%201\right)}$

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

$https://latex.codecogs.com/gif.latex?g_{\xi,Q(p)}(x)=\frac{\displaystyle{\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{\xi[Q(p)-u]}\left(1+\frac{\displaystyle{\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{[Q(p)-u]}\cdot%20x\right)^{-\frac{1}{\xi}-1}$

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

$https://latex.codecogs.com/gif.latex?\log\mathcal{L}(\xi,Q(p);\boldsymbol{x})=\sum_{i=0}^{N_u-1}%20\log%20g_{\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 $https://latex.codecogs.com/gif.latex?ES(p)$ instead of $https://latex.codecogs.com/gif.latex?\sigma$ . Here, we get

$https://latex.codecogs.com/gif.latex?g_{\xi,ES(p)}(x)=\frac{\displaystyle{\xi+\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{\xi(1-\xi)[ES(p)-u]}\left(1+\frac{\displaystyle{\left(\frac{n}{N_u}(1-p)\right)^{-\xi}-1}}{(1-\xi)[ES(p)-u]}\cdot%20x\right)^{-\frac{1}{\xi}-1}$

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

Extreme Value, MAT8595

Bias of Hill Estimators

28/01/2014 Arthur Charpentier 3 Comments

In the MAT8595 course, we’ve seen yesterday Hill estimator of the tail index. To be more specific, we did see see that if $https://latex.codecogs.com/gif.latex?\overline{F}(x)=C%20x^{-\alpha}$ , with $https://latex.codecogs.com/gif.latex?\alpha%3E0$ , then Hill estimators for $https://latex.codecogs.com/gif.latex?\alpha$ are given by

$https://latex.codecogs.com/gif.latex?\widehat{\alpha}_k%20=%20\left[\frac{1}{k}\sum_{i=0}^{k-1}%20\log%20X_{n-i,n}%20-\log%20X_{n-k,n}\right]^{-1}$
for $https://latex.codecogs.com/gif.latex?k\in\{1,2,\cdots,n\}$ . Then we did say that $https://latex.codecogs.com/gif.latex?\widehat{\alpha}_k$ satisfies some consistency in the sense that $https://latex.codecogs.com/gif.latex?\widehat{\alpha}_k%20\overset{\mathbb{P}}{\rightarrow}%20\alpha$ if $https://latex.codecogs.com/gif.latex?k\rightarrow\infty$ , but not too fast, i.e. $https://latex.codecogs.com/gif.latex?k/n\rightarrow0$ (under additional assumptions on the rate of convergence, it is possible to prove that $https://latex.codecogs.com/gif.latex?\widehat{\alpha}_k%20\overset{a.s.}{\rightarrow}%20\alpha$ ). Further, under additional technical conditions

$https://latex.codecogs.com/gif.latex?\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal%20L}{\rightarrow}\mathcal{N}(0,\alpha^2)$

In order to illustrate this point, consider the following code. First, let us consider a Pareto survival function, and the associated quantile function

> alpha=1.5
> S=function(x){ifelse(x>1,x^(-alpha),1)}
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}

The code here is obviously too complicated, since this power function can easily be inverted. But later on, we will consider a more complex survival function. Here are the survival function, and the quantile function,

> u=seq(0,5,by=.01)
> plot(u,Vectorize(S)(u),type="l",col="red")
> u=seq(0,99/100,by=.01)
> plot(u,Vectorize(Q)(u),type="l",col="blue",ylim=c(0,20))

Here, we need the quantile function to generate a random sample from this distribution,

> n=500
> set.seed(1)
> X=Vectorize(Q)(runif(n))

Hill plot is here

> library(evir)
> hill(X)
> abline(h=alpha,col="blue")

We can now generate thousands of random samples, and see how those estimators behave (for some specific $https://latex.codecogs.com/gif.latex?k$ ‘s).

> ns=10000
> HillK=matrix(NA,ns,10)
> for(s in 1:ns){
+ X=Vectorize(Q)(runif(n))
+ H=hill(X,plot=FALSE)
+ hillk=function(k) H$y[H$x==k]
+ HillK[s,]=Vectorize(hillk)(15*(1:10))
+ }

and if we compute the average,

> plot(15*(1:10),apply(HillK,2,mean)

we do get a series of estimators that can be considered as unbiased.

So far, so good. Now, recall that being in the max-domain of attraction of the Fréchet distribution does not mean that $https://latex.codecogs.com/gif.latex?\overline{F}(x)=C%20x^{-\alpha}$ , with $https://latex.codecogs.com/gif.latex?\alpha%3E0$ , but is means that

$https://latex.codecogs.com/gif.latex?\overline{F}(x)=%20x^{-\alpha}%20\mathcal{L}(x)$

for some slowly varying function $https://latex.codecogs.com/gif.latex?\mathcal{L}$ , not necessarily constant! In order to understand what could happen, we have to be slightly more specific. And this can be done only by looking at second order regular variation property of the survival function. Assume, here that there is some auxilary function $https://latex.codecogs.com/gif.latex?a$ such that

$https://latex.codecogs.com/gif.latex?\lim_{t\rightarrow\infty}\frac{\overline{F}(xt)/\overline{F}(t)-x^{-\alpha}}{a(t)}=x^{-\alpha}\frac{1-x^{-\beta}}{\beta}{}$

This (positive) constant $https://latex.codecogs.com/gif.latex?\beta$ is – somehow – related to the speed of convergence of the ratio of the survival functions to the power function (see e.g. Geluk et al. (2000) for some examples).

To be more specific, assume that

$https://latex.codecogs.com/gif.latex?\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot%20%20x^{-\alpha}$

then, the second order regular variation property is obtained using $https://latex.codecogs.com/gif.latex?a(t)=\beta%20t^{-\beta}$ , and then, if $https://latex.codecogs.com/gif.latex?k$ goes to infinity too fast, then the estimator will be biased. More precisely (see Chapter 6 in Embrechts et al. (1997)), if $https://latex.codecogs.com/gif.latex?k=O(n^{2\beta/(\alpha+2\beta)})$ , then, for some $https://latex.codecogs.com/gif.latex?\lambda%3E0$ ,

$https://latex.codecogs.com/gif.latex?\sqrt{k}\left(\widehat{\alpha}_k-\alpha\right)\overset{\mathcal%20L}{\rightarrow}\mathcal{N}\left(\frac{\alpha^3}{\beta-\alpha}\lambda,\alpha^2\right)$

The intuitive interpretation of this result is that if $https://latex.codecogs.com/gif.latex?k$ is too large, and if the underlying distribution is not exactly a Pareto distribution (and we do have this second order property), then Hill’s estimator is biased. This is what we mean when we say

if $https://latex.codecogs.com/gif.latex?k$ is too large, $https://latex.codecogs.com/gif.latex?\widehat{\alpha}_k$ is a biased estimator
if $https://latex.codecogs.com/gif.latex?k$ is too small, $https://latex.codecogs.com/gif.latex?\widehat{\alpha}_k$ is a volatile estimator

(the later comes from properties of a sample mean: the more observations, the less the volatility of the mean).

Let us run some simulations to get a better understanding of what’s going on. Using the previous code, it is actually extremly simple to generate a random sample with survival function

$https://latex.codecogs.com/gif.latex?\overline{F}(x)=\underbrace{C(1+x^{-\beta})}_{\mathcal{L}(x)}\cdot%20%20x^{-\alpha}$

> beta=.5
> S=function(x){+ ifelse(x>1,.5*x^(-alpha)*(1+x^(-beta)),1) }
> Q=function(p){uniroot(function(x) S(x)-(1-p),lower=1,upper=1e+9)$root}

If we use the code above. Here, with

> n=500
> set.seed(1)
> X=Vectorize(Q)(runif(n))

the Hill plot becomes

> library(evir)
> hill(X)
> abline(h=alpha,col="blue")

But it’s based on one sample, only. Again, consider thousands of samples, and let us see how Hill’s estimator is behaving,

so that the (empirical) mean of those estimator is

Extreme Value, Inference, MAT8595, Statistics

Likelihood Based Methods, for Extremes

22/01/2014 Arthur Charpentier 1 Comment

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 $https://latex.codecogs.com/gif.latex%20?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 $https://latex.codecogs.com/gif.latex?\boldsymbol{\theta}=(\xi,\sigma)$ , so that

$https://latex.codecogs.com/gif.latex?%20%20%20%20%20F_{(\xi,\sigma)}(x)%20=%20\begin{cases}%201%20-%20\left(1+%20\frac{\xi%20x}{\sigma}\right)^{-1/\xi}%20&,%20\xi%20\neq%200%20\\%201%20-%20\exp%20\left(-\frac{x}{\sigma}\right)%20&,%20\xi%20=%200%20\end{cases}$

If we solve (numerically) the first order condition of the maximum likelihood, we get an estimator $https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma%20}_n)$ which satisfies

$https://latex.codecogs.com/gif.latex?\sqrt{n}\left(\left[\begin{array}{c}\widehat{\xi}_n\\\widehat{\sigma%20}_n\end{array}\right]-\left[\begin{array}{c}\xi_0\\\sigma_0%20\end{array}\right]\right)\rightarrow%20\mathcal{N}\left(\left[\begin{array}{c}0\\end{array}\right],\left[\begin{array}{cc}(1+\xi_0)^2%20&%20\sigma_0[1+\xi_0]\\%20\sigma_0%20[1+\xi_0]%20&%202\sigma^2_0(1+\xi_0)%20\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 $https://latex.codecogs.com/gif.latex%20?n$ is large enough, then $https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma%20}_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 $https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\theta}}_n=(\widehat{\xi}_n,\widehat{\sigma%20}_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 $https://latex.codecogs.com/gif.latex%20?\widehat{\boldsymbol{\theta}}_n$ is asymptotically normal, and if $h:\Theta\subset\mathbb R^k \rightarrow \mathbb R^m$ is sufficently smooth, then $https://latex.codecogs.com/gif.latex%20?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 $https://latex.codecogs.com/gif.latex%20?\widehat{\alpha}_n=\widehat{\xi}_n^{-1}$ (which is another parametrization used in extreme value models), or on any quantile, $https://latex.codecogs.com/gif.latex%20?\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 $https://latex.codecogs.com/gif.latex%20?\widehat{\alpha}_n$ and $https://latex.codecogs.com/gif.latex%20?\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 $https://latex.codecogs.com/gif.latex%20?\xi$ , $https://latex.codecogs.com/gif.latex%20?\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 $https://latex.codecogs.com/gif.latex%20?\xi$ , we have to find an optimal $https://latex.codecogs.com/gif.latex%20?\sigma^\star(\xi)$ . We compute the (profile) likelihood i.e. $https://latex.codecogs.com/gif.latex%20?\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 $https://latex.codecogs.com/gif.latex%20?\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 $https://latex.codecogs.com/gif.latex%20?\xi$ .

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

MAT8595

Central Limit Theorem

06/01/2014 Arthur Charpentier Leave a comment

This week, in the MAT8595 course, before proving Fisher-Tippett theorem, we will get back on the proof of the Central Limit Theorem, and the class of stable distribution (in Lévy’s sense). In order to illustrate the problem of heavy tails, on the behavior of the mean, consider a sequence of i.i.d. Gaussian random variables $https://latex.codecogs.com/gif.latex?X_i$ ‘s. On top, we visualize the sequence, and below, we visualize the associate random walk

$https://latex.codecogs.com/gif.latex?S_n=\sum_{i=1}^n%20X_i$

(the central limit theorem will give a limiting distribution for $https://latex.codecogs.com/gif.latex?n^{-1}S_n$ in the case where the variance of the $https://latex.codecogs.com/gif.latex?X_i$ ‘s is finite)

If we consider a sequence of i.i.d. random variables $https://latex.codecogs.com/gif.latex?X_i$ ‘s whith heavier tails (possibly with infinite variance), we can still define $https://latex.codecogs.com/gif.latex?S_n$ , but as we can see below, $https://latex.codecogs.com/gif.latex?S_n$ can be quite heratic.

As we will see this Thursday, the key to derive stable distribution for the central limit theorem, or possible limiting distributions for the maximum is Cauchy’s function equation. I strongly recommand to look at the proof.

MAT8595

Copules et valeurs extrêmes, syllabus

05/01/2014 Arthur Charpentier Leave a comment

Le plan de cours pour le cours MAT8595 Copules et Valeurs Extremes est en ligne. L’entente d’évaluation sera signée au premier cours, ce lundi à 9:00 (salle SH-2140). D’autre billets seront mis en ligne dans les jours à venir, avec quelques exercices, et les articles qui serviront de base pour les projets, sur http://freakonometrics.hypotheses.org/courses/copulas-and-extremes.

Copulas, MAT8595

Graduate Course on Copulas and Extreme Values

18/11/2013 Arthur Charpentier 1 Comment

This Winter, I will be giving a (graduate) course on extreme values, and copulas (more generally multivariate models and dependence), MAT8595. It is an ISM course, and even if it will probably be given in French, I will upload information here, in English. I will upload the (detailed) syllabus of the course during the Christmas holidays. But to give an overview, for those willing to register, the first part of the course will focus on extreme value theory. The references will be

Beirlant, J, Goegebeur, Y., Segers, S. & Teugels, J. (2004) Statistics of Extremes. Wiley.
Embrechts, P. Kluppelberg, C. & Mikosch, T. (1997) Modelling Extremal Events. Springer Verlag.
Resnick, S. (1987). Extreme Value, Regular Variation and Point Processes. Springer Verlag.

The second part of the course will be on multivariate distributions. The references will be

Joe, H. (1997). Multivariate Models and Multivariate Dependence Concepts. Chapman & Hall.
Nelsen, R. (1999). An Introduction to Copulas. Springer Verlag.

Specific references and more details about the chapters will be given during the course. I will upload exercises this winter, as well as a list of articles that will be used for projects. Examples will be illustrated using R functions from dedicated packages.

Grades will be based on exercises (homework), report (based on a published paper) and final writen exam.

Freakonometrics