# Archimedean Crew

We have recently seen[1][2][3] how we can define dependencies between random variables with Archimedean copulas which calculate the probability that they each fall below given values by applying a generator function $$\phi$$ to the results of their cumulative distribution functions, or CDFs, for those values, and applying its inverse to their sum.
Like all copulas they are effectively the CDFs of vector valued random variables whose elements are uniformly distributed when considered independently. Whilst those Archimedean CDFs were relatively trivial to implement, we found that their probability density functions, or PDFs, were somewhat more difficult and that the random variables themselves required some not at all obvious mathematical manipulation to get right.
Having done all the hard work implementing the ak.archimedeanCopula, ak.archimedeanCopulaDensity and ak.archimedeanCopulaRnd functions we shall now use them to implement some specific families of Archimedean copulas.

Formally, given a continuous generator function $$\phi$$ and its inverse $$\phi^{-1}$$ that satisfy
$\phi : [0,1] \mapsto [0,\infty]\\ \phi(1) = 0\\ x > y \rightarrow \phi(x) < \phi(y)\\ (-1)^k \frac{\mathrm{d}^k}{\mathrm{d}t^k} \phi^{-1}(t) \geqslant 0$
where $$\mapsto$$ means maps to, $$\rightarrow$$ means implies and $$\frac{\mathrm{d}^k}{\mathrm{d}t^k}\phi^{-1}(t)$$ is the $$k^\mathrm{th}$$ derivative of $$\phi^{-1}$$ evaluated at $$t$$, for $$k$$ less than or equal to $$n$$, the function
$C_\phi\left(\mathbf{u}\right) = \phi^{-1}\left(\min\left(\sum_{i=0}^{n-1} \phi(u_i), \phi(0)\right)\right)$
where $$\sum$$ is the summation sign, is an $$n$$ dimensional Archimedean copula.

I've previously stated that those conditions for $$\phi$$ are a little bit stricter than necessary and in this post we'll be depending on the precise conditions; specifically that $$\phi$$ is a continuous convex function satisfying
$\phi : [0,1] \mapsto [0,\infty]\\ \phi(1) = 0\\ x > y \rightarrow \phi(x) < \phi(y)\\ (-1)^k \frac{\mathrm{d}^k}{\mathrm{d}t^k} \phi^{-1}(t) \geqslant 0$
for $$k$$ from zero to $$n-2$$ with $$(-1)^{n-2} \frac{\mathrm{d}^{n-2}}{\mathrm{d}t^{n-2}} \phi^{-1}(t)$$ being a non-increasing convex function, where a function $$f$$ is convex if, given a linear function $$g$$ connecting the points $$\left(x_0, f\left(x_0\right)\right)$$ and $$\left(x_1, f\left(x_1\right)\right)$$, it satisfies
$x_0 \leqslant x \leqslant x_1 \rightarrow f(x) \leqslant g(x)$
Note that these conditions are equivalent to the stricter ones if the $$(n-1)^\mathrm{th}$$ and $$n^\mathrm{th}$$ derivatives of $$\phi^{-1}$$ are well-defined. In fact, the whole point of the precise conditions is that they cover those Archimedean copulas whose generators don't have well-defined higher derivatives of their inverses.

We shall be using Nelson's An Introduction to Copulas[4] as the source for Archimedean generators but first there are a few special cases that we need to define.

### Bounds And Independence

