The COPULA Procedure

Archimedean Copulas

Subsections:

Overview of Archimedean Copulas

Let function $\phi : [0,1] \rightarrow [0, \infty )$ be a strict Archimedean copula generator function and suppose its inverse $\phi ^{-1}$ is completely monotonic on $[0, \infty )$. A strict generator is a decreasing function $\phi :[0,1]\rightarrow [0, \infty )$ that satisfies $\phi (0)=\infty $ and $\phi (1)=0$. A decreasing function $f(t):[a,b]\rightarrow (-\infty ,\infty )$ is completely monotonic if it satisfies

\[ (-1)^ k \frac{d^ k}{dt^ k} f(t)\ge 0, k\in \mathbb {N}, t\in (a,b) \]

An Archimedean copula is defined as follows:

\[ C(u_1, u_2,\ldots , u_ m) = \phi ^{-1}\Bigl ( \phi (u_1) + \cdots + \phi (u_ m) \Bigr ) \]

The Archimedean copulas available in the COPULA procedure are the Clayton copula, the Frank copula, and the Gumbel copula.

Clayton Copula

Let the generator function $\phi (u) = { \theta }^{-1} \left(u^{-\theta } -1\right)$. A Clayton copula is defined as

\[ C_{\theta }(u_1, u_2,{\ldots }, u_ m) = \left[ \sum _{i=1}^ m u_ i^{-\theta } - m+1\right] ^{-1/\theta } \]

with $\theta > 0$.

Frank Copula

Let the generator function be

\[ \phi (u) = - \log \left[ \frac{\exp (-{\theta }u)-1}{\exp (-{\theta })-1}\right] \]

A Frank copula is defined as

\[ C_{\theta }(u_1, u_2,{\ldots }, u_ m) = \frac{1}{\theta } \log \left\{ 1 + \frac{\prod _{i=1}^{m}[\exp (-{\theta }u_ i)-1]}{[\exp (-{\theta })-1]^{m-1}} \right\} \]

with $\theta \in (-\infty ,\infty ) \backslash \{ 0\} $ for $m=2$ and $\theta >0$ for $m\ge 3$.

Gumbel Copula

Let the generator function $\phi (u) = (-\log u) ^{\theta }$. A Gumbel copula is defined as

\[ C_{\theta }(u_1, u_2,{\ldots }, u_ m) = \exp \left\{ - \left[ \sum _{i=1}^{m}(-\log u_ i)^\theta \right] ^{1/{\theta }} \right\} \]

with $\theta > 1$.

Simulation

Suppose the generator of the Archimedean copula is $\phi $. Then the simulation method using Laplace-Stieltjes transformation of the distribution function is given by Marshall and Olkin (1988) where $\tilde{F}(t)= \int _0^\infty e^{-t x}dF(x) $:

  1. Generate a random variable $V$ with the distribution function $F$ such that $\tilde{F}(t)= \phi ^{-1}(t)$.

  2. Draw samples from independent uniform random variables $X_1,\ldots , X_ m$.

  3. Return $\bm U= (\tilde{F}(-\log (X_1)/V),\ldots \tilde{F}(-\log (X_ m)/V))^ T$.

The Laplace-Stieltjes transformations are as follows:

  • For the Clayton copula, $\tilde{F}= (1+t)^{-1/\theta }$, and the distribution function $F$ is associated with a Gamma random variable with shape parameter $\theta ^{-1}$ and scale parameter one.

  • For the Gumbel copula, $\tilde{F} = \exp (-t^{1/\theta })$, and $F$ is the distribution function of the stable variable $\textrm{St}(\theta ^{-1},1,\gamma ,0)$ with $\gamma = [\cos (\pi /(2\theta ))]^\theta $.

  • For the Frank copula with $\theta >0$, $\tilde{F}= - \log \{ 1-\exp (-t)[1- \exp (-\theta )]\} /\theta $, and ${F}$ is a discrete probability function $P(V=k)=(1-\exp (-\theta ))^ k/(k\theta )$. This probability function is related to a logarithmic random variable with parameter value $1-e^{-\theta }$.

For details about simulating a random variable from a stable distribution, see Theorem 1.19 in Nolan (2010). For details about simulating a random variable from a logarithmic series, see Chapter 10.5 in Devroye (1986).

For a Frank copula with $m=2$ and $\theta <0$, the simulation can be done through conditional distributions as follows:

  1. Draw independent $v_1,v_2$ from a uniform distribution.

  2. Let $u_1= v_1$.

  3. Let $u_2 = -\frac1\theta \log \left(1+\frac{v_2(1-e^{-\theta })}{v_2(e^{-\theta v_1}-1)-e^{-\theta v_1}}\right)$.

