# Freakonometrics

## Tag - copulas

Wednesday, October 17 2012

## Fractals and Kronecker product

A few years ago, I went to listen to Roger Nelsen who was giving a talk about copulas with fractal support. Roger is amazing when he gives a talk (I am also a huge fan of his books, and articles), and I really wanted to play with that concept (that he did publish later on, with Gregory Fredricks and José Antonio Rodriguez-Lallena). I did mention that idea in a paper, writen with Alessandro Juri, just to mention some cases where deriving fixed point theorems is not that simple (since the limit may not exist).

The idea in the initial article was to start with something quite simple, a the so-called transformation matrix, e.g.

$T=\frac{1}{8}\left(\begin{matrix}1& 0 & 1 \\ 0 & 4 & 0 \\ 1 & 0&1\end{matrix}\right)$

Here, in all areas with mass, we spread it uniformly (say), i.e. the support of $T(C^\perp)$ is the one below, i.e. $1/8$th of the mass is located in each corner, and $1/2$ is in the center. So if we spread the mass to have a copula (with uniform margin,)we have to consider squares on intervals $[0,1/4]$, $[1/4,3/4]$ and $[3/4,1]$,

Then the idea, then, is to consider $T^2=\otimes^2T$, where  $\otimes^2T$ is the tensor product (also called Kronecker product) of $T$ with itself. Here, the support of $T^2(C^\perp)$ is

Then, consider $T^3=\otimes^3T$, where $\otimes^3T$ is the tensor product of $T$ with itself, three times. And the support of $T^3(C^\perp)$ is

Etc. Here, it is computationally extremely simple to do it, using this Kronecker product. Recall that if $\mathbf{A}=(a_{i,j})$, then

$\mathbf{A}\otimes\mathbf{B} = \begin{pmatrix} a_{11} \mathbf{B} & \cdots & a_{1n}\mathbf{B} \\ \vdots & \ddots & \vdots \\ a_{m1} \mathbf{B} & \cdots & a_{mn} \mathbf{B} \end{pmatrix}$