The Fréchet-Hoeffding Theorem states that the functions of an $$n$$ dimensional vector $$\mathbf{u}$$ with elements between zero and one
\begin{align*} W(\mathbf{u}) &= \max\left(1-n + \sum_{i=1}^{n-1} u_i, 0\right)\\ M(\mathbf{u}) &= \min\left(u_0, \dots, u_{n-1}\right) \end{align*}
are the lower and upper bounds for any $$n$$ dimensional copula $$C$$
$W(\mathbf{u}) \leqslant C(\mathbf{u}) \leqslant M(\mathbf{u})$
Note that, for any number of dimensions, $$M$$ is the CDF of a random variable $$\mathbf{U}$$ that satisfies
\begin{align*} U_0 &\sim U(0,1)\\ U_i &= U_0 & \forall i>0 \end{align*}
where $$\sim$$ means drawn from, $$U(0,1)$$ is the standard uniform distribution and $$\forall$$ means for all, and, for two dimensions, $$W$$ is the CDF of a random variable $$\mathbf{U}^\prime$$ satisfying
\begin{align*} U^\prime_0 &\sim U(0,1)\\ U^\prime_1 &= 1-U^\prime_0 \end{align*}
so that $$M$$ is a copula in any number of dimensions and $$W$$ is only a copula in two.
Furthermore, since both result in observations that lie upon a line, their densities will be infinite along those lines and zero everywhere else.

The ak.lowerCopula function given in listing 1 implements $$W$$, treating zero and one dimensions as special cases.

Listing 1: ak.lowerCopula
ak.lowerCopula = function(n) {
var f;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0 || n>2) {
throw new Error('invalid dimension in ak.lowerCopula');
}

switch(n) {
case 0: f = function(u) {
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==0) {
throw new Error('invalid argument in ak.lowerCopula');
}
return 1;
};
break;

case 1: f = function(u) {
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==1) {
throw new Error('invalid argument in ak.lowerCopula');
}
return Math.max(0, Math.min(1, u.at(0)));
};
break;

case 2: f = function(u) {
var u0, u1;
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==2) {
throw new Error('invalid argument in ak.lowerCopula');
}
u0 = Math.max(0, Math.min(1, u.at(0)));
u1 = Math.max(0, Math.min(1, u.at(1)));
return Math.max(u0+u1-1, 0);
};
break;
}

f.dims = function() {return n;};

return Object.freeze(f);
};


This function, together with its associated ak.lowerCopulaDensity and ak.lowerCopulaRnd functions, are implemented in LowerCopula.js.
Similarly, $$M$$ is implemented by the ak.upperCopula function in listing 2

Listing 2: ak.upperCopula
ak.upperCopula = function(n) {
var f;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0) {
throw new Error('invalid dimension in ak.upperCopula');
}

switch(n) {
case 0:  f = function(u) {
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==0) {
throw new Error('invalid argument in ak.upperCopula');
}
return 1;
};
break;

case 1:  f = function(u) {
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==1) {
throw new Error('invalid argument in ak.upperCopula');
}
return Math.max(0, Math.min(1, u.at(0)));
};
break;

default: f = function(u) {
var ui, i;
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==n) {
throw new Error('invalid argument in ak.upperCopula');
}
ui = u.at(0);
for(i=1;i<n;++i) ui = Math.min(ui, u.at(i));
return Math.max(0, Math.min(1, ui));
};
break;
}

f.dims = function() {return n;};

return Object.freeze(f);
};


which can be found in UpperCopula.js along with the ak.upperCopulaDensity and ak.upperCopulaRnd functions.

Program 1 plots both of these copulas, with the upper bound on top and the lower bound beneath it.

Program 1: The Lower And Upper Bounds

Between these two extremes lie independent random variables, which are represented by the copula
$\Pi(\mathbf{u}) = \prod_{i=0}^{n-1} u_i$
where the big $$\prod$$ is the product sign, as implemented by ak.independentCopula in listing 3.

Listing 3: ak.independentCopula
ak.independentCopula = function(n) {
var f;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0) {
throw new Error('invalid dimension in ak.independentCopula');
}

f = function(u) {
var p = 1.0;
var i;
if(ak.type(u)!==ak.VECTOR_T || u.dims()!==n) {
throw new Error('invalid argument in ak.independentCopula');
}
for(i=0;i<n;++i) p *= Math.max(0, Math.min(1, u.at(i)));
return p;
};

f.dims = function() {return n;};

