The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
85,759 characters
myblue Cutting Feedback in Misspecified Copula Models
\setlength{\abovedisplayskip}{0.15cm}
\setlength{\belowdisplayskip}{0.15cm}
\pagestyle{empty}
\begin{titlepage}
\title{\bfseries\sffamily\color{myblue}
Cutting Feedback in Misspecified Copula Models}
\author{Michael Stanley Smith, Weichang Yu, David J. Nott and David T. Frazier}
\date{\today}
\maketitle
\noindent
{\small Michael Smith is Professor of Management (Econometrics) at Melbourne Business School in the University of Melbourne, Weichang Yu is Lecturer in the School of Mathematics and Statistics at the University of Melbourne, David Nott is an Associate Professor in the Department of Statistics and Data Science at the National University of Singapore and is affiliated with the NUS Institute of Operations Research and Analytics, and
David Frazier is Professor in Econometrics and Business Statistics at Monash University. David Frazier acknowledges funding from the Australian Research Council under Projects DE200101070 and DP200101414, while
David Nott and Michael Smith's research was supported by the Ministry of Education, Singapore, under the Academic Research Fund Tier 2 (MOE-T2EP20123-0009). The authors thank two anonymous referees whose insightful comments have improved the paper.
}
\newpage
\begin{center}
\mbox{}\vspace{2cm}\\
{\LARGE \title{\bfseries\sffamily\color{myblue} Cutting Feedback in Misspecified Copula Models}
}\\
\vspace{1cm}
{\Large Abstract}
\end{center}
\vspace{-1pt}
\onehalfspacing
\noindent
In copula models the marginal distributions and
copula function are specified separately. We treat these as two modules in a modular Bayesian inference framework, and propose conducting modified Bayesian inference by ``cutting feedback''.
Cutting feedback
limits the influence of potentially misspecified modules in posterior inference.
We consider two types of cuts. The first limits
the influence of a misspecified copula on inference for the marginals, which is
a Bayesian analogue of the popular Inference for Margins (IFM) estimator. The second limits the
influence of misspecified marginals on inference for the copula parameters by using
a pseudo likelihood of the ranks to define the cut model. We establish that if only one of the modules is misspecified, then the appropriate cut posterior gives accurate uncertainty quantification asymptotically for the parameters in the other module. Computation of the cut posteriors
is difficult, and new variational inference methods to do so are proposed.
The efficacy of the new methodology is demonstrated using both simulated data and a substantive multivariate time series
application from macroeconomic forecasting.
In the latter, cutting feedback from misspecified marginals to a 1096 dimension copula
improves posterior inference and predictive accuracy greatly, compared
to conventional Bayesian inference.
\vspace{20pt}
\noindent
{\bf Keywords}: Modular Inference, Posterior Consistency, Pseudo Rank Likelihood, Time Series Copula, Variational Inference.
\end{titlepage}
\newpage
\pagestyle{plain}
\setcounter{equation}{0}
\section{Introduction}\label{sec:intro}
A copula model
specifies a multivariate distribution using its marginal distributions and a copula function to capture the dependence structure~\citep{nelsen06,joe2014dependence}. This simplifies multivariate stochastic modelling, making copula models popular in hydrology~\citep{genest2007metaelliptical},
financial econometrics~\citep{patton2006},
transportation studies~\citep{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'' \citep{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~\citep{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~\citep{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 \citep{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~\citep{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~\citep{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~\cite{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{sec:bcm} gives some necessary background on copula models, variational inference and cutting feedback. Section~\ref{sec:cutting1} 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{sec:cutting2} 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{sec:macro}
contains the macroeconomic application, while
Section~\ref{sec:conc} concludes.
\section{Background}\label{sec:bcm}
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.
\subsection{Copula Models}
If $\bm{Y}=(Y_1,\ldots,Y_m)^\top\sim F_Y$ with marginals $Y_j\sim F_{j}$,
then the joint distribution function
\begin{equation}
F_Y(\text{\boldmath$y$})=C\left(F_1(y_1),\ldots,F_m(y_m)\right)\,,\label{eq:copmod}
\end{equation}
where $\text{\boldmath$y$}=(y_1,\ldots,y_m)^\top$ and $C:[0,1]^m \rightarrow \mathbb{R}^+$ is
a copula function; see~\citet[p.45]{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~\cite{nelsen06} and \cite{joe2014dependence} for some examples.
If $F_1,\ldots,F_m$ are all continuous distributions, the joint density of $\bm{Y}$ is
\begin{equation}
f_Y(\text{\boldmath$y$};\text{\boldmath$\theta$},\text{\boldmath$\psi$})=c\left(F_1(y_1;\text{\boldmath$\theta$}_1),\ldots,F_m(y_m;\text{\boldmath$\theta$}_m);\text{\boldmath$\psi$}\right)\prod_{j=1}^m f_j(y_j;\text{\boldmath$\theta$}_j)\,,\label{eq:copden}
\end{equation}
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~\cite{genest2007}.
Let ${\cal D}=\{\text{\boldmath$y$}_1,\ldots,\text{\boldmath$y$}_n\}$ be $n$ observations drawn independently from $F_Y$
at~\eqref{eq:copmod}. 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~\cite{PitChaKoh2006,silva2008copula,min2010,smithmin2010} and~\cite{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.
\subsection{Variational Inference}
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
\begin{equation}
q^*(\text{\boldmath$\theta$}, \text{\boldmath$\psi$}) = \operatorname*{{arg\,min}}_{q \in \mathcal F} \int \int q(\text{\boldmath$\theta$}, \text{\boldmath$\psi$}) \log \left \{ \frac{q(\text{\boldmath$\theta$}, \text{\boldmath$\psi$})}{p(\text{\boldmath$\theta$}, \text{\boldmath$\psi$} \, \vert \mathcal D)} \right \} \, \text{d} \text{\boldmath$\theta$} d \text{\boldmath$\psi$}.\label{eq:viprob}
\end{equation}
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~\citep{titsias2014doubly,kucukelbir2017automatic,tan2018gaussian}.
It is straightforward to show (e.g. see~\cite{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~\citep{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~\cite{zeiler12} that
we use here. The combination of stochastic optimization and generic approximations is often called black box VI~\citep{ranganath14,titsias2014doubly}.
VI has been used to estimate copula models
by~\cite{loaiza2019VBDA}, \cite{nguyen2020VI} and~\cite{smithklein2021}.
\subsection{Cutting feedback methods}\label{sec:cfm}
Cutting feedback is a form of Bayesian
modular inference~\citep{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
\cite{lunn+bsgn09}, \cite{plummer15}, \cite{jacob+mhr17} and~\cite{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
\begin{equation}
g({\cal D}|\bm{\eta}) = g_1({\cal D}|\text{\boldmath$\eta$}_1)g_2({\cal D}|\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2)\,.\label{eq:cflikefactor}
\end{equation}
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~\cite{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
\begin{align}
p(\text{\boldmath$\eta$}_1|{\cal D}) & =\frac{p(\text{\boldmath$\eta$}_1)g_1({\cal D}|\text{\boldmath$\eta$}_1)\bar{g}_2({\cal D}|\text{\boldmath$\eta$}_1)}{\bar{g}({\cal D})}, \label{theta-marginal}
\end{align}
and the conditional posterior density is
\begin{align}
p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1,{\cal D}) & =
\frac{p(\text{\boldmath$\eta$}_2|\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)}. \label{psi-conditional}
\end{align}
In~\eqref{theta-marginal}, $\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~\eqref{theta-marginal} to define the following
marginal cut posterior density
\begin{align}
p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D}) & = \frac{p(\text{\boldmath$\eta$}_1)g_1({\cal D}|\text{\boldmath$\eta$}_1)}{\int p(\text{\boldmath$\eta$}_1')g_1({\cal D}|\text{\boldmath$\eta$}_1')\,d\text{\boldmath$\eta$}_1'}. \label{cut-theta-marginal}
\end{align}
The joint cut posterior density is then defined as
\begin{align}
p_{\text{cut}}(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2|{\cal D}) & =
p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1,{\cal D}), \label{cut-joint}
\end{align}
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\,.
\]
\subsection{Exact cut posterior computation}
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~\eqref{cut-joint} 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~\cite{plummer15} suggested
using
``nested MCMC'' as in Algorithm~\ref{alg:nestedmcmc} below.
Other methods for cut posterior evaluation are discussed by \cite{liu+g20}, \cite{jacob2020unbiased} and
\cite{pompe+j21}.
\begin{algorithm}
\caption{Nested MCMC Sampler for Cut Posterior}
\label{alg:nestedmcmc}
\begin{algorithmic}
\State Generate sample $\text{\boldmath$\eta$}_1^{(1)},\ldots,\text{\boldmath$\eta$}_1^{(S)} \sim p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})\propto p(\text{\boldmath$\eta$}_1)g_1({\cal D}|\text{\boldmath$\eta$}_1)$ using an MCMC scheme
\For{$s=1,\ldots,S$}
\State Generate single value $\text{\boldmath$\eta$}_2^{(s)}\sim p(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1^{(s)},{\cal D})$ using an MCMC scheme
\EndFor
\end{algorithmic}
\end{algorithm}
\subsection{Variational cut posterior computation}
Given the difficulty of exact cut posterior computation,
variational inference methods to do so have been suggested
by \cite{yu+ns21} and~\cite{carmona+n22}. Lemma 1 of \cite{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~\eqref{eq:viprob} 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~\cite{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{alg:vigeneral} below outputs an approximation
to the joint cut posterior.
\begin{algorithm}
\caption{VI for Cut Posterior}
\label{alg:vigeneral}
\begin{algorithmic}
\State 1. Select a fixed form variational approximation $q_{\lambda}(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2)=q_{\widetilde{\lambda}}(\text{\boldmath$\eta$}_1)q_{\breve{\lambda}}(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)$
\State 2. Solve the optimization
$$\widetilde{\bm{\lambda}}^*=\arg \min_{\widetilde{\lambda}} D_{KL}\left( q_{\widetilde{\lambda}}(\text{\boldmath$\eta$}_1)||p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})\right)$$
to obtain an approximation $q_{\widetilde{\bm{\lambda}}^*}(\text{\boldmath$\eta$}_1)$ of $p_{\text{cut}}(\text{\boldmath$\eta$}_1|{\cal D})$
\State 3. Solve the optimization
$$\breve{\bm{\lambda}}^*=\arg \min_{\breve{\lambda}}
D_{KL}\left( q_{\widetilde{\lambda}^*}(\text{\boldmath$\eta$}_1)q_{\breve{\lambda}}(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)||
p(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2|{\cal D})\right)$$
to obtain an approximation $q^*(\text{\boldmath$\eta$})=q_{\widetilde{\lambda}^*}(\text{\boldmath$\eta$}_1)q_{\breve{\lambda}^*}(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)$ of $p_{\text{cut}}(\text{\boldmath$\eta$}_1,\text{\boldmath$\eta$}_2|{\cal D})$
\end{algorithmic}
\end{algorithm}
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[\begin{array}{c}
\bm{\mu}_{\eta_1} \\
\bm{\mu}_{\eta_2}
\end{array}\right], \;\;\;\;
L=\left[\begin{array}{cc}
L_{\eta_1} & \bm{0} \\
L_{\eta_1,\eta_2} & L_{\eta_2}
\end{array} \right],$$
then $q_{\widetilde{\lambda}}(\text{\boldmath$\eta$}_1)$ and $q_{\breve{\lambda}}(\text{\boldmath$\eta$}_2|\text{\boldmath$\eta$}_1)$
are both Gaussian densities with parameters
$\widetilde{\text{\boldmath$\lambda$}}=(\bm{\mu}_{\eta_1}^\top,
\text{vech}(L_{\eta_1})^\top)^\top$, and
$\breve{\text{\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. \citet{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.
\section{Cutting Feedback for Misspecified Copulas}\label{sec:cutting1}
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{sec:cutting2}.
\subsection{Type 1 cut posterior specification}
For the copula model with density at~\eqref{eq:copden}, 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 \eqref{cut-theta-marginal} 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
\begin{equation}
p_{\text{cut}}(\bm{\theta},\text{\boldmath$\psi$}|{\cal D})=p_{\text{cut}}(\bm{\theta}|{\cal D})p(\bm{\psi}|\bm{\theta},{\cal D})\,,\label{eq:jcutpost1}
\end{equation}
and we consider its computation using both the nested MCMC and variational approaches in Algorithms~\ref{alg:nestedmcmc} and~\ref{alg:vigeneral}.
\subsection{Theoretical equivalence with IFM}\label{sec:theory_cut1}
Among the most popular methods for estimating copula models is the ``inference for margins''(IFM) procedure of~\cite{joexu1996} and~\cite{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~\eqref{eq:jcutpost1}
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{lem:ifm1} below establishes that the IFM point estimator and the cut posterior mean are asymptotically equivalent.
\begin{lemma}\label{lem:ifm1}
If Assumptions~\ref{ass:cons} and~\ref{ass:dist2} in Part~\ref{app:dtf} of the Web Appendix are satisfied, then
$
\sqrt{n}(\overline{\bm{\eta}}-\widehat\bm{\eta})=o_p(1).
$
\end{lemma}
\noindent Assumptions \ref{ass:cons} and~\ref{ass:dist2} are similar to the standard regularity conditions employed in two-step copula modeling to deduce asymptotic normality of the IFM point estimator in \cite{joe2005}; see Part~\ref{app:dtf} of the Web Appendix for their specification and a
detailed discussion.
Lemma~\ref{lem:ifm1} 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{ass:cons} and \ref{ass:dist2} in Part~\ref{app:dtf} 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}$) :
\begin{flalign*}
\mathcal{I}=-\lim_{n\rightarrow+\infty}n^{-1}\operatorname{E}\left(\nabla^2_{\bm{\theta}\bm{\theta}}\log g_1(\mathcal{D}\mid\bm{\theta}_{0})\right)\,,
\end{flalign*}
where $\operatorname{E}$ is the expectation with respect to $P_0$.
Then Lemma~\ref{lem:two} below shows how the cut posterior for $\bm{\theta}$ quantifies uncertainty.
\begin{lemma}\label{lem:two}
If Assumptions \ref{ass:cons} and \ref{ass:dist2} in Part~\ref{app:dtf} of the Web Appendix are satisfied, then
$$
\int\left|p_{\mathrm{cut}}(\bm{\theta}|\mathcal{D})
-\phi_N\left(\bm{\theta};\widehat\bm{\theta},[n\mathcal{I}]^{-1}\right)\right| d\bm{\theta}
=o_p(1). $$
\end{lemma}
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:
\begin{flalign*}
\mathcal{M}=\begin{pmatrix} \mathcal{M}_{\bm{\theta}\bm{\theta}}&\mathcal{M}_{\bm{\theta}\bm{\psi}}\\\mathcal{M}_{\bm{\psi}\bm{\theta}}&\mathcal{M}_{\bm{\psi}\bm{\psi}}
\end{pmatrix},\text{ where }\mathcal{M}_{\bm{\theta}\bm{\psi}}=-\lim_{n\rightarrow+\infty}n^{-1}\operatorname{E}\left(\nabla^2_{\bm{\theta}\bm{\psi}}\log g_2(\mathcal{D}\mid\bm{\theta}_{0},\bm{\psi}_0)\right),
\end{flalign*}
with $\mathcal{M}_{\bm{\theta}\bm{\theta}}$ and $\mathcal{M}_{\bm{\psi}\bm{\psi}}$ defined analogously; then the following lemma holds.
\begin{lemma}\label{lem:ifm2}
If Assumptions \ref{ass:cons} and \ref{ass:dist2} in Part~\ref{app:dtf} of the Web Appendix are satisfied, then
$$
\int\left|p_{\mathrm{cut}}(\bm{\psi}|\mathcal{D})-
\phi_N\left(\bm{\psi};\widehat{\bm{\psi}},n^{-1}(\mathcal{M}_{\bm{\psi}\bm{\psi}}^{-1}+\mathcal{M}_{\bm{\psi}\bm{\psi}}^{-1}\mathcal{M}_{\bm{\psi}\bm{\theta}}\mathcal{I}^{-1}\mathcal{M}_{\bm{\theta}\bm{\psi}}\mathcal{M}_{\bm{\psi}\bm{\psi}}^{-1})\right)\right| d \bm{\psi}=o_p(1).
$$
\end{lemma}
Lemma \ref{lem:ifm2} 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.
\subsection{Simulation~1}
\label{sec::simExample1}
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 \citep{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.
\begin{table}[htbp]
\caption{Parameter Point Estimation Accuracy in Simulation 1 ($n=1000$)}
\label{tab:sim1biasrmse}
\begin{center}
\begin{tabular}{ccccccccc}
\hline\hline
& \multicolumn{5}{c}{Misspecified Copula Fit} & &\multicolumn{2}{c}{Correct Copula Fit} \\ \cline{2-6}\cline{8-9}
& Uncut/ & Cut/ & IFM & Uncut/ & Cut/ & &Uncut/ &Cut/\\
& MCMC & MCMC & &VI &VI & &MCMC &MCMC\\ \hline
Parameter &{\em Bias} & & & & & & &\\ \cline{2-6}\cline{8-9}
$\mu$ & 0.0102 & {\bf 0.0019} & 0.0039 & 0.0089 & {\bf 0.0019} & &0.0014 &0.0019\\
$\sigma^2$ & 0.0360 & 0.0069 & {\bf 0.0012} & 0.0386 & 0.0288 & &0.0057 &0.0069\\
$\alpha$ & -0.1803 & -0.0073 & 0.0137 & -0.1869 & {\bf -0.0033} & &-0.0035 &-0.0073\\
$\beta$ & -0.0924 & -0.0046 & 0.0045 & -0.0938 & {\bf -0.0028} & &-0.0019 &-0.0046\\
$\tau$ & 0.0192 & 0.0107 & 0.0139 & 0.0176 & {\bf 0.0090} & &0.0002 &-0.0101\\
&{\em RMSE} & & & & & & &\\ \cline{2-6}\cline{8-9}
$\mu$ & 0.0327 & {\bf 0.0295} & 0.0296 & 0.0298 & 0.0307 & & 0.0263 & 0.0295\\
$\sigma^2$ & 0.0612 & 0.0472 & {\bf 0.0333} & 0.0620 & 0.0560 & & 0.0447 & 0.0472\\
$\alpha$ & 0.3774 & {\bf 0.3114} & 0.3132 & 0.3580 & 0.3114 & & 0.3048 & 0.3114\\
$\beta$ & 0.1689 & {\bf 0.1406} & 0.1413 & 0.1644 & 0.1407 & & 0.1395 & 0.1406\\
$\tau$ & 0.0244 & 0.0181 & 0.0209 & 0.0232 & {\bf 0.0175} & & 0.0132 & 0.0182\\ \hline\hline
\end{tabular}
\end{center}
The bias and RMSE values of different point estimators computed over the $S=500$ simulation replicates, with the lowest values
in bold. Results on the left are where the misspecified copula is fit using the posterior mean from the conventional (i.e. uncut) and cut (type 1) posteriors, computed exactly using MCMC or approximately using VI. IFM is included for comparison. Results on the right are where the correct copula is fit using the
posterior mean from the conventional (i.e. uncut) and cut (type 1) posteriors
computed exactly using MCMC.
\end{table}
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{tab:sim1biasrmse} (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{tab:sim1biasrmse} (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{tab:sim1covkl} 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.
\begin{table}[htbp]
\caption{Posterior Coverage and Copula Model Accuracy in Simulation 1 ($n=1000$)} \label{tab:sim1covkl}
\begin{center}
\begin{tabular}{ccccccccc}
\hline\hline
& \multicolumn{5}{c}{Misspecified Copula Fit} & &\multicolumn{2}{c}{Correct Copula Fit} \\ \cline{2-6}\cline{8-9}
& Uncut/ & Cut/ & IFM & Uncut/ & Cut/ & &Uncut/ &Cut/\\
& MCMC & MCMC & &VI &VI & &MCMC &MCMC\\ \hline
Parameter &\multicolumn{3}{l}{\em Coverage Probabilities} & & & & &\\ \cline{2-6}\cline{8-9}
$\mu$ & 0.9600 & 0.9660 & - & 0.9760 & {\bf 0.9520} & &0.9510 &0.9660\\
$\sigma^2$ & 0.8660 & 0.9380 & - & 0.9780 & {\bf 0.9560} & &0.9460 & 0.9380\\
$\alpha$ & 0.8920 & {\bf 0.9560} & - & 0.8240 & 0.9760 & &0.9460 &0.9560\\
$\beta$ & 0.8660 & {\bf 0.9380} & - & 0.9840 & 0.9680 & &0.9460 &0.9380\\
$\tau$ & 0.5400 & 0.7440 & - & 0.7060 & {\bf 0.7780} & &0.9500 &0.9260\\
Component &\multicolumn{3}{l}{\em Mean Predictive KL Divergence} & & & & &\\ \cline{2-6}\cline{8-9}
$f_1$ & 0.0014 & {\bf 0.0010} & {\bf 0.0010} & 0.0014 & 0.0012 & & 0.0008 & 0.0010\\
$f_2$ & 0.0014 & {\bf 0.0010} & 0.0013 & 0.0012 & {\bf 0.0010} & & 0.0010 & 0.0010\\
$c$ & 0.1488 & 0.1424 & 0.1457 & 0.1499 & {\bf 0.1411} & & 0.0008 & 0.0012\\ \hline\hline
\end{tabular}
\end{center}
Top: coverage probabilities for 95\% credible intervals for each parameter, with
those closest to 0.95 in bold. Bottom: the mean predictive KL divergence for the two marginals and the copula density for their estimate, with the lowest values
in bold.
Results on the left are where the misspecified copula is fit using the conventional (i.e. uncut) and cut (type 1) posteriors, computed exactly using MCMC or approximately using VI. IFM is included for comparison. Results on the right are where the correct copula is fit using the conventional (i.e. uncut) and cut (type 1) posteriors
computed exactly using MCMC.
\end{table}
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$.
\section{Cutting Feedback for Misspecified Marginals}\label{sec:cutting2}
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{sec:cutting1}, and to do so in higher dimensions we introduce an
extension of this pseudo likelihood~\citep{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.
\subsection{Type 2 cut posterior specification}
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
\begin{equation}
p_{\text{PL}}(r({\cal D})|\text{\boldmath$\psi$}):=\prod_{i=1}^n \Delta_{a_{i1}}^{b_{i1}}\cdots \Delta_{a_{im}}^{b_{im}} C(\text{\boldmath$v$};\text{\boldmath$\psi$})\,,\label{eq:rlike}
\end{equation}
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
\begin{equation*}
\Delta_{a_{ij}}^{b_{ij}} C(\text{\boldmath$v$};\text{\boldmath$\psi$}):=C(v_1,\ldots,v_{j-1},b_{ij},v_{j+1},\ldots,v_m;\text{\boldmath$\psi$})-
C(v_1,\ldots,v_{j-1},a_{ij},v_{j+1},\ldots,v_m;\text{\boldmath$\psi$})\,,
\end{equation*}
is a differencing operator over element $j$~\citep[p.43]{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~\cite{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~\cite{genest95}, but corrects for the discrete nature of the rank data.
The pseudo rank likelihood at~\eqref{eq:rlike} 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
\begin{equation}
p_{\text{cut}}(\text{\boldmath$\psi$}|{\cal D})=\frac{p_{\text{PL}}(r(D)|\text{\boldmath$\psi$})p(\text{\boldmath$\psi$})}{\int p_{\text{PL}}(r(D)|\text{\boldmath$\psi$}')p(\text{\boldmath$\psi$}')d\text{\boldmath$\psi$}'}\,.\label{eq:mcutpsi}
\end{equation}
This definition fits into the two module system described in Section~\ref{sec:cfm} by considering
the factorization at~\eqref{eq:cflikefactor} 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~\eqref{eq:mcutpsi} 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
\begin{equation}
p_{\text{cut}}(\text{\boldmath$\psi$},\text{\boldmath$\theta$}|{\cal D})=p_{\text{cut}}(\text{\boldmath$\psi$}|{\cal D})p(\text{\boldmath$\theta$}|\text{\boldmath$\psi$},{\cal D})\,,
\label{eq:sec4jntcutpost}
\end{equation}
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{alg:vigeneral}.
\subsection{Theoretical behavior of the type 2 cut posterior}\label{sec:theory_cut2}
An advantage of the pseudo rank likelihood at~\eqref{eq:rlike} is that it is both computationally and
theoretically more tractable than the rank likelihood of~\cite{hoff07} and others.
As~\cite{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,~\cite{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~\eqref{eq:rlike} 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 \eqref{eq:mcutpsi} and~\eqref{eq:sec4jntcutpost}. 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{thm:ranks} below characterizes the behavior of $p_{\mathrm{cut}}(\bm{\psi}|\mathcal{D})$.
\begin{theorem}\label{thm:ranks}
If Assumptions \ref{ass:DGP1}-\ref{ass:crit1} in Part~\ref{app:rank} of the Web Appendix are satisfied, then
$$
\int\left|p_{\mathrm{cut}}(\bm{\psi}|\mathcal{D})-\phi_N\left(\bm{\psi};\widehat\bm{\psi}_{r},n^{-1}\mathcal{M}_{\bm{\psi}\bm{\psi}}(\bm{\psi}_\star)^{-1}\right)\right|d\bm{\psi}=o_p(1).
$$
\end{theorem}
Assumptions \ref{ass:DGP1}-\ref{ass:crit1} are given in Part~\ref{app:rank} 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{sec:discuss2} of the Web Appendix for further discussion.
Theorem~\ref{thm:ranks} 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{cor:two} in Web Appendix~\ref{app:dtf}.
In such cases, Theorem \ref{thm:ranks} 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{thm:cut2} below characterizes the
behavior of $p_\mathrm{cut}(\bm{\theta}|\mathcal{D})$.
\begin{theorem}\label{thm:cut2}
If Assumptions \ref{ass:cons} and \ref{ass:dist2} in Part~\ref{app:dtf} of the Web Appendix, and Assumptions \ref{ass:DGP1}--\ref{ass:crit1} in Part~\ref{app:rank} of the Web Appendix, and all satisfied, then
$$
\int\left|p_{\mathrm{cut}}(\bm{\theta}|\mathcal{D})-\phi_N\left(\bm{\theta};\widehat\bm{\theta}_{r},n^{-1}\Omega^{-1}\right)\right|d\bm{\theta} =o_p(1).
$$
\end{theorem}
Theorem~\ref{thm:cut2} 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{lem:two}) .
\subsection{Simulation 2}
\label{sec::simExample3}
The simulation study in Section~\ref{sec::simExample1} 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{Simulation3KL} (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{Simulation3KL} (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$.
\begin{table}[htbp]
\caption{Copula Model Estimation Accuracy in Simulation 2 ($n=1000$)}
\label{Simulation3KL}
\begin{center}
\begin{tabular}{lcccccccc}
\hline \hline
& \multicolumn{4}{c}{Misspecified Copula Fit} & &\multicolumn{2}{c}{Correct Copula Fit} \\ \cline{2-5}\cline{7-8}
& Uncut/ & Cut/ & Uncut/ & Cut/ & &Uncut/ &Cut/\\
& MCMC & MCMC &VI &VI & &MCMC &MCMC\\ \hline
Marginal $f_1$ & 0.8993 & {\bf 0.7945} & 0.8941 & 0.8090 & &0.0003 & 0.0003 \\
Marginal $f_2$ & 0.2103 & {\bf 0.1513} & 0.2034 & 0.1597 & & 0.0005 &0.0005 \\
Copula & 0.0166 & {\bf 0.0010} & 0.0084 & {\bf 0.0010} & & 0.0007 & 0.0010\\
\hline\hline
\end{tabular}
\end{center}
Mean predictive KL divergence metrics for the two marginals and the copula density
of the bivariate copula model estimate. The lowest values are in bold. Results are given for the type 2 cut
and conventional (i.e. uncut) posteriors, computed exactly using MCMC or approximately using variational inference. Results are given for both the misspecified copula model (left hand side) and the correctly
specified copula model (right hand side).
\end{table}
\subsection{Augmented type 2 cut posterior}
When $m$ is small, the cut posterior can be evaluated by direct application of Algorithm~\ref{alg:vigeneral}. However, for even moderate values of $m$,
the pseudo rank likelihood at~\eqref{eq:rlike} 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 \cite{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~\cite{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~\eqref{eq:mcutpsi}.
Again, this setup fits into the two module system discussed in Section~\ref{sec:cfm}, 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
\begin{equation}
p_{\text{cut}}(\text{\boldmath$\psi$},\text{\boldmath$u$},\text{\boldmath$\theta$}|{\cal D})=p_{\text{cut}}(\text{\boldmath$\psi$},\text{\boldmath$u$}|{\cal D})p(\text{\boldmath$\theta$}|\text{\boldmath$\psi$},{\cal D})\,,\label{eq:jaugcut}
\end{equation}
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~\eqref{eq:sec4jntcutpost}. We now discuss how to approximate~\eqref{eq:jaugcut} using recent developments in
variational inference methods.
\subsection{Variational inference for the augmented type 2 cut posterior}\label{sec:viaug}
The augmented cut posterior at~\eqref{eq:jaugcut} is estimated using Algorithm~\ref{alg:vigeneral} with approximation
\begin{equation} q_\lambda(\text{\boldmath$\theta$},\text{\boldmath$\psi$},\text{\boldmath$u$})=q_{\widetilde{\lambda}}(\text{\boldmath$\psi$},\text{\boldmath$u$})
q_{\breve{\lambda}}(\text{\boldmath$\theta$}|\text{\boldmath$\psi$})\,.\label{eq:cutva2}
\end{equation}
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, \cite{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 \cite{loaiza2019VBDA},
where further details can be found.
\section{Macroeconomic Example}\label{sec:macro}
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~\cite{smith2015}
and~\cite{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.
\subsection{Gaussian copula model}
We consider the Gaussian copula model of~\cite{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{app:a} 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
\begin{eqnarray}
\text{\boldmath$\phi$}(0) &= &\{\phi_{l_1,l_2}^0\} \mbox{ for } l_2=2,\ldots,N\;\text{and } l_2<l_1\,, \nonumber\\
\text{\boldmath$\phi$}(k) &= &\{\phi_{l_1,l_2}^k\} \mbox{ for } l_2=1,\ldots,N\;\text{and } l_2=1,\ldots,N\;\text{and }k=1,\ldots,p\,.
\label{eq:partials}
\end{eqnarray}
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
\cite{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{sec:viaug}. When solving the variational optimizations, a SGD algorithm with ADADELTA learning rate was used with 2000 steps.
\subsection{Prediction and log-score metric}
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
\begin{equation}
f_{t+h|t}(\text{\boldmath$y$}_{t+h}|\text{\boldmath$y$}_{t},\ldots,\text{\boldmath$y$}_{t-p+1})\equiv\int p(\text{\boldmath$y$}_{t+h}|\text{\boldmath$y$}_{t},\ldots,\text{\boldmath$y$}_{t-p+1},\text{\boldmath$\theta$},\text{\boldmath$\psi$})\pi_t(\text{\boldmath$\theta$},\text{\boldmath$\psi$})\mbox{d}\text{\boldmath$\theta$}\mbox{d}\text{\boldmath$\psi$}\,, \mbox{ for } h\geq 1\,.\label{eq:postpred}
\end{equation}
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~\cite{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~\cite{smithvahey2016}.
\begin{figure}[thb]
\centering
\includegraphics[width=0.7\textwidth]{LogScorePlots.pdf}
\caption{Plots of the log-score posterior predictive metric $LS_{j,h}$ for
the type~1 cut posterior (blue dashed line), type~2 cut posterior (red thick line) and
the conventional (i.e. uncut) posterior (black thin line). Panels (a-d) correspond to variables GDP Growth ($j=1$), Inflation ($j=2$), Interest
Rate ($j=3$) and Unemployment Rate ($j=4$), respectively. In each panel the metric values are plotted for predictions $h=1,\ldots,8$ quarters ahead. Higher values correspond to greater predictive accuracy.}
\label{fig:macroLS}
\end{figure}
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~\eqref{eq:postpred}, 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{fig:macroLS} 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~\cite{carmona+n22} in the context
of semi-modular inference.
\begin{figure}[thb]
\centering
\includegraphics[width=0.7\textwidth]{DensityPlotsMarginals.pdf}
\caption{Plots of the estimated skew-t marginal densities for
(a)~GDP Growth, (b)~Inflation, (c)~Interest Rate, and (d)~Unemployment Rate.
These estimates are evaluated at the posterior means of $\text{\boldmath$\theta$}$ for the type~1 cut posterior (blue dashed line), type~2 cut posterior (red thick line) and
the conventional (i.e. uncut) posterior (black thin line). Histograms of the data are also plotted.}
\label{fig:macromargins}
\end{figure}
\subsection{Estimates}
Figure~\ref{fig:macromargins}
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~\cite{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.
\begin{figure}[p]
\centering
\includegraphics[height=0.9\textheight]{TilePlotSpearmanCorCondensed.pdf}
\caption{Posterior means of the matrices of Spearman pairwise correlations $R(k)$ for $k=0,1,2,3$ from the fitted copula model.
The left hand panels (a,c,e,g) contain results for the conventional (i.e. uncut) joint posterior, while the right hand panels
(b,d,f,h) contain results for the type~2 cut posterior. For example, the estimated Spearman correlation between
the Interest Rate (IR) at time $t-3$ (variable $Y_{4,t-3}$) and the Unemployment Rate (UR) at time $t$ (variable $Y_{3,t}$) is $-0.008$
using the conventional posterior in panel~(g), and $0.370$ using the type~2 cut posterior in panel~(h).}
\label{fig:spearman}
\end{figure}
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~\cite{smithvahey2016}, and their estimates provide important macroeconomic insights.
Figure~\ref{fig:spearman} 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.
\section{Discussion}\label{sec:conc}
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~\citep{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~\cite{oakes1994} and~\cite{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~\citep{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. \cite{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{alg:nestedmcmc} 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.
\begin{table}[htbp]
\caption{Summary of Theoretical Results for Cut Posteriors}
\label{tab:theory}
\begin{center}
\begin{tabular}{lcccccc}
\toprule
\multirow{2}{*}{} &
\multicolumn{2}{c}{Both Correct} &
\multicolumn{2}{c}{Marginal Miss. } &
\multicolumn{2}{c}{Copula Miss. }
\\
& {$\bm{\theta}$} & {$\bm{\psi}$} & {$\bm{\theta}$}& {$\bm{\psi}$} & {$\bm{\theta}$}& {$\bm{\psi}$}\\
\midrule
Conventional Posterior& $\checkmark$ & $\checkmark$ & $\times$ & $\times$& $\times$ & $\times$\\
Type 1 Cut Posterior& $\checkmark$ & $\checkmark$ & $\times$ & $\times$& $\checkmark$ & $\times$\\
Type 2 Cut Posterior & $\checkmark$ & $\checkmark$ & $\times$ &$\checkmark$ & $\times$ & $\times$\\
\bottomrule
\end{tabular}
\end{center}
In the column headings ``Both Correct" indicates that both marginals and copula function are correctly specified; ``Marginal Miss." refers to the case where the marginals are misspecified, but the copula function is correct; and ``Copula Miss." refers to the case where the copula function is misspecified, but the marginals are correct. The parameters $\bm{\theta}$ and $\bm{\psi}$ refer to calibration of that specific parameter with ``$\checkmark$'' indicating (asymptotically) correct calibration, and ``$\times$'' denoting (possibly) inaccurate calibration.
\end{table}
Methods that control for misspecification of the marginals when estimating the copula function are rare, especially in high-dimensions. \cite{kim2007comparison} demonstrates that adopting nonparametric marginals as in~\cite{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{tab:theory} 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~\citep{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.
\newpage
\oldappendix
\newcommand{\appendixname~A\arabic{section}\quad}{\appendixname~A\arabic{section}\quad}
\setcounter{table}{0}
\setcounter{figure}{0}
\setcounter{algorithm}{0}
\section{Copula multivariate time series model}\label{app:a}
This appendix gives the likelihood for the Gaussian copula model of~\cite{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[ \begin{array}{ccc}
\Omega(0) &\cdots &\Omega(T-1) \\
\vdots &\ddots &\vdots \\
\Omega(T-1) &\cdots &\Omega(0)
\end{array} \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~\cite{daniels2009}. For a stationary VAR($p$) model, the
majority of the elements in $\text{\boldmath$\varphi$}$ are either exactly zero or replicated values.
\cite{smithvahey2016} show how to identify the unique values, which are denoted as $\text{\boldmath$\phi$}$ in Section~\ref{sec:macro}, and organize these
into the blocks at~\eqref{eq:partials} that capture serial dependence at different lags.
Our copula model in Section~\ref{sec:macro} 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~\citep{czado2019} with
density
\begin{equation}
c(\text{\boldmath$u$};\text{\boldmath$\phi$})=\prod_{i=2}^m \prod_{j=1}^{i-1} c_{i,j}(u_{i|j+1},u_{j|i-1};\varphi_{i,j})\,,
\label{eq:dvine}
\end{equation}
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~\cite{smithvahey2016}. This also gives an expression for
the product at~\eqref{eq:dvine} 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~\eqref{eq:copden}.
\setcounter{section}{0}