library(igcop)

The generic copula class (DJ and interpolated DJ) is characterized by a generating function \(\psi:[0, \infty) \rightarrow (0, 1]\), which needs to be a concave distribution function or convex survival function.

This manuscript provides results for a generic copula generated by some \(\psi\), in addition to the IG and IGL classes. Note that \(\psi\) need not have an inverse – the results still hold for the left-inverse, \(\psi^{\leftarrow}\).

Preliminary Functions

Incomplete Gamma Function

The (Upper) Incomplete Gamma function is defined as \[\Gamma^{*}(\alpha, x) = \int_x^{\infty} t ^ {\alpha - 1} e ^ {-t} dt\] for \(\alpha > 0\) and \(x \ge 0\), and the Gamma function is defined as \(\Gamma(\alpha) = \Gamma^{*}(\alpha, 0)\).

In R, the Gamma function is accessible through the function gamma(), and the igcop package makes the Incomplete Gamma function accessible through the igamma() function.

To make the igamma() function, the igcop package uses the Gamma distribution function pgamma() with scale or rate equal to 1, and shape equal to \(\alpha\) as an indirect way of accessing the Incomplete Gamma function. Denoting this Gamma cdf as \(F_{\alpha}\), the relationship is \[F_{\alpha}(x) = 1 - \frac{\Gamma^{*}(\alpha, x)}{\Gamma(\alpha)}.\]

Generating Function

The IG and IGL copula families result from a parametric generating function \(\Psi_k\) for \(k > 1\), defined below.