return Object.freeze(f);
};


This can be found with the ak.independentCopulaDensity and ak.independentCopulaRnd functions in IndependentCopula.js and is demonstrated by program 2.

Program 2: The Independence Copula

When it comes to implementing Archimedean copulas, these three copulas will frequently turn up as special cases.

### The Gumbel Copula

So far we've been using the generator and inverse
\begin{align*} \phi_\theta(u) &= (-\ln u)^\theta\\ \phi_\theta^{-1}(t) &= e^{-t^{\frac{1}{\theta}}} \end{align*}
as the example for our Archimedean copula functions, which I've asserted yield a valid copula for $$\theta$$ greater than or equal to one. That copula is known as the Gumbel copula.

The natural logarithm $$\ln$$ is a continuous function with the properties
$\ln(0) = -\infty\\ \ln(1) = 0\\ x > y \rightarrow \ln(x) > \ln(y)$
which imply that the Gumbel generator satisfies the first four of the stricter requirements for defining a valid copula whenever $$\theta$$ is greater than zero. For the last we trivially have
$(-1)^0 \frac{\mathrm{d}^0}{\mathrm{d}t^0}\phi_\theta^{-1}(t) = \phi_\theta^{-1}(t) = e^{-t^{\frac{1}{\theta}}} \geqslant 0$
and derivation 1 proves that it is satisfied for all higher derivatives if $$\theta$$ is greater than or equal to one.

Derivation 1: The Derivatives Of The Gumbel Generator's Inverse
Firstly, by the chain rule, we have
$\frac{\mathrm{d}}{\mathrm{d}t} \phi_\theta^{-1}(t) = e^{-t^{\frac{1}{\theta}}} \left(-\frac{1}{\theta} t^{\frac{1}{\theta}-1}\right) = -\phi_\theta^{-1}(t) \frac{1}{\theta} t^{\frac{1}{\theta}-1}$
which will be negative if $$\theta$$ is greater than zero.

Next, let use assume that
$\frac{\mathrm{d}^k}{\mathrm{d}t^k} \phi_\theta^{-1}(t) = (-1)^k \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \, t^{\frac{i}{\theta}-k}$
for some non-negative constants $$c_{k,i}$$, which is clearly the case for $$k$$ equal to one when $$\theta$$ is greater than zero. By the product rule
$\frac{\mathrm{d}^{k+1}}{\mathrm{d}t^{k+1}} \phi_\theta^{-1}(t) = (-1)^k \sum_{i=1}^k c_{k,i} \, t^{\frac{i}{\theta}-k} \frac{\mathrm{d}}{\mathrm{d}t} \phi_\theta^{-1}(t) + (-1)^k \phi_\theta^{-1}(t) \frac{\mathrm{d}}{\mathrm{d}t} \sum_{i=1}^k c_{k,i} \, t^{\frac{i}{\theta}-k}\\ \quad \quad = (-1)^k \sum_{i=1}^k c_{k,i} \, t^{\frac{i}{\theta}-k} \left(-\phi_\theta^{-1}(t) \left(\frac{1}{\theta} t^{\frac{1}{\theta}-1}\right)\right) + (-1)^k \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \left(\frac{i}{\theta}-k\right) t^{\frac{i}{\theta}-k-1}\\ \quad \quad = (-1)^{k+1} \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \frac{1}{\theta} t^{\frac{i+1}{\theta}-k-1} + (-1)^{k+1} \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \left(k-\frac{i}{\theta}\right) t^{\frac{i}{\theta}-k-1}\\ \quad \quad = (-1)^{k+1} \phi_\theta^{-1}(t) \sum_{i=2}^{k+1} c_{k,i-1} \frac{1}{\theta} t^{\frac{i}{\theta}-(k+1)} + (-1)^{k+1} \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \left(k-\frac{i}{\theta}\right) t^{\frac{i}{\theta}-(k+1)}$
Now, $$c_{k,i-1} \frac{1}{\theta}$$ is positive and finite for $$\theta$$ greater than zero and $$\left(k-\frac{i}{\theta}\right)$$ is non-negative for $$i$$ from one to $$k$$ if $$\theta$$ is greater than or equal to one, so that in the latter case
$\frac{\mathrm{d}^{k+1}}{\mathrm{d}t^{k+1}} \phi_\theta^{-1}(t) = (-1)^{k+1} \phi_\theta^{-1}(t) \sum_{i=1}^{k+1} c_{k+1,i} \, t^{\frac{i}{\theta}-(k+1)}$
for non-negative constants $$c_{k+1,i}$$ and, by induction, the assumption will be true for all $$k$$ greater than zero which means that
$(-1)^k \frac{\mathrm{d}^k}{\mathrm{d}t^k} \phi_\theta^{-1}(t) = (-1)^{2k} \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \, t^{\frac{i}{\theta}-k} = \phi_\theta^{-1}(t) \sum_{i=1}^k c_{k,i} \, t^{\frac{i}{\theta}-k} \geqslant 0$
as required.