Fitting

One method to estimate the parameters is to calibrate with Kendall’s tau. The relation between the parameter $\theta $ and Kendall’s tau is summarized in the following table for the three Archimedean copulas.

Table 10.2: Calibration Using Kendall’s Tau

Copula Type

$\tau $

Formula for $\theta $

Clayton

$\theta /(\theta +2)$

$2\tau /(1-\tau )$

Gumbel

$1-1/\theta $

$ 1/(1-\tau )$

Frank

$1-4\theta ^{-1}(1-D_1(\theta ))$

No closed form


In Table 10.2, $D_1(\theta )=\theta ^{-1}\int _0^\theta t/(\exp (t)-1)dt$ for $\theta >0$, and $D_1(\theta ) =D_1(\theta )+0.5\theta $ for $\theta <0$. In addition, for the Frank copula, the formula for $\theta $ has no closed form. The numerical algorithm for root finding can be used to invert the function $\tau (\theta )$ to obtain $\theta $ as a function of $\tau $.

Alternatively, you can use the MLE or the CMLE method to estimate the parameter $\theta $ given the data $\mathbf{u}=\{ u_{i,j}\} $ and $i=1,\ldots ,n,j=1,\ldots ,m$. The log-likelihood function for each type of Archimedean copula is provided in the following sections.

Fitting the Clayton Copula
For the Clayton copula, the log-likelihood function is as follows (Cherubini, Luciano, and Vecchiato, 2004, Chapter 7):

\begin{align*} l & = n\left[ m \log (\theta ) + \log \left(\Gamma \left(\frac{1}{\theta } +m \right) \right) -\log \left(\Gamma \left(\frac{1}{\theta }\right) \right) \right] - (\theta +1) \sum _{i,j} \log u_{ij} \\ & - \left(\frac{1}{\theta } +m \right) \sum _ i \log \left( \sum _ j u_{ij} ^{-\theta } -m +1 \right) \end{align*}

Let $g(\cdot )$ be the derivative of $\log (\Gamma (\cdot ))$. Then the first order derivative is

\begin{align*} \frac{d l}{d \theta } & = n \left[ \frac{m}{\theta }+ g\left(\frac{1}{\theta }+m\right) \frac{-1}{\theta ^2}- g\left(\frac{1}{\theta }\right)\frac{-1}{\theta ^2}\right] \\ & - \sum _{i,j} \log (u_{ij}) +\frac{1}{\theta ^2} \sum _ i \log \left( \sum _ j u_{ij} ^{-\theta } -m +1 \right) \\ & -\left(\frac{1}{\theta }+m \right) \sum _ i \frac{- \sum _ j u_{ij}^{-\theta }\log (u_{ij})}{\sum _ j u_{ij}^{-\theta } - m +1}\\ \end{align*}

The second order derivative is

\begin{align*} \frac{d^2 l}{d \theta ^2} & = n\left\{ \frac{-m}{\theta ^2}+ g’\left(\frac1\theta +m\right)\frac{1}{\theta ^4} + g\left(\frac1\theta +m\right) \frac{2}{\theta ^3} - g’\left(\frac{1}{\theta }\right)\frac{1}{\theta ^4} - g\left(\frac{1}{\theta }\right)\frac{2}{\theta ^3} \right\} \\ & -\frac{2}{\theta ^3}\sum _ i \log \left(\sum _ j u_{ij}^{-\theta } -m+1 \right) \\ & +\frac{2}{\theta ^2}\sum _ i \frac{-\sum _ j u_{ij}^{-\theta }\log u_{ij}}{\sum _ j u_{ij}^{-\theta }-m+1} \\ & - \left(\frac1\theta +m \right) \sum _ i \left\{ \frac{\sum _ j u_{ij}^{-\theta }(\log u_{ij})^2}{\sum _ j u_{ij}^{-\theta }-m+1 } -\left(\frac{\sum _ j u_{ij}^{-\theta }\log u_{ij}}{\sum _ j u_{ij}^{-\theta }-m+1 }\right)^2 \right\} \end{align*}

Fitting the Gumbel Copula

A different parameterization $\alpha =\theta ^{-1}$ is used for the following part, which is related to the fitting of the Gumbel copula. For Gumbel copula, you need to compute $\phi ^{-1\left( m\right) }$. It turns out that for $k=1,2,\ldots ,m$,

\[ \phi ^{-1(k)}(u)= (-1)^ k\alpha \exp (-u^\alpha ) u^{-k+\alpha }\Psi _{k-1}(u^\alpha ) \]

