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.
189,514 characters · 37 sections · 90 citation commands
Fast and Accurate Variational Inference for Models with Many Latent Variables
\pagestyle{empty}
\pagestyle{plain} \doublespacing {0.15cm} {0.15cm}
Models with large numbers of latent variables\footnote{In the machine learning literature latent variables are sometimes called “local variables” and model parameters are sometimes called “global variables”; for example, see hoffman2013.} are increasingly popular for capturing nuances in large or complex datasets. Examples include topic models hoffman2013, state space models durbin2012time, choice models train2009 and mixed models gelman2006, among others. In some cases it is possible to integrate out the latent variables, although frequently it is not. Instead, Bayesian estimation usually proceeds by considering the joint posterior distribution of the model parameters and latent variables, often called an “augmented posterior”. Conventional Markov chain Monte Carlo (MCMC) methods for evaluating the augmented posterior are computationally burdensome, and optimization-based variational inference ormerod2010explaining,blei2017 is a scalable alternative. However, existing variational methods make strong independence or other parametric assumptions which often result in a poor approximation to the augmented posterior. In this paper we suggest a general variational inference method which provides a more accurate approximation to the augmented posterior, while being scalable to the case where there are a large number of latent variables and global model parameters.
Key to effective variational inference (also called “variational Bayes”) is the selection of a suitable variational approximation (VA). We suggest a VA which uses an arbitrary parametric approximation for the model parameters, along with the exact conditional posterior distribution for the latent variables. We show that this approximation is tractable, and derive a method for estimating the re-parameterization gradient of the variational lower bound Kingma+w13,rezende+mw14, which is the main input to efficient stochastic gradient ascent (SGA) methods widely used for calibrating VAs bottou10,hoffman2013,salimans2013fixed. Importantly, our method does not require evaluation of the conditional posterior density of the latent variables, or its derivative, which are often unavailable or intractable. Instead, only a draw from the conditional posterior (either exactly or approximately) is needed, which is typically straightforward using a variety of well-explored methods for most models. Thus, our approach can combine the strengths of recent developments in MCMC and variational inference.
In some applications poor parameter uncertainty quantification may not harm predictive inference wang+b19. However, this is not the case for many latent variable models, including the two examples considered here. Although marginal posterior distributions of global model parameters are often well approximated as Gaussian for large data sets titsias2014doubly,archer2015,kucukelbir2016automatic,ong2017gaussian, observation-specific latent variables about which there is little information may exhibit highly non-Gaussian marginal posteriors. Furthermore, when the number of observation-specific latent variables grows with the sample size, poor uncertainty quantification for the latent variables can cause poor inference for the global model parameters, including inaccurate point estimates. A review of the existing literature on variational inference for complex latent variable models is given later in Section (ref).
To illustrate the advantages of our approach we employ it to estimate two contemporary and complex econometric models. The first is a time-varying parameter vector autoregression with stochastic volatility (TVP-VAR-SV) proposed by huber2020. This is a nonlinear state space model, for which mean field approximations are known to be poor wang2004, and Gaussian VAs with parsimonious structured covariance matrices are the most popular choice quiroz+nk18. Because this is a setting where the exact posterior can be computed using existing MCMC methods, the accuracy of our VA can be assessed. We find that our method provides an approximation that is close to the exact posterior and is much more accurate than even a well-structured Gaussian VA.
The second example is a random coefficient tobit (i.e. censored regression) model applied to a large consumer panel originating from danaher2020. Random coefficient models are widely used in marketing to capture consumer-level heterogeneity allenby1998 and in these models estimates of the random coefficients are central as they are employed to tailor advertising and promotions at the consumer level. It is impractical to employ MCMC methods to estimate such models for large datasets, so that recent marketing studies employ variational inference instead ansari2018,danaher2020. In their original study danaher2020 employ a structured Gaussian VA, and we show that this is less accurate than our proposed VA for a small sample of 1,000 consumers, which is a case where the exact posterior can be computed. However, evaluation of the exact posterior is impractical for a large sample of 20,000 individuals, yet we show that our approach is readily applied to this case and more effective than the Gaussian VA.
We also show how our proposed VA is suitable for sub-sampling variational inference for augmented posteriors that can be factorized. This is common in practice, and includes the random coefficient tobit model considered here. Application of sub-sampling further reduces the computational burden of calibrating our proposed VA for this example.
In the next section we introduce variational inference, outline our new methodology, and discuss related existing methods. Section (ref) applies our approach to the TVP-VAR-SV model, and Section (ref) to the mixed effects tobit model. Section (ref) outlines sub-sampling variational inference for our approximation, while Section (ref) concludes. Appendices (ref) and (ref) provide a proof and implementation details for our method and Appendix (ref) gives key derivations for the TVP-VAR-SV model. An extensive Online Appendix provides additional details and results for the two applications.
We first provide a short overview of variational inference, and then outline our new family of VAs for models with latent variables, along with how they can be used to provide variational inference.
We consider Bayesian inference with data $\text{\boldmath$y$}$ having density $p(\text{\boldmath$y$}|\text{\boldmath$\psi$})$, where in our paper, $\text{\boldmath$\psi$}$ contains the model parameters augmented with a potentially large number of latent variables. Assuming a prior density $p(\text{\boldmath$\psi$})$, Bayesian inference is based on the density $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})\propto p(\text{\boldmath$\psi$})p(\text{\boldmath$y$}|\text{\boldmath$\psi$})$, which in our paper is the augmented posterior. We will consider variational inference methods, in which a member $q_\lambda(\text{\boldmath$\psi$})$ of some parametric family of densities is used to approximate the target $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})$, where $\text{\boldmath$\lambda$}\in \Lambda$ is a vector of variational parameters. Approximate Bayesian inference is then formulated as an optimization problem, where a measure of divergence between $q_\lambda(\text{\boldmath$\psi$})$ and $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})$ is minimized with respect to $\text{\boldmath$\lambda$}$. The Kullback-Leibler divergence is typically used, and it is straightforward to show (see, for example, ormerod2010explaining) that this is equivalent to maximizing the variational lower bound (also called the evidence lower bound, or “ELBO”) given by \[ \mathcal L(\text{\boldmath$\lambda$})=\int \log \frac{p(\text{\boldmath$\psi$})p(\text{\boldmath$y$}|\text{\boldmath$\psi$})}{q_\lambda(\text{\boldmath$\psi$})} q_\lambda(\text{\boldmath$\psi$}) d\text{\boldmath$\psi$}\,. \]
The lower bound takes the form of an intractable integral, so it seems challenging to optimize. However, notice that it can be written as an expectation with respect to $q_\lambda$ as
where $g(\text{\boldmath$\psi$})=p(\text{\boldmath$\psi$})p(\text{\boldmath$y$}|\text{\boldmath$\psi$})$. This expression allows easy application of stochastic gradient ascent (SGA) methods bottou10. In SGA, given an initial value $\text{\boldmath$\lambda$}^{(1)}$, $\text{\boldmath$\lambda$}$ is updated recursively as
where $\text{\boldmath$\rho$}_i=(\rho_{i1},\dots, \rho_{im})^\top$ is a vector of step sizes, `$\circ$' denotes the element-wise product of two vectors, and $\widehat{\nabla_\lambda \mathcal L(\text{\boldmath$\lambda$}^{(i)})}$ is an unbiased estimate of the gradient of $\mathcal L(\text{\boldmath$\lambda$})$ evaluated at $\text{\boldmath$\lambda$}=\text{\boldmath$\lambda$}^{(i)}$. For appropriate step size choices this will converge to a local optimum of $\mathcal L(\text{\boldmath$\lambda$})$. Adaptive step size choices are often used in practice, and we use the popular automatic ADADELTA method of zeiler2012adadelta.
To implement SGA, unbiased estimates of the gradient of the lower bound are the key requirement. These can be obtained directly by differentiating (ref) with respect to $\text{\boldmath$\lambda$}$, and evaluating the expectation in a Monte Carlo fashion by simulating from $q_\lambda$. However, variance reduction methods for the gradient estimation are often also important for fast convergence and stability. One of the most useful is the “re-parametrization trick” Kingma+w13. In this approach, it is assumed that an iterate $\text{\boldmath$\psi$}$ can be generated from $q_\lambda$ by first drawing $\text{\boldmath$\varepsilon$}$ from a density $f_\varepsilon$ which does not depend on $\text{\boldmath$\lambda$}$, and then transforming $\text{\boldmath$\varepsilon$}$ by a deterministic function $\text{\boldmath$\psi$}=h(\text{\boldmath$\varepsilon$},\text{\boldmath$\lambda$})$. From (ref), the lower bound can be written as the following expectation with respect to $f_\varepsilon$:
Differentiating under the integral sign in ((ref)) gives the “re-parameterization gradient”
and approximating the expectation at ((ref)) by one or more draws from $f_\varepsilon$ gives an unbiased estimate of $\nabla_\lambda \mathcal L(\text{\boldmath$\lambda$})$. An intuitive reason for the success of the re-parameterization trick is that it allows gradient information from the log-posterior to be used, by moving the variational parameters inside $g(\text{\boldmath$\psi$})$ in ((ref)). For a well-chosen VA, only one draw from $f_\varepsilon$ is typically sufficient for the SGA to converge reliably. We employ the re-parameterization trick throughout, and specify $h$ in Sections (ref) and (ref).
In this paper we consider the case where $\text{\boldmath$\psi$}^\top=(\text{\boldmath$\theta$}^\top,\text{\boldmath$z$}^\top)$ is a parameter vector $\text{\boldmath$\theta$}=(\theta_1,\ldots,\theta_m)^\top$ augmented with additional latent variables $\text{\boldmath$z$}$. Examples include state space models where $\text{\boldmath$z$}$ are latent states durbin2012time, discrete choice models where $\text{\boldmath$z$}$ are latent utilities train2009 and mixed models where $\text{\boldmath$z$}$ are the realizations of random coefficients gelman2006. The prior $p(\text{\boldmath$\psi$})=p(\text{\boldmath$z$}|\text{\boldmath$\theta$})p(\text{\boldmath$\theta$})$, where $p(\text{\boldmath$z$}|\text{\boldmath$\theta$})$ is known from the model specification, and we approximate the augmented posterior density $p(\text{\boldmath$\psi$}|\text{\boldmath$y$})=p(\text{\boldmath$\theta$},\text{\boldmath$z$}|\text{\boldmath$y$})\propto p(\text{\boldmath$y$}|\text{\boldmath$\theta$},\text{\boldmath$z$})p(\text{\boldmath$z$}|\text{\boldmath$\theta$})p(\text{\boldmath$\theta$})\equiv g(\text{\boldmath$\theta$},\text{\boldmath$z$})$ with VAs of the form
As discussed below, it must be feasible to generate directly or approximately from $p(\text{\boldmath$z$}|\bm{y},\bm{\theta})$, although it is unnecessary to evaluate this density or its derivatives. The density $q^0_\lambda$ is chosen to be analytically tractable and from which it is convenient to sample, and we outline two effective choices in Section (ref). Calibration of $q_\lambda$ at (ref) has the potential of being much more efficient than exact sampling using MCMC, because at each VB step we draw the vector $\bm{\theta}$ jointly from $q^0_\lambda$, and do not require partitioning of $\text{\boldmath$\theta$}$ or Metropolis-Hastings steps as in many MCMC schemes.
A key reason why (ref) is an attractive choice can be seen by evaluating the lower bound (ref) as
From Bayes theorem $p(\text{\boldmath$z$}|\bm{y},\bm{\theta})= p(\bm{y}|\text{\boldmath$z$},\bm{\theta})p(\text{\boldmath$z$}|\bm{\theta})/p(\bm{y}|\bm{\theta})$, substituting this into (ref) gives
Here, $\mathcal L^0(\text{\boldmath$\lambda$})$ is the variational lower bound arising from approximating the posterior $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ directly by the VA $q^0_\lambda$. Thus, maximizing $\mathcal L(\text{\boldmath$\lambda$})$ using SGA methods is equivalent to maximizing $\mathcal L^0(\text{\boldmath$\lambda$})$ for the posterior of $\text{\boldmath$\theta$}$ with $\text{\boldmath$z$}$ marginalized out exactly, yet avoids the computation of the (often intractable) density $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ and its derivative.
A second major advantage of the VA at (ref) is that the gradient of the lower bound at (ref) has a simplified expression when using the re-parameterization trick, as summarized in the theorem below.
{\em Proof}: See Appendix (ref).
In (ref), the term $\nabla_\theta \log p(\text{\boldmath$z$}|\bm{\theta},\bm{y})$ is not needed, nor are derivatives with respect to $\text{\boldmath$z$}$, greatly simplifying calculation of the re-parameterization gradient. Instead, only a draw from the conditional posterior $p(\text{\boldmath$z$}|\text{\boldmath$\theta$},\text{\boldmath$y$})$ is required. There is a large literature on drawing either exactly, or approximately, from this distribution for a wide range of latent variable models using filtering, particle, MCMC, Hamiltonian Monte Carlo or other methods. In the applications considered here, we use either a single, or a small number, of sweeps from a Gibbs sampler initialized at the draw from the previous SGA step, which proves effective and simple to implement.
A third reason why (ref) is an attractive choice is because reducing the approximation error due to $\text{\boldmath$z$}$ also improves calibration of the marginal VA of $\text{\boldmath$\theta$}$, as summarized in the corollary below.
{\em Proof}: See Appendix (ref).
Algorithm (ref) calibrates our proposed VA to the augmented posterior using SGA with the re-parameterization trick and the ADADELTA learning rate.
The final ingredient of our VA at (ref) is $q_\lambda^0(\text{\boldmath$\theta$})$, along with a matching re-parameterization transformation $\text{\boldmath$\theta$}=h^0(\text{\boldmath$\varepsilon$}^0,\text{\boldmath$\lambda$})$. An advantage of our approach is that any existing variational family can be used for $q^0_\lambda$. Popular candidates include mean field Gaussian or elliptical distribution approximations, although we stress this is not the same as employing such approximations for the entire augmented posterior $p(\text{\boldmath$\theta$},\text{\boldmath$z$}|\text{\boldmath$y$})$. Below we discuss two flexible choices for $q^0_\lambda$ that are generic in that they do not exploit any model-specific properties of the posterior distributions. They are more accurate than mean field approximations, have generative representation allowing easy computation of the gradient at (ref), and involve computations that increase only linearly with the dimension of $\text{\boldmath$\theta$}$.
Gaussian VAs are popular, but can be computationally burdensome or inaccurate when an unrestricted covariance matrix is employed and the dimension $m$ is high. miller2017 and ong2017gaussian use a factor covariance structure to reduce the number of variational parameters. Here, $q^0_\lambda(\text{\boldmath$\theta$})=\phi_m(\text{\boldmath$\theta$};\text{\boldmath$\mu$}_\theta,B_\theta B_\theta^\top+D_\theta^2)$, where $B_\theta$ is an $(m\times k)$ matrix with $k<<m$ and zero upper triangular elements, $D_\theta=\mbox{diag}(\text{\boldmath$d$}_\theta)$ is a diagonal matrix, $\phi_m(\cdot;\text{\boldmath$a$},A)$ is the density of a $m$-dimensional $N(\text{\boldmath$a$},A)$ distribution, and $\text{\boldmath$\lambda$}=(\text{\boldmath$\mu$}_\theta^\top,\mbox{vech}(B_\theta)^\top,\text{\boldmath$d$}_\theta^\top)^\top$. The “vech” operator is the half-vectorization of a rectangular matrix (i.e. the vectorization of the non-zero elements). The distribution has generative representation $\text{\boldmath$\theta$}=\text{\boldmath$\mu$}_\theta+B_\theta \text{\boldmath$\zeta$}_1 + \text{\boldmath$d$}_\theta \circ \text{\boldmath$\zeta$}_2$, where $\text{\boldmath$\varepsilon$}^0=(\text{\boldmath$\zeta$}_1^\top,\text{\boldmath$\zeta$}_2^\top)^\top\sim N(\bm{0},I_{k+m})$, which defines the re-parameterization transformation $h^0$ and is used at step (a) of Algorithm (ref). ong2017gaussian derive analytical expressions for $\nabla_{\theta}\log q_\lambda^0(\text{\boldmath$\theta$})$ and $\frac{\partial \text{\boldmath$\theta$}}{\partial \text{\boldmath$\lambda$}}^\top$ that can be used for fast evaluation of the unbiased estimate of the re-parameterization gradient at step (c) of Algorithm (ref). When $k=0$ this approximation reduces to a product of independent Gaussians (i.e. a mean field approximation), while these authors found that setting $k=5$ balanced accuracy and computations for their examples.
Gaussian copula VAs are more flexible because they can accommodate different marginal distributions for each element of $\text{\boldmath$\theta$}$. han2016 suggest such an approximation, but with nonparametric or mixture density margins, which are difficult to calibrate quickly. Instead, we follow Smith2020 and employ a Gaussian copula VA constructed using an element-wise transformation\footnote{This is not to be confused with the transformation $h^0$ associated with the re-parameterization trick.} of $\text{\boldmath$\theta$}$ and a factor decomposition of the copula correlation matrix. We provide an outline here, but refer to the work of these authors for details.
Let $t_{\gamma}:\mathbb{R} \rightarrow \mathbb{R}$ be a family of one-to-one transformations with parameter vector $\text{\boldmath$\gamma$}$. If a parameter $\theta_i$ is constrained we first transform it to the real line; for example, with a scale or variance parameter we set $\theta_i$ to its logarithm. To construct $q^0_\lambda$, we transform each parameter as $\vartheta_i=t_{\gamma_i}(\theta_i)$ and adopt the Gaussian distribution $\text{\boldmath$\vartheta$}=(\vartheta_1,\ldots,\vartheta_m)^\top \sim N(\text{\boldmath$\mu$}_\vartheta,\Sigma_\vartheta)$. The density $q^0_\lambda$ can be recovered by computing the Jacobian of the transformation, so that
where the variational parameters are $\text{\boldmath$\lambda$}^\top=(\text{\boldmath$\gamma$}_1^\top,\ldots,\text{\boldmath$\gamma$}_m^\top,\text{\boldmath$\mu$}_\vartheta^\top, \mbox{vech}(\Sigma_\vartheta)^\top)$ and $t_{\gamma_i}'(\theta_i)=\frac{d\vartheta_i}{d\theta_i}$. Let $\text{\boldmath$\mu$}_{\vartheta}=(\mu_{\vartheta,1},\ldots,\mu_{\vartheta,m})^\top$ and $\sigma_{\vartheta,i}^2$ be the $i$th leading diagonal element of $\Sigma_{\vartheta}$, then the marginal densities of the approximation are
with $\text{\boldmath$\lambda$}_i^\top=(\text{\boldmath$\gamma$}_i^\top,\mu_{\vartheta,i},\sigma_{\vartheta,i}^2)$ a sub-vector of $\text{\boldmath$\lambda$}^\top$. It is straightforward to show that the distribution with density at (ref) has a Gaussian copula. Smith2020 point out that while $q^0_\lambda$ has a copula representation, it is more computationally efficient to utilize that at (ref), and we do so here.
To allow for large $m$, we adopt a factor structure for $\Sigma_\vartheta$ as above. Let $B_\vartheta$ be an $(m\times k)$ matrix with zeros in the upper triangle, $D_\vartheta=\mbox{diag}(\bm{d}_\vartheta)$ a diagonal matrix, then we assume that $\Sigma_\vartheta = B_\vartheta B_\vartheta^\top+D_\vartheta^2$. Thus, the variational parameters are $\text{\boldmath$\lambda$}^\top=(\text{\boldmath$\gamma$}_1^\top,\ldots,\text{\boldmath$\gamma$}_m^\top,\text{\boldmath$\mu$}_\vartheta^\top, \mbox{vech}(B_\vartheta)^\top,\text{\boldmath$d$}_\vartheta^\top)$. We note that this copula is equivalent to the Gaussian factor copula used to model data oh2017, although here it is a VA. The distribution has generative representation $\text{\boldmath$\theta$}=h^0(\text{\boldmath$\varepsilon$}^0,\text{\boldmath$\lambda$})=(t_{\gamma_1}^{-1}(\vartheta_1),\ldots, t_{\gamma_m}^{-1}(\vartheta_m))^\top$, where $\text{\boldmath$\vartheta$}=\text{\boldmath$\mu$}_\vartheta+B_\vartheta \text{\boldmath$\zeta$}_1 + \text{\boldmath$d$}_\vartheta \circ \text{\boldmath$\zeta$}_2$, and $\text{\boldmath$\varepsilon$}^0=(\text{\boldmath$\zeta$}_1^\top,\text{\boldmath$\zeta$}_2^\top)^\top\sim N(\bm{0},I_{k+m})$. Closed form expressions for $\frac{\partial\bm{\theta}}{\partial\bm{\lambda}}$ and $\nabla_\theta\log q^0_\lambda(\bm{\theta})$ required to compute (ref) can be derived as in Appendix (ref).
We employ a transformation suggested by yeojohnson2000 that has proven successful in transforming data to near normality. This extends the Box-Cox transformation to the entire real line and, for $0<\gamma<2$, is given by \[ t_\gamma(\theta)= \left\{
\right. \] When implementing SGA $t_\gamma$ is not evaluated, but its (closed form) inverse $t_\gamma^{-1}(\vartheta)$ is repeatedly. With this transformation, the margins at (ref) are flexible parametric densities that can exhibit skew and/or excess kurtosis, as illustrated in Smith2020. Other transformations may also be used for $t_\gamma$, producing alternative margins.
The alternative of applying generic approximations to the augmented posterior $p(\text{\boldmath$\theta$},\text{\boldmath$z$}|\text{\boldmath$y$})$, or the marginal posterior of the latent variables $p(\text{\boldmath$z$}|\text{\boldmath$y$})$, has been considered previously. For example, braun2010 do so for the posterior of a multinomial logistic regression augmented with random coefficient realizations, hui2017 for latent variables in a generalized linear model, LoaizaSmith2018 for the augmented posterior of a discrete-margined copula model, and archer2015 for the augmented posterior of a state space model, among others. Such approaches require a much larger number of variational parameters, resulting in three drawbacks compared to adopting the VA at (ref). First, additional error is introduced into the variational estimate of $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ through imprecision in the approximation of the posterior of $\text{\boldmath$z$}$. Second, there is an increased computational burden at each step of the SGA algorithm because the re-parameterization gradient is a much larger vector than that at (ref). Last, the SGA algorithm typically requires more steps, because additional noise is introduced into the Monte Carlo estimate of the gradient.
A seminal paper on variational inference for latent variable models is hoffman2013. These authors consider mean field approximations in models with global model parameters and latent variables (which are called “local variables”), and describe how sub-sampling methods can be used in the variational optimization. Their method requires the model to have a conjugate exponential family structure. hoffman+b15 consider structured stochastic variational inference methods where dependence between global parameters and local variables can be accommodated. They consider the possibility of using the exact conditional posterior for the latent variables as part of the approximation, similar to our approach, but unlike our method theirs also requires conjugacy. tan+n14 consider a stochastic variational inference implementation of non-conjugate variational message passing using sub-sampling which is useful for generalized linear mixed models, and some diagnostics for prior-data conflict checking. tan18 consider affine transformations for re-parametrizations of latent variable models which improve accuracy, and nolan+mw19 consider efficient variational inference in models with multi-level random effects structures.
tran2017variational consider variational inference in latent variable models when only an unbiased estimate of the likelihood is available. Their approach requires tuning log-likelihood estimates to achieve constant variance, and gunawan+tk17 overcome this disadvantage using sub-sampling methods and re-parametrization gradient estimates based on Fisher's identity. Their method is related to ours, but our work differs from theirs in two main ways. First, their applications to latent variable models are limited to random effects models for panel data, and we extend to more general latent variable models such as state space models. Second, gunawan+tk17 use importance sampling for integrating out any latent variables in gradient estimation, whereas here we use MCMC methods. The use of importance sampling limits applications to models with low-dimensional latent variables, unlike with our approach. For example, the panel data application in gunawan+tk17 uses a logistic model with a scalar intercept random effect, whereas later we consider a tobit model with a larger random effect vector. tomasetti2019 consider a VA similar to (ref), but for a partition of the parameter vector, not for the augmented posterior of a latent variable model. They do not use re-parameterization gradients, nor simulate from the conditional posterior, but instead use an importance sampler for sequential posterior inference.
There have been a number of general efforts to combine variational inference and MCMC. These include salimans2015markov, domke17, li+tl17, zhang2018ergodic and ye+bdh20 among others, although the methods do not focus specifically on computation for latent variable models. The method of ruiz+t19 does consider such models. The VA they consider corresponds to a parametrized approximation to an initial value of an MCMC algorithm which is run for a fixed small number of iterations. For the case of latent variable models, scalable amortized variational inference methods can be used in learning the variational parameters for the latent variables. hoffman17 considers SGA algorithms in which gradient estimates based on short MCMC runs starting from samples from a VA are considered. The approach is used for marginal maximum likelihood computations in deep Gaussian latent variable models.
To illustrate our method, we use it to estimate a time-varying parameter vector autoregression with stochastic volatility. This complex time series model is a leading model for forecasting macroeconomic variables clark2015,carriero2019 and we employ a variant suggested by huber2020 where a horseshoe prior is used to regularize the time-varying parameters. This is a challenging example for two reasons. First, the TVP-VAR-SV model is an example of a nonlinear state space model, for which mean field approximations of the augmented posterior are known to be poor wang2004,karl2016. To address this, for state space models previous authors have used structured Gaussian or other parametric VAs to the augmented posterior or marginal posterior of the latent states ghahramani2000,daunizeau2009,archer2015,naesseth2017. In contrast, our approach has a smaller number of variational parameters, involves less computation and provides greater accuracy than these approaches. The second challenge is that horseshoe regularization produces posterior densities that are funnel-shaped and difficult to approximate. To solve this problem we adopt the re-parameterization of ingraham2017.
Let $\{\text{\boldmath$y$}_t\}_{t=1}^T$ be a time series of $N$ macroeconomic variables with $\text{\boldmath$y$}_t=(y_{1,t},y_{2,t},\ldots,y_{N,t})^\top$. We employ the following conditionally Gaussian TVP-VAR-SV model:
Here, $\text{\boldmath$\epsilon$}_t=(\epsilon_{1,t},\ldots,\epsilon_{N,t})^\top\sim N(\bm{0},H_t)$, $H_t=\mbox{diag}(e^{h_{1,t}},\ldots,e^{h_{N,t}})$ is a diagonal matrix, $L_t$ is a lower triangular matrix with unit-valued leading diagonal (i.e. a “unitriangular” matrix), $\text{\boldmath$\beta$}_{0,t}$ is an intercept vector and $B_{1,t},\ldots,B_{p,t}$ are $(N\times N)$ autoregressive parameter matrices. The intercept and autoregressive parameters follow a random walk, with the $(pN^2+N)$-dimensional vector $\text{\boldmath$\beta$}_t^\top\equiv (\text{\boldmath$\beta$}_{0,t}^\top,\mbox{vec}(B_{1,t})^\top,\ldots,\mbox{vec}(B_{p,t})^\top)$ and $\text{\boldmath$w$}_t\sim N(0,V)$ with $V=\mbox{diag}(v_1,\ldots,v_K)$. The $N(N-1)/2$ free elements of $L_t$ are assumed to follow independent random walks. The terms $h_{i,1},\ldots,h_{i,T}$ are the logarithm of the volatilities for the $i$th variable, and these follow a stationary first order autoregression with mean $\bar{h}_i$, autoregressive parameter $-1<\rho^h_i<1$ and independent disturbances $e_{i,t}\sim N(0,\sigma_i^2)$.
Likelihood-based estimation using the parameterization in (ref) is difficult, and following huber2020 and others, it is convenient to both transform the system and adopt a non-centered parameterization. As outlined in Appendix (ref), this results in $i=1,\ldots,N$ unrelated regressions
where $\text{\boldmath$\alpha$}_i$ and $\text{\boldmath$x$}_{i,t}$ are vectors of dimension $J_i=2(pN+i)$. Each equation can then be estimated separately. The vector $\text{\boldmath$x$}_{i,t}$ is a function of both the observed times series values and time-varying latent variables $\widetilde \text{\boldmath$\eta$}_{i,t}$ specified at (ref) in Appendix (ref). huber2020 compare different priors for the regularization of $\text{\boldmath$\alpha$}_i=(\alpha_{i,1},\ldots,\alpha_{i,J_i})^\top$, including the horseshoe prior of carvalho2010 which we also use here. For each element $\alpha_{i,j}$, the horseshoe is specified by adopting the hyper-priors \[ \alpha_{i,j}|\xi_i,\chi_{i,j}~\sim N(0,\xi_i\chi_{i,j}),\ \ \chi_{i,j}|\nu_{i,j}\sim \mathcal{G}^{-1}\left(\frac{1}{2},\frac{1}{\nu_{i,j}}\right),\ \ \xi_i|\kappa_{i}\sim \mathcal{G}^{-1}\left(\frac{1}{2},\frac{1}{\kappa_{i}}\right),\ \ \nu_{i,1},\dots,\nu_{i,J_i},\kappa_{i}\sim \mathcal{G}^{-1}\left(\frac{1}{2},1\right), \] where $\mathcal{G}^{-1}(a,b)$ denotes an inverse gamma distribution with parameters $a,b$. This provides both global and local shrinkage of the coefficients in each vector $\text{\boldmath$\alpha$}_i$.
For the $i$th equation, let $\text{\boldmath$\chi$}_i=(\chi_{i,1},\ldots,\chi_{i,J_i})^\top$ and $\text{\boldmath$\nu$}_i=(\nu_1,\ldots,\nu_{J_i})^\top$. Then the parameter vector is $\bm{\theta}_i=\left(\bm{\alpha}_i^\top,\bm{\chi}_i^\top,\xi_i,\bm{\nu}_i^\top,\kappa_i,\bar{h}_i,\rho^h_i,\sigma_i^2\right)^\top$, which has a total of $3J_i+5$ elements. The latent variables consist of the $T$-dimensional vector of log-volatilities $\text{\boldmath$h$}_i=(h_{i,1},\ldots,h_{i,T})^\top$ and the $(TJ_i/2)$-dimensional vector $\widetilde \text{\boldmath$\eta$}_i=(\widetilde \text{\boldmath$\eta$}_{i,1}^\top,\ldots,\widetilde \text{\boldmath$\eta$}_{i,T}^\top)^\top$. While the likelihood is intractable because the latent variables cannot be integrated out analytically, the augmented posterior can be expressed in closed form. Let $\text{\boldmath$y$}_{(i)} \equiv (y_{i,1},\ldots,y_{i,T})^\top$ be the observations on the $i$th macroeconomic variable, and $\text{\boldmath$y$}_{(\backslash i)}$ be the observations on the other $N-1$ macroeconomic variables, then the augmented posterior is
where observe that $p(\text{\boldmath$h$}_i,\widetilde \text{\boldmath$\eta$}_i|\text{\boldmath$\theta$}_i)=p(\text{\boldmath$h$}_i|\text{\boldmath$\theta$}_i)p(\widetilde \text{\boldmath$\eta$}_i|\text{\boldmath$\theta$}_i)$. The remaining priors are as specified in huber2020, with $\bar{h}_i\sim N(0,10^2)$, the beta prior $\frac{\rho^h_i+1}{2}\sim\mathcal{B}(25,5)$ and $\sigma_i^2\sim\mathcal{G}\left(\frac{1}{2},\frac{1}{2}\right)$. These authors provide an MCMC scheme to compute the augmented posterior above.
The posterior densities of parameters that are regularized using the horseshoe prior have pathological funnel-shaped geometries, making them difficult to approximate betancourt2015,ghosh2019. ingraham2017 suggest a “non-centered” re-parameterization of $\text{\boldmath$\alpha$}_i$ that simplifies the posterior densities so that they are easier to approximate, where each element $\alpha_{i,j}=\tau_{i,j}\sqrt{\xi_i\chi_{i,j}}$ with $\tau_{i,j}\sim N(0,1)$. This re-parameterization is unnecessary when using MCMC, although doing so may affect sampling efficiency. Let $\text{\boldmath$\tau$}_{i}=(\tau_{i,1},\ldots,\tau_{i,J_i})^\top$, so that $\text{\boldmath$\alpha$}_i=\sqrt{\xi_i}(\text{\boldmath$\tau$}_i \circ \sqrt{\text{\boldmath$\chi$}_i})$, then the model parameters that we employ for variational inference are $\bm{\theta}_i=\left(\bm{\tau}_i^\top,\bm{\chi}_i^\top,\xi_i,\bm{\nu}_i^\top,\kappa_i,\bar{h}_i,\rho^h_i,\sigma_i^2\right)^\top$ which remains of length $3J_i+5$. Here, $\text{\boldmath$\chi$}_{i}$, $\xi_i$, $\text{\boldmath$\nu$}_{i}$, $\sigma^2_i$ and $\kappa_i$ are all transformed to their logarithms, and $\rho^h_i$ to $\Phi_1^{-1}\left(\frac{\rho^h_i+1}{2}\right)$, so that all parameters are unconstrained on the real line.
Let $\text{\boldmath$z$}\equiv(\text{\boldmath$h$}_i,\widetilde \text{\boldmath$\eta$}_i)^\top$, then we approximate the augmented posterior at (ref) using our suggested VA at (ref) with $q^0_\lambda$ the Gaussian factor model at Section (ref) with $k=5$ factors. The variational parameters $\text{\boldmath$\lambda$}$ are calibrated using Algorithm (ref).\footnote{We apply Algorithm (ref) to estimate each equation separately, so that the parameter vector is $\text{\boldmath$\theta$}_i$ in the variational approximation at (ref). However, we continue to denote the latent variable vector as $\text{\boldmath$z$}$ and variational parameters as $\text{\boldmath$\lambda$}$ to be consistent with the notation in Section (ref).} At step (b), drawing directly from $p(\text{\boldmath$z$}|\text{\boldmath$\theta$}_i^{(s)},\text{\boldmath$y$})$ is infeasible. Therefore, we draw $\text{\boldmath$z$}$ by generating between one and five sweeps from the densities $p(\text{\boldmath$h$}_i|\bm{y},\bm{\theta}_i^{(s)},\widetilde{\text{\boldmath$\eta$}}_i^{(s-1)})$ and $p(\widetilde{\text{\boldmath$\eta$}}_i|\bm{y},\bm{\theta}_i^{(s)},\bm{h}_i^{(s)})$ using standard state space methods detailed in Appendix (ref) and initialized at the value of $\text{\boldmath$z$}$ from the previous SGA step. While this provides only an approximate draw from $p(\text{\boldmath$z$}|\text{\boldmath$\theta$}_i^{(s)},\text{\boldmath$y$})$, it is fast and we find it provides accurate estimates as documented below, although other approaches to generating a draw may also be used. When computing the re-parameterization gradient at (ref), $\nabla_{\theta_i} \log g(\text{\boldmath$\theta$}_i,\text{\boldmath$z$})$ is available in closed form; see the Online Appendix.
The MCMC estimator of huber2020 is used with their recommended number of draws (15,000 burn-in and 15,000 Monte Carlo sample) to compute the “exact” posterior, and we employ it to judge the accuracy of our VA. For comparison we also calibrate the following Gaussian VA
where $K_i=J_i(3+T/2)+5$, $C_{G,h}$ is a band one lower triangular Cholesky factor, and a factor decomposition with five factors is used for $\Sigma_G$. The banded structure of the covariance matrix $C_{G,h}C_{G,h}^\top$ can better capture the Markov dependence in $\text{\boldmath$h$}_i$ than a factor decomposition; see quiroz+nk18 for a similar structured Gaussian VA.
The TVP-VAR-SV with $p=2$ lags is used to model the same data and eight macroeconomic variables used by huber2020 that these authors label the “medium case”. The data are quarterly observations from 1980:Q3 to 2017:Q4, so that $T=150$. In the eight equations, the number of parameters $|\text{\boldmath$\theta$}_i|$ is between 107 and 149, while the number of latent variables $|\text{\boldmath$z$}|$ is between 2700 and 3750. The purpose of this empirical analysis is to demonstrate the high degree of accuracy of our proposed VA.
Figure (ref) plots the exact posterior mean of the time-varying standard deviations $\exp(h_{i,t}/2)$ for each of the eight series, along with the mean of $\exp(h_{i,t}/2)$ for both calibrated VAs. Our proposed VA (labeled “Hybrid VA”) is more accurate than the structured Gaussian VA. Figure (ref) plots the exact posterior densities of $(\rho_i^h,\bar{h}_i,\sigma^2_i)$, along with those of the two VAs, for the variable Real GDP (equation $i=1$). Again, the densities from the structured Gaussian VA are poor, whereas those from our method are accurate. The Online Appendix provides the equivalent densities for the other seven variables, and very similar results are observed. Note that an accurate approximation for $\text{\boldmath$h$}_i$ is necessary to also obtain an accurate approximation of $(\rho_i^h,\bar{h}_i,\sigma^2_i)$, and vice-versa.
Figure (ref) plots the exact posterior of a subset of the time-varying autoregressive coefficients in black. These are the first rows of $B_{1,t}$ and $B_{2,t}$ which correspond to the first equation (for Real GDP). Both the posterior mean and the upper/lower 90% intervals are plotted. Some coefficients (e.g. those in panels (f,g,m,n,o)) are heavily regularized to zero, whereas others vary substantially over time (e.g. those in panels (b,c,i,j,k,p)). The equivalent mean and intervals are also plotted for the structured Gaussian VA (in red) and our hybrid VA (in yellow). The latter is much more accurate than the former. This posterior output for the VAs was computed using Monte Carlo draws of $(\text{\boldmath$\theta$}_i,\text{\boldmath$z$})$ from the calibrated VAs, and from which draws of $B_{1,t}$ and $B_{2,t}$ were evaluated using the relationships in Appendix (ref). Estimates of the time-varying lower unitriangular matrix $L_t$ (see the Online Appendix) also show that our hybrid VA is accurate, and the structured Gaussian VA much less so.
Last, we explore the impact of drawing $\text{\boldmath$z$}$ using multiple sweeps of a Gibbs sampler at step (b) of Algorithm (ref). Figure (ref) plots the means of the VAs of $\exp{(h_{1,t}/2)}$ against their exact posterior means when 1, 5 and 10 sweeps of a Gibbs sampler are used at step (b). While a greater number of sweeps results in increased accuracy, it is only a minor improvement. Similar results are observed for the other seven equations.
It is typical to judge the accuracy and calibration speed of different VAs using the lower bound ${\cal L}$. However, this is unavailable for the VA at (ref) because the density $p(\text{\boldmath$z$}|\text{\boldmath$\theta$}_i,\text{\boldmath$y$})$ is intractable. Therefore, we instead compute the accuracy of the one-step-ahead posterior predictive densities estimated using each VA. Let $\text{\boldmath$z$}_t\equiv (\widetilde \text{\boldmath$\eta$}_{i,t},h_{i,t})$ and $\text{\boldmath$y$}_{1:t}=(\text{\boldmath$y$}_1^\top,\ldots,\text{\boldmath$y$}_t^\top)^\top$, then the predictive density is computed for the TVP-VAR-SV model as
Here, the integral in $\widetilde{\text{\boldmath$\eta$}}_{i,t+1}$ is evaluated analytically by recognising a Gaussian density as in Appendix (ref), and the integral in $h_{i,t+1}$ is computed numerically. The accuracy of these predictive densities is measured using their average Kullback-Leibler divergence over time when computed from the VA and also from the exact posterior as follows. Let $\bar{\text{\boldmath$\theta$}}_i=E(\text{\boldmath$\theta$}_i|\text{\boldmath$y$})$ and $\bar{\text{\boldmath$z$}}_{t}=E(\text{\boldmath$z$}_{t}|\text{\boldmath$y$})$ be the exact posterior means (computed using MCMC), and $\widetilde{\text{\boldmath$\theta$}}_i,\widetilde{\text{\boldmath$z$}}_t$ be the means from the VA $q_\lambda$. The average Kullback-Leibler divergence is then $\overline{\mbox{KL}}(\text{\boldmath$\lambda$})=\frac{1}{T}\sum_{t=1}^T\mbox{KL}_{t+1|t}(\text{\boldmath$\lambda$})$, where \[ \mbox{KL}_{t+1|t}(\text{\boldmath$\lambda$})=\int p_{t+1|t}(y_{i,t+1}|\widetilde{\text{\boldmath$\theta$}}_i,\widetilde{\text{\boldmath$z$}}_{t}) \log\left(\frac{p_{t+1|t}(y_{i,t+1}|\widetilde{\text{\boldmath$\theta$}}_i,\widetilde{\text{\boldmath$z$}}_{t})}{p_{t+1|t}(y_{i,t+1}|\bar{\text{\boldmath$\theta$}}_i,\bar{\text{\boldmath$z$}}_{t})}\right)\mbox{d}y_{i,t+1}\,, \] is computed by numerical integration. Lower values of $\overline{\mbox{KL}}(\text{\boldmath$\lambda$})$ suggest the VA $q_\text{\boldmath$\lambda$}$ has increased accuracy, with $\overline{\mbox{KL}}(\text{\boldmath$\lambda$})=0$ when the VA is exact (i.e. when $q_\text{\boldmath$\lambda$}(\text{\boldmath$\theta$}_i,\text{\boldmath$z$})=p(\text{\boldmath$\theta$}_i,\text{\boldmath$z$}|\text{\boldmath$y$})$).
We compute this diagnostic for the variable Real GDP (i.e. equation $i=1$) although note that results for the other seven variables are similar. The minimum value of $\overline{\mbox{KL}}(\text{\boldmath$\lambda$})$ was 0.0578 and 0.0282 for the Gaussian and our hybrid VA, respectively, so that our VA is substantially more accurate than the benchmark by this metric, consistent with the other empirical results given above. Last, Figure (ref) plots $\overline{\mbox{KL}}(\text{\boldmath$\lambda$})$ for our hybrid VA with 1 and 5 sweeps of the sampler at step (b) of Algorithm (ref), and also the structured Gaussian VA. Panel (a) plots the divergence against walk clock time, and panel (b) against step number of the SGA, with all methods coded in MATLAB and run on a standard laptop. A single sweep of the sampler appears sufficient to calibrate the posterior predictive densities well by this metric. Calibration of the VA takes about 30s for one sweep and 180s for five sweeps, whereas in comparison the MCMC sampler took 246s using the same computing environment.
The second example applies our variational estimator to a mixed effects tobit model for $T=100$ weekly sales by members of a large consumer panel of U.S. customers from the field of marketing. The panel originates from danaher2020, who create a rich set of covariates that they model with both fixed and individual-level random coefficients. The authors point out that exact Bayesian estimation using MCMC is computationally infeasible for large numbers of customers from this panel, and use variational inference with a structured Gaussian VA as in ong2017gaussian to fit models for 4,000 individuals. Using this panel, we first show that for a small sample of 1,000 individuals (for which MCMC can be used), our proposed VA is very close to the exact posterior and also more accurate than either the Gaussian VA used by danaher2020 or a mean field Gaussian VA. We then apply our approach to a large sample of 20,000 customers (so that the number of weekly sales observations is two million) and show it improves upon both Gaussian VAs.
The response $y_{i,t}$ for individual $i = \{1,\dots,N\}$ at week $t = \{1,\dots,T\}$ in a tobit model is an observation of a latent variable $y^\star_{i,t}$ censored at zero, so that $y_{i,t}=y^{*}_{i,t}$ if $ y^{*}_{i,t}>0$, and $y_{i,t}=0$ if $y^{*}_{i,t}\leq 0$. The latent response follows a Gaussian mixed effects model \[ y^{*}_{i,t} = \bm{x}_{i,t}^\top\bm{\beta}+\text{\boldmath$w$}_{i,t}^\top\bm{\alpha}_i+\sigma\epsilon_{i,t}\,,\;\;\epsilon_{i,t}\sim N(0,1)\,,\;\; \text{\boldmath$\alpha$}_i\sim N(\bm{0},V_\alpha)\,, \] where $\text{\boldmath$x$}_{i,t}$ is a $(p\times 1)$ vector of fixed effect covariates, $\text{\boldmath$w$}_{i,t}$ is a $(r\times 1)$ vector of random effect covariates that is a sub-vector of $\text{\boldmath$x$}_{i,t}$. Bayesian analysis of panel data in marketing is popular because estimates of random coefficient values are often a key output allenby1998,manchanda2004response, as is the case here with $\text{\boldmath$\alpha$}=(\text{\boldmath$\alpha$}_1^\top,\ldots,\text{\boldmath$\alpha$}_N^\top)^\top$. A Bayesian analysis requires specification of priors. Here we set $V_\alpha=LL^\top+\Omega$, with $L$ a $r\times k_\alpha$ factor loading matrix (with zeros in the upper triangle and positive leading diagonal elements for identification) and the diagonal matrix $\Omega=\mbox{diag(\text{\boldmath$\omega$})}$, and then adopt independent uninformative priors for $\mbox{vech}(L),\text{\boldmath$\omega$},\text{\boldmath$\beta$},\sigma^2$ as outlined in the Online Appendix.
Let $\text{\boldmath$y$}$ be the vector of the $n=NT$ values of $y_{i,t}$, and $\text{\boldmath$y$}^\star_{{\tiny U}}$ be the vector of values of $y^\star_{i,t}$ which are unobserved (i.e. those where $y_{i,t}=0$), and set $\eta_{i,t}=\bm{x}_{i,t}^\top\bm{\beta}+\text{\boldmath$w$}_{i,t}^\top\bm{\alpha}_i$. The likelihood $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$ is intractable, so the focus is often on the posterior augmented with $\text{\boldmath$\alpha$}$, which has density
In many studies, including ours, weekly individual-level sales values $y_{i,t}$ are mostly zero, so that the first product at (ref) is over more terms than the second. For $r$ and/or $N$ large, integration over $\text{\boldmath$\alpha$}$ to obtain the likelihood is computationally difficult, and MCMC methods that simulate $\text{\boldmath$\alpha$}$ are popular.
It is often simpler and faster to consider the posterior augmented with both $\text{\boldmath$\alpha$}$ and $\text{\boldmath$y$}^\star_{\tiny U}$. This can be derived by first noting that $p(\bm{y}^*_{\tiny U}|\bm{y},\bm{\alpha},\bm{\theta})= \prod_{\{i,t|y_{i,t}=0\}}p(y_{i,t}^*|y_{i,t},\bm{\alpha},\bm{\theta})$, with
Then the augmented posterior is
which was computed by substituting in the expressions above. Using (ref) a Gibbs sampler that generates from both $\text{\boldmath$y$}^\star_{\tiny U}$ and $\text{\boldmath$\alpha$}$ and evaluates the augmented posterior exactly is given in the Online Appendix. However, even MCMC methods applied to either (ref) or (ref) become impractical for larger values of $r$ and/or $N$, and danaher2020 propose using variational inference instead.
The model is estimated using our approach with $\text{\boldmath$z$}=\left((\text{\boldmath$y$}^\star_{\tiny U})^\top,\text{\boldmath$\alpha$}^\top\right)^\top$ and parameter vector $\text{\boldmath$\theta$}=(\text{\boldmath$\beta$}^\top,\mbox{vech}(L)^\top,\text{\boldmath$\omega$}^\top,\sigma^2)^\top$, where some elements are transformed to the real line for ready application of the VA as outlined in the Online Appendix. Thus, the VA is to the augmented posterior at (ref), where the Gaussian copula at Section (ref) with $k=10$ factors is used for $q^0_\lambda$. When implementing Algorithm (ref), at step (b) $\text{\boldmath$z$}$ is generated from the conditionals $p(\text{\boldmath$\alpha$}|(\text{\boldmath$y$}^\star_{\tiny U})^{(s-1)},\text{\boldmath$\theta$}^{(s)},\text{\boldmath$y$})$ and $p(\text{\boldmath$y$}^\star_{\tiny U}|\text{\boldmath$\alpha$}^{(s)},\text{\boldmath$\theta$}^{(s)}),\text{\boldmath$y$})$ given in the Online Appendix, and initialized at the iterate from the previous SGA step. This forms a Gibbs sampler to produce a draw from $p(\text{\boldmath$z$}|\text{\boldmath$\theta$}^{(s)},\text{\boldmath$y$})$, and five sweeps are used to produce our results.
Two Gaussian VAs to (ref) are included as benchmarks, which are strong contenders because $\text{\boldmath$y$}_{\tiny U}^*$ (but not $\text{\boldmath$\alpha$}$) is integrated out from the posterior exactly. The first VA is that proposed by danaher2020 with density $q^{G}_\lambda(\text{\boldmath$\theta$},\text{\boldmath$\alpha$})=q^{G,a}_{\lambda^a}(\text{\boldmath$\theta$}) \prod_{i=1}^N q^{G,i}_{\lambda^i}(\text{\boldmath$\alpha$})$. Here, $q^{G,a}_{\lambda^a}$ and $q^{G,i}_{\lambda^i}$ are the densities of Gaussian factor VAs as in ong2017gaussian, and $\text{\boldmath$\lambda$}^\top=\left((\text{\boldmath$\lambda$}^a)^\top,(\text{\boldmath$\lambda$}^1)^\top, \ldots,(\text{\boldmath$\lambda$}^N)^\top\right)$. One factor is used for $q^{G,i}_{\lambda^i}$, while 10 factors are used for $q^{G,a}_{\lambda^a}$ to produce a richer approximation. The second benchmark is a Gaussian mean field VA (i.e. where the VA density is a product of univariate Gaussian densities), which is a popular choice in practice for models with many parameters. Both benchmarks are calibrated using SGA with the re-parameterization trick.
To assess the accuracy of our proposed VA we first consider weekly sales by a sample of $N=1000$ individuals, for which the exact posterior can be calculated in reasonable time (27 hours) using MCMC. There are $p=32$ covariates which are described in the Online Appendix, and include measurements of an individual's exposure to advertisements (ads) in three media (email, catalogs and paid search) to their in-store and online purchases for three retailer-brands (“B1”, “B2” and “B3”) in the clothing category. The response $y_{i,t}$ is the logarithm of in-store spend (plus one) for retailer-brand B1 (called the “focal brand” here) by individual $i$ in week $t$. We employ random coefficients at the individual consumer level for the three media ad exposure variables of the focal brand and the intercept, so that $r=4$.
Figure (ref) presents the posteriors of the elements of $V_\alpha$, where there is a striking difference between the three VAs. Both Gaussian VAs under-estimate the level of posterior variance, and are incorrectly located for some elements. In contrast, our proposed hybrid VA estimates the posteriors very well-- not only getting the correct location and variance, but also the skewness for all elements of $V_\alpha$. Part B.4 of the Online Appendix shows that the hybrid VA also approximates the posterior of $\text{\boldmath$\theta$}$ more accurately than the two Gaussian VAs.
A key output of the study by danaher2020 are estimates of the random coefficient values $\text{\boldmath$\alpha$}$. To judge their accuracy, Figure (ref) plots the $rN=4,000$ point estimates against their exact posterior means using the three different VAs. Scatters more closely aligned along the 45 degree line indicate increased accuracy, and the hybrid VA is much more accurate than both Gaussian VAs. Figure (ref) produces the same plot for posterior of $\text{\boldmath$\alpha$}$, but for the hybrid VA when 1, 5 and 10 sweeps of the Gibbs sampler are used at step (b) of Algorithm (ref). Increasing the number of sweeps improves the estimation accuracy of $\text{\boldmath$\alpha$}$, although a reasonable degree is obtained with only five sweeps. Additional results for this example are in the Online Appendix.
We now apply our method to the large dataset using the same response but for $N=20,000$ individuals and the same $r=10$ random coefficients considered by danaher2020, for which exact posterior inference cannot be easily computed. Table (ref) reports the estimates of $\text{\boldmath$\beta$}$ and $\sigma$ for the three VAs. The results are broadly consistent, although the choice of VA affects some key parameters. For example, the coefficient of “Log Price” (a key effect in the study) is $-0.271$ with a 95% posterior interval of $(-0.361,-0.178)$ using our VA, compared to $-0.169$ with a 95% posterior interval of $(-0.333, -0.004)$ using the Gaussian factor VA. Given the retail category is off-the-peg clothing, a sizable negative coefficient is to be expected. Similarly, the coefficient of “B1 Catalog” is twice as large in our analysis at $1.335$, compared to $0.751$ using the Gaussian VA, suggesting that focal brand advertising through catalogs is more effective.
Table (ref) reports the estimate of $V_\alpha$. In danaher2020 the correlations were estimated as close to zero, whereas with the larger sample and our VA we find many to be non-zero. In particular, the correlations between the intercept and the three focal brand advertising variables are all negative ($-0.62, -0.468, -0.248$) suggesting that individuals that are heavier buyers of the focal brand (i.e. with a larger intercept random effect value) are less affected by advertising; a marketing insight not made previously. A key objective of the original study was to measure the individual-level heterogeneity in the effects of focal and cross brand advertising, and we assess this here as a function of $V_\alpha$ as follows. For individual $i$, let $\text{\boldmath$w$}_{i,t}^{B1}$ be the three focal brand advertising spend covariates, and $\text{\boldmath$w$}_{i,t}^{B23}$ be the six other brand advertising spend covariates. Further, let $V_\alpha^{B1}$ and $V_\alpha^{B23}$ be the corresponding sub-matrices of $V_\alpha$. Then for each VA we compute the distribution of the following three measures of heterogeneity:
Table (ref) reports their estimates, and total heterogeneity is similar for both dependent VAs. However, the hybrid VA estimates a substantially higher level of heterogeneity in advertising effectiveness, particularly for cross-brand advertising.
We measure calibration speed using the point predictions $\hat{y}_{i,t}=E(y_{i,t}|\text{\boldmath$\alpha$}_i,\text{\boldmath$\theta$})$, which can be calculated analytically from the tobit model. From these we can compute the root mean squared error $\mbox{RMSE}(\text{\boldmath$\alpha$},\text{\boldmath$\theta$})=(\frac{1}{NT}(\sum_{i,t}(y_{i,t}-\hat{y}_{i,t})^2))^{1/2}$ for the values of $\text{\boldmath$\alpha$},\text{\boldmath$\theta$}$ obtained during the SGA. Figure (ref) plots the RMSE against both step number and clock time for all three VAs. By this metric, convergence of the SGA is several times faster for our proposed VA than for the two Gaussian benchmarks.
The estimated computation time for the MCMC sampler is between 13.6 and 54.4 days, depending on the number of draws required, as discussed in the Online Appendix (Part B.4). In contrast, 15,000 steps of Algorithm (ref) were sufficient to calibrate our hybrid VA in 15.8 hours using 5 sweeps at step (b), although less sweeps or sub-sampling decreases computation time further as discussed below.
Sub-sampling methods that use only mini-batches of data at each step of an SGA algorithm are common in machine learning for large data sets and have two main advantages. The first is that the need to read the whole data set into computer memory is avoided, which is important in some applications. A second advantage is that a favourable trade-off between computation time per iteration and gradient variance can be made, reducing overall computation time to reach convergence.
Many models with latent variables have an augmented posterior density that factors as \[ p(\text{\boldmath$\theta$},\text{\boldmath$z$}|\text{\boldmath$y$})\propto g(\text{\boldmath$\theta$},\text{\boldmath$z$}) = \prod_{i=1}^n p(\text{\boldmath$y$}_i,\text{\boldmath$z$}_i|\text{\boldmath$\theta$}) p(\text{\boldmath$\theta$})\,, \] where $\text{\boldmath$z$}^\top=(\text{\boldmath$z$}_1^\top,\ldots,\text{\boldmath$z$}_n^\top)$ and $\text{\boldmath$y$}^\top=(\text{\boldmath$y$}_1^\top,\ldots,\text{\boldmath$y$}_n^\top)$. For these models the VA at (ref) provides a convenient framework to implement sub-sampling. Let $\widehat{G}(\text{\boldmath$\theta$},\text{\boldmath$z$},\text{\boldmath$u$})$ be an unbiased estimator of $G(\text{\boldmath$\theta$},\text{\boldmath$z$}) \equiv \nabla_\theta \log g(\text{\boldmath$\theta$},\text{\boldmath$z$})$, where $\text{\boldmath$u$}\sim f_u$ are a set of random variables determining a sub-sampling mechanism. Then the re-parameterization gradient at Theorem (ref) can be re-written as the expectation
where $f_{\varepsilon,u}(\text{\boldmath$\varepsilon$},\text{\boldmath$u$})=f_\varepsilon(\text{\boldmath$\varepsilon$})f_u(\text{\boldmath$u$})$. An estimator with one draw $\text{\boldmath$\varepsilon$}^{(s)}=(\text{\boldmath$\varepsilon$}^{0,(s)},\text{\boldmath$z$}^{(s)})\sim f_\varepsilon$ (and thus also $\text{\boldmath$\theta$}^{(s)}=h^0(\text{\boldmath$\varepsilon$}^{0,(s)},\text{\boldmath$\lambda$})$) and one draw $\text{\boldmath$u$}^{(s)}\sim f_u$, is simply
The key to successful sub-sampling within variational inference is the choice of $\widehat G$. gunawan+tk17 propose an unbiased estimator for models without latent variables based on a Taylor series expansion around $\text{\boldmath$\theta$}=\bar{\text{\boldmath$\theta$}}$ with coefficients that only need to be computed once. However, this is not applicable to the VA at (ref) because these coefficients would be functions of $\text{\boldmath$z$}$ and have to be re-computed every step of the SGA, negating the computational improvements. For random coefficient models gunawan+tk17 also suggest integrating out $\text{\boldmath$z$}$ using importance sampling, but this is also prohibitively slow when the number and complexity of the random coefficients is high. We use a more conventional sub-sampling approach that is much faster here. Let $k_i(\text{\boldmath$\theta$},\text{\boldmath$z$}_i)=\nabla_\theta \log p(\text{\boldmath$y$}_i,\text{\boldmath$z$}_i|\text{\boldmath$\theta$})$, so that $G(\text{\boldmath$\theta$},\text{\boldmath$z$})=\sum_{i=1}^n k_i(\text{\boldmath$\theta$},\text{\boldmath$z$}_i)+\nabla_\theta \log p(\text{\boldmath$\theta$})$. Consider a sub-sample $S(\text{\boldmath$u$})\subset \{1,2,\ldots,n\}$ of size $|S(\text{\boldmath$u$})|$ with $\text{\boldmath$u$}\sim f_u$ (e.g. uniform random variables that correspond to a simple random sample without replacement), then we use the unbiased estimator of $G(\text{\boldmath$\theta$},\text{\boldmath$z$})$ below: \[\widehat{G}(\text{\boldmath$\theta$},\text{\boldmath$z$},\text{\boldmath$u$})=\frac{n}{|S(\text{\boldmath$u$})|}\sum_{j\in S(\text{\boldmath$u$})} k_j(\text{\boldmath$\theta$},\text{\boldmath$z$}_j) + \nabla_\theta \log p(\text{\boldmath$\theta$})\,. \] Employing (ref) with this approximation produces two immediate computational savings in Algorithm (ref): (i) only $\text{\boldmath$z$}_i$ for $i\in S(\text{\boldmath$u$})$ have to be generated at step (b), and (ii) the evaluation of $\widehat{G}(\text{\boldmath$\theta$},\text{\boldmath$z$},\text{\boldmath$u$})$ is much faster than $G$ at step (c). The trade-off is that more steps of the SGA may be necessary as the gradient approximation at (ref) is likely to have higher variance.
To illustrate the potential of sub-sampling, we employ it to compute variational inference for the tobit model with the augmented posterior at (ref) factored across individuals $N$. We consider the large dataset in Section (ref) and sample at random without replacement 5%, 10%, 25% and 50% of the individuals at each step of the SGA. Table (ref) reports the speed (time per SGA step) and the RMSE metric averaged over the last half of a 10,000 step run of Algorithm (ref). It does so for sub-samples of different size, as well as for 1 to 5 sweeps of the Gibbs sampler at step (b) of the algorithm. By the RMSE metric, sub-samples of size 25% give comparable estimates to the full sample, even when only one sweep is used at step (b), resulting in a substantial decrease in computation time.
Variational inference has great potential for big models; especially those that employ a large number of latent variables. However, popular mean field and structured Gaussian VAs for the augmented posterior can lack accuracy, and we present a more accurate VA for this case. It combines a parametric approximation to the marginal posterior of the global parameters with the exact conditional posterior of the latent variables. This removes the approximation error due to the latent variables, also improving calibration of the VA of the global parameter posterior as in Corollary (ref). Our VA admits an efficient representation of the re-parameterization gradient in Theorem (ref) that allows the SGA algorithm to be implemented without requiring computation of $p(\text{\boldmath$z$}|\text{\boldmath$\theta$},\text{\boldmath$y$})$ or its gradient. However, our approach requires a draw from this distribution for the estimate of $\nabla_\lambda {\cal L}(\text{\boldmath$\lambda$})$ to be unbiased. For more complex latent variable models this may be difficult or slow. This is the case in our two applications, although we find that even an approximate draw using a few sweeps of simple Gibbs samplers produces much more accurate VAs than Gaussian approximations to the augmented posterior. One major advantage of our method is that it provides a flexible framework for incorporating MCMC within variational inference. For augmented posteriors that can be factorized---a common case for latent variable models---our VA is well-suited to sub-sampling variational inference. The two applications demonstrate the generality and accuracy of our method. In particular, the tobit example shows how our approach can be applied to large data sets, and how sub-sampling is easily implemented.
We finish with some comments on using our approach in practice. First, in the computation of Bayesian inference, the parameters and latent variables in a model are often treated similarly. In such cases, we recommend selecting $\text{\boldmath$z$}$ based on the viability of sampling efficiently from their conditional posterior at step (b) of Algorithm (ref). Second, while any tractable parametric density can be used for $q_\lambda^0$, either a Gaussian or Gaussian copula with a factor correlation structure is a good default choice. ong2017gaussian and Smith2020 show they are scalable and well-suited to the re-parameterization trick, improving the performance of SGA optimization. The efficient selection of the number of factors $k$ is an open problem, but these authors found that between three and ten factors worked well in a wide range of applications. In our empirical work we show that our choices for $k$ produce VAs that are fast to calibrate and accurate. Finally, we list here the circumstances where our approach has potential to provide much faster Bayesian inference than MCMC, while retaining a high degree of accuracy. They include whenever: (i) generation from $p(\text{\boldmath$\theta$}|\text{\boldmath$z$},y)$ in MCMC is either slow or difficult; (ii) generating a draw from the conditional posterior at step (b) is fast and accurate; and, (iii) the augmented posterior is factorizable so as to allow sub-sampling.
\onehalfspacing
\oldappendix {\appendixname \thesection\quad}
First, because $\mathcal L(\text{\boldmath$\lambda$}) =\mathcal L^0(\text{\boldmath$\lambda$})=E_{q_\lambda}\left[\log p(\bm{y}|\bm{\theta})+\log p(\bm{\theta})-\log q^0_\lambda\left(\bm{\theta}\right)\right]$, the re-parameterization gradient of $\mathcal L$ is the same as that of $\mathcal L^0$, so that
Here, the random vector $\text{\boldmath$\varepsilon$}^0$ has density $f_{\text{\boldmath$\varepsilon$}^0}$ that does not depend on $\text{\boldmath$\lambda$}$, and $h^0$ is the one-to-one vector-valued re-parameterization transformation from $\text{\boldmath$\varepsilon$}^0$ to the parameter vector, such that $\text{\boldmath$\theta$}=h^0(\text{\boldmath$\varepsilon$}^0,\text{\boldmath$\lambda$})$.
Next, note that Fisher's identity gives (see, for example Equation (4) of poyiadjis2011)
Substituting this expression into Equation ((ref)), and writing $E_{f_\varepsilon}\left(.\right)$ for expectation with respect to $f_{\varepsilon}(\text{\boldmath$\varepsilon$})=f_{\varepsilon^0}\left(\text{\boldmath$\varepsilon$}^0\right)p\left(\text{\boldmath$z$}|\bm{\theta},\bm{y}\right)$, and because $g(\bm{\theta},\text{\boldmath$z$})=p(\text{\boldmath$y$}|\text{\boldmath$z$},\text{\boldmath$\theta$})p(\text{\boldmath$z$}|\text{\boldmath$\theta$})p(\text{\boldmath$\theta$})$, we get
which is the required result.
For any approximating density $q(\text{\boldmath$\theta$},\text{\boldmath$z$})=q(\text{\boldmath$\theta$})q(\text{\boldmath$z$}|\text{\boldmath$\theta$})$ of $p(\text{\boldmath$\theta$},\text{\boldmath$z$}|\text{\boldmath$y$})$, the Kullback-Leibler divergence is
For the VA at (ref), the second term on the right-hand side of (ref) is always zero because $q(\text{\boldmath$z$}|\text{\boldmath$\theta$})=p(\text{\boldmath$z$}|\text{\boldmath$\theta$},\text{\boldmath$y$})$. Therefore, minimizing the divergence between $q(\text{\boldmath$\theta$},\text{\boldmath$z$})$ and $p(\text{\boldmath$\theta$},\text{\boldmath$z$}|\text{\boldmath$y$})$ for $q(\text{\boldmath$\theta$},\text{\boldmath$z$})$ of the form at (ref) is equivalent to minimizing the KL divergence between $q^0_{\lambda}(\text{\boldmath$\theta$})$ and $p(\text{\boldmath$\theta$}|\text{\boldmath$y$})$. For any other VA of the form at (ref), the approximation to the marginal posterior distribution of $\text{\boldmath$\theta$}$ cannot improve on the KL-optimal approximation within the family $q_\lambda^0(\text{\boldmath$\theta$})$ for approximation (ref). This proves the result.
Here we derive closed form expressions for two terms required to compute the re-parameterization gradient at (ref) for the factor Gaussian copula for $q_\lambda^0$ in Section (ref). The necessary results for the factor Gaussian for $q_\lambda^0$ in Section (ref) are found in ong2016variational. The third term, $\nabla_\theta \log g(\text{\boldmath$\theta$},\text{\boldmath$z$})$, is computed from the augmented likelihood, and is therefore model specific.
From (ref) and the decomposition $\Sigma_\vartheta=B_\vartheta B_\vartheta^\top+D_\vartheta^2$, we get the gradient $$\nabla_\theta\log q^0_\lambda(\bm{\theta}) = -\left[\frac{\partial \text{\boldmath$\vartheta$}}{\partial\bm{\theta}}\right]^\top (B_\vartheta B_\vartheta^\top+D_\vartheta^2)^{-1}\left(\text{\boldmath$\vartheta$}-\bm{\mu}_{\vartheta}\right)+\left[\frac{\partial }{\partial\bm{\theta}}\sum_{i=1}^m \log t_{\gamma_i}'(\theta_i)\right]^\top\,,$$ where the diagonal matrix $\frac{\partial \text{\boldmath$\vartheta$}}{\partial\bm{\theta}} = \text{diag}\left(t_{\gamma_1}'(\theta_1),\dots,t_{\gamma_m}'(\theta_m)\right)$. For the Yeo-Johnson transformation, the derivative \[ t_\gamma'(\theta)= \left\{
\right. \] so that $\frac{\partial }{\partial\bm{\theta}}\sum_{i=1}^m \log t_{\gamma_i}'(\theta_i) = \left(\frac{\gamma_1-1}{|\theta_1|+1},\dots,\frac{\gamma_m-1}{|\theta_m|+1}\right)$. For large $m$ the inverse of the $(m\times m)$ matrix $(B_\vartheta B_\vartheta^\top+D_\vartheta^2)$ (or the solutions of linear systems in this matrix) can be computed efficiently using the Woodbury formula \[ (BB^\top+D^2)^{-1}=D^{-2}-D^{-2}B(I+B^\top D^{-2}B)^{-1}B^\top D^{-2}\,, \] when $(I+B^\top D^{-2}B)$ is a $(k\times k)$ matrix with $k<<m$.
The variational parameter vector is $\text{\boldmath$\lambda$}=(\text{\boldmath$\mu$}_\vartheta^\top,\mbox{vech}(B_\vartheta)^\top,\text{\boldmath$d$}_\vartheta^\top,\text{\boldmath$\gamma$}^\top)^\top$, so that we compute the derivative of $\text{\boldmath$\theta$}$ with respect to each of these four parameter vectors. To compute these first note that from the copula model and re-parameterization trick we can write $$\bm{\theta}=t^{-1}_{\bm{\gamma}}\left(\bm{\text{\boldmath$\vartheta$}}\right)=t^{-1}_{\bm{\gamma}}\left(\text{\boldmath$\mu$}_\vartheta+B_\vartheta\text{\boldmath$\zeta$}+D_\vartheta\text{\boldmath$\epsilon$}\right)\,,$$ where we denote $t^{-1}_{\gamma}\left(\bm{\text{\boldmath$\vartheta$}}\right)\equiv (t^{-1}_{\gamma_1}(\vartheta_1),\ldots,t^{-1}_{\gamma_m}(\vartheta_m))^\top$. Then by repeated application of the chain rule
with the diagonal matrices $\frac{\partial\bm{\theta}}{\partial\text{\boldmath$\vartheta$}} = \text{diag}\left(1/t_{\gamma_1}'(\theta_1),\dots,1/t_{\gamma_m}'(\theta_m)\right)$, $\frac{\partial\text{\boldmath$\vartheta$}}{\partial\bm{\gamma}} = \text{diag}\left(\frac{\partial}{\partial \gamma_1}t_{\gamma_1}(\theta_1),\dots,\frac{\partial}{\partial \gamma_m}t_{\gamma_m}(\theta_m)\right)$ and $\text{diag}(\text{\boldmath$\epsilon$})$, where \[ \frac{\partial}{\partial \gamma}t_{\gamma}(\theta)= \left\{
\right. \] Notice that evaluation of the derivatives at (i)--(iv) above only involves sparse matrix computations, which can be employed for larger values of $m$ in practice. Last, to obtain $\frac{\partial\bm{\theta}}{\partial \mbox{\footnotesize vech}(B_\vartheta)}$ simply extract the corresponding elements from $\frac{\partial\bm{\theta}}{\partial B_\vartheta}$.
In this appendix we provide further details on the TVP-VAR-SV model in Section (ref).
Following huber2020 and references therein, the response equation of the TVP-VAR-SV model at (ref) can be rewritten as the set of $N$ regressions at (ref) as follows. First, pre-multiply (ref) by $L_t^{-1}$, so that
Notice that the parameters of the TVP-VAR are easily recovered by setting $\text{\boldmath$\beta$}_{0,t}=L_t\text{\boldmath$\gamma$}_{0,t}$ and $B_{s,t}=L_t\Gamma_{s,t}$. Because $L_t$ is a lower unitriangular matrix, then $L_t^{-1}=\{l_{i,j,t}\}$ is also a lower unitriangular matrix with $l_{i,i,t}=1$ and $l_{i,j,t}=0$ for $j>i$. Denote the non-fixed elements of the $i$th row of $L^{-1}_t$ as $\bm{l}_{1:i-1,t} = \left(l_{i,1,t},\dots,l_{i,i-1,t}\right)^\top$ for $i\geq 2$, so that the entire $i$th row of $L_t^{-1}$ is $(\bm{l}_{1:i-1,t}^\top,1,\bm{0}_{N-i}^\top)$. Then each of the $i=1,\ldots,N$ individual equations of the TVP-VAR-SV can be written as
where $\bm{y}_{1:i-1,t} = \left(y_{1,t},\dots,y_{i-1,t}\right)^\top$, $\bm{\gamma}_{i,t} = \left(\Gamma_{i,1,t},\dots,\Gamma_{i,p,t},\gamma_{i,0,t}\right)^\top$, $\Gamma_{i,s,t}$ denotes the $i$th row of $\Gamma_{s,t}$, $\gamma_{i,0,t}$ is the $i$th element in $\bm{\gamma}_{0,t}$, and $\epsilon_{i,t}\sim N(0,\exp(h_{i,t}))$. Collecting together terms in (ref) and rearranging, the $i$th equation can be expressed as the regression with time-varying parameters \[ y_{i,t}=\widetilde\text{\boldmath$x$}_{i,t}^\top\text{\boldmath$\eta$}_{i,t}+\epsilon_{i,t}\,, \] where $\widetilde\text{\boldmath$x$}_{i,t}^\top=\left(\bm{y}_{t-1}^\top,\dots,\bm{y}_{t-p}^\top,1,-\bm{y}_{1:i-1,t}^\top\right)$ and $\text{\boldmath$\eta$}_{i,t}^\top = \left(\text{\boldmath$\gamma$}_{i,t}^\top,\bm{l}_{1:i-1,t}^\top\right)$. The coefficient vector is of dimension $Np+i=J_i/2$ and follows the random walk $\text{\boldmath$\eta$}_{i,t}=\text{\boldmath$\eta$}_{i,t-1}+\operatorname{diag}{(\sqrt{\text{\boldmath$v$}_i})}\text{\boldmath$\varepsilon$}_{i,t}$, with $\text{\boldmath$\varepsilon$}_{i,t}\sim N(\bm{0},I)$ independently. (Here, we use $\mbox{diag}(\text{\boldmath$v$})$ to denote a diagonal matrix with the leading diagonal $\text{\boldmath$v$}$.)
The coefficients $\text{\boldmath$\eta$}_{i,t}$ are transformed to the “non-centered” representation $\text{\boldmath$\eta$}_{i,t}=\text{\boldmath$\eta$}_{i,0}+\mbox{diag}(\sqrt{\text{\boldmath$v$}_i})\widetilde\text{\boldmath$\eta$}_{i,t}$ as a sum of a time-invariant term $\text{\boldmath$\eta$}_{i,0}$ and scaled time-varying deviations $\widetilde\text{\boldmath$\eta$}_{i,t}$. Substituting in this parameterization gives the regression representation at (ref)
with $\text{\boldmath$x$}_{i,t}^\top =(\widetilde \text{\boldmath$x$}_{i,t}^\top, \widetilde{\text{\boldmath$x$}}_{i,t}^\top\mbox{diag}(\widetilde \text{\boldmath$\eta$}_{i,t}))$ and $\text{\boldmath$\alpha$}_i^\top=(\text{\boldmath$\eta$}_{i,0}^\top,\sqrt{\text{\boldmath$v$}_i}^\top)$. Note that $\text{\boldmath$x$}_{i,t}$ is a function of the observed time series data and the time-varying latent variables $\widetilde \text{\boldmath$\eta$}_{i,t}$ only. The latter, along with the log-volatilties, are the latent variables used in our VA in Section (ref).
Another rearrangement of (ref) produces the linear state space model
Thus, the approach of carter1994gibbs can be used to obtain draws from the conditional posterior density $p(\widetilde \text{\boldmath$\eta$}_i|\text{\boldmath$h$}_i,\text{\boldmath$\theta$}_i,\text{\boldmath$y$})$. To generate the log-volatilities from $p(\text{\boldmath$h$}_i|\widetilde \text{\boldmath$\eta$}_i,\text{\boldmath$\theta$}_i,\text{\boldmath$y$})$, we follow kim1998 and use a mixture of seven normals to approximate the distribution of $\log(\epsilon_{i,t}^2)$ and generate $\text{\boldmath$h$}_i$ joint with mixture indicators.
\baselineskip=15pt
\baselineskip=20pt \setcounter{page}{1}
\setcounter{figure}{0} \setcounter{table}{0} This Online Appendix has two parts:
{\bf {Part A: Additional Details and Results for Section 3}}\\ \ \\ This appendix is split into two sub-sections. Further details on the priors and the required gradients and derivatives for the SGA algorithm are presented in Part A.1. Additional empirical results in Part A.2. \ \\ { A.1 Priors and gradients}\\ {Priors}\\ The priors for all the parameters in the model are $$\tau_{ij}~\sim N(0,1),\ \ \ \chi_{ij}|\nu_{ij}\sim \mathcal{G}^{-1}\left(\frac{1}{2},\frac{1}{\nu_{ij}}\right),\ \ \ \xi_i|\kappa_{i}\sim \mathcal{G}^{-1}\left(\frac{1}{2},\frac{1}{\kappa_{i}}\right)$$ $$ \nu_{i1},\dots,\nu_{i2J},\kappa_{i}\sim \mathcal{G}^{-1}\left(\frac{1}{2},1\right),\ \ \ \bar{h}_i\sim N(0,10^2),\ \ \ \frac{\rho_i^h+1}{2}\sim\mathcal{B}(25,5),\ \ \ \sigma_i^2\sim\mathcal{G}\left(\frac{1}{2},\frac{1}{2}\right)$$
To conduct variational inference we transform all parameters to the real line as follows: $$\text{(i) }\chi_{ij} \ \text{is transformed to } \tilde{\chi}_{ij} = \log \chi_{ij}; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{(ii) }\xi_{i} \ \text{is transformed to } \tilde{\xi}_{i} = \log \xi_{i};$$ $$\ \text{(iii) }\nu_{ij} \ \text{is transformed to } \tilde{\nu}_{ij} = \log \nu_{ij}; \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{(iv) }\kappa_{i} \ \text{is transformed to } \tilde{\kappa}_{i} = \log \kappa_{i};$$ $$\ \text{(v) }\rho_i^h \ \text{is transformed to } \tilde{\rho}_{i} = \Phi_1^{-1}\left(\frac{\rho_i^h+1}{2}\right); \ \ \ \ \ \ \text{(vi)} \sigma_i^2 \ \text{is transformed to } \tilde{\sigma}_i = \log \sigma_i^2;$$ After the transformations we obtain the following prior density functions (using the Jacobians of the change of variables)
{Gradient}\\ The model-specific gradient vector for the parameters of the TVP-VAR-SV is:
The different terms in this gradient can be computed as
where $s_i^2 = \frac{\sigma_i^2}{1-(\rho_i^h)^2}$, $\nabla_{\tilde{\chi}_i}\log p(\tilde{\bm{\chi}}_i|\tilde{\bm{\nu}}_i) = \left(-\frac{1}{2}+\frac{1}{\nu_{i1}\chi_{i1}},\dots,-\frac{1}{2}+\frac{1}{\nu_{iJ_i}\chi_{iJ_i}}\right)^\top$,\\ $\nabla_{\tilde{\nu}_i}\log p(\tilde{\bm{\chi}}_i|\tilde{\bm{\nu}}_i) = \left(-\frac{1}{2}+\frac{1}{\nu_{i1}\chi_{i1}},\dots,-\frac{1}{2}+\frac{1}{\nu_{iJ_i}\chi_{iJ_i}}\right)^\top$, $\nabla_{\tilde{\nu}_i}\log p(\tilde{\bm{\nu}}_i) = \left(-\frac{1}{2}+\frac{1}{\nu_{i1}},\dots,-\frac{1}{2}+\frac{1}{\nu_{iJ_i}}\right)^\top$.
{\bf {A.2 Supplemental figures}}\\
{\bf {Part B: Additional Details and Results for Section 4}}\\ \ \\ This appendix is split into four sub-sections. Part B.1 provides further details on the parameterization and choice of priors. In Part B.2 all the required gradients and derivatives for the SGA algorithm are derived. Part B.3 provides details on MCMC estimation of the tobit model, how to draw $\text{\boldmath$z$}$ in step (b) of Algorithm 1 and details on how long the MCMC sampler took to run for both the small and large tobit examples. Finally, Part B.4 provides additional results for the empirical application.\\ { B.1 Parameterization and Priors}\\ The parameters of the mixed effects tobit model are $\text{\boldmath$\theta$}=(\text{\boldmath$\beta$}^\top,\mbox{vech}(L)^\top,\text{\boldmath$\omega$}^\top,\sigma^2)^\top$. All the elements in $\bm{\theta}$ need to be transformed to the real line. To do this, we introduce $c = \log(1/\sigma)$, $\xi_i = \log(\omega_i)$, $\bm\xi = \left(\xi_1,\dots,\xi_r\right)^\top$, $\kappa_j = \log(L_{j,j})$, $\bm\kappa = \left(\kappa_1,\dots,\kappa_{k_\alpha}\right)^\top$, and $\bm{l}$ which denotes the elements in $\mbox{vech}(L)$ with the diagonal elements of $L$ excepted. We then re-express this parameter vector as $\text{\boldmath$\theta$}=(\text{\boldmath$\beta$}^\top,\text{\boldmath$\xi$}^\top,c,\bm{\kappa}^\top,\bm{l}^\top)^\top$.
The prior is defined as $p\left(\bm{\theta}\right) = p\left(\bm{\beta}\right)p\left(\bm{\xi}\right)p\left(c\right)p\left(\bm{\kappa}\right)p\left(\bm{l}\right)$, with
Here, $p(c)$ was constructed by considering the prior $p(\sigma^2)\propto\frac{1}{\sigma^2}$ and deriving the corresponding prior on $c$. The prior $p\left(\bm{\xi}\right)$ was constructed by using the prior $p\left(\bm{\omega}\right) = \prod_{i=}^{r}p\left(\omega_{i}\right), \text{ with } p_\omega\left(\omega_{i}\right) = \text{Inv-Gamma}\left(\omega_{i},1,1\right)$, and then using the Jacobian of the transformation to derive the corresponding prior on $\bm{\xi}$. Finally, $p\left(\bm{\kappa}\right)$ was constructed by considering the truncated normal prior $2\phi_1\left(L_{j,j};0,\sigma_l^2\right)$ $I\left(L_{j,j}>0\right)$ on $L_{j,j}$ and deriving the corresponding prior for $\kappa_j$ via the Jacobian of the transformation.
{ B.2 Computing gradients}\\ Given the expressions for the augmented posterior and prior densities above, the function $\log g(\bm\theta,\bm{z})$ can be written as:
where $\text{\boldmath$y$}^\star$ is the $Nr$ vector of the elements $y_{i,t}^\star$. The required gradients with respect to $c$ and $\bm{\beta}$ can be computed as: $$\nabla_c \log g(\bm\theta,\bm{z}) = n-e^{2c}\sum_{i,t}\left(y_{i,t}^*-\eta_{i,t}\right)^2$$ $$\nabla_\beta \log g(\bm\theta,\bm{z})= e^{2c}\sum_{i,t}\left(y_{i,t}^*-\eta_{i,t}\right)\bm{x}_{i,t}^\top-\bm{\beta}^\top\Sigma_{\beta}^{-1}$$ To construct the gradients with respect to $\bm{l}$ and $\bm{\kappa}$, we must first compute
The elements of the gradients $\nabla_{l} f(\bm{y}^*,\bm{\alpha},\bm{\theta})$ and $\nabla_{\kappa} f(\bm{y}^*,\bm{\alpha},\bm{\theta})$ can then be constructed using
Finally, the gradient with respect to $\bm\xi$ can be constructed as:
All the expression above can be computed by noting that $\frac{\partial V_{\alpha}}{\partial L} = \left(I_{r^2}+K_{r,r}\right)\left(L\otimes I_r\right)$ and $\frac{\partial V_{\alpha}^{-1}}{\partial V_{\alpha}} = -\left(V_{\alpha}^{-1}\otimes V_{\alpha}^{-1}\right)$. We can further simplify $\frac{\partial V_{\alpha}^{-1}}{\partial V_{\alpha}}\frac{\partial V_{\alpha}}{\partial L} =-\left(I_{r^2}+K_{r,r}\right)\left(V_{\alpha}^{-1}L\otimes V_{\alpha}^{-1}\right)$. For the gradeints with respect to $\bm{\xi}$ we can use $\frac{\partial V_{\alpha}}{\partial \bm{\omega}} = I_{r^2}P$, where $P$ is the matrix of ones and zeros that extract columns 1, $r+2$, $2r+3$,...$r^2$, and $\frac{\partial \bm{\omega}}{\partial \bm{\xi}}=\Omega$. \ \\ \ \\ { B.3 Exact Bayesian inference and drawing $\text{\boldmath$z$}$}\\ For exact Bayesian inference on the augmented posterior we employ the following MCMC sampling scheme:\\ \ \\ $\underline{\text{Sampling Scheme}}$\\ \ \ Step 1: Generate from $\bm{\alpha}|\text{\boldmath$y$}_{\tiny U}^\star,\bm{\theta},\bm{y}$.\\ \ \ Step 2: Generate from $\bm{y}^\star_{{\tiny U}}|\bm{\alpha},\bm{\theta},\bm{y}$.\\ \ \ Step 3: Generate from $\bm{\theta}|\text{\boldmath$y$}_{\tiny U}^\star,\bm{\alpha},\bm{y}$.\\ \ \\ To perform Step 1, note that $p(\bm{\alpha}|\text{\boldmath$y$}_{\tiny U}^\star,\bm{\theta},\bm{y})=\prod_{i=1}^Np\left(\bm{\alpha}_i|\bm{y},\text{\boldmath$y$}_{\tiny U}^\star,\bm{\beta}\right)$, where each density in the product is an $r$-dimensional Gaussian density, so that $p\left(\bm{\alpha}_i|\bm{y},\text{\boldmath$y$}_{\tiny U}^\star,\bm{\beta}\right)=\phi_r\left(\bm{\alpha}_i;A_i^{-1}M_i^\top,A_i^{-1}\right)$ with $M_i = \frac{1}{\sigma^2}\sum_{t=1}^{T}\left(y_{i,t}^*-\bm{x}_{i,t}^\top\bm{\beta}\right)\bm{w}_{i,t}^\top$ and $A_i = V_{\alpha}^{-1}+\frac{1}{\sigma^2}\sum_{t=1}^{T}\bm{w}_{i,t}\bm{w}_{i,t}^\top$. In Step 2 we draw from $p(\text{\boldmath$y$}_{\tiny U}^\star|\bm{\alpha},\bm{\theta},\bm{y})=\prod_{i,t|y_{i,t}=0}p(y_{i,t}^\star|y_{i,t},\bm{\alpha},\bm{\theta})$, with
In Step 3, generation is conducted via random walk Metropolis-Hastings. At the beginning of each iteration, the elements of $\bm{\theta}$ are randomly assigned to groups of $10$ elements. The groups are then sampled, one group conditional on the other, with the $10$-dimensional proposal density equal to the product of $10$ independent univariate normals. The variances of the univariate normals are set adaptively to target acceptance rates between $10\%$ and $20\%$.
Steps 1 and 2 in this sampling scheme are also used to generate from $\text{\boldmath$z$}$ in step (b) of Algorithm 1 used to calibrate the Hybrid VA. The latent vector $\text{\boldmath$z$}$ is initialized at its last value in the previous step of the SGA algorithm. If $\text{\boldmath$\theta$}$ was constant, this would define a Gibb sampler for $\text{\boldmath$z$}$ over the SGA steps. However, $\text{\boldmath$\theta$}$ is not constant over steps, so our approach only provides an approximate draw from $p(\text{\boldmath$z$}|\text{\boldmath$\theta$},\text{\boldmath$y$})$. Nevertheless, once the SGA algorithm converges the value of $\text{\boldmath$\theta$}$ does not change much over SGA step, so that our approach to generate $\text{\boldmath$z$}$ at step (b) is an approximate Gibbs sampler in this sense. Our empirical results suggest that this works well.
Sampling schemes for this tobit model can be slow to mix when there are large numbers of latent variables, as in our application. Computational details for the two data sets are given below.
\ \\ {\bf {B.4 Supplemental figures and tables}}\\