Note that for $$\theta$$ equal to one the Gumbel copula is the independent copula $$\Pi$$
$C_1(\mathbf{u}) = e^{-\sum_{i=0}^{n-1} -\ln u_i} = e^{\sum_{i=0}^{n-1} \ln u_i} = \prod_{i=0}^{n-1} u_i = \Pi(\mathbf{u})$
and that for $$\theta$$ equal to infinity, it's the upper bound $$M$$
\begin{align*} C_\infty(\mathbf{u}) &= e^{-\left(\sum_{i=0}^{n-1} \left(-\ln u_i\right)^\infty\right)^\frac{1}{\infty}} = e^{-\max\left(-\ln u_0, \dots, -\ln u_{n-1}\right)}\\ &= e^{\min\left(\ln u_0, \dots, \ln u_{n-1}\right)} = \min\left(u_0, \dots, u_{n-1}\right) = M(\mathbf{u}) \end{align*}
The Gumbel generator, its derivative, its inverse and its inverse's $$n^\mathrm{th}$$ derivative are implemented in listing 4 with p, dp, q and dnq respectively.

Listing 4: The Gumbel Generator, Inverse And Derivatives
function p(theta) {
return function(u) {return Math.pow(-Math.log(u), theta);};
}

function dp(theta) {
return function(u) {return -theta*Math.pow(-Math.log(u), theta-1)/u;};
}

function q(theta) {
return function(t) {return Math.exp(-Math.pow(t, 1/theta));};
}

function dnq(theta) {
return function(t, n) {return ak.exp(ak.neg(ak.pow(ak.surreal(n, t, 1), 1/theta)));};
}


Note that we're again using our ak.surreal[5] type to calculate all $$n$$ derivatives of the inverse with a single expression.

The ak.gumbelCopula function given in listing 5 switches between the independence copula, the upper bound and an ak.archimedeanCopula initialised with p and q depending on the dimension and the value of theta.

Listing 5: ak.gumbelCopula
ak.gumbelCopula = function(n, theta) {
var f, g;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0) {
throw new Error('invalid dimension in ak.gumbelCopula');
}

if(ak.nativeType(theta)!==ak.NUMBER_T || isNaN(theta) || (n>1 && theta<1)) {
throw new Error('invalid theta in ak.gumbelCopula');
}

if(n<2 || theta===1) g = ak.independentCopula(n);
else if(theta===ak.INFINITY) g = ak.upperCopula(n);
else g = ak.archimedeanCopula(n, p(theta), q(theta));

f = function(u) {return g(u);};
f.dims = function() {return n;};
f.theta = function() {return theta;};

return Object.freeze(f);
};


The ak.gumbelCopulaDensity function given in listing 6 takes the same approach

