In mlysy/LocalCop: Local Likelihood Inference for Conditional Copula Models

knitr::read_chunk("calculations.R")
library(LocalCop)
library(VineCopula)  # for simulating copula data
library(readr)       # for reading in the copula table
options(
  kableExtra.latex.load_packages = FALSE
)                    # for styling the copula table
library(kableExtra)  # add packages manually for joss submission
library(tidyr)       # for data manipulations
library(dplyr)       #
library(ggplot2)     # for plots
library(ggpubr)      # 
library(gamCopula)   # other packages for fitting
library(CondCopulas) # conditional copula models

\newcommand{\R}{{\textsf{R}}} \newcommand{\cpp}{{\textsf{C++}}}

Summary

Conditional copulas models allow the dependence structure between multiple response variables to be modelled as a function of covariates. LocalCop [@localcop] is an \R/\cpp package for computationally efficient semiparametric conditional copula modelling using a local likelihood inference framework developed in @ACY2011, @ACY2013 and @ACL2019.

Statement of Need

There are well-developed \R packages such as copula [@copula1; @copula2; @copula3; @copula4] and VineCopula [@vinecopula] for fitting copulas in various multivariate data settings. However, these software focus exclusively on unconditional dependence modelling and do not accommodate covariate information.

Aside from LocalCop, \R packages for fitting conditional copulas are gamCopula [@gamcopula] and CondCopulas [@condcopulas]. gamCopula estimates the covariate-dependent copula parameter using spline smoothing. While this typically has lower variance than the local likelihood estimate provided by LocalCop, it also tends to have lower accuracy [@ACL2019]. CondCopulas estimates the copula parameter using a semi-parametric maximum-likelihood method based on a kernel-weighted conditional concordance metric. LocalCop also uses kernel weighting, but it uses the full likelihood information of a given copula family rather than just that contained in the concordance metric, and is therefore more statistically efficient.

Local likelihood methods typically involve solving a large number of low-dimensional optimization problems and thus can be computationally intensive. To address this issue, LocalCop implements the local likelihood function in \cpp, using the \R/\cpp package TMB [@kristensen.etal16] to efficiently obtain the associated score function using automatic differentiation. Thus, LocalCop is able to solve each optimization problem very quickly using gradient-based algorithms. It also provides a means of easily parallelizing the optimization across multiple cores, rendering LocalCop competitive in terms of speed with other available software for conditional copula estimation.

Background

For any bivariate response vector $(Y_1, Y_2)$, the conditional joint distribution given a covariate $X$ is given by \begin{equation} F_X(y_1, y_2 \mid x) = C_X (F_{1\mid X} (y_1 \mid x),F_{2\mid X} (y_2 \mid x) \mid x ), \label{eq:fullmodel} \end{equation} where $F_{1\mid X}(y_1 \mid x)$ and $F_{2\mid X}(y_2 \mid x)$ are the conditional marginal distributions of $Y_1$ and $Y_2$ given $X$, and $C_X(u, v \mid x)$ is a conditional copula function. That is, for given $X = x$, the function $C_X(u, v \mid x)$ is a bivariate CDF with uniform margins.

The focus of LocalCop is on estimating the conditional copula function, which is modelled semi-parametrically as \begin{equation} C_X(u, v \mid x) = \mathcal{C}(u, v\mid \theta(x), \nu), \label{eq:copmod} \end{equation} where $\mathcal{C}(u, v \mid \theta, \nu)$ is a parametric copula family, the copula dependence parameter $\theta \in \Theta$ is an arbitrary function of $X$, and $\nu \in \Upsilon$ is an additional copula parameter present in some models. Since most parametric copula families have a restricted range $\Theta \subsetneq \mathbb{R}$, we describe the data generating model (DGM) in terms of the calibration function $\eta(x)$, such that \begin{equation} \theta(x) = g^{-1}(\eta(x)), \end{equation} where $g^{-1}: \mathbb{R} \to \Theta$ an inverse-link function which ensures that the copula parameter has the correct range. The choice of $g^{-1}(\eta)$ is not unique and depends on the copula family.

