Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
103,476 characters · 63 sections · 205 citation commands
Implicit Copulas: An Overview
Copulas are widely used to specify multivariate distributions for the statistical modeling of data. Fields where copula models have had a significant impact include (but are not limited to) actuarial science frees1998, finance cherubini2004,McNFreEmb2005,patton2006, hydrology favre2004,genest2007meta, climatology schoelzel2008, transportation bhat2009,smithbike2011 and marketing danaher2011,park2012. Copula models are popular because they simplify the specification of a distribution, allowing the marginals to be modeled arbitrarily, and then combined using a copula function. In practice, a major challenge is the selection and estimation of a copula function that captures the dependence structure well and is tractable. One choice are “implicit copulas”, which are copulas constructed from existing multivariate distributions by the inversion of Sklar's theorem as in nelsen06. This is a large and flexible family of copulas, which share an auxiliary representation that makes estimation tractable in high dimensions. Thus, they are suitable for modeling the large datasets that arise in many modern applications. The objective of this paper is to introduce and survey implicit copulas and their use in statistical modeling in an accessible manner.
Implicit copulas have a long history with key developments spread across multiple fields, including actuarial studies, econometrics, operations research, probability and statistics. Yet while there are many excellent existing monographs and surveys on copulas and copula models (see genest1986,joe97,McNFreEmb2005,nelsen06,genest2007,jaworski2010,patton2012,nikoloulopoulos2013,joe2014dependence and durante2015 for prominent examples) there does not appear to be a dedicated survey or overview on this important class of copulas. This paper aims to fill this gap and provides an overview that stresses common features of the implicit copula family, likelihood-based estimation, and the usefulness of implicit copulas in statistical modeling. Particular focus is given to recent developments on implicit copula processes for regression and time series data, along with Bayesian inference that extends the earlier overview by smithbcop2013 to these copula processes.
Two econometric applications illustrate the use of implicit copula models with non-Gaussian data. The first is a time-varying heteroscedastic time series model for U.S. inflation between 1954:Q1 and 2020:Q2. The implicit copula is a copula process constructed from a nonlinear state space model as in smithman2018. It is a “time series copula” that captures serial dependence. The second application is a five factor asset pricing regression model fama2015five with an asymmetric Laplace marginal distribution for monthly equity returns. The implicit copula here is a “regression copula” process with respect to the covariates as in KleSmi2019. The copula model forms a distributional regression KleKneLan2015,kneib2021, where the five factors affect the entire distribution of equity returns, not just its first or other moments. In both applications the implicit copulas are of dimension equal to the number of observations, so that they are high-dimensional. Nevertheless, their auxiliary representation allows for likelihood-based estimation of the copula parameters. In both examples the marginal distribution of the response variables exhibit strong asymmetries.
The overview is organized as follows. Section (ref) introduces general copula models, and then implicit copulas specifically. Their interpretation as transformations and specifications for variables that are continuous, discrete or mixed are also discussed. Section (ref) covers elliptical and skew-elliptical copulas, including the Gaussian, $t$, skew $t$ and factor copulas. Implicit copulas that capture serial dependence in time series data are covered in Section (ref). Section (ref) extends these to implicit copulas that capture both serial and cross-sectional dependence in multivariate time series. Section (ref) covers regression copula processes, with the implicit copula constructed from a regularized linear regression given in detail. It is shown that when this copula is combined with flexible marginals, it defines a promising new distributional regression model. Last, Section (ref) discusses the advantages of using implicit copula models for modeling data, and future directions.
All copula models are based on the theorem of sklar59 (i.e. “Sklar's theorem”), which states that for every random vector $\bm{Y}=(Y_1,\ldots,Y_m)^\top$ with distribution function $F_Y$ and marginals $F_{Y_1},\ldots,F_{Y_m}$, there exists a “copula function” $C:[0,1]^m \rightarrow [0,1]$, such that
where $\bm{y}=(y_1,\ldots,y_m)^\top$. The copula function $C$ is a well-defined distribution function for a random vector $\bm{U}=(U_1,\ldots,U_m)^\top$ on the unit cube with uniform marginal distributions. To construct a copula model, select $F_{Y_1},\ldots,F_{Y_m}$ (i.e. the “marginal models”) and a copula function $C$, to define $F_Y$ via (ref).
If all the elements of $\bm{Y}$ are continuous, then differentiating through (ref) gives the density
where $f_{Y_j}=\frac{\partial}{\partial y_j} F_{Y_j}$, and $c(\bm{u})=\frac{\partial^m}{\partial u_1\cdots \partial u_m} C(\bm{u})$ is widely called the “copula density” with $\bm{u}=(u_1,\ldots,u_m)^\top$. (Throughout this paper the notation $c(\text{\boldmath$u$})$ and $c(u_1,u_2,\ldots,u_m)$ are used interchangeably, as are $C(\text{\boldmath$u$})$ and $C(u_1,u_2,\ldots,u_m)$.) The decomposition at (ref) is used to specify the likelihood of a continuous response vector $\bm{Y}$ in a statistical model.
If all the elements of $\bm{Y}$ are discrete-valued (e.g. as with ordinal or binary data) the probability mass function is obtained by differencing over the elements of $\bm{Y}$ as follows. Let $b_j=F_{Y_j}(y_j)$ and $a_j=F_{Y_j}(y_j^-)$ be the left-hand limit of $F_j$ at $y_j$ (which is $a_j=F_{Y_j}(y_j-1)$ for ordinal $Y_j$). Then the mass function is
where $\text{\boldmath$v$}=(v_1,\ldots,v_m)^\top$ is a differencing vector, and the notation
Evaluating the mass function at (ref) is an $O(2^m)$ computation, so that its direct evaluation is impractical for high values of $m$ when undertaking likelihood-based estimation nikoloulopoulos2013. One solution suggested by smithkhaled2012 is to consider the joint distribution of $(\bm{Y},\bm{U})$. To do so, note that when $Y_j$ is discrete, $F_{Y_j}$ is a many-to-one function and $Y_j|U_j$ is a degenerate distribution with density $f(y_j|u_j)=\mathds{1}(a_j \leq u_j < b_j)$, where the indicator function $\mathds{1}(X)=1$ if $X$ is true, and zero otherwise. (An alternative notation is to use the Dirac delta function, with $f(y_j|u_j)=\delta_{y_j}(F_{Y_j}^-(u_j))$ where $F_{Y_j}^-$ is the quantile function of $Y_j$.) Then the mixed density of $(\bm{Y},\bm{U})$ is
Marginalizing over $\bm{U}$ gives the probability mass function at (ref) (i.e. $f_Y(\text{\boldmath$y$})=\int f_{Y,U}(\text{\boldmath$y$},\text{\boldmath$u$})\mbox{d}\text{\boldmath$u$}$); see Proposition 1 in smithkhaled2012.
Equation (ref) can be used to define an “extended likelihood” for estimation using computational methods for latent variables, where the observations on $\bm{U}$ are the latents. This has two advantages. First, the $O(2^m)$ computation at (ref) is avoided, allowing estimation for higher values of $m$. Second, only the copula density $c$ is required and not the copula function $C$, which is an advantage for some copulas where only $c$ can be computed, as is the case with most vine copulas joe1996,AasCzaFriBak2009. Bayesian data augmentation can be used based on (ref), and evaluated using Markov chain Monte Carlo (MCMC) as in smithkhaled2012 or variational Bayes methods as in loaiza2019VB. The latter is particularly attractive, because it allows for the estimation of discrete-margined copulas of very high dimensions, with examples up to $m=792$ presented by these authors.
If some elements of $\bm{Y}$ are continuous and others discrete, then $f_Y$ is often called a “mixed density”. In this case, an extended likelihood can be constructed from the distribution of $\bm{Y}$ joint with the elements of $\bm{U}$ that correspond only to the discrete variables; see smithkhaled2012. Similarly, if some individual elements $Y_j$ have distributions that are mixtures of continuous and discrete distributions (such as a zero-inflated continuous distribution) then an extended likelihood can also be constructed for this case; see gunawan2020 for how to do so.
McNFreEmb2005 use the term “implicit copula” for the copula that is implicit in the multivariate distribution of a continuous random vector $\bm{Z}=(Z_1,\ldots,Z_m)^\top$. It is obtained by inverting Sklar's theorem, which nelsen06 calls the “inversion method”, so that copulas derived in this fashion are also called “inversion copulas” (e.g. smithman2018). If $\bm{Z}$ has distribution function $F_Z$ with marginals $F_{Z_1},\ldots, F_{Z_m}$, then its implicit copula function is
Differentiating with respect to $\text{\boldmath$u$}$ gives the implicit copula density
where $\text{\boldmath$z$}=(z_1,\ldots,z_m)^\top$ is a function of $\text{\boldmath$u$}$ with elements $z_j=F_{Z_j}^{-1}(u_j)$ for $j=1,\ldots,m$. The implicit copula function $C_Z$ and density $c_Z$ above can be employed in (ref), (ref) and (ref). Thus, an implicit copula model uses Sklar's theorem twice: once to form the joint distribution $F_Y$ with arbitrary marginals, and a second time to construct the implicit copula from the joint distribution $F_Z$.
Because implicit copulas are an immediate consequence of Sklar's theorem, they have a long history. Early uses for modelling data include ruschendorf1976 and deheuvels1979, who both construct a non-parametric implicit copula from the empirical distribution function (although neither called it a copula). ruschendorf2009 gives an overview of the early developments of implicit copulas, pointing out that many transformation-based multivariate models---which themselves have a long history---are also copula models based on implicit copulas (although in the early literature this was often unrecognized and the term “copula” not used).
Note that only a continuous distribution $F_Z$ is used to construct an implicit copula here. This is because the implicit copula of a discrete distribution $F_Z$ is not unique genest2007.
One way to look at all copula models is that they are a transformation from $\bm{Y}$ to $\bm{U}=(U_1,\ldots,U_m)^\top \in[0,1]^m$. The key observation is that it is usually easier to capture multivariate dependence using $C$ on the vector space $[0,1]^m$, rather than directly on the domain of the original vector $\bm{Y}$. Implicit copulas go one step further, with a second transformation from $\bm{U}$ to $\bm{Z}=(F_{Z_1}^{-1}(U_1),\ldots,F_{Z_m}^{-1}(U_m))^\top$, and then capture the dependence structure using the distribution $F_Z$. Table (ref) provides a summary of these transformations, along with the marginal and joint distribution and density/mass functions of $\bm{Y}$. Throughout this paper, the vector $\bm{U}$ is referred to as the “copula vector” and $\bm{Z}$ as the “auxiliary vector” (the latter is also called a “pseudo vector” in smith+k19). Simulation from an implicit copula model is straightforward if $F_Z$ is tractable using Algorithm (ref), which produces a draw $\text{\boldmath$y$}\sim F_Y$.
Notice that the transformation $U_j=F_{Z_j}(Z_j)\sim \mbox{Uniform}[0,1]$ removes all features of the marginal distribution of $Z_j$. This becomes an important observation for establishing parameter identification when constructing implicit copulas, as discussed in Sections (ref), (ref) and (ref).
For the case where the elements of $\bm{Y}$ are discrete-valued, for an implicit copula model there exists an alternative extended likelihood based on the joint density of $(\bm{Y},\bm{Z})$, rather than that of $(\bm{Y},\bm{U})$ given previously at (ref). This alternative joint density is
with $a_j,b_j$ as defined above in Section (ref). Marginalizing over $\bm{Z}$ produces the probability mass function at (ref); i.e. $f_Y(\text{\boldmath$y$})=\int f_{Y,Z}(\text{\boldmath$y$},\text{\boldmath$z$})\mbox{d}\text{\boldmath$z$}$. An advantage is that it is often simpler to use computational methods to estimate an implicit copula using (ref) rather than (ref). Moreover, an extended likelihood is also easily defined for vectors $\bm{Y}$ with combinations of discrete, continuous or even mixed valued elements, by simplifying (ref) to only include elements of $\bm{Z}$ that correspond to the non-continuous valued variables.
Bayesian data augmentation is a suitable method for estimation using this extended likelihood. Here, values for $\text{\boldmath$z$}$ are generated in an MCMC sampling scheme to evaluate an “augmented posterior” proportional to the extended likelihood multiplied by a parameter prior. This has been used to estimate the elliptical and skew elliptical copulas discussed in Section (ref) below. For example, pitt2006 do so for a Gaussian copula, while danaher2011 do so for the $t$ copula, and smith2012 for the skew $t$ copula. hoff2007 considered the extended likelihood above using empirical marginals and rank data, danaher2011 and dobra2011 provide early applications to higher dimensional Gaussian $F_Z$. Last, the multivariate probit model is a Gaussian copula model, and the popular approach of chib1998 is a special case of these data augmentation algorithms.
In practice, parametric copulas $C(\text{\boldmath$u$};\text{\boldmath$\theta$})$ with parameter vector $\text{\boldmath$\theta$}$ are almost always used in statistical modelling, with McNFreEmb2005, nelsen06 and joe2014dependence giving overviews of choices. However, the implicit copulas of elliptical distributions, and more recently skew elliptical distributions, are common choices for capturing dependence in many applications. An attractive feature is that because elliptical and skew elliptical distributions are closed under marginalization, so are their implicit copulas.
The simplest and most popular elliptical copula is the “Gaussian copula”, which is constructed from $\bm{Z}\sim N_m(\bm{0},\Omega)$ with $\Omega$ an $m\times m$ correlation matrix. If $\Phi_m(\cdot;\text{\boldmath$a$},\Omega)$ denotes an $N_m(\text{\boldmath$a$},\Omega)$ distribution function, and $\Phi(\cdot)$ a $N(0,1)$ distribution function, then from (ref) the Gaussian copula function is
If $\phi_m(\cdot;\text{\boldmath$a$},\Omega)$ is a $N_m(\text{\boldmath$a$},\Omega)$ density, and $\phi$ is a standard normal density, then plugging the Gaussian densities into (ref) gives the Gaussian copula density
with $\text{\boldmath$z$}=(\Phi^{-1}(u_1),\ldots,\Phi^{-1}(u_m))^\top$.
There are a number of immediate observations on the Gaussian copula. First, the auxiliary vector $\bm{Z}$ has a distribution with a zero mean and unit marginal variances. This is because information about the first two marginal moments of $Z_j$ are lost in the transformation $U_j=F_{Z_j}(Z_j)$ and are unidentified in the copula density. Second, adopting any constant mean value (other than zero) and marginal variances (other than unit values) for $\bm{Z}$ produces the same Gaussian copula $C_{\tiny{\mbox Ga}}$. Third, closure under marginalization means that if $\bm{U}$ has distribution function $C_{\tiny{\mbox Ga}}(\text{\boldmath$u$};\Omega)$, then any subset $\bm{U}_0$ of elements of $\bm{U}$ has distribution function $C_{\tiny{\mbox Ga}}(\text{\boldmath$u$}^0;\Omega^0)$, where $\Omega^0$ is a correlation matrix made up of the corresponding rows and columns of $\Omega$.
A fourth observation is that any parametric correlation structure for $\bm{Z}$ is inherited by the Gaussian copula. It is this property that has led the widespread adoption of Gaussian copula models for modeling time series cario1996, longitudinal lambert2002, cross-sectional murray2013 and spatial bai2014,hughes2015 data. The Gaussian copula has a long history, particularly when formed implicitly via transformation (e.g. LiHammond), although some early and influential mentions include joe1993, clemen1999 and wang1999, while li2000 popularized its use in finance. A comprehensive overview of the Gaussian copula and its properties is given by Song2000.
fang2002 and embrechts2002 use an elliptical distribution for $\bm{Z}$, and study the resulting class of “elliptical copulas”. When combined with choices for the marginals of $\bm{Y}$ in a copula model, fang2002 call the distribution $F_Y$ “meta-elliptical”, and an overview of their dependence properties is given by abdous2005. After the Gaussian copula, the most popular elliptical copula is the $t$ copula, where a multivariate $t$ distribution with degrees of freedom $\nu>0$ is adopted for $\bm{Z}$. embrechts2002 and venter2003 study this copula, and the main advantage is that it can capture higher dependence in extreme values, which is important for financial and actuarial variables. A lesser known property is that values of $\nu$ close to zero allow for positive dependence between squared elements of $\bm{Y}$. This is useful for capturing the serial dependence in heteroscedastic time series, such as equity returns in finance; see loaiza2018hetero and bladt2021.
Elliptical copulas exhibit radial symmetry, where the distributions of $(U_i,U_j)$ and $(1-U_i,1-U_j)$ are the same. Yet there are applications where this is unrealistic, including for the dependence between equity returns longin2001,ang2002asymmetric and regional electricity spot prices smith2012. The implicit copulas of skew elliptical distributions genton2004 allow for asymmetric pairwise dependence, with the most common being those constructed from the differing skew $t$ distributions. demarta2005 were the first to construct an implicit copula from a skew $t$ distribution (i.e. a “skew $t$ copula”), for which they used a special case of the generalized hyperbolic distribution, and chan2010 do so for an adjustment of the skew normal distribution of azzalini1996. The most popular variants of the skew $t$ distribution are those of Azz2003 and sahu2003, which share a similar conditionally Gaussian representation. smith2012 show how to construct implicit copulas from these latter two skew $t$ distributions, and estimate them using MCMC. yoshiba2018 considers maximum likelihood estimation for the skew $t$ copula constructed from the distribution of Azz2003, and oh2020dynamic consider a dynamic extension of the skew $t$ copula of demarta2005 for high dimensions.
Write $t_d(\text{\boldmath$a$},\Omega,\nu)$ for a $d$-dimensional $t$ distribution with location $\text{\boldmath$a$}$, scale matrix $\Omega$ and degrees of freedom $\nu$, with density $f_t(\cdot;\text{\boldmath$a$},\Omega,\nu)$. Let $\bm{X}$ and $\bm{Q}$ be $(m\times 1)$ vectors with joint distribution
Here, $D=\mbox{diag}(\delta_1,\ldots,\delta_m)$ is a diagonal matrix and $\Gamma$ is positive definite. Then the skew $t$ distribution of sahu2003 (with location parameter equal to zero) is given by $\bm{Z}=(\bm{X}|\bm{Q}>\bm{0})$, which has density
where $\bm{V}\sim t_m\left(D(\Gamma+D^2)^{-1}\text{\boldmath$z$},\frac{S(\text{\boldmath$z$})+\nu}{m+\nu}(I-D(\Gamma+D^2)^{-1}D,m+\nu\right)$ and $S(\text{\boldmath$z$})=\text{\boldmath$z$}'(\Gamma+D^2)^{-1}\text{\boldmath$z$}$. This manner of constructing a skew $t$ distribution is called “hidden conditioning” because $\bm{Q}$ is latent. The skew $t$ distribution of Azz2003 is constructed in a similar way, but where $\bm{Q}$ is a scalar.
The parameter $\text{\boldmath$\delta$}=(\delta_1,\ldots,\delta_m)^\top$ controls the level of asymmetry in the distribution of $\bm{Z}$, but in the implicit copula it controls the level of asymmetric dependence. This is a key observation as to why skew $t$ copulas have strong potential for applied modeling. To construct this copula, first fix the leading diagonal elements of $\Gamma$ to ones (i.e. restrict $\Gamma$ to be a correlation matrix), and note that the marginal of $Z_j$ is also a skew $t$ distribution with density $f_{\mbox{\tiny St}}(z_j;1,\delta_j,\nu)$. Then, the copula function and density are given by (ref) and (ref), respectively. These require computation of the distribution function $F_{Z_j}(z_j)=\int_{-\infty}^{z_j}f_{\mbox{\tiny St}}(z_j';1,\delta_j,\nu)\mbox{d}z_j'$ and its inverse (i.e. the quantile function) which can either be undertaken numerically using standard methods, or using the interpolation approach outlined in (ref) for large datasets. Simulation from a skew $t$ copula model is straightforward using (ref) and the representation of a $t$ distribution as Gaussian conditional on a Gamma variate. To do so, at Step 1 of Algorithm 1 generate a draw $\text{\boldmath$z$}\sim F_Z$ by drawing sequentially as follows:
A computational bottleneck for the evaluation of the skew $t$ copula density is the evaluation of the multivariate integral $\mbox{Pr}(\bm{V}>\bm{0};\text{\boldmath$z$})$ at (ref).
However, this can be avoided in likelihood-based estimation by considering the tractable conditionally Gaussian representation motivated by (ref). Let $W\sim \mbox{Gamma}(\nu/2, \nu/2)$, then consider the joint distribution of $(\bm{X},\bm{Q},W|\bm{Q}>\bm{0})$ with density
where $(\bm{X}|\bm{Q}=\text{\boldmath$q$},W=w)\sim N_m(D\text{\boldmath$q$},\frac{1}{w}\Gamma)$ and $(\bm{Q}|W=w)\sim N_m(\bm{0},\frac{1}{w}I_m)$. Marginalizing out $(\text{\boldmath$q$},w)$ gives the skew $t$ density at (ref) in $\text{\boldmath$x$}$. smith2012 use this feature to design Bayesian data augmentation algorithms for the skew $t$ copula that generate $(\text{\boldmath$q$},w)$ as latent variables in Markov chain Monte Carlo (MCMC) sampling schemes for both continuous-valued and discrete-valued $\bm{Y}$.
The density of the Azz2003 skew $t$ distribution does not feature the multivariate probability term $\mbox{Pr}(\bm{V}>0)$, so that it is easier to evaluate its implicit copula density, as in yoshiba2018. But when computing the Bayesian posterior using data augmentation it makes little difference, because the copula density is never evaluated directly.
To capture dependence in high dimensions, “factor copulas” are increasingly popular, and there are two main types in the literature. The first links a small number of independent factors by a pair-copula construction to produce a higher dimensional copula, as proposed by krupskii2013. Flexibility is obtained by using different bivariate copulas for the pair-copulas and a different number of factors, with applications and extensions found in nikoloulopoulos2015factor,mazo2016,schamberger2017,tan2019 and krupskii2020. In general, this type of factor copula is not an implicit copula. The second type of factor copula is the implicit copula of a traditional elliptical or skew-elliptical factor model. This type of copula emerged in the finance literature for low-dimensional applications laurent2005, but is increasingly used to model dynamic dependence in high dimensions; see creal2015,oh2017,oh2018 and oh2020dynamic. Estimation issues grow with the dimension and complexity of the copula, and this remains an active field of research.
One of the simplest factor copulas is a Gaussian static factor copula, which laurent2005 suggest for a single factor, and murray2013 consider for a larger number of factors. The multiple factor copula can be defined as follows. Let $\widetilde{\bm{Z}} \sim N_m(\bm{0},\Lambda \Lambda^\top +D)$, where $\Lambda=\{\lambda_{j,k}\}$ is an $m\times p$ matrix of factor loadings, $D=\mbox{diag}(d_1,\ldots,d_m)$ is a diagonal matrix of idiosyncratic variations, and typically $p<<m$. The implicit copula of $\widetilde{\bm{Z}}$ is a Gaussian copula, as outlined in Section (ref). To derive the parameter matrix $\Omega$, set the diagonal matrix \[ S=\mbox{diag}(\Lambda \Lambda^\top +D)=\mbox{diag}\left(\sum_{k=1}^p \lambda_{1,k}^2+d_1,\ldots,\sum_{k=1}^p \lambda_{m,k}^2+d_m\right)\,, \] then $\bm{Z}=S^{-1/2}\widetilde{\bm{Z}}$, so that $\Omega=S^{-1/2}(\Lambda \Lambda^\top + D)S^{-1/2}$.
murray2013 identify the loadings and idiosyncratic variations by setting $D=I$, the upper triangular elements of $\Lambda$ to zero and the leading diagonal elements to positive values $\lambda_{i,i}>0$. The copula parameters are then $\text{\boldmath$\theta$}=(\mbox{vecl}(\Lambda),d_1,\ldots,d_m)$, where $\mbox{vecl}(\Lambda)$ is the half-vectorization operator applied to the lower triangle of the rectangular matrix $\Lambda$. In the non-copula factor model literature, there are alternative ways to identify $\Lambda$ and $D$ kaufmann2017,fruhwirth2018, and similar restrictions may be adapted for the correlation matrix $\Omega$ as well. In a Bayesian framework, priors also have to be adopted for $\Lambda$ and $D$, and these can be used to provide further regularization as in murray2013 and elsewhere.
Simulation from this factor copula model is fast using the latent variable representation of the factor structure given by $\text{\boldmath$\eta$}\sim N_p(\bm{0},I)$ and $\widetilde{\bm{Z}}|\text{\boldmath$\eta$} \sim N_m(\Lambda \text{\boldmath$\eta$},D)$. To do so, at Step 1 of Algorithm (ref) generate a draw $\text{\boldmath$z$}\sim F_Z$ by drawing sequentially as follows:
Copulas have been used extensively to capture the cross-sectional dependence in multivariate time series; see patton2012 for a review. However, they can also be used to capture the serial dependence in a univariate series. The resulting time series models are extremely flexible, and there are many potential applications to continuous, discrete or mixed data.
If $\bm{Y}=(Y_1,\ldots,Y_T)^\top$ is a time series vector, then the copula $C$ at (ref) with $m=T$ captures the serial dependence in the series and is called a “time series copula”. While there has been less work on time series copulas than those used to capture cross-sectional dependence, they are increasingly being used for both time series data (where there is a single observation on the vector $\bm{Y}$) and longitudinal data (where there are multiple observations on the vector $\bm{Y}$). Early contributions include darsow1992, joe97, freeswang2005, chen2006tscop, ibragimov2009 and beare2010 for Markov processes, Wilson2010 for the implicit copulas of Gaussian processes popular in machine learning, and smith2010vine for vine copulas that exploit the time ordering of the elements of $\bm{Y}$.
For a continuous-valued stochastic process $\{Y_t\}$, denote the copula model for the joint density of time series variables $\bm{Y}_{1:t}=(Y_1,\ldots,Y_t)^\top$ as \[ f_{Y_{1:t}}(y_1,\ldots,y_t)=c_{1:t}(u_1,\ldots,u_t)\prod_{s=1}^t f_{Y_s}(y_s)\,, \] where $c_{1:t}$ is a $t$-dimensional copula density that defines a {\em copula process} for stochastic process $\{U_t\}$, with $U_t=F_{Y_t}(Y_t)$. Then the conditional distribution $Y_{t+1}|\bm{Y}_{1:t}$ has density
Here, $f_{U_{t+1|1:t}}$ is the density of $(U_{t+1}|U_1,\ldots,U_t)$, which is not uniform on $[0,1]$ (whereas the marginal distribution of $U_{t+1}$ is uniform on $[0,1]$). This conditional density can be used to form predictions from the copula model. It can also be used in likelihood-based estimation because $f_Y(\text{\boldmath$y$})= \prod_{t=2}^T\left\{ f_{U_{t|1:t-1}}(u_{t}|u_1,\ldots,u_{t-1})f_{Y_{t}}(y_{t})\right\}f_{Y_1}(y_1)$, with $\text{\boldmath$y$}=(y_1,\ldots,y_T)^\top$, which can be computed efficiently for many choices of copula $c_{1:T}$. In drawable vine copulas (D-vines) $f_{U_{t+1|1:t}}$ is further decomposed into a product of bivariate copulas called “pair-copulas” AasCzaFriBak2009, allowing for a flexible representation of the serial dependence structure, as discussed by smith2010vine, beare2015, smith2015, loaiza2018hetero, bladt2021 and others.
For longitudinal data with a sufficient number of observations on $\bm{Y}$, it is possible to estimate the marginal distribution functions $F_{Y_1},\ldots, F_{Y_T}$ at (ref) separately as in smith2010vine. But for time series data it is necessary to impose some structure on these marginal densities. For example, freeswang2005,freeswang2006 employ generalized linear regression models with time-based covariates in an actuarial setting. In the absence of common covariates, the marginals may be assumed time-invariant, so that $F_{Y_t}\equiv G$ for all $t$ as in chen2006tscop and smith2015. Flexible marginals, such as a skew $t$ distribution, or non-parametric estimators such as smoothed empirical distribution functions or kernel density estimators, can be used.
Time series copulas can also be used for discrete-valued data; see joe97 for an early exploration of such models. smithkhaled2012 do so for longitudinal data using the extended likelihood at (ref) and the copula decomposition above, so that for $\text{\boldmath$y$}=(y_1,\ldots,y_T)^\top$ and $\text{\boldmath$u$}=(u_1,\ldots,u_T)^\top$,
with $U_1$ marginally uniform on $[0,1]$. These authors employ a D-vine copula, and show how estimation using this extended likelihood can be undertaken by Bayesian data augmentation, where the values of $\text{\boldmath$u$}$ are generated in an MCMC sampling scheme. Alternatively, loaiza2019VB show how to estimate the copula parameters using variational Bayes methods blei2017. These calibrate tractable approximations to the augmented posterior obtained from the extended likelihood above. They call this approach “variational Bayes data augmentation” (VBDA) and show it is faster than MCMC and can be employed for much larger $T$ for many choices of copula.
In early work, lambert2002 and freeswang2005,freeswang2006 suggested adopting the implicit copula of an auxiliary stochastic process $\{Z_t\}$. In this case, the copula density $c_{1:t}$ has the form at (ref), so that for $t\geq 2$ \[ c_{1:t}(u_1,\ldots,u_t)=f_{Z_{1:t}}(z_1,\ldots,z_t)/\prod_{s=1}^t f_{Z_s}(z_s). \] The conditional density at (ref) is therefore
A major advantage of an implicit time series copula is that for many processes $\{Z_t\}$, the densities $f_{Z_{t+1|1:t}}$ and $f_{Z_{t+1}}$ are straightforward to compute and simulate from, simplifying parameter estimation and evaluation of predictive distributions. It is straightforward to show (e.g. see chen2006tscop,smith2015) that if $\{Z_t\}$ is a (strongly) stationary stochastic process, then $F_{Z_t}$ is time invariant and $\{U_t\}$ is also stationary because $U_t=F_{Z_t}(Z_t)$ is a monotonic transformation. In addition, if the marginal distribution $F_{Y_t}$ is also time invariant, the process $\{Y_t\}$ is also stationary.
For an implicit copula, the extended likelihood at (ref) based on $\bm{Z}=(Z_1,\ldots,Z_T)^\top$ with realization $\text{\boldmath$z$}=(z_1,\ldots,z_T)^\top$ can be used instead of that at (ref), which is
The simplest implicit time series copulas are those based on stationary Gaussian time series models. cario1996 and joe97 suggest using a zero mean stationary autoregression of lag length $p$, so that \[ Z_s = \sum_{k=1}^p \rho_k Z_{s-k}+ e_s\,,\mbox{ for } s=1,2,\ldots\,, \] with $e_s\sim N(0,\sigma^2)$ an independent disturbance, and parameters $\{\rho_1,\ldots,\rho_p,\sigma^2\}$. Then $\bm{Z}_{1:t}=(Z_1,\ldots,Z_t)^\top \sim N_t(\bm{0},\sigma^2\Sigma_{1:t})$, with $\sigma^2\Sigma_{1:t}$ the usual full rank autocovariance matrix with $\Sigma^{-1}_{1:t}$ a band $p$ matrix that is a function of $\text{\boldmath$\rho$}=(\rho_1,\ldots,\rho_p)^\top$ only.
Therefore, the implicit copula of $\bm{Z}_{1:t}$ is the Gaussian copula $C_{\mbox{\tiny Ga}}(\text{\boldmath$u$};\Omega_{1:t})$ with the autocorrelation matrix $\Omega_{1:t}=\mbox{diag}(\Sigma_{1:t})^{-1/2}\, \Sigma_{1:t} \,\mbox{diag}(\Sigma_{1:t})^{-1/2}$. The parameter $\sigma$ does not feature in $\Omega_{1:t}$ (i.e. it is unidentified in the copula), so that it is sufficient to fix it to an arbitrary value such as $\sigma^2=1$, as is done here. Thus, $\Omega_{1:t}$ is only a function of $\text{\boldmath$\rho$}$, so that $\text{\boldmath$\theta$}=\text{\boldmath$\rho$}$ are the copula parameters. The marginal distribution $Z_t\sim N(0,\gamma_0)$, with variance $\gamma_0$ computed from $\text{\boldmath$\rho$}$. Denoting the density of a standard normal as $\phi(\cdot)$, and that of a $N(\mu,\sigma^2)$ as $\phi_1(\cdot;\mu,\sigma^2)$, the conditional density
with $z_{t}=\Phi^{-1}_1(u_{t};0,\gamma_0)$ a $N(0,\gamma_0)$ distribution function evaluated at $u_t$. (The dependence of this conditional density on $\text{\boldmath$\theta$}$ is tacit here.) Thus, the likelihood of a continuous-valued series, or the extended likelihood of a discrete-valued series, can be expressed in terms of the copula parameters $\text{\boldmath$\rho$}$ and the marginals $F_{Y_1},\ldots,F_{Y_T}$. A variety of estimation methods, including standard maximum likelihood, can then be used to estimate the time series copula parameters.
This copula model extends the stationary autoregression from a marginally Gaussian process to one with any other marginal distribution. This is why cario1996 originally labeled it an “autoregression-to-anything” transformation, although these authors did not recognize it as a Gaussian copula. Interestingly, even though the auxiliary stochastic process $\{Z_t\}$ is conditionally homoscedastic (i.e. $\mbox{Var}(Z_{t+1}|Z_{1:t})=1$) the process $\{Y_t\}$ need not be so (i.e. it can be heteroscedastic). To see this, notice that even when $f_{Y_t}=g$ is time invariant, the conditional density of $Y_{t+1}|Y_{1:t}$ is \[ f_{Y_{t+1|1:t}}(y_{t+1}|y_{1},\ldots,y_t)= \phi\left(z_{t+1}-\sum_{k=1}^p \rho_k z_{t-k+1} \right) \frac{g(y_{t+1})}{\phi_1\left(z_{t+1};0,\gamma_0\right)}\,. \] The second moment of this density is not necessarily a constant with respect to time, as demonstrated in smith+vahey2016.
The usual measures of serial dependence for an autoregression (e.g. autocorrelation or partial autocorrelation matrices) can be computed for $\{Z_t\}$. Spearman correlations, which are unaffected by the choice of continuous margin(s) $F_{Y_t}$, provide equivalent metrics for $\{Y_t\}$. For example, the Spearman autocorrelation at lag $h$ is \[ \rho^S_{h}=\frac{6}{\pi} \mbox{arcsin}\left(\frac{\gamma_h}{2\gamma_0 }\right)\,, \] where $\gamma_h\equiv \mbox{Cov}(Z_{t+h},Z_t)$ is the autocovariance at lag $h$ for the auxiliary stochastic process, and is a function of $\text{\boldmath$\rho$}$. Other popular measures of concordance, can also be computed easily for different values of $h$.
Last, while the Gaussian autoregression copula---or indeed other Gaussian time series copulas, such as those based on Gaussian processes Wilson2010---produces a flexible family of time series models, the form of serial dependence is still limited. For example, serial dependence is both symmetric and has zero tail dependence, which are properties of the Gaussian copula. This motivates the construction of more flexible time series copulas, as now discussed.
A wide array of time series and other statistical models can be written in state space form; see durbin2012 for an overview of this extensive class. smithman2018 outline how to construct and estimate the implicit time series copulas of such models, as is now outlined.
A nonlinear state space model for $\{Z_t\}$ is given by the observation and transition equations
Here, $H_t$ is the distribution function of $Z_t$, conditional on an $r$-dimensional state vector $\bm{X}_t$. The states follow a Markov process, with conditional distribution function $K_t$. Typically, tractable parametric distributions are adopted for $H_t$ and $K_t$, with the parameters denoted collectively as $\bm{\theta}$.
A key requirement in evaluating (ref) and (ref) is the computation of the marginal distribution and density functions of $Z_t$. Marginalizing over $\bm{X}_t$ gives these as
where the dependence on $\bm{\theta}$ is denoted explicitly here. The density $h_t(z_t|\bm{x}_t;\bm{\theta}) =\frac{d}{d z_t} H_t(z_t|\bm{x}_t;\bm{\theta})$, and $f_{X_t}(\bm{x}_t|\bm{\theta})$ is the marginal density of the state variable $\bm{X}_t$. Evaluation of the integrals in ((ref)) is straightforward either analytically or numerically for many choices of state space model used in practice. Note that the quantile function $z_t=F^{-1}_{Z_t}(u_t|\bm{\theta})$ is a function of $\text{\boldmath$\theta$}$, which can be computed quickly using the interpolation method outlined in (ref) when $F_{Z_t}$ is time invariant.
A more challenging problem is the evaluation of the numerator in (ref). To compute this, the state vector $\bm{x}=(\bm{x}_1^\top,\ldots,\bm{x}_T^\top)^\top$ with $Tr$-dimensional joint density $f_X$ needs to be integrated out, with
where $k_t(\bm{x}_t|\bm{x}_{t-1};\bm{\theta})=\frac{d}{d \bm{x}_t} K_t(\bm{x}_t|\bm{x}_{t-1};\bm{\theta})$. While there a number of existing methods in the state space literature to evaluate $f_Z(\bm{z}|\bm{\theta})$ above, robust Bayesian MCMC methods that generate the states $\text{\boldmath$x$}$ are very popular. The same methods can also be employed estimate the implicit copula as outlined below.
Conditional on the states, a continuous time series copula model likelihood is
where all components on the right-hand side of (ref) are known densities. Computationally, it is much easier to work with (ref), rather than with the decomposition (ref) and copula density (ref). Adopting the prior $\pi_\theta(\bm{\theta})$, Bayesian estimation and inference of the copula parameters $\text{\boldmath$\theta$}$ can be based on the MCMC sampler at Algorithm (ref) below, which produces Monte Carlo draws from the posterior of $\text{\boldmath$\theta$}$ augmented with the latent states $\text{\boldmath$x$}$.
Unlike the states $\text{\boldmath$x$}$, the values $\bm{z}=(z_1,\ldots,z_T)^\top$ are not generated in the sampling scheme, but instead are computed as $z_t=F_{Z_t}^{-1}(u_t|\text{\boldmath$\theta$})$ for each draw of the parameters $\bm{\theta}$. Crucially, Step 1 is exactly the same as that for the underlying state space model, so that any of the wide range of existing procedures for generating $\text{\boldmath$x$}$ can be employed. Step 2 can be undertaken using a Metropolis-Hastings step, with a proposal based on a numerical or other approximation to the conditional posterior. In Algorithm (ref) the marginal distributions $F_{Y_1},\ldots,F_{Y_T}$ are assumed known. It is common to estimate these prior to estimating the copula parameters joe2005, although joint estimation of the marginals and copula parameters may also be considered. In a Bayesian analysis the prior $\pi_\theta(\bm{\theta})$ reflects any constraints required to identify $\text{\boldmath$\theta$}$.
smithman2018 constructed the implicit copulas of three specific state space models, and estimated their parameters for U.S. inflation between 1954:Q1 and 2013:Q4. These included an unobserved component stochastic volatility (UCSV) model, as is now outlined. To illustrate, it is then applied to the same quarterly U.S. inflation series used by these authors, but updated to include all observations up to 2020:Q2. This includes the impact of the Covid-19 pandemic, which this flexible copula time series model is well-suited to capture.
The UCSV model is specified for bivariate state vector $\bm{x}_t=(\mu_t,\zeta_t)^\top$ as
The parameters $|\rho_\mu|<1$ and $|\rho_\zeta|<1$, which ensures $\{Z_t\}$ is a (strongly) stationary first order Markov process. The mean $E(Z_t)=\bar \mu$, which is unidentified in the implicit copula at (ref), and set $\bar \mu=0$ here. The marginal variance $\mbox{Var}(Z_t)=s^2_\mu+\exp(\bar \zeta + s^2_\zeta/2)$, where $s^2_\mu=\sigma^2_\mu/(1-\rho_\mu^2)$ and $s^2_\zeta=\sigma^2_\zeta/(1-\rho_\zeta^2)$. The variance $\mbox{Var}(Z_t)$ is unidentified in the copula, and setting this equal to one provides an equality constraint on $\bar \zeta = \log(1-s^2_\mu)-\frac{s_\zeta^2}{2}$. In addition, $\exp(\bar \zeta+s^2_\zeta/2)\geq 0$, giving the inequality constraint $0<\sigma^2_\mu\leq (1-\rho_\mu^2)$. With these identifying constraints, the dependence parameters of the resulting implicit copula are $\bm{\theta}=\{\rho_\mu,\rho_\zeta,\sigma^2_\mu,\sigma^2_\zeta\}$.
Because $\{Z_t\}$ is stationary, the marginal density $f_{Z_t}$ at ((ref)) is time-variant and given by \[ f_{Z_1}(z;\bm{\theta})=\int\int \phi_1\left( z;\mu,\exp(\zeta)\right) \phi_1(\zeta;\bar \zeta, s_\zeta^2) \phi_1(\mu;0,s^2_\mu)d\mu d\zeta\,. \] The integral in $\mu$ can be recognized as that of a Gaussian density to give
with $w(\zeta)^2=s^2_\mu+\exp(\zeta)$. Computing the (log) copula density at ((ref)) requires evaluating $\log(f_{Z_1})$ and the quantile function $F_{Z_1}^{-1}$ at all $T$ observations. To do so, the accurate and fast numerical method described in (ref) is used.
The parameters of this time series copula model are estimated using their Bayesian posterior with the prior $\pi_{\theta}(\text{\boldmath$\theta$})\propto \frac{1}{\sigma^2_\mu \sigma^2_\zeta}\mathds{1}(\text{\boldmath$\theta$} \in R_\theta)$, where $R_\theta$ is the region of parameter values that correspond to the constraints outlined above. Algorithm (ref) can be used to estimate the copula parameters, where at Step 1 the state vector $\text{\boldmath$x$}$ is partitioned into $\bm{\mu}=(\mu_1,\ldots,\mu_T)^\top$ and $\bm{\zeta}=(\zeta_1,\ldots,\zeta_T)^\top$, and generated using the two separate steps:
The posterior of $\bm{\mu}$ in Step 1a can be recognized as normal with zero mean and a band one precision matrix, so that generation is both straightforward and fast. There are a number of efficient methods to generate $\bm{\zeta}$ in Step 1b in the literature, and the fast “precision sampler” for the latent states outlined in chan2009 is used here. In Step 2 of the sampler, a normal approximation is used as a proposal density for the Metropolis-Hastings step, which has high acceptance rates in practice.
The adaptive kernel density estimator (AKDE) of shimazaki2010 is used to estimate a time-invariant marginal distribution $G$ of $Y_t$, and is presented in Figure (ref). The estimated density is smooth, positively skewed, and heavy-tailed; it accounts for both high (e.g. $2.9\%$ in 1974:Q3) and low (e.g. $-0.529\%$ in 2020:Q1) values. Figure (ref) plots the time series, plus the copula data $u_t=G(y_t)$ for $t=1,\ldots,T$.
To summarize the posterior estimate of the implicit copula, Figure (ref) plots the posterior means and 90% posterior intervals for $\text{\boldmath$\mu$}$ and $\exp(\text{\boldmath$\zeta$}/2)$, which are the mean and standard deviation of the auxiliary vector $\bm{Z}$. While these are not the mean and standard deviation of $\bm{Y}$, they do account for movements in the moments of this variable, and the impact of the Covid-19 pandemic on 2020 can be seen as a sharp jump in $\exp(\zeta_t/2)$ in panel (a), while the inflationary period of the 1970's can be see in high values of $\mu_t$ in panel (b).
There is high serial dependence in both state variables. One way to show how this affects the time series copula is to consider the bivariate margin $c_{t-1:t}(u_{t-1},u_t|\text{\boldmath$\theta$})$ of the copula density, which is time invariant. It is given by $c_{1:2}(u_1,u_2|\text{\boldmath$\theta$})=f_{Z_{1:2}}(z_1,z_2|\text{\boldmath$\theta$})/f_{Z_1}(z_1|\text{\boldmath$\theta$})f_{Z_1}(z_2|\text{\boldmath$\theta$})$, where the numerator is computed by numerical integration. Figure (ref) plots this density at the posterior mean of the copula parameters $\text{\boldmath$\theta$}$, and two interesting features can be seen. First, “spikes” at the corners (i.e. near (0,0), (0,1), (1,0) and (1,1)) are indicative of strong dependence in the {\em volatility} of the series; see loaiza2018hetero and bladt2021 for a discussion of such a pattern in a time series copula. Second, the positive “ridge” running from (0,0) to (1,1) is indicative of positive dependence in the {\em level} of the series. Both these features are well-known aspects of inflation time series, and the implicit copula captures them both while also allowing for the asymmetric marginal distribution in Figure (ref). The time series copula model therefore allows for more realistic modeling of tail risk than the standard UCSV model, and improves density forecast accuracy, including in the tails.
Copulas have been used extensively to capture cross-sectional dependence in a multivariate stochastic process $\{\bm{Y}_t\}$, where $\bm{Y}_t=(Y_{1,t},\ldots,Y_{d,t})^\top$; see patton2006, rodriguez2007, hafner2012 and creal2015 for just some examples. These models typically capture serial dependence through existing marginal time series models; for example, heteroscedastic models are normally used for financial returns. An alternative is to use a single high-dimensional---but parsimonious---copula to capture both serial and cross-sectional dependence jointly. An advantage of this approach is that the marginal distribution of each variable can be modeled directly, including as non-parametric. It is this type of time series copula that is the focus of this section.
If the random vector $\bm{Y}=(\bm{Y}_1^\top,\ldots,\bm{Y}_T^\top)^\top$, then $F_Y$ is given by (ref) with $m=Td$ and the order of the elements of $\bm{Y}$ determines the interpretation of $C$. If all variables are continuous, $f_Y$ is given by (ref), so that
where $\text{\boldmath$y$}=(\text{\boldmath$y$}_1^\top,\ldots,\text{\boldmath$y$}_T^\top)^\top$, $\text{\boldmath$y$}_t=(y_{1,t},\ldots,y_{d,t})^\top$, $\text{\boldmath$u$}=(\text{\boldmath$u$}_1^\top,\ldots,\text{\boldmath$u$}_T^\top)^\top$ and $\text{\boldmath$u$}_t=(u_{1,t},\ldots,u_{d,t})^\top$. For discrete-valued variables the mass function is given by (ref), and an extended likelihood for $(\bm{Y},\bm{U})$ is given by (ref); see loaiza2019VB. Marginal models for each of the $d$ series are required, and one option is to assume they are time-invariant with distribution functions $G_1,\ldots,G_d$, which can be estimated separately.
Selecting an appropriate $Td$-dimensional copula with density $c$ at (ref) is difficult because it needs to capture three forms of dependence: (i) cross-sectional contemporaneous, (ii) within-series serial, and (iii) cross-series serial. One solution is to use a vine copula; for example see brechmann2015 for the two-dimensional case, smith2015 and loaiza2018hetero for D-vine copulas, beare2015 for an M-vine and zhao2020 for an alternative vine-based copula; see also remillard2012,nagler2020. However, implicit copulas constructed from existing multivariate time series models offer a tractable alternative to vines, particularly for series where $T$ and/or $d$ are large.
The most popular implicit copula for multivariate time series is that of a Gaussian vector autoregression (VAR) for $\{\bm{Z}_t\}$. This is an extension of the autoregression copula in Section (ref). Consider the following VAR with lag $p$,
The mean is set to zero because it is unidentified in the copula, and the variances are fixed so that $\mbox{Var}(Z_{j,t})=1$. Then $\bm{Z}=(\bm{Z}_1^\top,\ldots, \bm{Z}_T^\top)^\top\sim N_{Td}(\bm{0},\Omega)$, where $\Omega$ is the block Toeplitz correlation matrix of this process. This is a matrix of $(T\times T)$ blocks, with the $(s,t)$th block being given by $\Omega_h \equiv \mbox{Corr}(\bm{Z}_{t+h},\bm{Z}_t)$ for $h=|t-s|$ and $t\geq s$; for example, see lutkepohl2005. Because the Gaussian copula is closed under marginalization, the $d$-dimensional marginal distribution in $\bm{Y}_t$ also has a Gaussian copula function $C_{\mbox{\tiny Ga}}(\text{\boldmath$u$}_t;\Omega_0)$. For example, for continuous data $f_{Y_t}(\text{\boldmath$y$}_t)=c_{\mbox{\tiny Ga}}(\text{\boldmath$u$}_t;\Omega_0)\prod_{j=1}^d f_{Y_{j,t}}(y_{j,t})$.
For continuous time series, a straightforward approach to estimate the model is to first estimate appropriate marginal distributions $F_{Y_{j,t}}$; for example, by assuming time-invariance in the marginals and applying a kernel density estimator to each of the $d$ series. Second, compute the auxiliary data $z_{j,t}=\Phi^{-1}(F_{Y_{j,t}}(y_{j,t}))$ for $i=1,\ldots,d$ and $t=1,\ldots,T$. Third, apply standard likelihood-based methods for Gaussian VARs directly to this auxiliary data to estimate the unknown parameters $B_1,\ldots,B_p,\Sigma$. From these $\Omega$ can be computed, although this can be impractical to evaluate when $m=Td$ is large and there is often no need to do so.
The conditional density for $\bm{Y}_{t+1}|\bm{Y}_{1:t}$ can be derived in a similar manner as for the univariate autoregression copula model. This is given by
where time-invariant marginal densities $g_1,\ldots,g_d$ are assumed. Drawing from this conditional distribution is straightforward by first simulating $\bm{Z}_{t+1}$ directly from (ref), and then transforming to a draw $\bm{Y}_{t+1}=(G_1^{-1}(\Phi(Z_{1,t+1})),\ldots,G_d^{-1}(\Phi(Z_{d,t+1})))^\top$, which can be used to compute the predictive distribution.
biller2003 were the first to construct the Gaussian VAR copula via transformation, but did not recognize it as a Gaussian copula and called it a “Vector-Autoregressive-To-Anything” distribution. smith2015 and smith+vahey2016 also consider this Gaussian copula model, its D-vine representation and apply it to multivariate macroeconomic and financial forecasting. Similar to the univariate case, the Gaussian VAR copula model can capture a degree of heteroscedasticity in the time series given a suitable choice of marginal distributions $G_1,\ldots,G_d$ for the $d$ time series; see smith+vahey2016 for a demonstration. There is also a growing interest in multivariate times series copulas in machine learning. For example, salinas2019high construct a Gaussian copula from low rank factor decomposition where the small number of factors follow a Gaussian process with recurrent neural network (RNN) dynamics. klein2020deep propose constructing a Gaussian copula process that is constructed as the implicit copula of an RNN with Gaussian errors.
Existing econometric applications of multivariate time series often have parameters that vary over time (widely called a “dynamic” model) along with substantial regularization; see bitto2019, huber2020 and references therein. These features can also be employed for the parameters of implicit copulas. For example, smith+vahey2016 use Bayesian selection on the D-vine representation of the Gaussian VAR copula for regularization, creal2015, oh2017 and opschoor2020 allow the parameters of elliptical copulas to vary over time, and oh2020dynamic consider a dynamic skew $t$ copula. In another approach loaiza2018hetero extend the UCSV model in Section (ref) to the multivariate case and show how to construct its implicit copula. In all these studies, the copula models are more accurate than non-copula benchmarks, and the implicit copulas used are scalable to high dimensions.
Copula models with regression margins have been used widely; for examples, see pitt2006, song2009, masarotto2012 and klein2016. However, another usage of a copula with regression data is to capture the dependence between multiple observations on a single dependent variable $Y$, conditional on the covariate values. This defines a copula process Wilson2010 on the covariate space, which smith+k19 call a “regression copula”. When combined with a flexible marginal distribution for $Y$, it specifies a new distributional regression model. This is where the covariates affect the entire distribution of $Y$. KleSmi2019 and smith+k19 consider a regression copula that is the implicit copula of the joint distribution of observations in an auxiliary regression model. They are inherently high dimensional, yet can be estimated in reasonable time using Bayesian methods. The idea is outlined in this section for continuous $Y$, and greater detail can be found in these papers.
Consider $N>1$ realizations $\bm{Y}_{1:N}=(Y_1,\ldots,Y_N)^\top$ of a dependent variable with corresponding values $\text{\boldmath$x$}_{1:N}=\{\text{\boldmath$x$}_1,\ldots,\text{\boldmath$x$}_N\}$ for $p$ covariates, with $\text{\boldmath$x$}_i=(x_{i,1},\ldots,x_{i,p})^\top$. Then application of Sklar's theorem to the distribution of $\bm{Y}_{1:N}|\text{\boldmath$x$}_{1:N}$ gives \[ F_{Y_{1:N}}(\text{\boldmath$y$}_{1:N}|\text{\boldmath$x$}_{1:N})=C^\dagger_{1:N}\left(F_{Y_1}(y_1|\text{\boldmath$x$}_1),\ldots, F_{Y_N}(y_N|\text{\boldmath$x$}_{N})\,;\,\text{\boldmath$x$}_{1:N}\right)\,. \] The $N$-dimensional copula function $C^\dagger_{1:N}(\cdot\,;\,\text{\boldmath$x$}_{1:N})$ is a copula process on the covariate space, and $F_{Y_i}(y_i|\text{\boldmath$x$}_i)$ is the distribution function of $Y_i|\text{\boldmath$x$}_i$. Both are typically unknown, and in a copula model these are selected to define the distribution. One tractable but effective simplification is to allow the covariates to only affect the dependent variable through the copula function, so that $Y_i$ is {\em marginally} independent of $\text{\boldmath$x$}_i$. In this case,
with $\text{\boldmath$\theta$}$ unknown copula parameters that are unaffected by the dimension $N$ and require estimation. Here, the {\em joint} distribution of $\bm{Y}_{1:N}$ is dependent on $\text{\boldmath$x$}_{1:N}$ via the copula, so that the conditional distribution $Y_{N}|(\bm{Y}_{1:N-1}=\text{\boldmath$y$}_{1:N-1}),\text{\boldmath$x$}_{1:N}$ is also. The latter is employed as the predictive distribution of the regression model, as discussed further below.
When the dependent variable is continuous, the joint density is
with $u_i=F_{Y_i}(y_i)$. An advantage (that is also in common with the time series copula models discussed in Section (ref)) is that if $F_{Y_i}(y_i)\equiv G(y_i)$ is assumed to be invariant with respect to the index $i$, then $G$ can be estimated using non-parametric or other flexible estimators. The remaining component of the copula model at (ref) is the choice of copula process, which is aptly called a regression copula because it is a function of $\text{\boldmath$x$}_{1:N}$.
To see how (ref) defines a distributional regression model, consider the predictive density for a continuous-valued dependent variable. For a sample of size $n$ with covariate values $\text{\boldmath$x$}_{1:n}$ and dependent variable values $\bm{Y}_{1:n}=\text{\boldmath$y$}_{1:n}$ arising from (ref), the predictive distribution of the subsequent value $Y_{n+1}$ with observed covariates $\text{\boldmath$x$}_{n+1}$ is defined to be that of $Y_{n+1}|\text{\boldmath$x$}_{1:n+1},\text{\boldmath$y$}_{1:n}$, which has density
Thus, the predictive density is a function of the covariate vector $\text{\boldmath$x$}_{n+1}$, as well as those of the sample $\text{\boldmath$x$}_{1:n}$. Moreover, the entire distribution (not just the first or other moments of $Y_{n+1}$) is a function of $\text{\boldmath$x$}_{n+1}$ as illustrated empirically in Section (ref).
One regression copula process $C_{1:N}$ that can be used at (ref) is an implicit copula derived from an existing regression model, as now discussed.
Implicit regression copulas are constructed as in Section (ref), but when also conditioning on the covariate values; i.e. from an “auxiliary regression” model. Consider a regression model for the auxiliary vector $\bm{Z}_{1:N}=(Z_1,\ldots,Z_N)^\top$ with covariate values $\text{\boldmath$x$}_{1:N}$ and parameter vector $\text{\boldmath$\theta$}$. Denote the joint distribution function of $\bm{Z}_{1:N}|\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$}$ as $F_{Z_{1:N}}(\cdot|\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$})$, with $i$th marginal $F_{Z_i}(\cdot|\text{\boldmath$x$}_i,\text{\boldmath$\theta$})$. Then, extending the definition in Table (ref), the following transformations define a regression copula model \[ U_i=F_{Z_i}(Z_i|\text{\boldmath$x$}_i,\text{\boldmath$\theta$})\,,\mbox{ and } Y_i=F_{Y_i}^{-1}(U_i)\,. \] If $\text{\boldmath$z$}_{1:N}=(z_1,\ldots,z_N)^\top$, $z_i=F_{Z_i}^{-1}(u_i|\text{\boldmath$x$}_i,\text{\boldmath$\theta$})$ and $m=N$, then the implicit copula function at (ref) and density at (ref) for this model are given by
In (ref) $f_{Z_i}(z_i|\text{\boldmath$x$}_i,\text{\boldmath$\theta$})$ is the density function of the auxiliary variable $Z_i$, conditional on the covariates $\text{\boldmath$x$}_i$. These expressions for $C_{Z_{1:N}}$ and $c_{Z_{1:N}}$ can then be used in (ref) and (ref) to specify a distributional regression.
If in the auxiliary regression $f_{Z_{1:N}}(\text{\boldmath$z$}_{1:N}|\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$})=\prod_{i=1}^N f_{Z_i}(z_i|\text{\boldmath$x$}_i,\text{\boldmath$\theta$})$, then from (ref) the implicit copula is the trivial independence copula. Thus, only distributions where $\bm{Z}_{1:N}$ are dependent are useful for constructing an implicit regression copula, as in Section (ref) below.
Employing the copula density at (ref) for that in the predictive density at (ref), gives the following:
To evaluate (ref) in practice, a point estimate of $\text{\boldmath$\theta$}$ can be used. In a Bayesian analysis another option exists, where $\text{\boldmath$\theta$}$ is integrated out with respect to its posterior density $f(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ to obtain \[ f_{\mbox{\tiny pred}}^{\mbox{\tiny Bayes}}(y_{n+1}|\text{\boldmath$x$}_{n+1})= \int f_{\mbox{\tiny pred}}(y_{n+1}|\text{\boldmath$x$}_{n+1},\text{\boldmath$\theta$})f(\text{\boldmath$\theta$}|\text{\boldmath$y$})\mbox{d}\text{\boldmath$\theta$}\,. \] This is called the “posterior predictive density”, and evaluation of the integral is usually undertaken using draws obtained from an MCMC sampling scheme.
In principle, implicit copula processes outlined above can be constructed from a wide range of different regression models. KleSmi2019 suggest doing so for a Gaussian linear regression, as now outlined.
For a dependent variable $\widetilde{Z}_i$, consider the linear regression \[ \widetilde{Z}_i=\text{\boldmath$x$}_i^\top \text{\boldmath$\beta$} + \sigma e_i\,, \] with $e_i$ distributed independently $N(0,1)$. Conditional on both $\text{\boldmath$x$}_i$ and the parameters $\text{\boldmath$\beta$},\sigma^2$, the elements of $\widetilde{\bm{Z}}_{1:N}=(\widetilde{Z}_1,\ldots,\widetilde{Z}_N)^\top$ are distributed independently, so that their joint distribution cannot be used directly to specify a useful regression copula with density at (ref). However, a Bayesian framework can be employed where $\text{\boldmath$\beta$}$ is treated as random and marginalized out of the distribution for $\widetilde{\bm{Z}}_{1:N}$, the elements of which are then dependent. From this distribution a useful implicit regression copula can be formed as below.
If $B =[\text{\boldmath$x$}_1|\text{\boldmath$x$}_2|\cdots|\text{\boldmath$x$}_N]^\top$ is the $(N \times p)$ regression design matrix, then the regression can be written as the linear model
The conjugate proper prior
is used, where the precision matrix $P(\text{\boldmath$\theta$})$ is of full rank $p$ and a function of $\text{\boldmath$\theta$}$. It is necessary to assume a proper prior for $\text{\boldmath$\beta$}$, because it ensures that the distribution with $\text{\boldmath$\beta$}$ integrated out is also proper. Doing so (by recognizing a normal in $\text{\boldmath$\beta$}$) gives
with $\Sigma=(B ^\top B + P(\text{\boldmath$\theta$}))^{-1}$. Application of the Woodbury formula further simplifies the variance matrix at (ref) as \[ \sigma^2(I - B \Sigma B ^\top)^{-1}=\sigma^2 \left(I+B P(\text{\boldmath$\theta$})^{-1} B ^\top \right)\,. \] The variance of an individual observation $i$ is the $i$th leading diagonal element of this matrix, so that $\mbox{Var}(\widetilde{Z}_{i}|\text{\boldmath$x$}_i,\text{\boldmath$\theta$},\sigma^2) =\sigma^2(1+\text{\boldmath$x$}_i^\top P(\text{\boldmath$\theta$})^{-1} \text{\boldmath$x$}_i)$.
The copula of any normal distribution is the Gaussian copula discussed in Section (ref). The parameter matrix $R$ is the correlation matrix of (ref), and it is obtained by standardizing $\widetilde{Z}_i$ to have unit variance as follows. Let $s_i=(1+\text{\boldmath$x$}_i^\top P(\text{\boldmath$\theta$})^{-1} \text{\boldmath$x$}_i)^{-1/2}$, then define the auxiliary variable of the implicit copula as $Z_i=\frac{s_i}{\sigma}\widetilde{Z}_i$. Thus, if the diagonal matrix $S(\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$})=\mbox{diag}(s_1,\ldots,s_N)$, then from (ref) the conditional distribution of $\bm{Z}_{1:N}=(Z_1,\ldots,Z_N)^\top$ is $\bm{Z}_{1:N}|\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$},\sigma^2 \sim N\left(\bm{0},R\right)$ with correlation matrix
and has copula function $C_{\mbox{\tiny Ga}}(\text{\boldmath$u$};R)$. This is a copula process on the covariate space because $R$ is a function of the covariate vector $\text{\boldmath$x$}_{1:N}$ (the notation $R$ and $R(\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$})$ is used interchangeably here.) The parameter $\sigma^2$ does not feature in $R$ and is unidentified in the copula, so that $\sigma^2=1$ can be assumed throughout.
Different implicit copulas can be constructed by using different conditionally Gaussian priors for $\text{\boldmath$\beta$}$ at (ref). KleSmi2019 explore three different choices, including the horseshoe prior of CarPol2010 which is outlined here. This prior provides regularization of $\text{\boldmath$\beta$}$ in the auxiliary regression. The prior is given by
see polsonscott12. The hyper-parameters of this prior are the parameters of the implicit copula $\text{\boldmath$\theta$}=(\text{\boldmath$\lambda$}^\top,\tau)^\top$, while the precision matrix $P(\text{\boldmath$\theta$})=\mbox{diag}(\text{\boldmath$\lambda$})^{-2}$ is diagonal.
For a sample of $n$ observations, from (ref) and (ref), the likelihood is \[ f_{Y_{1:n}}(\text{\boldmath$y$}_{1:n}|\text{\boldmath$x$}_{1:n},\text{\boldmath$\theta$})=f_{Z_{1:n}}(\text{\boldmath$z$}_{1:n}|\text{\boldmath$x$}_{1:n},\text{\boldmath$\theta$})\prod_{i=1}^n \frac{f_{Y_i}(y_i)}{f_{Z_i}(z_i|\text{\boldmath$x$}_i,\text{\boldmath$\theta$})} = \phi_n\left(\text{\boldmath$z$}_{1:n};\bm{0},R(\text{\boldmath$x$}_{1:n},\text{\boldmath$\theta$})\right) \prod_{i=1}^n \frac{g(y_i)}{\phi(z_i)}\,, \] for an invariant marginal distribution with density $f_{Y_i}=g$. However, even though the likelihood is available in closed form, for large $n$ evaluating and inverting the $(n\times n)$ matrix $R(\text{\boldmath$x$}_{1:n},\text{\boldmath$\theta$})$ to compute the likelihood is computationally demanding.
Instead, it is more efficient to use the likelihood also conditional on $\text{\boldmath$\beta$}$, and integrate out $\text{\boldmath$\beta$}$ using an MCMC scheme. (It is stressed here that doing so does not change the implicit copula specification.) First, note that from (ref) when also conditioning on $\text{\boldmath$\beta$}$ the vector $\bm{Z}_{1:n}=S\widetilde{\bm{Z}}_{1:n}\sim N(SB\text{\boldmath$\beta$},S^2)$ (with $\sigma^2=1$ and $S=S(\text{\boldmath$x$}_{1:n},\text{\boldmath$\theta$})$). Also, the Jacobian of the transformation from $\bm{Z}_{1:n}$ to $\bm{Y}_{1:n}$ is $J_{Z_{1:n}\rightarrow Y_{1:n}}=\prod_{i=1}^n g(y_i)/\phi(z_i)$. Then, by a change of variables, the likelihood also conditional on $\text{\boldmath$\beta$}$ is
which can be evaluated in $O(n)$ operations because $S$ is a diagonal matrix. A Bayesian approach that employs this conditional likelihood, evaluates the augmented posterior $f(\text{\boldmath$\beta$},\text{\boldmath$\theta$}|\text{\boldmath$y$}_{1:n})$ using the sampler at (ref). Implementation details for this sampler are given in KleSmi2019.
One way to compute the predictive density that avoids computing $R(\text{\boldmath$x$}_{1:n},\text{\boldmath$\theta$})$ or its inverse (and is therefore faster than alternatives), is to also condition on $\text{\boldmath$\beta$}$. By a change of variables from $Y_{n+1}$ to $Z_{n+1}$,
where $s_{n+1}=(1+\text{\boldmath$x$}_{n+1}^\top P(\text{\boldmath$\theta$})^{-1} \text{\boldmath$x$}_{n+1})^{-1/2}$ and $z_{n+1}=\Phi^{-1}(F_{Y_{n+1}}(y_{n+1}))$.
The draws for $\text{\boldmath$\beta$},\text{\boldmath$\theta$}$ from Algorithm (ref) can be used to either integrate out $\text{\boldmath$\beta$},\text{\boldmath$\theta$}$ with respect to the augmented posterior, or to compute plug-in point estimates for $\text{\boldmath$\beta$}$ and also $s_{n+1}$. If $F_{Y_{n+1}}$ is fixed to its estimate, and $\{\text{\boldmath$\beta$}^{[1]},\text{\boldmath$\theta$}^{[1]},\ldots,\text{\boldmath$\beta$}^{[J]},\text{\boldmath$\theta$}^{[J]}\}$ are the Monte Carlo draws, then the two Bayesian posterior estimators for the predictive density are:
with \[ \hat{\text{\boldmath$\beta$}} = \frac{1}{J}\sum_{j=1}^J \text{\boldmath$\beta$}^{[j]}\,,\mbox{ and } \hat{s}_{n+1} = \frac{1}{J}\sum_{j=1}^J \left(1+\text{\boldmath$x$}_{n+1}^\top P(\text{\boldmath$\theta$}^{[j]})^{-1} \text{\boldmath$x$}_{n+1}\right)^{-1/2}\,. \] In their empirical work, KleSmi2019 and smith+k19 found that estimates from these two estimators were very similar.
Linear regression is widely used to estimate financial asset pricing models, where the dependent variable is the (excess) return on a stock. Yet stock returns are distributed far from Gaussian, so that a Gaussian regression model is mis-specified. To illustrate the regression copula model it is used to model monthly excess returns on American Express Company (which has NYSE ticker symbol “AXP”) using data from 07/1972 to 10/2020. The marginal distribution is assumed invariant with respect to observation, so that $F_{Y_i}(y_i)=G(y_i)$. A three parameter asymmetric Laplace distribution is fit, which better accounts for the distribution of returns as highlighted by chen2012 and taylor2019. Figure (ref) plots the density of the fitted margin, which is both asymmetric and has very heavy tails.
The monthly values of the five factors suggested and described by fama2015five were used as covariates. These are market risk (MktRf), size (SMB), value (HML), profit (RMW) and investment (CMA) factors, with data obtained from Kenneth French's website. The first three factors are widely employed, while the inclusion of the additional two factors RMW and CMA is more controversial. The linear regression copula constructed using the horseshoe prior for $\text{\boldmath$\beta$}$ was estimated using Algorithm (ref). Even though the implicit copula is of dimension $n=580$, employing the conditional likelihood at (ref) means that estimation is tractable, with a computation time of only 32s to draw 10,000 iterates on a standard laptop.
Table (ref) summarizes the posterior estimates of the coefficients $\text{\boldmath$\beta$}$ and the regularization parameters $\text{\boldmath$\lambda$}$. Of the five covariates, only the traditional three (MktRt, SMB and HML) were significant (“significant” here refers to the whether, or not, zero falls into the 95% posterior intervals for each coefficient $\beta_i$.) Thus, evidence for the inclusion of the two new factors RMW and CMA is weak. Also reported are the Metropolis-Hastings acceptance rates for the parameters. These are all high, suggesting Step 2 in Algorithm (ref) is effective.
To illustrate the effect of the three significant covariates on the distribution of excess AXP returns, Figure (ref) plots the predictive density (estimated using $\hat{f}^{\mbox{\tiny Bayes}}_{\mbox{\tiny pred}}$) for different values of each covariate, setting the other four covariates equal to their median values. For example, in panel (a) which focuses on variation in MktRf (the excess market return), the distribution is very different in location, spread, and shape for a typical month (MktRf=0.98) in comparison to a poor month ($\mbox{MktRf}=-9.35$) or a strong month ($\mbox{MktRf}=8.42$). This highlights that the regression copula process combined with $G$ defines a distributional regression model, where each covariate affects the entire distribution of $Y$.
While there are only $p=5$ covariates in the example here, the regularization provided by the horseshoe prior allows $\text{\boldmath$x$}_i$ to be of much higher dimension $p$. In particular, KleSmi2019 suggest forming $\text{\boldmath$x$}_i$ using a large number of functional basis terms, such as radial or p-spline bases. This produces a semiparametric distributional regression model that these authors call a “copula smoother”. KleNotSmi2019 instead suggest using the large number of terms from the output layer of a deep neural network (DNN) to form $\text{\boldmath$x$}_i$. The result is a “deep distributional regression” method that combines the flexibility of a DNN with the probabilistic calibration of the copula model.
The implicit copula process described in Section (ref) can also be extended in several directions. klein2020bvs derive the implicit copula for a linear regression with spike-and-slab priors for $\text{\boldmath$\beta$}$. This extends popular Bayesian variable selection methods to a dependent variable with an arbitrary marginal distribution. smith+k19 extend the homoscedastic regression for the auxiliary response $\tilde{Z}_i$ to a heteroscedastic regression. The resulting implicit copula is a mixture of Gaussian copulas, and more flexible than the linear regression copula outlined here. Implicit copula processes constructed from other regression models are also possible.
At (ref) the marginals $F_{Y_j}$ are assumed independent of the covariates, while the copula is not. In contrast, the copula can be assumed to be independent of the covariates, while the marginals are not; for examples, see oakes2000, pitt2006 and song2009. The first approach defines a copula process for a univariate response, whereas the latter approach defines a multivariate regression model for multiple response variables.
The implicit regression copulas outlined in Section (ref) have a dependence structure that is a parametric function of the covariates through the inversion of Sklar's theorem. For a linear regression, this is given by the expression for the correlation matrix $R(\text{\boldmath$x$}_{1:N},\text{\boldmath$\theta$})$ at (ref). An alternative is to make either the parameters or dependence metrics of a copula $C$ smooth functions of the covariates without directly using inversion. Such models are called “conditional copula models”, and there is a extensive literature dealing with this case; for example, see gijbels2011,veraverbeke2011,craiu2012,acar2013,sabeti2014,klein2016 and vatter2018. Another approach is to treat the covariates as a random vector $\bm{X}$ and model it jointly with $Y$ in a copula model, from which the conditional distribution $Y|\bm{X}=\text{\boldmath$x$}$ can be derived. This approach can be easily extended to multivariate responses, as in zhao2019 who employ an elliptical copula with regularization of the parameter space provided by penalization of the coefficients of $\bm{X}$.
While every copula $C$ has one or more implicit representations, it is often infeasible to derive the distribution $F_Z$ of the auxiliary variables. Instead, in this paper we consider implicit copulas to be those derived from a given parametric continuous distribution $F_Z$. Knowledge of $F_Z$ makes estimation of these implicit copulas tractable when using likelihood-based estimation methods, such as Bayesian MCMC or variational inference. This includes high-dimensional cases where $C$ or $c$ cannot be computed in reasonable time, such as the time series and regression copula processes discussed in this paper that have dimension equal to the number of observations.
Another copula family that can be employed in high dimensions are vine copulas joe1996,bedford2002. These are constructed from bivariate copula building blocks called “pair-copulas” by AasCzaFriBak2009. By selecting different pair-copulas, vines can be constructed with a wide range of dependence structures; see czado2019 for an overview. Vines are based on a decomposition into conditional distributions. In some applications an appropriate decomposition arises naturally, such as with time series as in smith2010vine, or when conditioning on latent factors as in krupskii2013,krupskii2020. But, in general, there are many different possibilities morales2010 and it can be difficult to select an appropriate choice, although there have been advances in approaches to do so czado2019. Another challenge in high dimensions is that it can be slow to evaluate the copula density and simulate from the vine, both of which are necessary for parameter estimation and inference. However, truncation as in brechmann2012truncated or other simplifications, such as for stationary Markov time series smith2010vine,smith2015,beare2015, can alleviate these problems. In contrast, it is often unnecessary to evaluate the implicit copula density when estimating the copula model, and simulation from high dimensional implicit copulas is typically fast and stable using Algorithm (ref).
Copulas with discrete and mixed marginals are very increasingly popular genest2007, although parameter estimation is challenging. Several approximate likelihood approaches based on the continuous extension of discrete random variables studied by denuit2005 have been suggested, although these typically exhibit significant bias in the copula parameter estimates; see niko2013JSPI and nikoloulopoulos2016 for demonstrations using the Gaussian copula. In contrast, Bayesian data augmentation approaches discussed here can evaluate the posterior of such parametric copula models in high dimensions without resorting to approximating the likelihood. For implicit copulas, estimation using data augmentation based on the extended likelihood in Section (ref) is popular in practice due to its simplicity and robustness. When using MCMC sampling, as in pitt2006 for the Gaussian copula, the posterior is evaluated exactly (up to Monte Carlo error) and can be used in high dimensions as demonstrated by danaher2011 and dobra2011. An alternative approach that can be used in even higher dimensions is variational inference, as loaiza2019VB outline, although this is an approximate estimation method.
While the Gaussian copula is by far the most popular choice when modeling the dependence of discrete data, it is not clear that it is always the best choice. For example, smith2012 found that a skew $t$ copula provided a substantial improvement over a symmetric $t$ copula for 15-dimensional discrete data. While not explored here, implicit copula processes can also be used for discrete time series data. Doing so provides an alternative to the Markov vine copula models currently popular, as in loaiza2019VB and emura2021.
Finally, this article aims to highlight the potential of time series and regression implicit copula processes. In machine learning they offer a computationally convenient avenue to extend existing deep models to allow for uncertainty quantification. Examples include salinas2019high who do so using Gaussian copula processes and klein2020deep who using the regression copulas in Section (ref) with deep basis functions. The state space copula proposed by smithman2018 and outlined in Section (ref) also has substantial potential. Many existing statistical and econometric models can be written in state space form, from simple time series models to smoothing splines. Combining their implicit copulas with flexible marginals extends these models to more complex data distributions in a straightforward fashion.