ak.gumbelCopulaDensity = function(n, theta) {
var f, g;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0) {
throw new Error('invalid dimension in ak.gumbelCopulaDensity');
}

if(ak.nativeType(theta)!==ak.NUMBER_T || isNaN(theta) || (n>1 && theta<1)) {
throw new Error('invalid theta in ak.gumbelCopulaDensity');
}

if(n<2 || theta===1) g = ak.independentCopulaDensity(n);
else g = ak.archimedeanCopulaDensity(n, p(theta), dp(theta), dnq(theta));

f = function(u) {return g(u);};
f.dims = function() {return n;};
f.theta = function() {return theta;};

return Object.freeze(f);
};


as does the ak.gumbelCopulaRnd function implemented in listing 7.

Listing 7: ak.gumbelCopulaRnd
ak.gumbelCopulaRnd = function(n, theta, arg3, arg4, arg5) {
var m, threshold, rnd, f, g;

if(ak.nativeType(arg3)!==ak.NUMBER_T) rnd = arg3;
else if(ak.nativeType(arg4)!==ak.NUMBER_T) rnd = arg4;
else rnd = arg5;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0) {
throw new Error('invalid dimension in ak.gumbelCopulaRnd');
}

if(ak.nativeType(theta)!==ak.NUMBER_T || isNaN(theta) || (n>1 && theta<1)) {
throw new Error('invalid theta in ak.gumbelCopulaRnd');
}

if(ak.nativeType(rnd)===ak.UNDEFINED_T) {
rnd = Math.random;
}
else if(ak.nativeType(rnd)!==ak.FUNCTION_T) {
throw new Error('invalid random number generator in ak.gumbelCopulaRnd');
}

if(n<2 || theta===1) g = ak.independentCopulaRnd(n, rnd);
else if(theta===ak.INFINITY) g = ak.upperCopulaRnd(n, rnd);
else g = ak.archimedeanCopulaRnd(n, p(theta), q(theta), dnq(theta), arg3, arg4, arg5);

f = function() {return g();};
f.dims = function() {return n;};
f.theta = function() {return theta;};
f.rnd = function() {return rnd;};

return Object.freeze(f);
};


These functions are defined in GumbelCopula.js and program 3 uses this implementation to plot the, by now rather familiar, density of the Gumbel copula.

Program 3: The Gumbel Copula Density

### The Clayton Copula

The next Archimedean copula that we'll implement is the Clayton copula, which is defined by the generator and inverse
\begin{align*} \phi_\theta(u) &= \frac{1}{\theta}\left(u^{-\theta}-1\right)\\ \phi_\theta^{-1}(t) &= (1+\theta t)^{-\frac{1}{\theta}} \end{align*}
Derivation 2: The Lower Bound Of Clayton's θ
The straight line connecting $$\phi_\theta(0)$$ and $$\phi_\theta(1)$$ is trivially given by
$f(u) = \frac{1}{\theta}(u-1)$
If $$\theta$$ is negative then this is equal to
$f(u) = \frac{1}{|\theta|}(1-u)$
in which case we also have
$\phi_\theta(u) = \frac{1}{|\theta|}(1-u^{|\theta|})$
Now, $$\phi_\theta$$ will be concave if
\begin{align*} \phi_\theta(u) &> f(u)\\ \frac{1}{|\theta|}\left(1-u^{|\theta|}\right) &> \frac{1}{|\theta|}(1-u)\\ 1-u^{|\theta|} &> 1-u\\ u^{|\theta|} &< u\\ \ln u^{|\theta|} &< \ln u\\ |\theta| \ln u &< \ln u\\ |\theta| &> 1 \end{align*}
for $$u$$ less than one, since that implies that $$\ln u$$ is negative.
For negative $$\theta$$ we have
\begin{align*} \phi_\theta(u) &= \frac{1}{\theta}\left(u^{|\theta|}-1\right)\\ \phi_\theta(0) &= \frac{1}{\theta}\left(0^{|\theta|}-1\right) = -\frac{1}{\theta} \end{align*}
where the vertical bars stand for the absolute value of the term between them, and so the inverse will instead be given by
$\phi_\theta^{-1}(t) = \begin{cases} (1+\theta t)^{-\frac{1}{\theta}} & t \leqslant -\frac{1}{\theta}\\ 0 & t \geqslant -\frac{1}{\theta} \end{cases}$
In particular, when $$\theta$$ is equal to minus one the inverse and its derivative take the forms
\begin{align*} \phi_{-1}^{-1}(t) &= \begin{cases} 1-t & t \leqslant 1\\ 0 & t \geqslant 1 \end{cases}\\ \frac{\mathrm{d}}{\mathrm{d}t}\phi_{-1}^{-1}(t) &= \begin{cases} -1 & t \leqslant 1\\ \phantom{-}0 & t \geqslant 1 \end{cases} \end{align*}
and so the latter is undefined at one. In fact this yields the lower bound copula and requires the weaker conditions on the generator to demonstrate that it's valid in two dimensions.
Note that $$\theta$$ cannot be less than minus one because that would imply that $$\phi_\theta$$ is concave, as proven by derivation 2.