Local likelihood estimation of the conditional copula parameter $\theta(x)$ uses Taylor expansions to approximate the calibration function $\eta(x)$ at an observed covariate value $X = x$ near a fixed point $X = x_0$, i.e., $$ \eta(x)\approx \eta(x_0) + \eta^{(1)}(x_0) (x - x_0) + \ldots + \frac{\eta^{(p)}(x_0)}{p!} (x - x_0)^{p}. $$ One then estimates $\beta_k = \eta^{(k)}(x_0)/k!$ for $k = 0,\ldots,p$ using a kernel-weighted local likelihood function \begin{equation} \ell(\boldsymbol{\beta}) = \sum_{i=1}^n \log\left{ c\left(u_i, v_i \mid g^{-1}( \boldsymbol{x}_{i}^T \boldsymbol{\beta}), \nu \right)\right} K_h\left(\frac{x_i-x_0}{h}\right), \label{eq:locallik} \end{equation} where $(u_i, v_i, x_i)$ is the data for observation $i$, $\boldsymbol{x}_i = (1, x_i - x_0, (x_i - x_0)^2, \ldots, (x_i - x_0)^p)$, $\boldsymbol{\beta}= (\beta_0, \beta_1, \ldots, \beta_p)$, and $K_h(z)$ is a kernel function with bandwidth parameter $h > 0$. Having maximized $\ell(\boldsymbol{\beta})$ in \autoref{eq:locallik}, one estimates $\eta(x_0)$ by $\hat \eta(x_0) = \hat \beta_0$. Usually, a linear fit with $p=1$ suffices to obtain a good estimate in practice.

Usage

LocalCop is available on CRAN and GitHub. The two main package functions are:

CondiCopLocFit(): For estimating the calibration function at a sequence of values $\boldsymbol{x}0 = (x{01}, \ldots, x_{0m})$.
CondiCopSelect(): For selecting a copula family and bandwidth parameter using leave-one-out cross-validation (LOO-CV) with subsampling as described in @ACL2019.

In the following example, we illustrate the model selection/tuning and fitting steps for data generated from a Clayton copula with conditional Kendall $\tau$ displayed in \autoref{fig:copcomp}. The CV metric for each combination of family and bandwidth are displayed in \autoref{fig:select1-plot}.

<<dgm>>

<<select-precalc>>

<<select-calc>>

saveRDS(cvselect, file = "cvselect1.rds")

cvselect <- readRDS("cvselect1.rds")

<<select-plot>>

<<locfit>>

<<gamfit>>

<<condfit>>

sumfit1 <- tibble(
  x = x0,
  True = BiCopEta2Tau(family = family, eta = eta_fun(x0)),
  LocalCop = tau_loc,
  gamCopula = tau_gam,
  CondCopulas = tau_cond
)

<<dgm>>

<<select-precalc>>

<<select-calc>>

saveRDS(cvselect, file = "cvselect2.rds")

cvselect <- readRDS("cvselect2.rds")

<<locfit>>

<<gamfit>>

<<condfit>>

sumfit2 <- tibble(
  x = x0,
  True = BiCopEta2Tau(family = family, eta = eta_fun(x0)),
  LocalCop = tau_loc,
  gamCopula = tau_gam,
  CondCopulas = tau_cond
)

<<plotfit>>

plot1 <- plotfit(sumfit1) + ggtitle("(a) n = 300")
plot2 <- plotfit(sumfit2) + ggtitle("(b) n = 1000")

ggarrange(plot1, plot2, nrow = 2, common.legend = TRUE)

In \autoref{fig:copcomp}, we compare the true conditional Kendall $\tau$ to estimates using each of the three conditional copula fitting packages LocalCop, gamCopula, and CondCopulas, for sample sizes $n = 300$ and $n = 1000$. In gamCopula, selection of the copula family smoothing splines is done using the generalized CV framework provided by the \R package mgcv [@wood17]. In CondCopulas, selection of the bandwidth parameter is done using LOO-CV. In this particular example, the sample size of $n = 300$ is not large enough for gamCopula to pick a sufficiently flexible spline basis, and CondCopulas picks a large bandwidth which oversmooths the data. For the larger sample size $n = 1000$, the three methods exhibit similar accuracy.

Acknowledgements

We acknowledge funding support from the Natural Sciences and Engineering Research Council of Canada Discovery Grants RGPIN-2020-06753 (Acar) and RGPIN-2020-04364 (Lysy).