---
title: Key characteristics of the Gumbel copula
author: "Anne-lise Caillat, Christophe Dutang, Marie Veronique Larrieu, Triet NGuyen"
date: '`r Sys.Date()`'
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    fig_caption: yes
    toc: true
    number_sections: yes
bibliography: gumbel.bib
link-citations: true
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Key characteristics of the Gumbel copula}
  %!\VignetteEncoding{UTF-8}
  \usepackage[utf8]{inputenc}
  \usepackage{amsmath}
---

<!-- rmarkdown::render("GumbelDoc.Rmd") -->
 
# Introduction

Copulas constitute a statistical tool used to model the dependence between random variables. The copula function indeed links the joint density to the marginal densities and thus contains all the information about the dependence structure of the model.

In this document, we will focus on a particular copula, namely the Gumbel copula. In the literature, it is sometimes referred to as the Gumbel-Hougaard copula. For clarity, we will consider the bivariate case here and will conclude with a generalization to the multivariate case. Most of the information is drawn from @nelsen.

# Characteristics

## Definition

The Gumbel copula is defined as
$$
C_\alpha\left(u,v\right)=\exp\left[-\left(\left(-\ln u\right)^\alpha+\left(-\ln v\right)^\alpha\right)^{\frac{1}{\alpha}}\right],
$$
where $\alpha\geq1$ is the parameter.

## Archimedean copula
The generator is $\phi_\alpha(t)=\left(-\ln t\right)^\alpha$ with $\alpha\geq1$ and $t\in [0,1]$.

\begin{itemize}
- the inverse is $\phi_\alpha^{-1}(t)=e^{-t^{1/\alpha}}$,
- $\phi_\alpha(0)=+\infty$, 
- $\phi_\alpha(1)=0$,
- $\phi_\alpha^\prime(t)=-\frac{\alpha}{t}(-\ln t)^{\alpha-1}$ et $\phi_\alpha^\prime(t)<0$ pour $\alpha\geq1$,
- $\phi_\alpha^{\prime\prime}(t)=\frac{\alpha}{t^2}(-\ln t)^{\alpha-2}\left[\alpha-1-\ln t\right]$ et $\phi_\alpha^{\prime\prime}(t)\geq0$ pour $\alpha\geq1$.
\end{itemize}
Thus, $\phi$ is a continuous, strictly decreasing function from $[0,1]$ to $[0,+\infty]$, convex, and it is a strict generator.

We retrieve the Gumbel family through the relation $\phi\alpha^{-1}\left(\phi\alpha(u)+\phi_\alpha(v)\right)$. Archimedean copulas form families of copulas that possess several interesting properties.


## Density
We know that the density of an Archimedean copula, with a twice-differentiable generator is given by
$$
c_\alpha(u,v)=-\frac{ \phi_{\alpha}^\prime\left(C_\alpha(u,v)\right) \phi_{\alpha}^\prime(u) \phi_{\alpha}^\prime(v)}{ \phi_{\alpha}^\prime\left(C_\alpha(u,v)\right)^3}
$$
This leads to the following density
\begin{center}
\fbox{$\forall u, v \in ]0,1[^2, c_\alpha(u,v)=C_\alpha(u,v) \left[ \phi_{\alpha}(u)+\phi_{\alpha}(u)  \right]^{\frac{1}{\alpha}-2}  \left[ \alpha-1+\left( \phi_{\alpha}(u)+\phi_{\alpha}(u)\right)^{\frac{1}{\alpha}} \right] \frac{\phi_{\alpha-1}(u)\phi_{\alpha-1}(v)}{uv}.
$}
\end{center}
The supra expression of the density is not valid on the edges of the square $[0,1] \times [0,1]$. 
On the edge of the domain, we have the following expressions
$$
\forall u>0,v<1, c_\alpha(u,0)=c_\alpha(0,u) = c_\alpha(1,v) =c_\alpha(v,1)=0,
$$
but the density is infinite at $(0,0)$ and $(1,1)$.

 