It's fairly easy to calculate that, for $$\theta$$ greater than minus one, the higher derivatives of $$\phi_\theta^{-1}$$ are given by
$\frac{\mathrm{d}^k}{\mathrm{d}t^k} \phi_\theta^{-1}(t) = (-1)^k \times \left(1+\theta t\right)^{-\frac{1}{\theta}-k} \times \prod_{i=1}^k \left(1+(i-1)\theta\right)$
If $$\theta$$ is positive then this trivially satisfies the stricter conditions for the definition of an Archimedean copula.

Derivation 3: The Concavity Of f
If $$-\frac{1}{\theta}-k$$ is positive, we have
\begin{align*} f_\theta(0) &= 1\\ f_\theta\left(-\tfrac{1}{\theta}\right) &= 0 \end{align*}
which are connected by the linear function
$g_\theta(t) = 1 + \theta t$
Now, $$f_\theta$$ is concave if
\begin{align*} \left(1+\theta t\right)^{-\frac{1}{\theta}-k} &> 1 + \theta t\\ \left(-\frac{1}{\theta}-k\right) \times \ln \left(1+\theta t\right) &> \ln \left(1 + \theta t\right)\\ -\frac{1}{\theta}-k &< 1\\ \theta &< -\frac{1}{k+1} \end{align*}
Finally, we have
$\theta = -\frac{1}{k+1} \rightarrow -\frac{1}{\theta}-k = 1$
satisfying the initial assumption.
If it's negative, derivation 3 proves that the function
$f_\theta(t) = \left(1+\theta t\right)^{-\frac{1}{\theta}-k}$
is concave for $$\theta$$ less than $$-\frac{1}{k+1}$$ and, furthermore, for $$\theta$$ greater than or equal to $$-\frac{1}{k+1}$$, we have
$1 + (i-1) \theta \geqslant 1 - \frac{i-1}{k+1}$
which is clearly positive for $$i$$ between one and $$k$$.

When $$\theta$$ equals zero, the generator's inverse can be expressed as
\begin{align*} \theta &= -\frac{1}{n}\\ \phi_0^{-1}(t) &= \lim_{n \rightarrow \infty} \left(1 - \frac{t}{n}\right)^n \end{align*}
which is the original definition of $$e^{-t}$$, so that
\begin{align*} \phi_0(u) &= -\ln u\\ \phi_0^{-1}(t) &= e^{-t} \end{align*}
and the Clayton copula is equivalent to the independent copula.