where $\Psi _{k-1}$ is a function that is described later. The copula density is given by

\begin{eqnarray*} c& =& \phi ^{-1\left( m\right) }\left( x\right) \prod \limits _{k}\phi ^{\prime }\left( u_{k}\right)\\ & =& \left( -1\right) ^{m}\alpha \exp \left( -x^{\alpha }\right) x^{ -k+\alpha }\Psi _{m-1}\left( x^{\alpha }\right) \prod \limits _{k}\phi ^{\prime }\left( u_{k}\right) \\ & =& \left( -1\right) ^{m}f_{1}f_{2}f_{3}f_{4}f_{5} \end{eqnarray*}

where $ x=\sum _{k}\phi \left( u_{k}\right) $, $f_1=\alpha $, $f_2=\exp (-x^\alpha )$,$f_3=x^{-k+\alpha } $,$f_4=\Psi _{m-1}(x^\alpha )$, and $f_5=(-1)^ m\prod _ k\phi ’(u_ k)$.

The log density is

\begin{eqnarray*} l & =& \log (c) \\ & =& \log \left( f_{1}\right) +\log \left( f_{2}\right) +\log \left( f_{3}\right) +\log \left( f_{4}\right) +\log \left( \left( -1\right) ^{m}f_{5}\right) \end{eqnarray*}

Now the first order derivative of the log density has the decomposition

\begin{eqnarray*} \frac{dl}{d\alpha } & =& \frac{1}{c}\frac{dc}{d\alpha }=\sum _{j=1}^{4}\frac{1}{f_{j}}\frac{df_{j}}{d\alpha }+\frac{d\sum _{k}\log \left( -\phi ^{\prime }\left( u_{k}\right) \right) }{d\alpha } \end{eqnarray*}

Some of the terms are given by

\begin{eqnarray*} \frac{1}{f_{1}}\frac{df_{1}}{d\alpha } & =& \frac{1}{\alpha } \\ \frac{1}{f_{2}}\frac{df_{2}}{d\alpha } & =& -x^{\alpha }\log \left( x\right) -\alpha x^{\alpha -1}\frac{dx}{d\alpha } \\ \frac{1}{f_{3}}\frac{df_{3}}{d\alpha } & =& \log \left( x\right) +\left( -k+\alpha \right) x^{-1}\frac{dx}{d\alpha } \end{eqnarray*}

where

\[ \frac{dx}{d\alpha }=\sum \left( -\log u_{k}\right) ^{1/\alpha }\log \left( -\log u_{k}\right) \left( \frac{-1}{\alpha ^{2}}\right) \]

The last term in the derivative of the $dl/d\alpha $ is

\begin{eqnarray*} \log \left( -\phi ^{\prime }\left( u_{k}\right) \right) & =& \log \left( \frac{1}{\alpha }\left( -\log u_{k}\right) ^{\frac{1}{\alpha }-1}\frac{1}{u_{k}}\right) \\ & =& -\log \alpha -\log \left( u_{k}\right) +\left( \frac{1}{\alpha }-1\right) \log \left( -\log \left( u_{k}\right) \right) \\ \frac{d\sum _{k}\log \left( -\phi ^{\prime }\left( u_{k}\right) \right) }{d\alpha } & =& \sum _{k=1}^{m}-\frac{1}{\alpha }-\frac{1}{\alpha ^{2}}\log \left( -\log \left( u_{k}\right) \right) \\ & =& -\frac{m}{\alpha }-\frac{1}{\alpha ^{2}}\sum _{k=1}^{m}\log \left( -\log \left( u_{k}\right) \right) \end{eqnarray*}

Now the only remaining term is $f_4$, which is related to $\Psi _{m-1}$. Wu, Valdez, and Sherris (2007) show that $\Psi _ k(x)$ satisfies a recursive equation

\[ \Psi _ k(x) =[\alpha (x-1)+k]\Psi _{k-1}(x)-\alpha x \Psi _{k-1}’(x) \]

with $\Psi _0(x)=1$.

The preceding equation implies that $\Psi _{k-1}(x)$ is a polynomial of x and therefore can be represented as

\[ \Psi _{k-1}\left( x\right) =\sum _{j=0}^{k-1}a_{j}\left( k-1,\alpha \right) x^{j} \]

In addition, its coefficient, denoted by $a_{j}(k-1,\alpha )$, is a polynomial of $\alpha $. For simplicity, use the notation $a_ j(\alpha )\equiv a_ j(m-1,\alpha ) $. Therefore,

