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.
85,759 characters · 23 sections · 77 citation commands
myblue Cutting Feedback in Misspecified Copula Models
{0.15cm} {0.15cm} \pagestyle{empty}
\pagestyle{plain} \setcounter{equation}{0}
A copula model specifies a multivariate distribution using its marginal distributions and a copula function to capture the dependence structure nelsen06,joe2014dependence. This simplifies multivariate stochastic modelling, making copula models popular in hydrology genest2007metaelliptical, financial econometrics patton2006, transportation studies bhat2009 and elsewhere. However, exploiting this modularity of copula models to improve the accuracy of statistical inference has been less explored. The goal of the present work is to do so using Bayesian modular inference methods which, to the best of our knowledge, have not been considered previously for copula models. We use a technique called “cutting feedback” liu+bb09,jacob+mhr17, which is applicable to models that comprise multiple components labelled modules. If some modules are misspecified, cutting feedback modifies conventional Bayesian inference to limit the influence of the unreliable modules on the others, producing a “cut posterior” that is more accurate than a conventional posterior for this case. By treating the marginals of a copula model as one module, and the copula function as a second module, we specify two types of cut posterior and develop methods for their evaluation. We establish both theoretically and empirically that the cut posteriors are more accurate than the conventional posterior under the given module misspecification.
Conventional Bayesian inference for copula models using the joint posterior can be unreliable when either the copula function or the marginals are misspecified, and we consider both scenarios. In the first scenario, the goal is to prevent misspecification of the copula fuction with unknown parameters $\text{\boldmath$\psi$}$ from contaminating inference about the marginals with unknown parameters $\text{\boldmath$\theta$}$. We construct a joint cut posterior for $(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$ by cutting feedback from $\text{\boldmath$\psi$}$ to $\text{\boldmath$\theta$}$. The approach is a Bayesian analogue of the IFM method joexu1996,joe2005, but with Bayesian propagation of uncertainty. We prove that this cut posterior asymptotically quantifies uncertainty about $\text{\boldmath$\theta$}$ accurately when the marginal distributions are well-specified, even if the copula function is misspecified, and that the cut posterior mean is first-order equivalent to the IFM estimator.
The second scenario is where the goal is to prevent misspecification of the marginals from contaminating inference about the copula function. This is the more challenging case. To cut feedback from $\text{\boldmath$\theta$}$ to $\text{\boldmath$\psi$}$ we define a novel marginal\footnote{The term “marginal” has two usages here. The first is to refer to the marginal distributions of the copula model, while the second is to Bayesian marginal cut or conventional posterior distributions. Care is taken to clarify between the two throughout.} cut posterior of $\text{\boldmath$\psi$}$ using a pseudo likelihood of the rank data. This is then combined with the conventional conditional posterior for $\text{\boldmath$\theta$}$, given $\text{\boldmath$\psi$}$, to define a joint cut posterior. We prove that this cut posterior asymptotically quantifies uncertainty accurately about the copula parameters if the copula is correctly specified, even if the marginal distributions are misspecified. We are unaware of any existing analog to this approach.
Computation of cut posteriors using Markov chain Monte Carlo (MCMC) or importance sampling methods is difficult because a cut posterior features an unignorable intractable term. Nested MCMC methods are a common approach to circumvent this problem plummer15, but they do not scale well. A major contribution of this paper is the development of efficient and scalable variational methods for the evaluation of the cut posteriors outlined above. Variational inference ormerod2010,blei2017 formulates the approximation of a Bayesian posterior distribution as an optimization problem. It is particularly attractive for evaluating cut posterior distributions because the problematic intractable term does not need to be computed during the optimization yu+ns21,carmona+n22. We show how to implement variational inference for both cut posteriors of the copula model.
In the scenario where feedback from $\text{\boldmath$\theta$}$ to $\text{\boldmath$\psi$}$ is cut, the introduction of the pseudo likelihood of the ranks adds an additional computational bottleneck because it is difficult to evaluate or optimize directly in even moderate dimensions. To solve this problem we introduce an extended likelihood PitChaKoh2006,hoff07,smith2012estimation and then define a cut version of the resulting augmented posterior which is both tractable and has the desired cut posterior as its marginal in $(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$. We develop variational approximations to this augmented cut posterior that are both accurate and allow for fast solution of the variational optimization.
We demonstrate the efficacy of the proposed methodology in the presence of both forms of module misspecification using two simulation studies. These are low dimensional to allow evaluation of the exact cut posteriors, which is difficult otherwise. Here, the cut posteriors provide substantially improved statistical inference in comparison to the conventional posterior, while their variational approximations are also shown to be accurate. The effectiveness of the variational inference methodology is demonstrated in a substantive macroeconomic application where it is infeasible to use exact methods. The example updates the four-dimensional multivariate time series analysis of smithvahey2016 to contemporary data. A copula model is used with four unique skew-t marginals and a Gaussian copula of dimension 1096. This copula arises from a four-dimensional Gaussian copula process with 72 unique parameters observed at 274 time points. The primary objective is density and tail forecasting, and parametric marginals and copula function are necessary to do so, but both are difficult to select. We show that cutting feedback from the marginals to the copula (i.e. from $\text{\boldmath$\theta$}$ to $\text{\boldmath$\psi$}$) improves statistical inference and forecasting accuracy substantially, compared to the conventional posterior.
The rest of the paper is structured as follows. Section (ref) gives some necessary background on copula models, variational inference and cutting feedback. Section (ref) considers cutting feedback when the copula is misspecified, but the marginals are not. Both the theoretical behaviour and computation of the cut posterior are discussed and demonstrated in a simulation study. Section (ref) considers cutting feedback when the marginals are misspecified, but the copula is not. The theoretical behaviour of the cut posterior is established, while an augmented posterior and an appropriate variational approximation is proposed for its computation, and a simulation study demonstrates. Section (ref) contains the macroeconomic application, while Section (ref) concludes.
We begin with a brief outline of copula models and their estimation using the conventional joint posterior. This is followed by an introduction to cutting feedback methods and the computation of cut posteriors for a model with two modules.
If $\bm{Y}=(Y_1,\ldots,Y_m)^\top\sim F_Y$ with marginals $Y_j\sim F_{j}$, then the joint distribution function
where $\text{\boldmath$y$}=(y_1,\ldots,y_m)^\top$ and $C:[0,1]^m \rightarrow \mathbb{R}^+$ is a copula function; see nelsen06. This decomposition provides a convenient modular way to construct a multivariate distribution, where the marginals $F_1,\ldots,F_m$ and copula function $C$ can be selected separately. Parametric marginals $F_j(y_j;\text{\boldmath$\theta$}_j)$ and copula function $C(\text{\boldmath$u$};\text{\boldmath$\psi$})$ are often used,\footnote{The notation $C(\text{\boldmath$u$};\text{\boldmath$\psi$})$ and $C(u_1,\ldots,u_m;\text{\boldmath$\psi$})$ are used interchangably throughout the paper, as are $c(\text{\boldmath$u$};\text{\boldmath$\psi$})$ and $c(u_1,\ldots,u_m;\text{\boldmath$\psi$})$ for the copula density.} with $\text{\boldmath$u$}=(u_1,\ldots,u_m)^\top$ and parameters $\text{\boldmath$\theta$}=(\text{\boldmath$\theta$}_1^\top,\ldots, \text{\boldmath$\theta$}_m^\top)^\top$ and $\text{\boldmath$\psi$}$. The resulting distribution $F_Y$ is commonly called a “copula model” and is employed widely. Many copula functions with different dependence properties have been studied previously; see nelsen06 and joe2014dependence for some examples.
If $F_1,\ldots,F_m$ are all continuous distributions, the joint density of $\bm{Y}$ is
where $c(\text{\boldmath$u$};\text{\boldmath$\psi$})=\frac{\partial}{\partial\text{\boldmath$u$}}C(\text{\boldmath$u$};\text{\boldmath$\psi$})$ is called the copula density, and $f_j(y_j;\text{\boldmath$\theta$}_j)=\frac{\partial}{\partial y_j} F_j(y_j;\text{\boldmath$\theta$}_j)$ is the marginal density of $Y_j$. When one or more $F_j$ is discrete or mixed, the joint mixed density function involves differencing over those dimensions; see genest2007.
Let ${\cal D}=\{\text{\boldmath$y$}_1,\ldots,\text{\boldmath$y$}_n\}$ be $n$ observations drawn independently from $F_Y$ at (ref). For continuous marginals, the joint posterior density is $p(\text{\boldmath$\theta$},\text{\boldmath$\psi$}|{\cal D})\propto \prod_{i=1}^n f_Y(\text{\boldmath$y$}_i|\text{\boldmath$\theta$},\text{\boldmath$\psi$})p(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$, where $p(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$ is the prior. Evaluation of the posterior using Markov chain Monte Carlo (MCMC) methods has been discussed previously by PitChaKoh2006,silva2008copula,min2010,smithmin2010 and murray2013 among others. However, evaluation of the posterior using MCMC methods can be slow for large $m$, and variational inference (VI) is a faster and more scalable alternative.
MCMC methods evaluate the posterior exactly (up to a controllable level of Monte Carlo error), whereas VI approximates the posterior by a density chosen from a family of tractable distributions with densities $q\in\mathcal F$. The density is chosen to minimize the distance between the two, with the Kullback-Leibler (KL) divergence the most commonly used measure, so that for a copula model
Many families $\mathcal F$ have been considered in the literature, but a Gaussian with density $q_\lambda(\text{\boldmath$x$})=\phi_N(\text{\boldmath$x$};\text{\boldmath$\mu$},\Sigma)$ indexed by its unique parameters $\text{\boldmath$\lambda$}=(\text{\boldmath$\mu$}^\top,\mbox{vech}(\Sigma)^\top)^\top$ is one of the most popular titsias2014doubly,kucukelbir2017automatic,tan2018gaussian. It is straightforward to show (e.g. see ormerod2010) that $q^*(\text{\boldmath$\theta$}, \text{\boldmath$\psi$}) = \operatorname*{{arg\,max}}_{q \in \mathcal F}\mathcal L(\text{\boldmath$\lambda$})$, where the function \[ \mathcal L(\text{\boldmath$\lambda$})=E_q\left(\log h(\text{\boldmath$\theta$},\text{\boldmath$\psi$})-\log q_\lambda(\text{\boldmath$\theta$},\text{\boldmath$\psi$})\right)\,, \] is called the Evidence Lower Bound (ELBO) and $h(\text{\boldmath$\theta$},\text{\boldmath$\psi$})=p(\mathcal D|\text{\boldmath$\theta$},\text{\boldmath$\psi$})p(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$.
A popular way to solve this problem is to use stochastic gradient optimization bottou10. This employs an unbiased approximation of the gradient $\nabla_\lambda \mathcal L(\text{\boldmath$\lambda$})$ along with automatic adaptive step sizes for the updates of $\text{\boldmath$\lambda$}$, such as the ADADELTA method of zeiler12 that we use here. The combination of stochastic optimization and generic approximations is often called black box VI ranganath14,titsias2014doubly. VI has been used to estimate copula models by loaiza2019VBDA, nguyen2020VI and smithklein2021.
Cutting feedback is a form of Bayesian modular inference liu+bb09 that removes the impact of mis-specifying one or more model components on inference for the other components. Comprehensive overviews of cutting feedback methods are provided by lunn+bsgn09, plummer15, jacob+mhr17 and yu+ns21. A short introduction is given here for a two module system because the methods developed later for copula models are two module systems.
Consider a model for data ${\cal D}$ with density $g({\cal D}|\bm{\eta})$ and parameter vector $\bm{\eta}$. Consider the partition $\bm{\eta}=(\text{\boldmath$\eta$}_1^\top,\text{\boldmath$\eta$}_2^\top)^\top$ and assume the density can be factorized as
Often in a two module system the data consists of two sources ${\cal D}_1$ and ${\cal D}_2$, with $g_1({\cal D}|\text{\boldmath$\eta$}_1)=g_1({\cal D}_1|\text{\boldmath$\eta$}_1)$ and $g_2({\cal D}|\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2)=g_2({\cal D}_2|\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2)$; e.g. see plummer15. However, we do not assume this simplification here because a more general perspective, where $g_1$ and $g_2$ represent different terms in a decomposition of the likelihood, is needed for the cut methods for copulas.
Denoting the prior density as $p(\text{\boldmath$\eta$})=p(\text{\boldmath$\eta$}_1)p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)$, we define two “modules”, with Module 1 consisting of $g_1({\cal D}|\text{\boldmath$\eta$}_1)$ and $p(\text{\boldmath$\eta$}_1)$, and Module 2 consisting of $g_2({\cal D}|\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2$) and $p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)$. The conventional joint posterior density is $$p(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2|{\cal D})= p(\text{\boldmath$\eta$}_1|{\cal D})\times p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1,{\cal D}).$$ Writing $\bar{g}({\cal D})=\int p(\text{\boldmath$\eta$}_1)p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1) g({\cal D}|\bm{\eta})d\bm{\eta},$ and $\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1)=\int p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)g_2({\cal D}|\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2)d\text{\boldmath$\eta$}_2,$ a simple derivation shows that the marginal posterior density is
and the conditional posterior density is
In (ref), $\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1)$ is called the “feedback” term, because it captures the effect of Module 2 on inference for $\text{\boldmath$\eta$}_1$. If Module 2 is misspecified the influence of the feedback term can result in misleading marginal inference for $\text{\boldmath$\eta$}_1$. Hence in joint Bayesian inference, even if Module 1 is correctly specified, misspecification of Module 2 can result in misleading inference about parameters appearing in both modules.
To eliminate the impact of a misspecification of Module 2 on inference for $\text{\boldmath$\eta$}_1$, the feedback term can be removed from (ref) to define the following marginal cut posterior density
The joint cut posterior density is then defined as
where $p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1,{\cal D})$ is the same conditional posterior for the cut and uncut cases. A key observation is that uncertainty about $\text{\boldmath$\eta$}_1$ is still propagated when computing marginal cut posterior inference for $\text{\boldmath$\eta$}_2$ with \[p_{\text{cut}}(\text{\boldmath$\eta$}_2|{\cal D})=\int p_{\text{cut}}(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2|{\cal D}) d\text{\boldmath$\eta$}_1\,. \]
Cut posterior computation is difficult. The joint cut posterior density is $$p_{\text{cut}}(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2|{\cal D}) \propto \frac{p(\text{\boldmath$\eta$}_1)p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)g_1({\cal D}|\text{\boldmath$\eta$}_1)g_2({\cal D}|\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2)}{\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1)},$$ where $\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1)$ is usually intractable. This makes it hard to implement MCMC or importance sampling methods to evaluate the cut posterior in many models. One approach is to draw samples from (ref) by first drawing $\text{\boldmath$\eta$}_1'\sim p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})$, and then $\text{\boldmath$\eta$}_2'|\text{\boldmath$\eta$}_1'\sim p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1',{\cal D})$. Because $\text{\boldmath$\eta$}_1'$ is fixed in the second stage, the intractable term $\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1')$ is not computed. However, direct generation from these distributions is often difficult, and plummer15 suggested using “nested MCMC” as in Algorithm (ref) below. Other methods for cut posterior evaluation are discussed by liu+g20, jacob2020unbiased and pompe+j21.
Given the difficulty of exact cut posterior computation, variational inference methods to do so have been suggested by yu+ns21 and carmona+n22. Lemma 1 of yu+ns21 establishes that the cut posterior distribution is closest in Kullback-Leibler divergence to the true posterior amongst distributions that have $\text{\boldmath$\eta$}_1$ marginal density $p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})$. Therefore, if the family of approximations ${\cal F}$ is restricted to those that have $\text{\boldmath$\eta$}_1$ marginal density $p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})$, solving the conventional variational optimization problem at (ref) will also provide the optimal variational approximation to the cut posterior. Crucially, solving this optimization does not require computation of the intractable term $\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1)$ which creates the computational bottleneck in MCMC.
This observation motivates a sequential VI procedure suggested by yu+ns21. In a first stage an approximation of $p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})$ is computed, which is then kept fixed in a second stage. Consider a family of densities of the form $q_{\lambda}(\text{\boldmath$\eta$})=q_{\widetilde{\lambda}}(\text{\boldmath$\eta$}_1)q_{\breve{\lambda}}(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)$, where $\bm{\lambda}=(\widetilde{\bm{\lambda}}^\top,\breve{\bm{\lambda}}^\top)^\top$ are variational parameters partitioned into two sets. The first set $\widetilde{\bm{\lambda}}$ parametrize the $\text{\boldmath$\eta$}_1$ marginal density, and the second set $\breve{\bm{\lambda}}$ parametrize the conditional density for $\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1$. If $D_{KL}\left(q\,||\,p\right)$ denotes the KL divergence of $q$ from $p$, then Algorithm (ref) below outputs an approximation to the joint cut posterior.
In our empirical work, Gaussian variational approximations are used with Gaussian density $q_\lambda(\bm{\eta})=\phi_N(\bm{\eta};\bm{\mu},\Sigma)$, with mean $\bm{\mu}$ and variance $\Sigma=LL^\top$, where is $L$ a lower triangular Cholesky factor. Partitioning $\text{\boldmath$\mu$}$ and $L$ to be conformable with $\bm{\eta}=(\text{\boldmath$\eta$}_1^\top,\text{\boldmath$\eta$}_2^\top)^\top$, so that $$\bm{\mu}=\left[
\right], \;\;\;\; L=\left[
\right],$$ then $q_{\widetilde{\lambda}}(\boldmath$\eta$_1)$ and $q_{\breve{\lambda}}(\boldmath$\eta$_2|\boldmath$\eta$_1)$ are both Gaussian densities with parameters $\widetilde{\boldmath$\lambda$}=(\bm{\mu}_{\eta_1}^\top, vech(L_{\eta_1})^\top)^\top$, and $\breve{\boldmath$\lambda$}=(\bm{\mu}_{\eta_2}^\top,\text{vec}(L_{\eta_1,\eta_2})^\top,\text{vech}(L_{\eta_2})^\top)^\top$, where `vec' and `vech' are the vectorization and half-vectorization matrix operators, respectively. Methods for optimizing a Gaussian variational density parametrized by a Cholesky factor are well-known in the literature (e.g. titsias2014doubly, among many others) and we do not describe this in detail here. Other fixed form approximating families can also be used in this framework.
A copula model can be viewed as a two module system, where the marginals $F_1(\cdot;\text{\boldmath$\theta$}_1),\ldots,F_m(\cdot;\text{\boldmath$\theta$}_m)$ form one module, and the copula $C(\cdot;\text{\boldmath$\psi$})$ is a second module. This section discusses cutting feedback when the marginals are thought to be adequate, but the copula function may be misspecified. We label the cut posterior for this case “type 1” to distinguish it from that in Section (ref).
For the copula model with density at (ref), we set $\text{\boldmath$\eta$}_1=\text{\boldmath$\theta$}$ and $\text{\boldmath$\eta$}_2=\text{\boldmath$\psi$}$ and factor the likelihood as $g({\cal D}|\bm{\theta},\bm{\psi})=g_1({\cal D}|\bm{\theta})g_2({\cal D}|\bm{\theta},\bm{\psi})$, where $$g_1({\cal D}|\bm{\theta})=\prod_{i=1}^n \prod_{j=1}^m f_j(y_{ij};\bm{\theta}_j),\;\; \mbox{ and }\;\; g_2({\cal D}|\bm{\theta},\bm{\psi})=\prod_{i=1}^n c(F_1(y_{i1};\bm{\theta}_1),\dots, F_m(y_{im};\bm{\theta}_m);\bm{\psi}).$$ Assuming prior density $p(\bm{\theta},\bm{\psi})=p(\bm{\theta})p(\bm{\psi})$, with $p(\bm{\theta})=\prod_{j=1}^m p(\bm{\theta}_j)$, then the marginal cut posterior at (ref) simplifies to $$p_{\text{cut}}(\bm{\theta}|{\cal D})=\prod_{j=1}^m p_j(\bm{\theta}_j|\bm{y}_{(j)}),$$ where $\bm{y}_{(j)}=(y_{1j},\dots, y_{nj})^\top$ denotes the data for the $j$th marginal and $p_j(\bm{\theta}_j|\bm{y}_{(j)})\propto p(\bm{\theta}_j)\prod_{i=1}^n f_j(y_{ij};\bm{\theta}_j).$
The ordinary and cut conditional posterior density is $$p(\bm{\psi}|\bm{\theta},{\cal D})= \frac{p(\bm{\psi})\prod_{i=1}^n c(F_1(y_{i1};\bm{\theta}_1),\dots, F_m(y_{im};\bm{\theta}_m);\bm{\psi})}{\bar{g}_2({\cal D}|\bm{\theta})},$$ where $\bar{g}_2({\cal D}|\bm{\theta})= \int p(\bm{\psi})\prod_{i=1}^n c(F_1(y_{i1};\bm{\theta}_1),\dots, F_m(y_{im};\bm{\theta}_m);\bm{\psi}) \,d\bm{\psi}$. The joint cut posterior is
and we consider its computation using both the nested MCMC and variational approaches in Algorithms (ref) and (ref).
Among the most popular methods for estimating copula models is the “inference for margins”(IFM) procedure of joexu1996 and joe2005. In IFM, each $\bm{\theta}_j$ is estimated by maximizing the likelihood of the $j$th marginal model, and then the copula parameters $\bm{\psi}$ are estimated by maximizing the likelihood conditional on these estimates. We now show for large $n$ the cut posterior at (ref) resembles a Bayesian version of IFM.
We first establish that the posterior mean of $p_{\text{cut}}(\bm{\theta},\bm{\psi}|{\cal D})$, denoted as $\bar\bm{\eta}:=\int \text{\boldmath$\eta$} p_{\text{cut}}(\bm{\eta}|{\cal D}) d\bm{\eta}$, is asymptotically equivalent to the IFM point estimator. To this end, define $\widehat{\bm{\theta}}$ as the IFM estimator obtained by first maximizing $\log g_1(\mathcal{D}|\bm{\theta})$, define $\widehat{\bm{\psi}}$ as the IFM estimator obtained by maximizing $\log g_2(\mathcal{D}\mid \bm{\psi},\widehat{\bm{\theta}})$ over $\bm{\psi}$, and set $\widehat{\bm{\eta}}=(\widehat{\bm{\theta}}^\top,\widehat{\bm{\psi}}^\top)^\top$. Lemma (ref) below establishes that the IFM point estimator and the cut posterior mean are asymptotically equivalent.
Assumptions (ref) and (ref) are similar to the standard regularity conditions employed in two-step copula modeling to deduce asymptotic normality of the IFM point estimator in joe2005; see Part (ref) of the Web Appendix for their specification and a detailed discussion.
Lemma (ref) does not address the accuracy with which the cut posterior quantifies uncertainty. To establish this we require the following additional definitions and observations. Let $P_0$ denote the true data generating process (DGP) for the observed data, and $p_0$ its density, then under Assumptions (ref) and (ref) in Part (ref) of the Web Appendix, it can be shown that both $\bar{\bm{\theta}}$ and $\widehat{\bm{\theta}}$ are consistent estimators of $$\bm{\theta}_0=\operatorname*{{arg\,min}}_{\bm{\theta}}D_{\text{KL}}\left(p_0\,||\,g_1(\cdot\mid\bm{\theta})\right)\,. $$ That is, $g_1(\cdot\mid\bm{\theta}_0)$ is the closest element of the class $\{\bm{\theta}: g_1(\cdot\mid\bm{\theta})\}$ to $P_0$ in terms of KL divergence, and $\bm{\theta}_0$ is the corresponding pseudo-true value. Further, define the following matrix of second derivatives for the marginal model parameters (i.e., $\bm{\theta}$) :
where $\operatorname{E}$ is the expectation with respect to $P_0$. Then Lemma (ref) below shows how the cut posterior for $\bm{\theta}$ quantifies uncertainty.
This result shows that asymptotically the cut posterior for $\text{\boldmath$\theta$}$ resembles a Gaussian distribution centred at the IFM $\widehat\bm{\theta}$, and with variance $\mathcal{I}^{-1}/n$. Therefore, when the marginals are correctly specified, the type 1 cut posterior for $\bm{\theta}$ correctly quantifies uncertainty\footnote{By this we mean that a level $(1-\alpha)$ credible set asymptotically has frequentist coverage at the $(1-\alpha)$ level under $P_0$; i.e., Bayesian credible sets agree asymptotically with frequentist confidence sets.} for the unknown parameter value $\bm{\theta}_0$, even if the copula function is misspecified. That is, $p_{\mathrm{cut}}(\bm{\theta}|\mathcal{D})$ delivers inferences that are asymptotically the same as IFM and also correctly quantifies uncertainty.
To understand how the marginal cut posterior for $\bm{\psi}$, $$p_\mathrm{cut}(\bm{\psi}|\mathcal{D})=\int p_{\text{cut}}(\bm{\theta},\text{\boldmath$\psi$}|{\cal D})d\bm{\theta}=\int p_{\text{cut}}(\bm{\theta}|{\cal D})p(\bm{\psi}|\bm{\theta},{\cal D})d\bm{\theta},$$ quantifies uncertainty, define $$\bm{\psi}_0=\operatorname*{{arg\,min}}_{\bm{\psi}}D_{\text{KL}}\left({p^{}_0}\,||\,g_2(\cdot\mid\bm{\psi},\bm{\theta}_0)\right), $$which is the pseudo-true value for the copula parameters $\bm{\psi}$ when the unknown $\bm{\theta}$ is replaced by $\bm{\theta}_0$, and define the following matrix of second derivatives:
with $\mathcal{M}_{\bm{\theta}\bm{\theta}}$ and $\mathcal{M}_{\bm{\psi}\bm{\psi}}$ defined analogously; then the following lemma holds.
Lemma (ref) demonstrates that the variability for the marginal cut posterior of $\bm{\psi}$ depends on the variability of the cut posterior for $\bm{\theta}$. Hence, uncertainty flows from $p_\mathrm{cut}(\bm{\theta}|\mathcal{D})$ to $p_\mathrm{cut}(\bm{\psi}|\mathcal{D})$, but not the other way. This implies that if the cut posterior for $\bm{\psi}$ is to correctly quantify uncertainty, then both the marginal components and the copula function must be well-specified.
A simulation study compares the accuracy of the type 1 cut posterior to that of the conventional (i.e. uncut) posterior and IFM. For each sample size $n \in \{100, 500, 1000\}$ a total of $S=500$ datasets are generated from a bivariate copula model. The marginal $f_1$ is a log-normal distribution, with mean and variance parameters $\mu=1$ and $\sigma^2=1$, while the marginal $f_2$ is a gamma distribution with shape and rate parameters $\alpha=7$ and $\beta=3$. A t-copula demarta2005 is used with Kendall's tau $\tau=0.7$ and unity degrees of freedom parameter.
For each dataset, we fit a copula model with the correct marginal distributional forms, along with a bivariate Gumbel copula. Thus, the copula is misspecified, but can still capture correlation, as measured by Kendall's tau $\tau$, equal to that of the DGP. We assign vague proper priors $\mu \sim N(0,100^2)$, $\sigma^2 \sim \text{Half-Normal}(0,100^2)$, $\alpha \sim \text{Half-Cauchy}(0,5)$, $\beta \sim \text{Half-Cauchy}(0,5)$, and $\tau \sim \text{Uniform}(0,1)$, where $\text{Half-Cauchy}(m,s)$ is a half Cauchy distribution with location $m$ and scale $s$. Both the outlined variational methodology and MCMC algorithms are used, with details given in Part A1 of the Web Appendix. This results in four Bayesian posteriors, the means of which are used as point estimators. IFM is also used for comparison.
To measure estimation accuracy of the true parameter values in the DGP, the bias and root mean square error (RMSE) is evaluated over the $S$ replicates. Table (ref) (left-hand side) reports these for the case where $n = 1000$, and we make four observations. First, the cut posterior has lower bias and RMSE than the conventional posterior for all parameters, so that cutting feedback from the misspecified copula improves estimation accuracy. Second, $\tau$ is estimated more accurately using its cut posterior than IFM. Third, the variational and exact posterior results are similar, suggesting the former is an accurate approximation. (Although, if a Gaussian VA underestimates uncertainty for other target posteriors, then a richer variational family can also be used.) Last, IFM provides a more accurate estimate of $\sigma^2$, but is less accurate than the cut posterior for all other parameters.
We also consider accuracy when the correctly specified copula model (i.e. a t-copula with the correct degrees of freedom) is fit. Table (ref) (right-hand side) reports the bias and RMSE for both the cut and conventional posterior in this case, both estimated exactly using MCMC and the same uniform prior on $-1<\tau<1$. The cut posterior is only slightly less accurate than the conventional (i.e. uncut) posterior.
To assess the accuracy of the marginal posterior distributions we compute the coverage of their 95% credible intervals. To assess the accuracy of the point estimates of the copula model components (in addition to their parameter values) we compute their predictive KL divergences. The latter is defined for marginal $j=1,2$ as $$ \text{KL}_j = \int f_j (y ; \widehat{\text{\boldmath$\theta$}}_j) \left [ \log f_j (y ; \widehat{\text{\boldmath$\theta$}}_j) - \log f_j^\star (y) \right ] \; \text{d} y, $$ where $\widehat{\text{\boldmath$\theta$}}_j$ is a point estimate of $\text{\boldmath$\theta$}_j$, and $f_j^\star$ is the true marginal density of the DGP. The predictive KL divergence for the copula is defined as $$ \text{KL}_{cop} = \int \int c (u, v ; \widehat{\tau}) \left [ \log c (u, v ; \widehat{\tau}) - \log c^\star (u, v) \right ] \; \text{d} u \text{d} v\,, $$ where $\widehat{\tau}$ is a point estimate of $\tau$, and $c^\star$ is the true copula density for the DGP. The integrals above are computed numerically.
Table (ref) reports the coverage probability of the credible intervals, along with the mean of the KL divergence metrics over the $S$ replicates. The results further confirm that under misspecification of the copula function (left-hand side of the table), the cut posterior is substantially more accurate than the conventional (uncut) posterior, and that the variational and exact posteriors are very similar. Despite IFM estimating $\sigma^2$ slightly more accurately than the cut posterior, the cut posterior either equals or out-performs estimation accuracy of all model components as measured by the KL divergences. Again, we see that when the correct model is fit (right-hand side of the table), the cut posterior is only slightly less accurate than the conventional posterior.
Results for the cases where $n = 100$ and $n=500$ are reported in Part A4 of the Web Appendix, and are very similar to those for $n=1000$.
This section discusses cutting feedback when the copula function $C(\cdot;\text{\boldmath$\psi$})$ is adequate, but the marginals $F_1(\cdot;\text{\boldmath$\theta$}_1),\ldots,F_m(\cdot;\text{\boldmath$\theta$}_m)$ are misspecified. We label the cut posterior for this case “type 2” and use a pseudo likelihood of the rank data for its specification. Evaluation of this cut posterior is more challenging than that in Section (ref), and to do so in higher dimensions we introduce an extension of this pseudo likelihood PitChaKoh2006,hoff07,smith2012estimation and then define a cut version of the resulting augmented posterior which is both tractable and has the desired type 2 cut posterior as its marginal.
Setting $\text{\boldmath$\eta$}_1=\text{\boldmath$\psi$}$ and $\text{\boldmath$\eta$}_2=\text{\boldmath$\theta$}$, to define the marginal cut posterior for $\text{\boldmath$\psi$}$ we use a pseudo likelihood based on the rank data. For each $y_{ij}$ define its rank within marginal $j$ as $r(y_{ij})$ \footnote{For example, in the absence of ties this is $r(y_{ij})=\sum_{k=1}^n \mathds{1}(y_{kj}\leq y_{ij})$.} and denote all the rank data as $r({\cal D})=\{r(y_{ij});i=1,\ldots,n, j=1,\ldots,m\}$. We employ the following probability mass function for the (discrete-valued) ranks
where $a_{ij}=(r(y_{ij})-1)/(n+1)$, $b_{ij}=r(y_{ij})/(n+1)$, $\text{\boldmath$v$}=(v_1,\ldots,v_m)^\top$\,, and where
is a differencing operator over element $j$ nelsen06. This is the likelihood under the assumption that each marginal is an empirical distribution function. It is related to the “rank likelihood” that is obtained from the exact distribution of the ranks; for example, see hoff07 for specification of the rank likelihood of a Gaussian copula. However, as we discuss later, it is more tractable than a rank likelihood. It is also related to the popular pseudo-likelihood in genest95, but corrects for the discrete nature of the rank data.
The pseudo rank likelihood at (ref) does not depend on the marginal parameters $\bm{\theta}$ because the ranks are a strictly increasing transformation of $\mathcal{D}$ and are unaffected by the marginal distributions. Therefore it can be used to define a marginal cut posterior for $\text{\boldmath$\psi$}$ with density
This definition fits into the two module system described in Section (ref) by considering the factorization at (ref) with $g_1({\cal D}|\text{\boldmath$\psi$})=p_{\text{PL}}(r(D)|\text{\boldmath$\psi$})$ and $g_2({\cal D}|\text{\boldmath$\psi$},\text{\boldmath$\theta$})=p({\cal D}|\text{\boldmath$\psi$},\text{\boldmath$\theta$})/p_{\text{PL}}(r(D)|\text{\boldmath$\psi$})$. With these definitions, the feedback term is \[ \bar{g}_2({\cal D}|\text{\boldmath$\psi$})=\int g_2({\cal D}|\text{\boldmath$\psi$},\text{\boldmath$\theta$})p(\text{\boldmath$\theta$})d\text{\boldmath$\theta$}= \int \frac{p({\cal D}|\text{\boldmath$\psi$},\text{\boldmath$\theta$})}{p_{\text{PL}}(r(D)|\text{\boldmath$\psi$})}p(\text{\boldmath$\theta$})d\text{\boldmath$\theta$}\,. \] In the two module system, the cut posterior at (ref) is obtained by removing this feedback term. Notice that if the likelihood $p({\cal D}|\text{\boldmath$\psi$},\text{\boldmath$\theta$})$ is close to the pseudo rank likelihood $p_{\text{PL}}(r(D)|\text{\boldmath$\psi$})$, then $\bar{g}_2({\cal D}|\text{\boldmath$\psi$})\approx 1$ and the cut and ordinary posteriors for $\text{\boldmath$\psi$}$ will also be close. Conversely, if the likelihood and the pseudo rank likelihood deviate, the cut and ordinary posteriors will differ.
The joint cut posterior is defined as
where the conditional $p(\text{\boldmath$\theta$}|\text{\boldmath$\psi$},{\cal D})=p({\cal D}|\text{\boldmath$\theta$},\text{\boldmath$\psi$})p(\text{\boldmath$\theta$})/\int p({\cal D}|\text{\boldmath$\theta$}',\text{\boldmath$\psi$}) p(\text{\boldmath$\theta$}') d\text{\boldmath$\theta$}'$. The normalizing constant of this conditional is not computed when implementing Algorithm (ref).
An advantage of the pseudo rank likelihood at (ref) is that it is both computationally and theoretically more tractable than the rank likelihood of hoff07 and others. As hoff2014information state, the rank likelihood “is the integral of a copula density over a complicated set defined by multivariate order constraints”, making it intractable and complicating the derivation of theoretical results for parameter inference. For example, hoff2014information control an accurate approximation of the rank likelihood in order to deduce their theoretical results.
In contrast, the type 2 cut posterior depends on (ref) and the parametric likelihood $p(\mathcal{D}|\bm{\theta},\bm{\psi})$, both of which are tractable. This allows direct analysis of the behavior of the cut posteriors $p_{\text{cut}}(\text{\boldmath$\psi$}|{\cal D})$ and $p_{\text{cut}}(\text{\boldmath$\psi$},\bm{\theta}|{\cal D})$ in (ref) and (ref). To this end, let $M_n(\bm{\psi}):=\log p_{\text{PL}}(r(\mathcal{D})|\bm{\psi})$, with $\mathcal{M}(\bm{\psi})=\lim_{n \rightarrow \infty}M_n(\bm{\psi})/(1+n)$; further define $\widehat\bm{\psi}_r=\operatorname*{{arg\,max}}_{\bm{\psi}} M_n(\bm{\psi})$, $\bm{\psi}_\star=\operatorname*{{arg\,max}}_{\bm{\psi}}\mathcal{M}(\bm{\psi})$, and $\mathcal{M}_{\bm{\psi}\bm{\psi}}(\bm{\psi}_\star)=-\nabla_{\bm{\psi}\bm{\psi}}^2\mathcal{M}(\bm{\psi}_\star)$. Theorem (ref) below characterizes the behavior of $p_{\mathrm{cut}}(\bm{\psi}|\mathcal{D})$.
Assumptions (ref)-(ref) are given in Part (ref) of the Web Appendix, and they ensure that $M_n(\bm{\psi})$ admits enough regularity so that the cut marginal posterior $p_{\mathrm{cut}}(\bm{\psi}|\mathcal{D})$ satisfies a Berstein-von Mises result. While these assumptions are specific to the copula function, they are satisfied for popular choices, including elliptical copulas, such as the student-t and Gaussian copulas, and key Archimedean copulas such as the Gumbel and Clayton copulas; see Part (ref) of the Web Appendix for further discussion.
Theorem (ref) implies that in large samples the cut posterior $p_{\mathrm{cut}}(\bm{\psi}|\mathcal{D})$ based on the pseudo rank likelihood resembles a Gaussian density centered at $\widehat{\bm{\psi}}_r$. If the copula is correctly specified, then the information matrix equality is satisfied and we have that $\mathcal{M}_{\bm{\psi}\bm{\psi}}(\bm{\psi}_\star)^{-1}\equiv \mathrm{var}\left\{\nabla_{\bm{\psi}} M_n(\bm{\psi}_\star)/\sqrt{n}\right\}$, which has a particular form given in Corollary (ref) in Web Appendix (ref). In such cases, Theorem (ref) implies that the cut posterior based on the pseudo rank likelihood correctly quantifies uncertainty.
To state the behavior of $p_{\mathrm{cut}}(\bm{\theta}|\mathcal{D})=\int p_{\text{cut}}(\text{\boldmath$\psi$},\text{\boldmath$\theta$}|{\cal D})d\bm{\psi}$, let ${Q}_n(\bm{\theta},\bm{\psi})= \log p(\mathcal{D}|\bm{\theta},\bm{\psi})$, and write $\mathcal{Q}(\bm{\theta},\bm{\psi}):=\lim_{n\rightarrow \infty}n^{-1}\operatorname{E}\left(\log p(\mathcal{D}|\bm{\theta},\bm{\psi}) \right)$, with derivatives of $\mathcal{Q}(\bm{\theta},\bm{\psi})$ denoted as $\mathcal{Q}_{ij}(\bm{\theta},\bm{\psi})=\nabla^2_{ij}\mathcal{Q}(\bm{\theta},\bm{\psi})$ for $i,j\in\{\bm{\theta},\bm{\psi}\}$. Further define $\widehat\bm{\theta}_r:=\operatorname*{{arg\,max}}_{\bm{\theta}} Q_n(\bm{\theta},\widehat\bm{\psi}_r)$, $\bm{\theta}_\star:=\operatorname*{{arg\,max}}_{\bm{\theta}}\mathcal{Q}(\bm{\theta},\bm{\psi}_\star)$, $$ \Omega^{-1}=\mathcal{Q}_{\bm{\theta}\bm{\theta}}(\bm{\eta}_\star)^{-1}+\mathcal{Q}_{\bm{\theta}\bm{\theta}}(\bm{\eta}_\star)^{-1}\mathcal{Q}_{\bm{\theta}\bm{\psi}}(\bm{\eta}_\star)\mathcal{M}_{\bm{\psi}\bm{\psi}}(\bm{\psi}_\star)^{-1}\mathcal{Q}_{\bm{\psi}\bm{\theta}}(\bm{\eta}_\star)\mathcal{Q}_{\bm{\theta}\bm{\theta}}(\bm{\eta}_\star)^{-1}\,, $$ and $\bm{\eta}_\star=(\bm{\theta}_\star^\top,\bm{\psi}_\star^\top)^\top$. Then Theorem (ref) below characterizes the behavior of $p_\mathrm{cut}(\bm{\theta}|\mathcal{D})$.
Theorem (ref) shows that the uncertainty for the cut posterior of $\bm{\theta}$ depends on the uncertainty in the cut posterior for $\bm{\psi}$ through the term $\mathcal{M}_{\bm{\psi}\bm{\psi}}(\bm{\psi}_\star)^{-1}$. Therefore, the cut posterior of $\bm{\theta}$ will only quantify uncertainty correctly if the copula model marginals and copula function are both well-specified. This is in contrast to the type 1 cut posterior where the marginal parameter posteriors delivered reliable uncertainty quantification as long as the marginal models were well-specified (see Lemma (ref)) .
The simulation study in Section (ref) is extended to compare the accuracy of the type 2 cut posterior to that of the conventional posterior. Data is generated from a bivariate copula model with the similar marginals as in Simulation 1 (except that $\sigma^2 = 0.25$), but using a Gumbel copula with Kendall's tau $\tau=0.7$. For each dataset we fit a copula model with the correct copula family (i.e. a Gumbel), along with normal marginals with mean and variance parameters $\mu_j,\sigma^2_j$ for $j=1,2$ and constrained to be positive. Thus, the marginals are misspecified but the distribution has the same support as the DGP. We employ the vague proper priors $\mu_j, \sim N(0,100^2)$, $\sigma_j^2 \sim \text{Half-Normal}(0,100^2)$, and $\tau \sim U(0,1)$. Both the outlined variational methodology and MCMC algorithms are used to evaluate the type 2 cut posterior, along with the conventional posteriors, resulting in four Bayesian estimators. Details are given in Part A1 of the Web Appendix.
The accuracy of each posterior is measured using the predictive KL divergence metrics. Table (ref) (left hand side) reports their mean values over the $S=500$ replicates for the case where $n=1000$. The cut posterior provides much more accurate estimates of both the marginal and copula components, compared to the conventional posteriors. Moreover, MCMC and variational estimates provide very similar levels of accuracy. Table (ref) (right hand side) reports the accuracy when the correctly specified copula model (i.e. with the correct forms for the marginals) is fit using both the conventional and cut posteriors computed using MCMC. The same vague proper priors are used for the marginal parameters, and the accuracy of the cut posterior is almost identical to the conventional posterior. Results for the cases where $n=100$ and $n=500$ are reported in Part A5 of the Web Appendix, and are very similar to those for $n=1000$.
When $m$ is small, the cut posterior can be evaluated by direct application of Algorithm (ref). However, for even moderate values of $m$, the pseudo rank likelihood at (ref) forms a computational bottleneck because it requires $O(n2^m)$ evaluations of $C$. In this case the computation can be avoided by employing the extended likelihood in smith2012estimation which is tractable for higher values of $m$.
Let $\text{\boldmath$u$}_i=(u_{i1},\ldots,u_{im})^\top\sim C(\cdot;\text{\boldmath$\psi$})$ and $\text{\boldmath$u$}=(\text{\boldmath$u$}_1^\top,\ldots,\text{\boldmath$u$}_n^\top)^\top$ be auxiliary variables, such that $p(r({\cal D})|\text{\boldmath$u$})=\prod_{ij}p(r(y_{ij})|u_{ij})=\prod_{ij}\mathds{1}(a_{ij}\leq u_{ij}<b_{ij})$. Then define an extended likelihood as \[ p(r({\cal D}),\text{\boldmath$u$}|\text{\boldmath$\psi$}):=p(r({\cal D})|\text{\boldmath$u$})p(\text{\boldmath$u$}|\text{\boldmath$\psi$})= \prod_{ij}\mathds{1}(a_{ij}\leq u_{ij}<b_{ij})\prod_{i=1}^n c(\text{\boldmath$u$}_i|\text{\boldmath$\psi$})\,. \] Theorem 1 in smith2012estimation shows that integrating over $\text{\boldmath$u$}$ retrieves the pseudo rank likelihood; i.e. $p_{\text{PL}}(r({\cal D})|\text{\boldmath$\psi$})=\int p(r({\cal D}),\text{\boldmath$u$}|\text{\boldmath$\psi$})d\text{\boldmath$u$}$. Using this extended likelihood, we define the marginal cut posterior of $\text{\boldmath$\psi$}$ augmented with $\text{\boldmath$u$}$ as \[ p_{\text{cut}}(\text{\boldmath$\psi$},\text{\boldmath$u$}|{\cal D})=\frac{p(r(D),\text{\boldmath$u$}|\text{\boldmath$\psi$})p(\text{\boldmath$\psi$})}{\int p_{\text PL}(r(D)|\text{\boldmath$\psi$}')p(\text{\boldmath$\psi$}')d\text{\boldmath$\psi$}'}\,. \] Integrating the density above over $\text{\boldmath$u$}$ gives the required cut posterior at (ref).
Again, this setup fits into the two module system discussed in Section (ref), but with $\text{\boldmath$\eta$}_1=(\text{\boldmath$\psi$}^\top,\text{\boldmath$u$}^\top)^\top$ and $\text{\boldmath$\eta$}_2=\text{\boldmath$\theta$}$, so that the cut posterior
which we call the “augmented cut posterior” (i.e. the joint cut posterior augmented with $\text{\boldmath$u$}$). In this augmented cut posterior, $p(\text{\boldmath$\theta$}|\text{\boldmath$\psi$},\text{\boldmath$u$},{\cal D})=p(\text{\boldmath$\theta$}|\text{\boldmath$\psi$},{\cal D})$ and the marginal in $(\text{\boldmath$\psi$}^\top,\text{\boldmath$\theta$}^\top)^\top$ is the required cut posterior at (ref). We now discuss how to approximate (ref) using recent developments in variational inference methods.
The augmented cut posterior at (ref) is estimated using Algorithm (ref) with approximation
As before, a $N(\text{\boldmath$\mu$},LL^\top)$ approximation is used in $(\text{\boldmath$\psi$}^\top,\text{\boldmath$\theta$}^\top)^\top$, which has marginal in $\text{\boldmath$\psi$}$ with density $q_{\widetilde{\lambda}_a}(\text{\boldmath$\psi$})=\phi_N(\text{\boldmath$\psi$};\text{\boldmath$\mu$}_\psi,L_{\psi}L_{\psi}^\top)$ and parameters $\widetilde{\text{\boldmath$\lambda$}}_a=(\text{\boldmath$\mu$}_\psi^\top,\text{vech}(L_\psi)^\top)^\top$. In Step 2 of the algorithm, $p_{\text{cut}}(\text{\boldmath$\psi$},\text{\boldmath$u$}|{\cal D})$ is approximated by $q_{\widetilde \lambda}$, for which we consider the family discussed below.
For copula models with discrete-valued marginals, loaiza2019VBDA study approximations to a posterior augmented by latents $\text{\boldmath$u$}$. They consider VAs of the form $q_{\widetilde{\lambda}}(\text{\boldmath$\psi$},\text{\boldmath$u$})= q_{\widetilde{\lambda}_a}(\text{\boldmath$\psi$})q_{\widetilde{\lambda}_b}(\text{\boldmath$u$})$ with $\widetilde{\text{\boldmath$\lambda$}}=(\widetilde{\text{\boldmath$\lambda$}}_a^\top,\widetilde{\text{\boldmath$\lambda$}}_b^\top)^\top$. They found approximations with marginal density in $\text{\boldmath$u$}$ given by \[ q_{\widetilde{\lambda}_b}(\text{\boldmath$u$})=\mathop{\prod_{i=1:n}}_{j=1:m} \frac{\phi_N(\zeta_{ij};\delta_{ij},\omega_{ij})}{(b_{ij}-a_{ij})\phi_N(\zeta_{ij};0,1)}\,,\;\; \zeta_{ij}=\Phi^{-1}\left(\frac{u_{ij}-a_{ij}}{b_{ij}-a_{ij}}\right)\,, \] provide a balance between scalability and accuracy. With this approximation, the parameters $\widetilde{\text{\boldmath$\lambda$}}_b$ consist of the $2nm$ mean and log-variance values $\{\delta_{ij},\log \omega_{ij};i=1,\ldots,n;\, j=1,\dots,m\}$.
This approximation is derived from adopting a normal distribution for a transformation of $u_{ij}\in(a_{ij},b_{ij}]$ to the real line. An advantageous property is that $q_{\widetilde{\lambda}_b}$ can be shown to converge to the exact marginal cut posterior in $\text{\boldmath$u$}$ as $n\rightarrow \infty$, so that for larger datasets it is a very accurate approximation. Another advantage of this approximation is that it is tractable, and fast to learn when combined with stochastic gradient descent (SGD). We implement this optimization with control variates as outlined in loaiza2019VBDA, where further details can be found.
Recent studies have applied high-dimensional copula models to multivariate economic and financial time series to capture both cross-sectional and serial dependence jointly; see smith2015 and nagler2022 for examples. Copula models are attractive because when the marginals are asymmetric, the predictive distributions exhibit time-varying asymmetry, which is an important feature of such data. Out-of-sample density and tail forecasting are the primary objectives of these studies, for which heavy-tailed parametric marginals are preferred. To illustrate the impact of cutting feedback, we use it to account for misspecification of either the marginals or copula function in such a model.
We consider the Gaussian copula model of smithvahey2016, who apply it to $N=4$ U.S. macroeconomic time series observed quarterly, which are $Y_{1,t}$ (Output Growth), $Y_{2,t}$ (Inflation), $Y_{3,t}$ (Unemployment Rate), and $Y_{4,t}$ (Interest Rate). These four variables are observed at times $t=1,\ldots,T$, so that the copula is of dimension $m=NT$, although $n=1$ because this is a single time series. The implicit copula of an $N$-dimensional stochastic process $\{\bm{W}_t\}_{t=1}^T$ that follows a stationary lag $p=4$ Gaussian vector autoregression (VAR) is used. It is a large Gaussian copula with parameter matrix $\Omega$ that is a correlation matrix with a sparse block Toeplitz structure.
Rather than define the copula model likelihood directly in terms of $\Omega$, these authors express it more efficiently in terms of the unique semi-partial correlations. Appendix (ref) shows how to do so, where the unique semi-partial correlations associated with each lag $k=0,1,\ldots,p$ are grouped together and denoted as
This is achieved by writing the Gaussian copula as a sparse D-vine where many of the component pair-copulas have density exactly equal to unity. Denote $\text{\boldmath$\phi$}=\{\text{\boldmath$\phi$}(0),\ldots,\text{\boldmath$\phi$}(p)\}$ as the set of unique semi-partial correlations, then there is a one-to-one relationship between $\text{\boldmath$\phi$}$ and $\Omega$.
The original study considered quarterly data from 1954:Q1 until 2011:Q1. The data were sourced from the Federal Reserve Economic Database and the 2022:Q3 vintage from Real-Time Dataset for Macroeconomists hosted by the Philadelphia Federal Reserve. In our analysis we extend the same economic time series to 2022:Q2, so that $T=274$. The matrix $\Omega$ is of dimension $m=1096$, although it is parsimonious because the underlying copula process has only 72 unique semi-partial correlations $\text{\boldmath$\phi$}$. Regularization is known to improve the predictive performance of standard VAR models, so that smithvahey2016 use a spike-and-slab prior on $\text{\boldmath$\phi$}$ for their copula model. In the current analysis, ridge priors with different levels of regularization at each lag are used. If $\widetilde{\phi}^k_{l_1,l_2}=\Phi^{-1}\left((\phi_{l_1,l_2}^k+1)/2\right)$ is a transformation of $\phi_{l_1,l_2}^k$ to the real line, then the prior $\widetilde{\phi}^k_{l_1,l_2}\sim N(0,\tau^2_k)$ with $\tau^2_k \sim C^+(0,1)$ a half-Cauchy distribution. The unconstrained copula and regularization parameters are therefore $\text{\boldmath$\psi$}=\{\widetilde{\text{\boldmath$\phi$}},\log \tau_0^2,\log \tau_1^2,\ldots,\log \tau_p^2\}$.
In this application, prediction of the distributional tails is necessary to quantify macroeconomic risk. A heavy-tailed parametric model is usually preferred to a non- or semi-parametric one because the latter tends to under-weight the possibility of extreme events, such as that observed during the recent pandemic. We follow the original study where time invariant skew-t marginals were used for each variable (truncated to positive values for the Interest Rate variable). However, given the economic shocks since 2011:Q1, it is uncertain whether or not this choice of marginals or Gaussian copula remain suitable for the extended dataset used here. Therefore, we consider cutting feedback first from the copula parameters $\text{\boldmath$\psi$}$ to the marginal parameters $\text{\boldmath$\theta$}$ (the type 1 cut posterior), and then also from $\text{\boldmath$\theta$}$ to $\text{\boldmath$\psi$}$ (the type 2 cut posterior). Because $m$ is large it is infeasible to compute the cut posteriors exactly, and VI was used. For the type 2 cut posterior, the variational approximation to the augmented posterior was employed as outlined in Section (ref). When solving the variational optimizations, a SGD algorithm with ADADELTA learning rate was used with 2000 steps.
To judge the accuracy of the different posteriors, we calculate a log-score metric using the posterior predictive distribution as follows. If $\text{\boldmath$y$}_t=(y_{1,t},\ldots,y_{N,t})^\top$ is the observed value of the $N=4$ variables $\bm{Y}_t=(Y_{1,t},\ldots,Y_{N,t})^\top$ at time $t$, then the posterior predictive density $h$ steps ahead is
Here, $\pi_t(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$ is a posterior density based on the data $y_1,\ldots,y_t$, for which we consider both variational cut posteriors and also the joint posterior. The integral is evaluated by averaging over 5000 draws from $\pi_t$. For the conventional posterior these are obtained using an MCMC scheme as in smithvahey2016 but where the regularization parameters $\tau_0^2,\ldots,\tau_p^2$ are also drawn. Drawing from the variational cut posteriors is straightforward because they are fixed form Gaussian approximations. Conditional on $(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$, draws from the predictive density $p(\text{\boldmath$y$}_{t+h}|\text{\boldmath$y$}_{t},\ldots,\text{\boldmath$y$}_{t-p+1},\text{\boldmath$\theta$},\text{\boldmath$\psi$})$ can be obtained using the sparse D-vine representation of the Gaussian copula as outlined in smithvahey2016.
A log-score metric for variable $j$ predicted $h$ steps ahead can be computed as \[ LS_{j,h}=\sum_{t=p}^{T-h} {\log \widehat{f_{t+h|t}}}(y_{j,t+h}|\text{\boldmath$y$}_{t},\ldots,\text{\boldmath$y$}_{t-p+1})\,. \] Here, $\log\widehat{f_{t+h|t}}(y_{j,t+h}|\text{\boldmath$y$}_{t},\ldots,\text{\boldmath$y$}_{t-p+1})$ is a kernel density estimate of the logarithm of draws from (ref), evaluated at the observed value $y_{j,t+h}$. Higher values of this log-score indicate better calibrated posterior distributions $\left\{\pi_p,\ldots,\pi_{T-h}\right\}$ for predictive purposes.
Figure (ref) plots $LS_{j,h}$ for each variable $h=1,\ldots,8$ quarters ahead, which matches the typical macroeconomic forecast horizon. By this metric, the type 2 cut posterior is a substantial improvement over the conventional posterior and type 1 cut posterior for GDP Growth (the main forecast variable), Inflation and the Interest Rate. This suggests that misspecification of the marginals impacts posterior inference, much more than any potential misspecification of the copula function. The approach of using predictive performance to select between cut and conventional posteriors has been discussed previously by carmona+n22 in the context of semi-modular inference.
Figure (ref) plots the estimated skew-t marginal densities for the four macroeconomic variables, along with histograms of the data. The three posterior estimates differ substantially, highlighting the impact of cutting feedback in this model. The histograms show that skew-t distributions are likely to be a misspecification for the copula model marginals in our extended dataset. For example, between 2011:Q1 and 2022, the Federal Reserve set interest rates to historical near-zero lows, corresponding to a mode at these values in the histogram in panel (c). While a truncated skew-t was an appropriate marginal for the pre-2011 data studied by smithvahey2016, it is inappropriate for the extended dataset that has a bimodal marginal in Interest Rate. For this reason, the type 2 cut posterior correctly cuts feedback from the misspecified marginals when computing inference about the $\text{\boldmath$\psi$}$. This increases the overall accuracy of inference, as measured by the log-score metrics, relative to the conventional posterior.
Finally, we consider the matrices $R(k)\equiv \{r_{i,j}(k)\}$ of pairwise Spearman's rho values $r_{i,j}(k)=\rho(Y_{i,t},Y_{j,t-k})$. These are a function of the posterior of $\text{\boldmath$\phi$}$ as outlined in smithvahey2016, and their estimates provide important macroeconomic insights. Figure (ref) plots mean estimates of $R(0)$, $R(1)$, $R(2)$ and $R(3)$ using the conventional posterior (left hand panels), and using the type 2 cut posterior (right hand side). Cutting feedback perturbs these Spearman correlation estimates. For example, the pairwise correlation between the Interest Rate at time $t-3$ and Inflation at time $t$ is estimated to be $0.092$ in the conventional posterior, whereas in the type 2 cut posterior it is $-0.076$. The latter is more consistent with monetary policy, where interest rate increases are often aimed at reducing future inflation.
The modular nature of copula models can greatly simplify the specification of many multivariate stochastic models. It can also be used to improve the accuracy of statistical inference under potential model misspecification. As far as we are aware, this is the first paper to propose cutting feedback methods to do so. We show theoretically and empirically that these methods can be more accurate in misspecified models than the conventional Bayesian posterior.
Previous inference methods that control for misspecification of the copula function when estimating the marginals include IFM joexu1996,joe2005. For parametric marginals this is usually implemented using a two-stage maximum likelihood procedure, to which we show the type 1 cut posterior mean is asympototically equivalent. For nonparametric marginals, a well-established approach is to estimate the marginals using their empirical distribution functions, followed by estimating the copula parameters using pseudo-maximum likelihood; see oakes1994 and genest95. This can be numerically unstable in higher dimensions, in which case kernel density estimators may be adopted for the marginals. If Bayesian nonparametric distributions hjort2010bayesian are used to model the marginals, then estimation using our proposed type 1 cut posterior provides a Bayesian equivalent which can be used in high dimensions when evaluated by variational methods. grazianliseo17 also suggest Bayesian estimation of a copula model by generating each $\text{\boldmath$\theta$}_j$ from their marginal posteriors, as at the first step of Algorithm (ref) when evaluating the type 1 cut posterior. However, they employ these draws to evaluate an approximate posterior of a dependence parameter based on an exponentially tilted likelihood, rather than a cut posterior.
Methods that control for misspecification of the marginals when estimating the copula function are rare, especially in high-dimensions. kim2007comparison demonstrates that adopting nonparametric marginals as in genest95 can guard against this, but this will be at the cost of reduced statistical efficiency when the marginals are in fact well-specified. Our type 2 cut posterior guards against this type of misspecification while attempting to limit any loss in statistical efficiency. Table (ref) summarizes our theoretical results. Along with the conventional posterior, both types of cut posterior are correctly calibrated asymptotically when the copula model is well-specified. However, unlike the conventional posterior, the type 2 cut posterior is also correctly calibrated under misspecification of the marginal models, and the type 1 cut posterior under misspecification of the copula function.
Evaluation of cut posteriors is difficult, and another contribution of our paper is the development of variational methods to do so for copula models. The definition of the type 2 cut posterior using a pseudo rank likelihood complicates computation, although this can be overcome by considering an augmented posterior with the cut posterior as its marginal in $(\text{\boldmath$\theta$},\text{\boldmath$\psi$})$. Application of the variational methods to a 1096 dimension Gaussian copula for a macroeconomic forecasting application demonstrates their speed and efficiency in high dimensions.
Finally, we note that macroeconomic example is also interesting in itself. Copula time series models have strong potential smithmin2010,smith2015,smithman2018,nagler2022, but selection of an appropriate copula function or parametric marginals can be difficult. In this case, guarding against misspecification is valuable, and our empirical work shows a cut posterior can increase density forecasting accuracy relative to the conventional posterior. Further useful applications of our new Bayesian methodology for cutting feedback in copula modeling await.
\oldappendix {\appendixname A\arabic{section}\quad}
\setcounter{table}{0} \setcounter{figure}{0} \setcounter{algorithm}{0}
This appendix gives the likelihood for the Gaussian copula model of smithvahey2016, to which we refer for full details. Let the $T$ values of the VAR($p$) process be stacked into vector $\bm{W}=(\bm{W}_1^\top,\ldots,\bm{W}_T^\top)^\top=(W_1,W_2,\ldots,W_m)^\top \sim N(0,\Omega)$. The VAR is stationary and constrained to have unit marginal variances, so that \[ \Omega= \left[
\right] \] is a block Toeplitz correlation matrix with $\text{corr}(\bm{W}_t,\bm{W}_s)=\Omega(|t-s|)$. For $i>j+1$, define the semi-partial correlation $\varphi_{i,j}=\mbox{corr}(W_i,W_j|W_{j+1},\ldots,W_{i-1})$ and $\varphi_{i+1,i}=\text{corr}(W_{i+1},W_{i})$. There is a one-to-one transformation between $\Omega$ and the semi-partial correlations $\text{\boldmath$\varphi$}=\{\varphi_{i,j}\}_{i=1:N,j<i}$ due to Yule; e.g. see daniels2009. For a stationary VAR($p$) model, the majority of the elements in $\text{\boldmath$\varphi$}$ are either exactly zero or replicated values. smithvahey2016 show how to identify the unique values, which are denoted as $\text{\boldmath$\phi$}$ in Section (ref), and organize these into the blocks at (ref) that capture serial dependence at different lags.
Our copula model in Section (ref) uses the implicit copula of $\bm{W}$, which is a Gaussian copula with parameter matrix $\Omega$. It is well-known that a Gaussian copula can be written as a D-vine czado2019 with density
where $c_{i,j}(\cdot,\cdot;\varphi_{i,j})$ is a bivariate Gaussian copula density with parameter $\varphi_{i,j}$ given by the semi-partial correlation defined above. When $\varphi_{i,j}=0$ the pair-copula is the independence copula with density $c_{i,j}(\cdot,\cdot;0)=1$. The arguments of each pair-copula, $u_{i|j+1}$ and $u_{j|i+1}$, can be computed from $\text{\boldmath$u$}$ and $\text{\boldmath$\phi$}$ efficiently using the recursive algorithm outlined in Appendix A of smithvahey2016. This also gives an expression for the product at (ref) in terms of only the non-independence pair-copula densities (i.e. those pair-copula densities which are not equal to unity). Finally, because this model is for a single time series, the likelihood is simply given by (ref).
\setcounter{section}{0}