Finally, for $$\theta$$ equal to infinity it defines the copula
\begin{align*} C_\infty(\mathbf{u}) &= \lim_{\theta \rightarrow \infty}\left(1+\theta \sum_{i=0}^{n-1} \frac{1}{\theta}\left(u_i^{-\theta}-1\right) \right)^{-\frac{1}{\theta}} = \lim_{\theta \rightarrow \infty}\left(1-n + \sum_{i=0}^{n-1} u_i^{-\theta} \right)^{-\frac{1}{\theta}}\\ &= \lim_{\theta \rightarrow \infty}\left(\max\left(u_0^{-\theta},\dots,u_{n-1}^{-\theta}\right)\right)^{-\frac{1}{\theta}} = \min\left(u_0,\dots,u_{n-1}\right) = M(\mathbf{u}) \end{align*}
since as $$\theta$$ tends to infinity the term involving the smallest $$u_i$$ will grow infinitely larger than all of the other terms. Unfortunately, if you look closely then you'll notice that I've not taken into account the possibility that two or more of the $$u_i$$ have the same smallest value. To justify this apparent oversight I shall start waving my hands vigorously, mumble something about infinitesimal perturbations and ask you to accept that it doesn't actually matter.

Putting all of these cases together means that the Clayton generator and its inverse satisfy the conditions for defining a valid copula in $$n$$ dimensions for $$\theta$$ greater than or equal to $$-\frac{1}{n-1}$$.

Listing 8 gives the Clayton generator p, derivative dp, inverse q, inverse derivative dnq and the ak.claytonCopula implementation of the copula which defers to our implementations of $$\Pi$$, $$W$$, $$M$$ and the general purpose Archimedean copula as appropriate.

Listing 8: ak.claytonCopula
function p(theta) {
return function(u) {return (Math.pow(u, -theta)-1)/theta;};
}

function dp(theta) {
return function(u) {return -Math.pow(u, -theta-1);};
}

function q(theta) {
return function(t) {return Math.pow(1+theta*t, -1/theta);};
}

function dnq(theta) {
return function(t, n) {
return ak.pow(ak.surreal(n, 1+theta*t, theta), -1/theta);
};
}

ak.claytonCopula = function(n, theta) {
var f, g;

if(ak.nativeType(n)!==ak.NUMBER_T || !isFinite(n) || n!==ak.floor(n) || n<0) {
throw new Error('invalid dimension in ak.claytonCopula');
}

if(ak.nativeType(theta)!==ak.NUMBER_T || isNaN(theta) || (n>1 && theta<-1/(n-1))) {
throw new Error('invalid theta in ak.claytonCopula');
}

if(n<2 || theta===0) g = ak.independentCopula(n);
else if(theta===-1) g = ak.lowerCopula(n);
else if(theta===ak.INFINITY) g = ak.upperCopula(n);
else g = ak.archimedeanCopula(n, p(theta), q(theta));

f = function(u) {return g(u);};
f.dims = function() {return n;};
f.theta = function() {return theta;};

return Object.freeze(f);
};


You can find these in ClaytonCopula.js, together with the random number generator ak.claytonCopulaRnd and the density ak.claytonCopulaDensity which is plotted by program 4.

Program 4: The Clayton Copula Density

### The Ali-Mikhail-Haq Copula

The Ali-Mikhail-Haq copula is defined by the generator and inverse
\begin{align*} \phi_\theta(u) &= \ln \frac{1-\theta(1-u)}{u}\\ \phi^{-1}_\theta(t) &= \frac{1-\theta}{e^t-\theta} \end{align*}
for $$\theta$$ in the interval $$[-1,1]$$ with the special cases
\begin{align*} \phi_0(u) &= -\ln u\\ \phi^{-1}_0(t) &= e^{-t} \end{align*}
which result in the independent copula, and
\begin{align*} \phi_1(u) &= \frac{1}{u} - 1\\ \phi^{-1}_1(t) &= \frac{1}{t+1} \end{align*}
Since we've already seen how to prove that a generator yields a valid copula, I shan't bother you with the maths from this point on. The implementations of further Archimedean copulas, their densities and their random vector generators aren't going to be much different from what we've seen so far either and so I'll just tell you what they're called and point you at the files containing them, which in this case are the ak.aliCopula, ak.aliCopulaDensity and ak.aliCopulaRnd functions in AliCopula.js (with all due apologies to Mikhail and Haq).
Program 5 plots the density of the Ali-Mikhail-Haq copula.