Function Formula igcop function
\(\Psi_k(y)\) \[ \begin{cases} y \frac{\Gamma(k) - \Gamma^{*}(k, y ^ {-1})}{\Gamma(k - 1)} + \frac{\Gamma^{*}(k - 1,y ^ {-1})}{\Gamma(k - 1)}, & y > 0;\\ 0, & y = 0 \end{cases} \] igl_gen()
\(\Psi_k'(y)\) \[ \frac{\Gamma(k) - \Gamma^{*}(k, 1 / y)} {\Gamma(k - 1)} \] igl_gen_D()
\(\Psi_k''(y)\) \[ -\frac{y ^ {-k - 1} \exp(-1 / y)} {\Gamma(k - 1)} \] igl_gen_DD()

The “kappa” version of the generating function is also important.

General definition.

Function Formula
\(\kappa_{\psi}(y)\) \(\psi(y) - y \psi'(y)\)
\(\kappa_{\psi}'(y)\) \(-y \psi''(y)\)

Taking our specific generating function in place of \(\psi\):

Function Formula igcop function
\(\kappa_k(y) := \kappa_{\Psi_k}(y)\) \[ \frac{\Gamma^{*}(k - 1, 1 / y)} {\Gamma(k - 1)} = 1 - \text{pgamma}(1 / y, k - 1) \] igl_kappa()
\(\kappa_k'(y)\) \[ \frac{\exp(-1/y)}{y ^ k \Gamma(k - 1)} = \frac{\text{dgamma}(1 / y, k - 1)}{y ^ 2} \] NA

Interpolating Function

The copula arising from the generating function \(\psi\) gives rise to a parametric copula family governed by a parameter \(\theta > 0\), where \(\theta = 0\) is the independence copula, and \(\theta \rightarrow \infty\) results in the original copula. As such, this copula family interpolates a DJ copula, giving rise to an even wider class of copulas.

The interpolated class of copulas relies on an interpolating function function \(H_\psi(\cdot; \theta): [1,\infty) \rightarrow (0, 1]\), defined below.

Generic definition of the interpolating function:

Function Formula
\(H_{\psi}(y; \eta)\) \[ \begin{cases} \frac{1}{y}\psi\left(\frac{1}{\eta \log y}\right), & y > 1;\\ 1, & y = 1 \end{cases} \]
\(\text{D}_1 H_{\psi}(t; \eta)\) \[ -\frac{1}{t ^ 2} \left[ \psi\left(\frac{1}{\eta \log t}\right) + \frac{1}{\eta (\log t) ^ 2} \psi'\left(\frac{1}{\eta \log t}\right) \right] \]
\(\text{D}_2 H_k(t; \eta)\) \[ -\frac{1}{t \eta ^ 2 \log t} \psi'\left(\frac{1}{\eta \log t}\right) \]

Definition of the interpolating function using the IGL generating function \(\Psi_k\) (taking \(H_k := H_{\Psi_k}\)):

Function Formula igcop function
\(H_k(t; \eta)\) \[ \begin{cases} \frac{1}{t}\Psi_k\left(\frac{1}{\eta \log t}\right), & t > 1;\\ 1, & t = 1 \end{cases} \] interp_gen()
\(\text{D}_1 H_k(t; \eta)\) \[ -\frac{ (\log t + 1) (\Gamma(k) - \Gamma^{*}(k, \eta \log t)) + \eta (\log t) ^ 2 \Gamma^{*}(k - 1, \eta \log t) }{ \eta (t \log t) ^ 2 \Gamma(k - 1) } \] interp_gen_D1()
\(\text{D}_2 H_k(t; \eta)\) \[ -\frac{1}{t \eta ^ 2 \log t} \frac{\Gamma(k) - \Gamma^{*}(k,\eta \log t)}{\Gamma(k)} \] NA

For the IG copula, \(H_k := H_{\Psi_k}\) (function interp_gen() in the igcop package). Note that \(H_k(\cdot; \theta)\) is strictly decreasing for all \(\theta > 0\) and \(k > 1\), so \(H_k^{-1}(\cdot; \theta)\) is the unique inverse function of \(H_k(\cdot; \theta)\).

OLD PLOTS (OF SOMETHING LIKE H_KAPPA)

Copula Formulas

DJ Copula Class

The Durante-Jaworski copula class:

Quantity Formula
\(C(u, v; \psi)\) \[ u + v - 1 + (1 - u) \psi\left((1 - u) ^ {-1} \psi^{\leftarrow}(1 - v)\right) \]
\(C_{2 | 1}(v | u; \psi)\) \[ 1 - \kappa_{\psi}\left((1 - u) ^ {-1} \psi ^ {\leftarrow}(1 - v) \right) \]
\(C_{1 | 2}(u | v; \psi)\) \[ 1 - \frac{\psi'\left((1 - u) ^ {-1} \psi ^ {\leftarrow}(1 - v)\right)} {\psi'\left(\psi ^ {\leftarrow}(1 - v)\right)} \]
\(c(u, v; \psi)\) \[ -\frac{ \psi ^ {\leftarrow}(1 - v) \psi''\left((1 - u)^{-1} \psi^{\leftarrow}(1 - v)\right) }{ (1 - u) ^ 2 \psi'\left(\psi^{\leftarrow}(1 - v)\right) } \]
\(C_{2 | 1}^{-1}(\tau | u; \psi)\) \[ 1 - \psi((1 - u) \kappa_{\psi}^{\leftarrow}(1 - \tau)) \]
\(C_{1 | 2}^{-1}(\tau | v; \psi)\)

Note that Coia (2017) in E.1.2 writes these distributional quantities as if the family was defined as reflections.

DJ Copula Class with \(\Psi_k\): IGL Copula Class

The parameter space is \(k > 1\). Note: Values of \(k < 2\) lead to unstable computations – especially with \(k\) close to 1.

Quantity Formula distplyr function
\(C(u, v; k)\) \[ u + v - 1 + (1 - u) \Psi_k\left((1 - u) ^ {-1} \Psi_k^{-1}(1 - v)\right) \] piglcop()
\(C_{2 | 1}(v | u; k)\) \[ 1- \frac{\Gamma^{*}\left(k - 1, (1 - u) / \Psi_k ^ {-1}(1 - v)\right)} {\Gamma\left(k-1\right)} \] pcondiglcop(), pcondiglcop21()
\(C_{1 | 2}(u | v; k)\) \[ 1 - \frac{\Gamma(k) - \Gamma^{*}\left(k, (1 - u) / \Psi_k ^ {-1}(1 - v)\right)} {\Gamma(k) - \Gamma^{*}\left(k, 1 / \Psi_k ^ {-1}(1 - v)\right)} \] pcondiglcop12()
\(c(u, v; k)\) \[\frac{ (1 - u)^{k-1} \left[\Psi_k^{-1}(1 - v)\right] ^ {-k} \exp\left\{-(1 - u) / \Psi_k^{-1}(1 - v)\right\} }{ \Gamma(k) - \Gamma^{*}\left(k, 1 / \Psi_k^{-1}(1 - v)\right) }\] diglcop()
\(C_{2 | 1}^{-1}(\tau | u; k)\) \[ 1 - \Psi_k\left(\frac{1 - u}{\text{qgamma}(\tau, k - 1)}\right) \] qcondiglcop(), qcondiglcop21()
\(C_{1 | 2}^{-1}(\tau | v; k)\) qcondiglcop12()

Interpolated DJ Copula Class

Quantity Formula
\(C(u, v; \theta, \psi)\) \[ u + v - 1 + (1 - u) H_{\psi}\left(H_{\psi} ^ {\leftarrow}(1 - v; \theta); \theta (1 - u)\right) \]
\(C_{2 | 1}(v | u; \theta, \psi)\) \[ 1 - H_{\kappa_{\psi}}\left(H_{\psi} ^ {\leftarrow}(1 - v; \theta); \theta (1 - u)\right) \]
\(C_{1 | 2}(u | v; \theta, \psi)\) \[ 1 - (1 - u) \frac{ \psi\left(\frac{1}{\theta (1 - u) \log t}\right) + \frac{1}{\theta (1 - u) (\log t) ^ 2} \psi'\left(\frac{1}{\theta (1 - u) \log t}\right) } { \psi\left(\frac{1}{\theta \log t}\right) + \frac{1}{\theta (\log t) ^ 2} \psi'\left(\frac{1}{\theta \log t}\right) } \] where \(t = H_{\psi} ^ {\leftarrow}(1 - v; \theta)\)
\(c(u, v; \theta, \psi)\) \[ \frac{ \kappa_{\psi}\left(\frac{1}{\theta (1 - u) \log t}\right) + \frac{1}{\theta (1 - u) (\log t) ^ 2} \kappa_{\psi}'\left(\frac{1}{\theta (1 - u) \log t}\right) } { \psi\left(\frac{1}{\theta \log t}\right) + \frac{1}{\theta (\log t) ^ 2} \psi'\left(\frac{1}{\theta \log t}\right) } \] where \(t = H_{\psi} ^ {\leftarrow}(1 - v; \theta)\)
\(C_{2 | 1}^{-1}(\tau | u; \theta, \psi)\) \[ 1 - H_{\psi}\left(H_{\kappa_{\psi}}^{\leftarrow}(1 - \tau; \theta (1 - u)); \theta \right) \]
\(C_{1 | 2}^{-1}(\tau | v; \theta, \psi)\)

Note that Coia (2017) in E.2.2 writes these distributional quantities in terms of the reflection copula.

Interpolated DJ Copula Class with \(\Psi_k\): IG copula

The parameter space is \(\theta \ge 0\), \(k > 1\). Note: Values of \(k < 2\) lead to unstable computations – especially with \(k\) close to 1.

Quantity Formula distplyr function
\(C(u, v; \theta, k)\) \[ u + v - 1 + (1 - u) H_k\left(H_k ^ {-1}(1 - v; \theta); \theta (1 - u)\right) \] pigcop()
\(C_{2 | 1}(v | u; \theta, k)\) \[ 1 - \frac{ \Gamma^{*}\left(k - 1, \theta (1 - u) \log\left(H_k^{-1}(1 - v; \theta)\right)\right) }{ H_k^{-1}(1 - v; \theta) \Gamma(k - 1) } \] pcondigcop(), pcondigcop21()
\(C_{1 | 2}(u | v; \theta, k)\) \[ 1 - (1 - u) \frac{ \text{D}_1 H_k\left(H_k^{-1}(1 - v; \theta); \theta(1 - u)\right) }{ \text{D}_1 H_k\left(H_k^{-1}(1 - v; \theta); \theta\right) } \] pcondigcop12()
\(c(u, v; \theta, k)\) \[ -\frac{ \theta ^ {k - 1} (1 - u) ^ {k-1} (\log t_v) ^ {k - 2} t_v ^ {-\theta (1 - u)} + \Gamma^{*}(k - 1, \theta (1 - u) \log t_v) }{ \Gamma(k - 1) t_v ^ 2 \text{D}_1 H_k(t_v; \theta) } \] digcop()
\(C_{2 | 1}^{-1}(\tau | u; \theta, k)\) qcondigcop(), qcondigcop21()
\(C_{1 | 2}^{-1}(\tau | v; \theta, k)\) qcondigcop12()

(Note that Coia (2017) mistakenly writes \(H_{\psi}\) instead of \(H_k\) in the formula for \(C_{\text{IG}, 1 | 2}\).)

The copula family is named “Integrated Gamma” because \(\Psi_k'(1 / \cdot)\) is proportional to the Gamma distribution cdf with scale parameter 1 and shape parameter \(k\).

These two families are related, in that \[\lim_{\theta \rightarrow \infty} C_{IG}(u, v; \theta, k) = C_{IGL}(u, v; k)\] for all \(k\). Further, taking \(\theta = 0\) results in the independence copula, so that the IG copula family is an “interpolation” between the independence copula and each IGL copula.

Translating formulas to code

Generating Function \(\Psi_k\)

Noting that \(\Gamma(k - 1) = \Gamma(k) / (k - 1)\),

\[ \Psi_k(t) = y (k - 1) \left(1 - \frac{\Gamma^{*}(k, y ^ {-1})}{\Gamma(k)} \right) + \left(\frac{\Gamma^{*}(k - 1,y ^ {-1})}{\Gamma(k - 1)}\right) \]

The two major parentheses can be calculated from the Gamma cdf: the first set being pgamma(1 / y, shape = k), and the second set being pgamma(1 / y, shape = k - 1, lower.tail = FALSE). Here is the definition of igl_gen:

igl_gen
## function (t, k) 
## {
##     tinv <- 1/t
##     (k - 1) * t * stats::pgamma(tinv, k) + stats::pgamma(tinv, 
##         k - 1, lower.tail = FALSE)
## }
## <bytecode: 0x7f8cb47903f0>
## <environment: namespace:igcop>

The derivative is proportional to the Gamma density. Using the same Gamma function identity:

\[ \Psi_k'(y) = (k - 1) \left(1 - \frac{\Gamma^{*}(k, y ^ {-1})}{\Gamma(k)}\right) \]

Computing inverses

To compute by solving \(\varphi_k(x; \eta) = p\) for \(x\), the Newton-Raphson algorithm is used to solve \(g(x) = 0\) for \(x\), where \[ g\left(x\right)=xp-\bar{F}_{k-1}\left(\eta\log x\right), \] with derivative \[ g'\left(x\right)=p+\frac{\eta}{x}f_{k-1}\left(\eta\log x\right). \] A starting value is used by noting that \(\varphi_k\) is the product of two invertible survival functions, and is therefore smaller than those two survival functions. The smaller of the two roots of these survival functions is therefore an upper bound for the root of \(\varphi_k\). The starting point is taken to be immediately to the left of this upper bound.

Computing interp_gen_inv

To compute \(H_k^{\leftarrow}\left(p,\theta\right)\) for \(p\in\left(0,1\right)\) and some \(\theta>0\), Newton-Raphson can be used to solve \(h\left(x\right)=0\) for \(x\), where \[ h\left(x\right)=xp-\Psi_k\left(\frac{1}{\theta\log x}\right). \] This function has derivative \[ h'\left(x\right)=p+\frac{1}{\theta x\left(\log x\right)^2}\Psi_k'\left(\frac{1}{\theta\log x}\right). \]

To get a good starting value for the algorithm, note that the function \(H_k\) is a product of two survival functions, just like \(\varphi_k\) is. And like inverting \(\varphi_k\), the inverse \(H_k^{\leftarrow}\left(p;\theta\right)\) for some \(\theta>0\) has an upper bound that’s the minimum of the roots of the individual survival functions: \[ \min\left\{ \frac{1}{p},\exp\left(\frac{1}{\theta\Psi_k^{\leftarrow}\left(p\right)}\right)\right\} , \] where \(\Psi_k\) is defined in Equation~(@ref(eq:cnstr-Psi)). Computing the inverse \(\Psi_k^{\leftarrow}\) should not be problematic in terms of accuracy of the final result, since this computation is only used to determine a starting value.

References

Coia, Vincenzo. 2017. “Forecasting of Nonlinear Extreme Quantiles using Copula Models.” PhD Dissertation, The University of British Columbia.