\[ f_4 = \Psi _{m-1}\left( x^{\alpha }\right) =\sum _{j=0}^{m-1}a_{j}\left( \alpha \right) x^{j\alpha } \]
\begin{align*} \frac{df_{4}}{d\alpha }& =\frac{d\Psi _{m-1}\left( x^{\alpha }\right) }{d\alpha }\\ & =\sum _{j=0}^{m-1}\left[ \frac{da_{j}\left( \alpha \right) }{d\alpha }x^{j\alpha }+a_{j}\left( \alpha \right) x^{j\alpha }\log \left( x\right) j+a_{j}\left( \alpha \right) \left( j\alpha \right) x^{j\alpha -1}\frac{dx}{d\alpha }\right] \end{align*}

Fitting the Frank copula


For the Frank copula,

\[ \phi ^{-1(k)}(u)= -\frac{1}{\theta }\Psi _{k-1}\left( (1+e^{-u}(e^{-\theta }-1))^{-1}\right) \]

When $\theta >0$, a Frank copula has a probability density function

\begin{eqnarray*} c & =& \varphi ^{-1\left( m\right) }\left( x\right) \prod \limits _{k}\varphi ^{\prime }\left( u_{k}\right) \\ & =& \frac{-1}{\theta }\Psi _{m-1}\left( \frac{1}{1+e^{-x}\left( e^{-\theta }-1\right) }\right) \prod \limits _{k}\varphi ^{\prime }\left( u_{k}\right) \end{eqnarray*}

where $x=\sum _{k}\varphi \left( u_{k}\right) $.

The log likelihood is

\[ \log c=-\log \left( \theta \right) +\log \left( \Psi _{m-1}\left( \frac{1}{1+e^{-x}\left( e^{-\theta }-1\right) }\right) \right) +\sum \log \left( \varphi ^{\prime }\left( u_{k}\right) \right) \]

Denote

\[ y=\frac{1}{1+e^{-x}\left( e^{-\theta }-1\right) } \]

Then the derivative of the log likelihood is

\begin{eqnarray*} \frac{d\log c}{d\theta } & =& -\frac{1}{\theta }+\frac{1}{ \Psi _{m-1}\left(y\right) }\frac{d\Psi _{m-1}}{d\theta }+\sum _ k \frac{1}{\varphi ^{\prime }\left( u_{k}\right) }\frac{d\varphi ^{\prime }\left( u_{k}\right) }{d\theta } \end{eqnarray*}

The term in the last summation is

\begin{eqnarray*} \frac{1}{\varphi ^{\prime }\left( u_{k}\right) }\frac{d\varphi ^{\prime }\left( u_{k}\right) }{d\theta } & =& \frac{1}{\theta \left( 1-e^{\theta u_{k}}\right) }\left[ 1-e^{\theta u_{k}}+\theta ue^{\theta u_{k}}\right] \end{eqnarray*}

The function $\Psi _{m-1}$ satisfies a recursive relation

\[ \Psi _ k(x)= x(x-1)\Psi _{k-1}’(x) \]

with $\Psi _0(x)=x-1$. Note that $\Psi _{m-1}$ is a polynomial whose coefficients do not depend on $\theta $; therefore,

\begin{eqnarray*} \frac{d\Psi _{m-1}}{d\theta } & =& \frac{d\Psi _{m-1}}{dy}\frac{dy}{d\theta } \\ & =& \frac{d\Psi _{m-1}}{dy}\left[ \frac{dy}{d\theta }+\frac{dy}{dx}\frac{dx}{d\theta }\right] \\ & =& \frac{d\Psi _{m-1}}{dy}\left[ \frac{e^{-x}e^{-\theta }}{\left[ 1+e^{-x}\left( e^{-\theta }-1\right) \right] ^{2}}+\frac{e^{-x}\left( e^{-\theta }-1\right) }{\left[ 1+e^{-x}\left( e^{-\theta }-1\right) \right] ^{2}}\frac{dx}{d\theta }\right] \end{eqnarray*}

where

\begin{align*} \frac{dx}{d\theta } & =\sum _ k \frac{d\varphi \left( u_{k}\right) }{d\theta }=\sum _ k \left[ -\frac{u_ k e^{-\theta u_ k}}{1-e^{-\theta u_ k}}+\frac{e^{-\theta }}{1-e^{\theta }}\right] \\ & =\sum _ k \left[ -\frac{u_{k}}{e^{\theta u_{k}}-1}+\frac{1}{e^{\theta }-1}\right] \end{align*}

For the case of $m=2$ and $\theta <0$, the bivariate density is

\[ \log c = \log (\theta (1-e^{-\theta }))-\theta (u_1+u_2) -\log ((1-e^{-\theta }-(1-e^{-\theta u_1})(1-e^{-\theta u_2}))^2) \]