Program 5: The Ali-Mikhail-Haq Copula Density

### The Frank Copula

Next up is the Frank copula which is defined by
\begin{align*} \phi_\theta(u) &= -\ln \frac{e^{-\theta u}-1}{e^{-\theta}-1}\\ \phi^{-1}_\theta(t) &= -\frac{1}{\theta} \ln \left(1 + \left(e^{-\theta}-1\right)e^{-t}\right) \end{align*}
for any value of $$\theta$$ and has the special cases
\begin{align*} C_{-\infty} &= W\\ C_0 &= \Pi\\ C_\infty &= M \end{align*}
Note that it's valid in any number of dimensions except when $$\theta$$ equals minus infinity, when it's only valid in two.

The ak.frankCopula, ak.frankCopulaDensity and ak.frankCopulaRnd functions are defined in FrankCopula.js, and program 6 plots the density.

Program 6: The Frank Copula Density

### The Joe Copula

The Joe copula is defined by the generator and inverse
\begin{align*} \phi_\theta(u) &= -\ln \left(1-\left(1-u\right)^\theta\right)\\ \phi^{-1}_\theta(t) &= 1 - \left(1-e^{-t}\right)^{\frac{1}{\theta}} \end{align*}
and is valid in any number of dimensions for $$\theta$$ greater than or equal to one with the limiting cases
\begin{align*} C_1 &= \Pi\\ C_\infty &= M \end{align*}
Program 7 plots the ak.joeCopulaDensity implementation of its density which can be found, together with ak.joeCopula and ak.joeCopulaRnd, in JoeCopula.js.

Program 7: The Joe Copula Density

### The Gumbel-Barnett Copula

The generator and inverse
\begin{align*} \phi_\theta(u) &= \ln \left(1 - \theta \ln u\right)\\ \phi^{-1}_\theta(t) &= e^\tfrac{1-e^t}{\theta} \end{align*}
define the Gumbel-Barnett copula, which is valid for $$\theta$$ between zero and one and is equal to the independent copula at the lower bound.

This copula, its density and its random vector generator are implemented by the ak.barnettCopula, ak.barnettDensity and ak.barnettRnd functions in BarnettCopula.js and program 8 plots the second of those.

Program 8: The Gumbel-Barnett Copula Density

### The Genest-Ghoudi Copula

Finally, we have the Genest-Ghoudi copula which is defined by
\begin{align*} \phi_\theta(u) &= \left(1-u^\tfrac{1}{\theta}\right)^\theta\\ \phi^{-1}_\theta(t) &= \left(1-t^\tfrac{1}{\theta}\right)^\theta \end{align*}
and is valid in $$n$$ dimensions if $$\theta$$ is an integer greater than one or if it's greater than or equal to $$n-1$$ with the limiting cases
\begin{align*} C_1 &= W\\ C_\infty &= M \end{align*}
Program 9 plots its density using the ak.genestCopulaDensity function which is implemented in GenestCopula.js along with the ak.genestCopula and ak.genestCopulaRnd functions.

Program 9: The Genest-Ghoudi Copula Density

That's enough about copulas for the time being, although there's more to be said and we shall return to the subject in the future.

### References

[1] Archimedean Skew, www.thusspakeak.com, 2018.

[2] Archimedean View, www.thusspakeak.com, 2018.

[3] Archimedean Review, www.thusspakeak.com, 2019.

[4] Nelsen, R. An Introduction to Copulas, Springer, 1998.

[5] You're Going To Have To Think! Why Automatic Differentiation Won't Cure Your Calculus Blues, www.thusspakeak.com, 2014.