So, we need a transformation matrix: consider the following $4\times4$ matrix,
> k=4
> M=matrix(c(1,0,0,1,
+            0,1,1,0,
+            0,1,1,0,
+            1,0,0,1),k,k)
> M
[,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    1    0
[3,]    0    1    1    0
[4,]    1    0    0    1
Once we have it, we just consider the Kronecker product of this matrix with itself, which yields a $4^2\times4^2$ matrix,
> N=kronecker(M,M)
> N[,1:4]
[,1]  [,2] [,3] [,4]
[1,]     1    0    0    1
[2,]     0    1    1    0
[3,]     0    1    1    0
[4,]     1    0    0    1
[5,]     0    0    0    0
[6,]     0    0    0    0
[7,]     0    0    0    0
[8,]     0    0    0    0
[9,]     0    0    0    0
[10,]    0    0    0    0
[11,]    0    0    0    0
[12,]    0    0    0    0
[13,]    1    0    0    1
[14,]    0    1    1    0
[15,]    0    1    1    0
[16,]    1    0    0    1
And then, we continue,
> for(s in 1:3){N=kronecker(N,M)}
After only a couple of loops, we have a $4^5\times4^5$ matrix. And we can plot it simply to visualize the support,
> image(N,col=c("white","blue"))
As we zoom in, we can visualize this fractal property,

Thursday, September 27 2012

## Bounding sums of random variables, part 1

For the last course MAT8886 of this (long) winter session, on copulas (and extremes), we will discuss risk aggregation. The course will be mainly on the problem of bounding  the distribution (or some risk measure, say the Value-at-Risk) for two random variables with given marginal distribution. For instance, we have two Gaussian risks. What could be be worst-case scenario for the 99% quantile of the sum ? Note that I mention implications in terms of risk management, but of course, those questions are extremely important in terms of statistical inference, see e.g. Fan & Park (2006).

This problem, is sometimes related to some question asked by Kolmogorov almost one hundred years ago, as mentioned in Makarov (1981). One year after, Rüschendorf (1982) also suggested a proof of bounds calculation. Here, we focus in dimension 2. As usual, it is the simple case. But as mentioned recently, in Kreinovich & Ferson (2005), in dimension 3 (or higher), "computing the best-possible bounds for arbitrary n is an NP-hard (computationally intractable) problem". So let us focus on the case where we sum (only) two random variable (for those interested in higher dimension, Puccetti & Rüschendorf (2012) provided interesting results for a dual version of those optimal bounds).

Let $\Delta$ denote the set of univariate continuous distribution function, left-continuous, on $\mathbb{R}$. And $\Delta^+$ the set of distributions on $\mathbb{R}^+$. Thus, $F\in\Delta^+$ if $F\in\Delta$ and $F(0)=0$. Consider now two distributions $F,G\in\Delta^+$. In a very general setting, it is possible to consider operators on $\Delta^+\times \Delta^+$. Thus, let $T:[0,1]\times[0,1]\rightarrow[0,1]$ denote an operator, increasing in each component, thus that $T(1,1)=1$. And consider some function $L:\mathbb{R}^+\times\mathbb{R}^+\rightarrow\mathbb{R}^+$ assumed to be also increasing in each component (and continuous). For such functions $T$ and $L$, define the following (general) operator, $\tau_{T,L}(F,G)$ as

$\tau_{T,L}(F,G)(x)=\sup_{L(u,v)=x}\{T(F(u),G(v))\}$

One interesting case can be obtained when $T$is a copula, $C$. In that case,

$\tau_{C,L}(F,G):\Delta^+\times\Delta^+\rightarrow\Delta^+$

and further, it is possible to write

$\tau_{C,L}(F,G)(x)=\sup_{(u,v)\in L^{-1}(x)}\{C(F(u),G(v))\}$

It is also possible to consider other (general) operators, e.g. based on the sum

$\sigma_{C,L}(F,G)(x)=\int_{(u,v)\in L^{-1}(x)} dC(F(u),G(v))$

or on the minimum,

$\rho_{C,L}(F,G)(x)=\inf_{(u,v)\in L^{-1}(x)}\{C^\star(F(u),G(v))\}$

where $C^\star$ is the survival copula associated with $C$, i.e. $C^\star(u,v)=u+v-C(u,v)$. Note that those operators can be used to define distribution functions, i.e.

$\sigma_{C,L}(F,G):\Delta^+\times\Delta^+\rightarrow\Delta^+$

and similarly

$\rho_{C,L}(F,G):\Delta^+\times\Delta^+\rightarrow\Delta^+$

All that seems too theoretical ? An application can be the case of the sum, i.e. $L(x,y)=x+y$, in that case $\sigma_{C,+}(F,G)$ is the distribution of sum of two random variables with marginal distributions $F$ and $G$, and copula $C$. Thus, $\sigma_{C^\perp,+}(F,G)$ is simply the convolution of two distributions,

$\sigma_{C^\perp,+}(F,G)(x)=\int_{u+v=x} dC^\perp(F(u),G(v))$

The important result (that can be found in Chapter 7, in Schweizer and Sklar (1983)) is that given an operator $L$, then, for any copula $C$, one can find a lower bound for $\sigma_{C,L}(F,G)$

$\tau_{C^-,L}(F,G)\leq \tau_{C,L}(F,G)\leq\sigma_{C,L}(F,G)$

as well as an upper bound

$\sigma_{C,L}(F,G)\leq \rho_{C,L}(F,G)\leq\rho_{C^-,L}(F,G)$

Those inequalities come from the fact that for all copula $C$, $C\geq C^-$, where $C^-$ is a copula. Since this function is not copula in higher dimension, one can easily imagine that get those bounds in higher dimension will be much more complicated...

In the case of the sum of two random variables, with marginal distributions $F$ and $G$, bounds for the distribution of the sum $H(x)=\mathbb{P}(X+Y\leq x)$, where $X\sim F$ and $Y\sim G$, can be written

$H^-(x)=\tau_{C^- ,+}(F,G)(x)=\sup_{u+v=x}\{ \max\{F(u)+G(v)-1,0\} \}$

for the lower bound, and

$H^+(x)=\rho_{C^- ,+}(F,G)(x)=\inf_{u+v=x}\{ \min\{F(u)+G(v),1\} \}$

for the upper bound. And those bounds are sharp, in the sense that, for all $t\in(0,1)$, there is a copula $C_t$ such that

$\tau_{C_t,+}(F,G)(x)=\tau_{C^- ,+}(F,G)(x)=t$

and there is (another) copula $C_t$ such that

$\sigma_{C_t,+}(F,G)(x)=\tau_{C^- ,+}(F,G)(x)=t$

Thus, using those results, it is possible to bound cumulative distribution function. But actually, all that can be done also on quantiles (see Frank, Nelsen & Schweizer (1987)). For all $F\in\Delta^+$ let $F^{-1}$ denotes its generalized inverse, left continuous, and let $\nabla^+$ denote the set of those quantile functions. Define then the dual versions of our operators,

$\tau^{-1}_{T,L}(F^{-1},G^{-1})(x)=\inf_{(u,v)\in T^{-1}(x)}\{L(F^{-1}(u),G^{-1}(v))\}$

and

$\rho^{-1}_{T,L}(F^{-1},G^{-1})(x)=\sup_{(u,v)\in T^\star^{-1}(x)}\{L(F^{-1}(u),G^{-1}(v))\}$

Those definitions are really dual versions of the previous ones, in the sense that $\tau^{-1}_{T,L}(F^{-1},G^{-1})=[\tau_{T,L}(F,G)]^{-1}$ and $\rho^{-1}_{T,L}(F^{-1},G^{-1})=[\rho_{T,L}(F,G)]^{-1}$.

Note that if we focus on sums of bivariate distributions, the lower bound for the quantile of the sum is

$\tau^{-1}_{C^{-},+}(F^{-1},G^{-1})(x)=\inf_{\max\{u+v-1,0\}=x}\{F^{-1}(u)+G^{-1}(v)\}$

while the upper bound is

$\rho^{-1}_{C^{-},+}(F^{-1},G^{-1})(x)=\sup_{\min\{u+v,1\}=x}\{F^{-1}(u)+G^{-1}(v)\}$

A great thing is that it should not be too difficult to compute numerically those quantities. Perhaps a little bit more for cumulative distribution functions, since they are not defined on a bounded support. But still, if the goal is to plot those bounds on $[0,10]$, for instance. The code is the following, for the sum of two lognormal distributions $LN(0,1)$

> F=function(x) plnorm(x,0,1)
> G=function(x) plnorm(x,0,1)
> n=100
> X=seq(0,10,by=.05)
> Hinf=Hsup=rep(NA,length(X))
> for(i in 1:length(X)){
+ x=X[i]
+ U=seq(0,x,by=1/n); V=x-U
+ Hinf[i]=max(pmax(F(U)+G(V)-1,0))
+ Hsup[i]=min(pmin(F(U)+G(V),1))}

If we plot those bounds, we obtain

> plot(X,Hinf,ylim=c(0,1),type="s",col="red")
> lines(X,Hsup,type="s",col="red")

But somehow, it is even more simple to work with quantiles since they are defined on a finite support. Quantiles are here

> Finv=function(u) qlnorm(u,0,1)
> Ginv=function(u) qlnorm(u,0,1)

The idea will be to consider a discretized version of the unit interval as discussed in Williamson (1989), in a much more general setting. Again the idea is to compute, for instance

$\sup_{u\in[0,x]}\{F^{-1}(u)+G^{-1}(x-u)\}$

The idea is to consider $x=i/n$ and $u=j/n$, and the bound for the quantile function at point $i/n$ is then

$\sup_{j\in\{0,1,\cdots,i\}}\left\{F^{-1}\left(\frac{j}{n}\right)+G^{-1}\left(\frac{i-j}{n}\right)\right\}$

The code to compute those bounds, for a given $n$ is here

> n=1000
> Qinf=Qsup=rep(NA,n-1)
> for(i in 1:(n-1)){
+ J=0:i
+ Qinf[i]=max(Finv(J/n)+Ginv((i-J)/n))
+ J=(i-1):(n-1)
+ Qsup[i]=min(Finv((J+1)/n)+Ginv((i-1-J+n)/n))
+ }

Here we have (several $n$s were considered, so that we can visualize the convergence of that numerical algorithm),

Here, we have a simple code to visualize bounds for quantiles for the sum of two risks. But it is possible to go further...

Thursday, September 20 2012

## (nonparametric) Copula density estimation

Today, we will go further on the inference of copula functions. Some codes (and references) can be found on a previous post, on nonparametric estimators of copula densities (among other related things).  Consider (as before) the loss-ALAE dataset (since we've been working a lot on that dataset)

> library(MASS)
> library(evd)
> X=lossalae
> U=cbind(rank(X[,1])/(nrow(X)+1),rank(X[,2])/(nrow(X)+1))

The standard tool to plot nonparametric estimators of densities is to use multivariate kernels. We can look at the density using

> mat1=kde2d(U[,1],U[,2],n=35)
> persp(mat1$x,mat1$y,mat1$z,col="green", + shade=TRUE,theta=s*5, + xlab="",ylab="",zlab="",zlim=c(0,7)) or level curves (isodensity curves) with more detailed estimators (on grids with shorter steps) > mat1=kde2d(U[,1],U[,2],n=101) > image(mat1$x,mat1$y,mat1$z,col=
+ rev(heat.colors(100)),xlab="",ylab="")
> contour(mat1$x,mat1$y,mat1$z,add= + TRUE,levels = pretty(c(0,4), 11)) Kernels are nice, but we clearly observe some border bias, extremely strong in corners (the estimator is 1/4th of what it should be, see another post for more details). Instead of working on sample $(U_i,V_i)$ on the unit square, consider some transformed sample $(Q(U_i),Q(V_i))$, where $Q:(0,1)\rightarrow\mathbb{R}$ is a given function. E.g. a quantile function of an unbounded distribution, for instance the quantile function of the $\mathcal{N}(0,1)$ distribution. Then, we can estimate the density of the transformed sample, and using the inversion technique, derive an estimator of the density of the initial sample. Since the inverse of a (general) function is not that simple to compute, the code might be a bit slow. But it does work, > gaussian.kernel.copula.surface <- function (u,v,n) { + s=seq(1/(n+1), length=n, by=1/(n+1)) + mat=matrix(NA,nrow = n, ncol = n) + sur=kde2d(qnorm(u),qnorm(v),n=1000, + lims = c(-4, 4, -4, 4)) + su<-sur$z
+ for (i in 1:n) {
+     for (j in 1:n) {
+ 	Xi<-round((qnorm(s[i])+4)*1000/8)+1;
+ 	Yj<-round((qnorm(s[j])+4)*1000/8)+1
+ 	mat[i,j]<-su[Xi,Yj]/(dnorm(qnorm(s[i]))*
+ 	dnorm(qnorm(s[j])))
+     }
+ }
+ return(list(x=s,y=s,z=data.matrix(mat)))
+ }

Here, we get

Note that it is possible to consider another transformation, e.g. the quantile function of a Student-t distribution.

> student.kernel.copula.surface =
+  function (u,v,n,d=4) {
+  s <- seq(1/(n+1), length=n, by=1/(n+1))
+  mat <- matrix(NA,nrow = n, ncol = n)
+ sur<-kde2d(qt(u,df=d),qt(v,df=d),n=5000,
+ lims = c(-8, 8, -8, 8))
+ su<-sur$z + for (i in 1:n) { + for (j in 1:n) { + Xi<-round((qt(s[i],df=d)+8)*5000/16)+1; + Yj<-round((qt(s[j],df=d)+8)*5000/16)+1 + mat[i,j]<-su[Xi,Yj]/(dt(qt(s[i],df=d),df=d)* + dt(qt(s[j],df=d),df=d)) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + } Another strategy is to consider kernel that have precisely the unit interval as support. The idea is here to consider the product of Beta kernels, where parameters depend on the location > beta.kernel.copula.surface= + function (u,v,bx=.025,by=.025,n) { + s <- seq(1/(n+1), length=n, by=1/(n+1)) + mat <- matrix(0,nrow = n, ncol = n) + for (i in 1:n) { + a <- s[i] + for (j in 1:n) { + b <- s[j] + mat[i,j] <- sum(dbeta(a,u/bx,(1-u)/bx) * + dbeta(b,v/by,(1-v)/by)) / length(u) + } + } + return(list(x=s,y=s,z=data.matrix(mat))) + } On those two graphs, we can clearly observe strong tail dependence in the upper (right) corner, that cannot be intuited using a standard kernel estimator... Tuesday, September 18 2012 ## Copulas and tail dependence, part 3 We have seen extreme value copulas in the section where we did consider general families of copulas. In the bivariate case, an extreme value can be written where $A(\cdot)$ is Pickands dependence function, which is a convex function satisfying Observe that in this case, where is Kendall'tau, and can be written For instance, if then, we obtain Gumbel copula. This is what we've seen in the section where we introduced this family. Now, let us talk about (nonparametric) inference, and more precisely the estimation of the dependence function. The starting point of the most standard estimator is to observe that if has copula $C$, then has distribution function And conversely, Pickands dependence function can be written Thus, a natural estimator for Pickands function is where $\widehat{H}_n$ is the empirical cumulative distribution function of This is the estimator proposed in Capéràa, Fougères & Genest (1997). Here, we can compute everything here using > library(evd) > X=lossalae > U=cbind(rank(X[,1])/(nrow(X)+1),rank(X[,2])/(nrow(X)+1)) > Z=log(U[,1])/log(U[,1]*U[,2]) > h=function(t) mean(Z<=t) > H=Vectorize(h) > a=function(t){ + f=function(t) (H(t)-t)/(t*(1-t)) + return(exp(integrate(f,lower=0,upper=t, + subdivisions=10000)$value))
+ }
> A=Vectorize(a)
> u=seq(.01,.99,by=.01)
> plot(c(0,u,1),c(1,A(u),1),type="l",col="red",
+ ylim=c(.5,1))
Even integrate to get an estimator of Pickands' dependence function. Note that an interesting point is that the upper tail dependence index can be visualized on the graph, above,

> A(.5)/2
[1] 0.4055346

Monday, September 17 2012

## Copulas and tail dependence, part 1

As mentioned in the course last week Venter (2003) suggested nice functions to illustrate tail dependence (see also some slides used in Berlin a few years ago).
• Joe (1990)'s lambda
Joe (1990) suggested a (strong) tail dependence index. For lower tails, for instance, consider
i.e

• Upper and lower strong tail (empirical) dependence functions

The idea is to plot the function above, in order to visualize limiting behavior. Define

for the lower tail, and
for the upper tail, where is the survival copula associated with , in the sense that
while

Now, one can easily derive empirical conterparts of those function, i.e.

and
Thus, for upper tail, on the right, we have the following graph

and for the lower tail, on the left, we have

For the code, consider some real data, like the loss-ALAE dataset.

> library(evd)
> X=lossalae

The idea is to plot, on the left, the lower tail concentration function, and on the right, the upper tail function.

> U=rank(X[,1])/(nrow(X)+1)
> V=rank(X[,2])/(nrow(X)+1)
> Lemp=function(z) sum((U<=z)&(V<=z))/sum(U<=z)
> Remp=function(z) sum((U>=1-z)&(V>=1-z))/sum(U>=1-z)
> u=seq(.001,.5,by=.001)
> L=Vectorize(Lemp)(u)
> R=Vectorize(Remp)(rev(u))
> plot(c(u,u+.5-u[1]),c(L,R),type="l",ylim=0:1,
+ xlab="LOWER TAIL          UPPER TAIL")
> abline(v=.5,col="grey")

Now, we can compare this graph, with what should be obtained for some parametric copulas that have the same Kendall's tau (e.g.). For instance, if we consider a Gaussian copula,

> tau=cor(lossalae,method="kendall")[1,2]
> library(copula)
> paramgauss=sin(tau*pi/2)
> copgauss=normalCopula(paramgauss)
> Lgaussian=function(z) pCopula(c(z,z),copgauss)/z
> Rgaussian=function(z) (1-2*z+pCopula(c(z,z),copgauss))/(1-z)
> u=seq(.001,.5,by=.001)
> Lgs=Vectorize(Lgaussian)(u)
> Rgs=Vectorize(Rgaussian)(1-rev(u))> lines(c(u,u+.5-u[1]),c(Lgs,Rgs),col="red")

or Gumbel's copula,

> paramgumbel=1/(1-tau)
> copgumbel=gumbelCopula(paramgumbel, dim = 2)
> Lgumbel=function(z) pCopula(c(z,z),copgumbel)/z
> Rgumbel=function(z) (1-2*z+pCopula(c(z,z),copgumbel))/(1-z)
> u=seq(.001,.5,by=.001)
> Lgl=Vectorize(Lgumbel)(u)
> Rgl=Vectorize(Rgumbel)(1-rev(u))
> lines(c(u,u+.5-u[1]),c(Lgl,Rgl),col="blue")

That's nice (isn't it?), but since we do not have any confidence interval, it is still hard to conclude (even if it looks like Gumbel copula has a much better fit than the Gaussian one). A strategy can be to generate samples from those copulas, and to visualize what we had. With a Gaussian copula, the graph looks like
> u=seq(.0025,.5,by=.0025); nu=length(u)
> nsimul=500
> MGS=matrix(NA,nsimul,2*nu)
> for(s in 1:nsimul){
+ Xs=rCopula(nrow(X),copgauss)
+ Us=rank(Xs[,1])/(nrow(Xs)+1)
+ Vs=rank(Xs[,2])/(nrow(Xs)+1)
+ Lemp=function(z) sum((Us<=z)&(Vs<=z))/sum(Us<=z)
+ Remp=function(z) sum((Us>=1-z)&(Vs>=1-z))/sum(Us>=1-z)
+ MGS[s,1:nu]=Vectorize(Lemp)(u)
+ MGS[s,(nu+1):(2*nu)]=Vectorize(Remp)(rev(u))
+ lines(c(u,u+.5-u[1]),MGS[s,],col="red")
+ }

(including - pointwise - 90% confidence bands)

> Q95=function(x) quantile(x,.95)
> V95=apply(MGS,2,Q95)
> lines(c(u,u+.5-u[1]),V95,col="red",lwd=2)
> Q05=function(x) quantile(x,.05)
> V05=apply(MGS,2,Q05)
> lines(c(u,u+.5-u[1]),V05,col="red",lwd=2)

while it is

with Gumbel copula. Isn't it a nice (graphical) tool ?

But as mentioned in the course, the statistical convergence can be slow. Extremely slow. So assessing if the underlying copula has tail dependence, or not, it now that simple. Especially if the copula exhibits tail independence. Like the Gaussian copula. Consider a sample of size 1,000. This is what we obtain if we generate random scenarios,

or we look at the left tail (with a log-scale)

Now, consider a 10,000 sample,

or with a log-scale

We can even consider a 100,000 sample,

or with a log-scale

On those graphs, it is rather difficult to conclude if the limit is 0, or some strictly positive value (again, it is a classical statistical problem when the value of interest is at the border of the support of the parameter). So, a natural idea is to consider a weaker tail dependence index. Unless you have something like 100,000 observations...

Wednesday, September 12 2012

## Kendall's function for copulas

As mentioned in the course on copulas, a nice tool to describe dependence it Kendall's cumulative function. Given a random pair with distribution  , define random variable . Then Kendall's cumulative function is

Genest and Rivest (1993) introduced that function to choose among Archimedean copulas (we'll get back to this point below).

From a computational point of view, computing such a function can be done as follows,

• for all , compute as the proportion of observation in the lower quadrant, with upper corner , i.e.

• then compute the cumulative distribution function of 's.

To visualize the construction of that cumulative distribution function, consider the following animation

Thus, here the code to compute simply that cumulative distribution function is

n=nrow(X)
i=rep(1:n,each=n)
j=rep(1:n,n)
S=((X[i,1]>X[j,1])&(X[i,2]>X[j,2]))
Z=tapply(S,i,sum)/(n-1)

The graph can be obtain either using

plot(ecdf(Z))

or

plot(sort(Z),(1:n)/n,type="s",col="red")

The interesting point is that for an Archimedean copula with generator , then Kendall's function is simply

If we're too lazy to do the maths, at least, it is possible to compute those functions numerically. For instance, for Clayton copula,

h=.001
phi=function(t){(t^(-alpha)-1)}
dphi=function(t){(phi(t+h)-phi(t-h))/2/h}
k=function(t){t-phi(t)/dphi(t)}
Kc=Vectorize(k)

Similarly, let us consider Gumbel copula,

phi=function(t){(-log(t))^(theta)}
dphi=function(t){(phi(t+h)-phi(t-h))/2/h}
k=function(t){t-phi(t)/dphi(t)}
Kg=Vectorize(k)

If we plot the empirical Kendall's function (obtained from the sample), with different theoretical ones, derived from Clayton copulas (on the left, in blue) or Gumbel copula (on the right, in purple), we have the following,

Note that the different curves were obtained when Clayton copula has Kendall's tau equal to 0, .1, .2, .3, ..., .9, 1, and similarly for Gumbel copula (so that Figures can be compared). The following table gives a correspondence, from Kendall's tau to the underlying parameter of a copula (for different families)

as well as Spearman's rho,

To conclude, observe that there are two important particular cases that can be identified here: the case of perfect dependent, on the first diagonal when , and the case of independence, the upper green curve, . It should also be mentioned that it is also common to plot not function , but function , defined as ,

Friday, February 24 2012

## MAT8886 copula dataset

To illustrate inference issues on the copula section, we will use the so-called loss-alae dataset. This dataset contains

• LOSS : the loss amount, up to the limit.
• ALAE : the associated expense,
• LIMIT : the limit (-99 means no limit), and
• CENSOR (1=censored/limit reached, 0=otherwise).

The dataset is the following

> lossalae=read.table(
+ "http://freakonometrics.free.fr/lossalae.txt",
> X=lossalae$LOSS > Y=lossalae$ALAE
> plot(X,Y)

It is possible to visualize that dataset, on a log-scale

> xhist <- hist(log(X), plot=FALSE)
> yhist <- hist(log(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="",log="xy",col="red")
> 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")

or the copula-type scatterplot (i.e. ranks of the observations),

> n=length(X)+1
> U=rank(X)/n
> V=rank(Y)/n
> xhist <- hist((U), plot=FALSE)
> yhist <- hist((V), 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(U,V, xlab="", ylab="",col="red")
> 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")

Thursday, December 1 2011

## MAT8886 Extreme values and copulas: assignment (part 1)

This winter, I will give a course on extremes and copulas, at UQAM (every Thursday and Friday, from 13:00 till 14:30). More details will come about the syllabus and the evaluation, but a large part of the evaluation will be based on a report on one of the following topics. On extreme values
1. power law and application, Katz & Katz (1999) applications in athletics, Savaglio & Carbone (2000) applications in sports, Gabaix, Gopikrishnan Plerou & Stanley (2003) applications in finance, see also Mandelbrot & Taleb. (2005), or for more theoretical material, Clauset, Shalizi & Newman (2007) for power law in empirical data
2. sport records and extreme value, Einmahl & Smeets (2009) ultimate 100m, Einmahl & Magnus (2008), records in athletics
3. smoothing techniques and extremes, Chavez-Desmoulins & Davison (2005) generalized additive models, Hall & Tajvidi (2000) temporal trend, Padoan & Wand (2008) , Yee & Stephensen (2007) vgam
4. extremes and exponential regression, Beirlant, Dierckx, Goegebeur & Matthys (1999) exponential regression, Beirlant, Dierckx, Guillou & Starica (2002) log-spacing, Matthys & Beirlant (2003), high quantile regression, Beirlant, Dierckx & Guillou (2005), generalized quantile plot
5. extremes with random sample size, Vorn (1989) stability, Barakat (1990), Silvestrov & Teugels (1998) limit theorem
6. sea and wind: multivariate extremes, de Haan & de Ronde (1998), and the Neptune project,
7. subexponential distributions Teugels (1975), Pitman (1980), Klüppelberg (2006) encyclopedia of actuarial sciences,
8. ECOMOR reinsurance treaty, Kremer (1982), Beirlant (2006)
On copulas
1. NLCA nonlinear canonical analysis, Lancaster (1958) bivariate distribution Dauxois & Pousse (1975), extension of canonical analysis,
2. maximum correlation, Gebelin (1941) Korrelation als variation, Lancaster (1957), normal distribution and contingency tables Rényi (1959), measures of dependence, Sethuraman (1990) properties of estimator, Buja (1990) functional canonical variates, Yu (2008) properties,
3. copulas and Markov chains, Darsow, Nguyen & Olsen (1992), Schmitz (2003) chapter 5
4. Levy copulas, Kallsen & Tankov (2006), introduction to Levy copulas, Barndorff-Nielsen & Lidner (2004) aspects of Levy copulas,
5. vine copulas, Bedford (2002) graphical model for dependence,
7. dependent mortality, Frees, Carrière & Valdez (1997), Ji, Hardy & Li (2011) Markovian approach
8. Levy-frailty copulas,  Marshall & Olkin (1967) multivariate exponential distributions, Mai & Scherer (2009), Gnedin & Pitman (2008)

and references therein...

Wednesday, November 16 2011

## PhD defense on copulas

This Wednesday I will be at Université Paris 1 Sorbonne as a member of the jury of the PhD thesis of Pierre-André Maugis, on conditional correlation and vine copula.

Vine copulas were born in 2002 with the paper of Tim Bedford and Roger M. Cooke Vines--a new graphical model for dependent random variables. The idea is to use the following decomposition for a multivariate density

(from Bayes formula, with synthetic notations). Then using the relationship between a bivariate density and its copula (density)
thus
Using again Bayes formula,
and we can write
Since and , the previous expression becomes
or to stress on the most important part (as I see it)
It is common then to assume that this conditional copula does not depend on the conditioning parameter. The more detailed expression of that joint trivariate density is
The (parametric) inference algorithm is defined in Cooke, Joe and Aas (2010) as follows

The important assumption in vine copula models is that conditional copulas are constant. And this assumption might be relevant in some cases. For instance, in the Gaussian case (the observations have a Gaussian joint distribution - or at least copula - and we fit a vine model with Gaussian bivariate copulas).

The code to fit a vine copula is the following,
> library(CDVine)> library(mnormt)> SIGMA=matrix(c(1,.6,.7,.6,1,.8,.7,.8,1),3,3)> X=rmnorm(n=100000,varcov=SIGMA)> CDVineSeqEst(dat=X, family = c(1,1,1),+ type = 1, method = "mle")$par[1] 0.6001505 0.7023699 0.6698215$par2[1] 0 0 0
Note that it is consistent with the following algorithm where conditional copulas are fitted. In the following, for all values of the given component, we wit a Gaussian copula for the conditional remaining pair,
> U=pnorm(X)> U1U2=U[,1:2]> U1U3=U[,c(1,3)]> GaussCop = normalCopula(param=.5, dim = 2)> U1U2=U[,1:2]> U1U3=U[,c(1,3)]> fit12.mpl = fitCopula(GaussCop, U1U2, method="mpl")@estimate> fit13.mpl = fitCopula(GaussCop, U1U3, method="mpl")@estimate> fit12.mpl[1] 0.5984932> fit13.mpl[1] 0.7005185> fit23a=fit23b=rep(NA,99)> for(i in 4:96){+ x=i/100+ C12=pcopula(normalCopula(param=fit12.mpl, dim = 2),U1U2)+ C13=pcopula(normalCopula(param=fit13.mpl, dim = 2),U1U3)+ U12=rank(C12)/(nrow(U)+1)+ U13=rank(C13)/(nrow(U)+1)+ U23=cbind(U12[abs(U[,1]-x)<.02],U13[abs(U[,1]-x)<.02])+ V23=cbind(rank(U23[,1])/(nrow(U23)+1),+ rank(U23[,2])/(nrow(U23)+1))+ fit23.mpl = fitCopula(GaussCop, V23, method="mpl")@estimate+ fit23a[i]=fit23.mpl+ }> plot(X,fit23a,col="red")
It looks like assuming the conditional copula as constant was a valid assumption here

But note that if the true distribution is not Gaussian, then assuming the conditional copula as constant is not valid anymore (here a trivariate Clayton copula was generated)

Wednesday, April 6 2011

## Talk at Laval University at the Actuarial Seminar

I was last Friday at Laval University for a conference by David Cummins and Mary Weiss (here). I will be back tomorrow, this time to give a talk, on "distorting probabilities in actuarial science" (the talk will be extremely close to the one I gave at McGill in November). "In this talk, we will first get back on properties of distortion operators for pricing financial and insurance risks. Based on the dual version of the expected utility framework, we will see how distorted risk measures have been introduced, from VaR and TVaR, to Esscher premium and Wang's measures. Then we will discuss extensions in higher dimension. We will discuss tail properties of distorted copulas (in the particular case of Archimedean copulas). A natural application will be aging problems (in survival analysis or in credit risk)." Slides can be downloaded from here.

This talk can be seen as a first part, the second one behing the talk I will give in 15 days, again at Laval University, but this time for the Seminar of Statistics. The talk will be on "Beta kernel and transformed kernel : applications to quantile estimation, and copula density estimation".

Tuesday, February 15 2011

## In statistics, having too much information might not be a good thing

A common idea in statistics is that if we don't know something, and we use an estimator of that something (instead of the true value) then there will be some additional uncertainty.

For instance, consider a random sample, i.i.d., from a Gaussian distribution. Then, a confidence interval for the mean is

where is the quantile of probability level of the standard normal distribution . But usually, standard deviation (the something is was talking about earlier) is usually unknown. So we substitute an estimation of the standard deviation, e.g.
and the cost we have to pay is that the new confidence interval is
where now is the quantile of the Student distribution, of probability level , with degrees of freedom.
We call it a cost since the new confidence interval is now larger (the Student distribution has higher upper-quantiles than the Gaussian distribution).
So usually, if we substitute an estimation to the true value, there is a price to pay.
A few years ago, with Jean David Fermanian and Olivier Scaillet, we were writing a survey on copula density estimation (using kernels,  here). At the end, we wanted to add a small paragraph on the fact that we assumed that we wanted to fit a copula on a sample i.i.d. with distribution , a copula, but in practice, we start from a sample with joint distribution (assumed to have continuous margins, and - unique - copula ). But since margins are usually unknown, there should be a price for not observing them.
To be more formal, in a perfect wold, we would consider
but in the real world, we have to consider
where it is standard to consider ranks, i.e. are empirical cumulative distribution functions.
My point is that when I ran simulations for the survey (the idea was more to give illustrations of several techniques of estimation, rather than proofs of technical theorems) we observed that the price to pay... was negative ! I.e. the variance of the estimator of the density (wherever on the unit square) was smaller on the pseudo sample than on perfect sample .
By that time, we could not understand why we got that counter-intuitive result: even if we do know the true distribution, it is better not to use it, and to use instead a nonparametric estimator. Our interpretation was based on the discrepancy concept and was related to the latin hypercube construction:

With ranks, the data are more regular, and marginal distributions are exactly uniform on the unit interval. So there is less variance.
This was our heuristic interpretation.
A couple of weeks ago, Christian Genest and Johan Segers proved that intuition in an article published in JMVA,

Well, we observed something for finite , but Christian and Johan obtained an analytical result. Hence, if we denote
the empirical copula in the perfect world (with known margins) and
the one constructed from the pseudo sample, they obtained that, everywhere
with nice graphs of ,
So I was very happy last week when Christian show me their results, to learn that our intuition was correct. Nevertheless, it is still a very counter-intuitive result.... If anyone has seen similar things, I'd be glad to hear about it !

Friday, November 26 2010

## Talk at McGill

I will give a talk at McGill university this afternoon, on "distorting probabilities in actuarial science". Note that Louis Paul Rivest will give a talk just after, but at UQAM, at the statistical seminar (here)

Tuesday, November 16 2010

## Course on copulas and correlated risks (in French, still)

The course on copulas, in Luminy, starts at 8.30 on Wednesday (here). The slides can be found here.

Saturday, October 16 2010

## Mandelbrot, fractals and counterexamples in applied probability

Benoît Mandelbrot died yesterday. Like most of the blogs dealing with applied mathematics, it looks like I have to mention this event. Unfortunately, I don't know much about fractals...
The first time I heard about Mandelbrot and chaos was when I have been working on fractional time series (see eg here). Murad Taqqu gave a very interesting short course in Paris, and I have been using it in two papers (actually one more should appear soon in Climate Change).
The second time was in Québec (city), five years ago, when Roger Nelsen gave a talk on copulas with fractal suppport (here). By that time, we were finishing our paper with Alessandro (in mathematical finance, and limiting theorems). I remember adding a Remark in the paper (that can be found here), since using that kind of copulas was a nice way to show that, without sufficient regularity conditions, the limit we were looking for had no sense.

A few months after, with Johan, in a paper on pitfalls on lower tail dependence for Archimedean copulas (here), we used again this fractal construction (here applied to Archimedean copulas) to find a nice counterexample to a (false) theorem on regular variation. Again, the goal was to understand how the dependence structure of  given  and  changed when  goes to 0. In the case of Archimedean copulas, the copula of the conditional pair is still Archimedean (with another generator, except for Clayton copula). The graph below show how  changes, as  decreases... Actually, I draw  since Archimedean generators are not unique.

I have also decided to plot .

Here, we see that there is no way of talk about a possible limit for the conditional copula because of a fractal behavior in 0 of the generator (even if my fractal are not as nice as the one you can find on the internet....). So thanks Benoît for giving us a nice toy to build interesting counterexamples !

Thursday, October 14 2010

## Copules et risques corrélés

J'avais promis dans un commentaire que je mettrais bientôt en ligne un survey sur les copules.... après plusieurs mois de retard, le document est en ligne ici, et il est sobrement intitulé copules et risques multiples. Toutes les remarques et critiques sont les bienvenues ! Il s'agit d'un chapitre pour un livre dont je mettrais la référence en ligne ultérieurement. A priori, ce document devrait servir de base pour le cours qui sera donné dans un mois au CIRM (mentionné ici, dans le cadre des Journées d'Études Statistiques).

- page 1 of 2