```{r, fig.cap="Density of the Gumbel copula", fig.label="density", fig.width=8, fig.height=4, echo=FALSE, message=FALSE, warning=FALSE}
require(gumbel)         
x <- y <- seq(.1, .9, length = 30)
nrz <- length(x)
ncz <- length(y)
# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("white", "green", "yellow", "orange", "red") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)

par(mfrow=1:2, mar=c(1,2,1,1))

z<-outer(x, y, dgumbel, alpha=1.5)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = color[facetcol],
       ltheta = 100, shade = 0.25, ticktype = "detailed",
       xlab = "u", ylab = "v", zlab = "density", main=expression(alpha==1.5))
         
z <- outer(x, y, dgumbel, alpha=4)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
                 
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = color[facetcol],
       ltheta = 100, shade = 0.25, ticktype = "detailed",
       xlab = "u", ylab = "v", zlab = "density", main=expression(alpha==4))
```





## Dependence measure


### Kendall's tau
Kendall's tau is
$$
\tau = \frac{\alpha-1}{\alpha}.
$$



### Spearman's rho

Spearman's rho has no explicit form
$$\rho_S = 12\int\!\!\!\!\int_{[0,1]^2}C(u,v)dudv-3$$



## Extreme copula

The concept of tail dependence provides information about the "amount" of dependence at the distribution tails. It is a relevant tool for studying the simultaneous occurrence of extreme values. This is a local measure, unlike Kendall's tau and Spearman's rho, which measure dependence across the entire distribution. The coefficients of tail dependence, on the left and right, for a pair $(X,Y)$, are defined by
$$
\lambda_L = \underset{t\rightarrow0^+}{\lim} P(Y > F_Y^{-1}(t) / X > F_X^{-1}(t)$ = 0
\text{and} 
\lambda_U = \underset{t\rightarrow1^-}{\lim} P(Y > F_Y^{-1}(t) / X > F_X^{-1}(t)$ =  2-2^{\frac{1}{\alpha}}.
$$

The Gumbel copula is max-stable, i.e. 
$$
C_\alpha\left(u^{\frac{1}{n}},v^{\frac{1}{n}}\right)^n = C_\alpha(u,v).
$$


## Multivariate copula


Using @nelsen, the Gumbel copula is
\begin{center}
\fbox{$C_\alpha\left(u_1,\dots,u_n\right)=\exp\left(-\left[\sum\limits_{i=1}^n\left(-\ln u_i\right)^\alpha\right]^{1/\alpha}\right)$}
\end{center}



## Simulation


### Algorithm using the Kendall function ``$K_C$''

Using @nelsen2, the algorihm is

-  simulate $(y,t) \stackrel{iid}{\sim} \mathcal U(0,1)$
-  $x:=K_C^{-1}(t)$ o\`u $K_C(t)=t-\frac{t\ln t}{\alpha}$
-  $u := \phi^{-1}\left(\phi(x)y \right)$ et $v := \phi^{-1}\left(\phi(x)(1-y) \right)$


## Frailty approach by @marshall

Marshall's approach involves using a random variable $\Theta$ for a random vector $(X_1, \dots, X_d)$ such that the components $X_i$ are conditionally independent given $\Theta$. The joint distribution of the vector $(X_1, \dots, X_d)$ is given by
$$
F_{X_1, \dots, X_d}(x_1, \dots, x_d)=L_\Theta\left( \sum_{i=1}^d L^{-1}_\Theta\left(F_{X_i}(x_i)\right) \right)
$$
where $L(\Theta)$ denotes the Laplace transform of the random variable $\Theta$. Moreover, the following algorithm allows for the simulation of random variables from an Archimedean copula $L_\Theta^{-1}$:

- Simulate $(x_1, \dots, x_d) \stackrel{iid}{\sim} \text{Exp}(1)$,
- Simulate $\theta$ following the distribution with Laplace transform $\phi^{-1}$,
- For each $i$, set $u_i := \phi^{-1}\left(\frac{x_i}{\theta}\right)$.

In the case of the Gumbel copula, $\phi^{-1}(t) = e^{-t^{1/\alpha}}$, which is the Laplace transform of a stable law with parameters $(1/\alpha, 0, 1, 0)$. For more details on stable laws, refer to @nolan. Finally, to simulate random variables with a stable distribution, we used @chambers's algorithm.


# Fitting methods

In this section, we briefly present the four usual fitting methods: 

- Method of Moments provided in `gumbel.MBE()`
- Exact max_imum likelihood provided in `gumbel.EML()`
- Inference for margins provided in `gumbel.IFM()`
- Canonical max_imum likelihood provided in `gumbel.CML()`


## Method of Moments - `gumbel.MBE`

This method involves estimating the parameters $\theta$ of the marginal laws and the parameter $\alpha$ of the copula using the method of moments, i.e.,

-  solve
$$
\left\{ \begin{array}{c}
\overline X_n = f(\theta_1,\dots,\theta_d)\\
S_n^2 = g(\theta_1,\dots,\theta_d)\\
\mu_{3,n} = h(\theta_1,\dots,\theta_d)\\
\vdots\\
\end{array}
\right. ,
$$
where $d$ is the length of $\theta$.

repeat this step over all marginal and then
-  Invert Kendall's tau or Spearman's rho to obtain the parameter $\alpha$ of the copula.


Example with exponential marginal, we get
$$
\hat \lambda_n = \frac{1}{\overline X_n}
\textrm{~~and~~}
\hat \alpha_n = \frac{1}{1-\tau_n}.
$$

## Exact max_imum likelihood - `gumbel.EML`

When the density of the copula ex_ists, max_imum likelihood estimators can be used. To simplify, we assume that a bivariate copula with a density is used, and that the marginal distributions have densities. 
Let $\theta_1$ and $\theta_2$ denote the parameters of the marginal laws. The log-likelihood is written as:
\begin{multline*}
\ln \mathcal{L}(\alpha, \theta_1, \theta_2, x_1, \dots, x_n, y_1, \dots, y_n) 
= \sum_{i=1}^n \ln\left( c \left( F_1(x_i, \theta_1), F_2(y_i, \theta_2), \alpha \right) \right) \\
    \sum_{i=1}^n \ln\left( f_1(x_i, \theta_1) \right) + \sum_{i=1}^n \ln\left( f_2(y_i, \theta_2) \right).
    \end{multline*}
    Often, explicit expressions for the estimators that max_imize $\ln \mathcal{L}$ do not ex_ist, and a numerical max_imization is performed instead.

## Inference for margins - `gumbel.IFM`

Under the assumption that the copula has a density, it is possible to combine the two previous approaches by first estimating the parameters of the marginal distributions, then estimating the copula parameter. This involves:

- Estimating the parameters $\theta_1$ and $\theta_2$ via max_imum likelihood
- Constructing pseudo-data $\forall 1 \leq i \leq n,; u_i = F_1(x_i, \hat{\theta}_1)$ and $v_i = F_2(y_i, \hat{\theta}_2)$
- Estimating the parameter $\alpha$ by maximizing the log-likelihood
    $$
    \ln \mathcal{L}(\alpha, u1, \dots, u_n, v_1, \dots, v_n) = \sum_{i=1}^n \ln\left( c \left( u_i, v_i, \alpha \right) \right).
    $$

This method has the advantage of using the "classical" max_imum likelihood estimators for the marginals.

## Canonical max_imum likelihood - `gumbel.CML`

It is a semi-parametric method based on the previous approach:

- Calculate the empirical distribution functions $F_{1,n}$ and $F_{2,n}$
- Construct the pseudo-data $\forall 1 \leq i \leq n,; u_i = F_{1,n}(x_i)$ and $v_i = F_{2,n}(y_i)$
- Estimate the parameter $\alpha$ by maximizing the log-likelihood
    $$
    \ln \mathcal{L}(\alpha, u1, \dots, u_n, v_1, \dots, v_n) = \sum_{i=1}^n \ln\left( c \left( u_i, v_i, \alpha \right) \right).
    $$




# Bibliography