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.
68,238 characters · 16 sections · 64 citation commands
Variational Bayes in State Space Models: Inferential and Predictive Accuracy
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} State space models; Variational inference; Probabilistic forecasting; Bayesian consistency; Scoring rules.
\spacingset{1.5}
A common class of models used for time series modelling and prediction is the class of {state space models (SSMs). This class includes nonlinear structures, like stochastic volatility models, regime switching models, mixture models, and models with random dynamic jumps; plus linear structures, such as linear Gaussian unobserved component models. (See { durbin2001, harvey2004,} and {giordani2011,} for extensive reviews).}
The key feature of SSMs is their dependence on hidden, {or latent,} `local'\ variables, or states, which govern the dependence of the observed{\ data, in conjunction} with a vector of unknown `global'\ parameters. {This feature leads to inferential challenges with, for example, the likelihood function { for the global parameters }being }analytically unavailable, except in {special cases. Whilst frequentist methods have certainly been adopted (see {Daniel1993, RUIZ1994, anderson1996, gallant1996, SANDMANN1998, Bates2006, AITSAHALIA2007,} and {AITSAHALIA2020}, amongst others), it is arguable that Bayesian Markov chain Monte Carlo (MCMC) methods have become the most {common tool} for analysing general SSMs, with such techniques expanded in more recent times to accommodate particle filtering, via pseudo-marginal variants such as particle MCMC (PMCMC) ({ andrieu:doucet:holenstein:2010}; {Flury2011). }See giordani2011 and fearnhead2011 for a detailed coverage of this literature, including the variety of MCMC-based algorithms adopted therein.}
Whilst (P)MCMC methods have been transformative in the SSM field, {they do suffer from certain well-known limitations. }Most notably, they require that either the {(complete) }likelihood function is\ available in closed form or that an unbiased estimator of it is available. Such methods also do not necessarily scale well to high-dimensional problems;{\ that is, to models with multiple observed and/or state processes. }If the{\ assumed data generating process (DGP)} is intractable, inference can proceed using { approximate Bayesian computation (ABC) ( dean:singh:jasra:peters:2011; CREEL2015; Martin2019), since ABC requires only simulation - not evaluation } - of the DGP. However, ABC also does not scale well to problems with a large number of parameters (see, e.g., Corollary 1 in frazier2018asymptotic for details).
Variational Bayes (VB) methods (see {blei2017variational} for a review) can be seen as a potential class of alternatives to either (P)MCMC- or ABC-based inference in SSMs. In particular, and in contrast to these methods, VB scales well to high-dimensional problems, using optimization-based techniques to effectively manage a large number of unknowns (tran2017variational; quiroz2018gaussian; koop2018variational; chan2020fast; loaiza2020fast).
{In this paper}, we make three contributions to the literature on {the application of VB to SSMs}. The first contribution is to {highlight the fundamental issue that lies at the heart of the use of VB in an SSM setting, linking this to an existing issue identified }in the literature as the `incidental parameter problem' (Neyman1948; lancaster2000incidental; westling2019beyond).{\ In brief, without due care, the application of VB to the local parameters in an SSM }leads to a lack of Bayesian consistency for the global parameters. { Moreover}{, in a class of common SSMs, we {demonstrate analytically} the impact of this inconsistency on the resulting state inference, and show that even in idealized settings} inconsistent inference for {the global parameters can }lead to{\ highly inaccurate inferences about the local parameters. }The second contribution is to review some existing variational methods, and to link their prospects for consistency to the manner in which they do, or do not, circumvent th{e incidental parameter problem. Thirdly, } we {undertake a numerical} comparison of several competing variational methods, in terms of {both }{inferential and predictive accuracy}. The key findings are that: i) correct management of the local variables leads to inferential accuracy that closely {matches that of exact (MCMC-based) Bayes; ii) }inadequate treatment of the local variables leads, in contrast, to noticeably less accurate inference; iii) predictive accuracy {shows some robustness to} inferential inaccuracy, {but only for small sample sizes. Once the size of the sample is very large, the consistency (or otherwise) of a VB }method impinges on{\ predictive accuracy, }with a clear ranking becoming evident across the methods for some {DGPs}; {with certain VB methods unable to produce similar out-of-sample accuracy results to exact Bayes {in some settings}. }
{We believe {that all three contributions serve as} novel }insights{ \ into the role of VB in SSMs, {which may} lead to best practice, if heeded. }
Throughout the remainder, we make use of the following notational conventions. Generic $p,g$ are used to denote densities, and $\pi $ is used to denote posteriors conditioned only on data, and where the conditioning is made explicit depending on the situation. For any arbitrary collection of data $(z_{1},\dots ,z_{n})$, we abbreviate this collection as $z_{1}^{n}$. For a sequence $a_{n}$, the terms $O_{p}(a_{n})$ , $o_{p}(a_{n})$ and $ \rightarrow _{p}$ have their usual meaning. Similarly, we let $ \operatornamewithlimits{plim\,}_{n}X_{n}=c$ denote $X_{n}\rightarrow _{p}c$. We let $d(\cdot ,\cdot )$ denote a metric on $\Theta \subseteq \mathbb{R} ^{d_{\theta }}$. {The proofs of all theoretical results, certain definitions, {plus additional tables and figures}, are included in the Supplementary Appendix.}
An SSM is a stochastic process consisting of the pair $\{(X_{t},Y_{t})\}$, where $\{X_{t}\}$ is a Markov chain taking values in the measurable space $( \mathcal{X},\mathcal{F}_{X},\mu )$, and $\{Y_{t}\}$ is a process taking values in a measure space $(\mathcal{Y},\mathcal{F}_{Y},\chi )$, such that, conditional on $\{X_{t}\}$, the sequence $\{Y_{t}\}$ is independent. The model is formulated through the following conditional and transition densities: for a vector of unknown random parameters $\theta $ taking values in the probability space $(\Theta ,\mathcal{F}_{\theta },P_{\theta })$, where $P_{\theta }$ admits the density function $p_{\theta }$,
where $\chi _{\theta }(\cdot ,\cdot )$ denotes the transition kernel with respect to the measure $\mu $. For simplicity, throughout the remainder we disregard {the }terms dependence on the initial measure $\nu $ and the invariant measure $\mu $, when no confusion will result. The order-one Markov assumption for $X_{t}$ is innocuous, and any finite (and known) Markov order can be accommodated via a redefinition of the state variables.
Given the independence of $Y_{t}$ conditional on $X_{t}$, and the Markovian nature of $X_{t}|X_{t-1}$, the complete data likelihood is
The (average) observed data log-likelihood is thus
and the maximum likelihood estimator (MLE) of $\theta $ is $\widehat{\theta } _{n}^{MLE}:=\operatornamewithlimits{argmax\,}_{\theta \in \Theta }\ell _{n}(\theta )$. {As is standard knowledge,} $\ell _{n}(\theta )$ {is available in closed form only for particular forms of }$g_{\theta }(y_{t}|x_{t})$ and $\chi _{\theta }(x_{t+1},x_{t});$ the canonical example being when ((ref)) and ((ref)) define a linear Gaussian state space model (LGSSM). Similarly, for $p(\theta )$ denoting the prior density, the exact (marginal) posterior for $\theta $, defined as
{{is available (e.g. via straightforward MCMC methods) only in limited cases, the LGSSM being one such case. In more complex settings and/or settings where either }$\theta $ {or }$\{(X_{t},Y_{t})\}$, {or both, are high-dimensional, accessing ((ref)) can be difficult, with standard MCMC methods leading to slow mixing, and thus potentially unreliable inferences (betancourt:2018).}}
To circumvent these issues, recent research has suggested the use of variational methods for SSMs: these methods can be used to approximate either the log-likelihood function in ((ref)) or the marginal posterior in ((ref)), depending on the mode of inference being adopted. The focus of this paper, as {already highlighted}, is on variational Bayes and, in particular, on the accuracy of such methods in SSMs. However, as part of the following section we also demonstrate the asymptotic behaviour of frequentist variational {point estimators of }$\theta $, as this result will ultimately help us interpret the behavior of the variational posterior in SSMs.
{The idea of VB is }to produce an approximation to the joint posterior $\pi (x_{1}^{n},\theta |y_{1}^{n})$ in (ref) by searching over a given family of distributions for the member that minimizes a user-chosen divergence measure between the posterior of interest and the family. This { replaces} the posterior sampling problem {with} one of optimization over the family of densities used to implement the approximation. We now review the use of variational methods in SSMs, paying particular attention to the Markovian nature of the states.
{VB} approximates the posterior $\pi (x_{1}^{n},\theta |y_{1}^{n})$ by minimizing the KL divergence between a family of densities {$\mathcal{Q}$}, with generic element $q(x_{1}^{n},\theta )$, and $\pi $:
Optimizing the KL divergence directly is not feasible {since it depends} on the unknown $\pi (x_{1}^{n},\theta |y_{1}^{n})$; the very quantity we are trying to approximate. However, minimizing the KL divergence between $q$ and $\pi $ is equivalent to maximizing the {so-called }variational evidence lower bound (ELBO):
which we can access. {Hence}, for a given class $\mathcal{Q}$, we may define the variational approximation as
The standard approach to obtaining $\widehat{q}$ is to consider a class of product distributions
with $Q$ often restricted to be mean-field, i.e., $\theta _{i}$ independent of $\theta _{j}$, $i\neq j$, and $x_{1}^{n}$ independent of $\theta $.
Regardless of the variational family adopted, $\text{KL}(q||\pi )$, and hence $\text{ELBO}(q||\pi )$, involve both $\theta $ and $x_{1}^{n}$. {The product form of }{$\mathcal{Q}$}{\ allows us to write:}
{where the last line follows from Fubini's theorem and the fact that }$ q_{x}(x_{1}^{n}|\theta )$, {by assumption, is a proper density function, for all }$\theta .$ {Further, defining}
by Jensen's inequality
Thus $\mathcal{L}_{n}(\theta )$ can be viewed as an approximation (from below) to the observed data log-likelihood. Defining
the $\text{ELBO}(q||\pi )$ can then be expressed as
This representation {decomposes }$\text{ELBO}(q||\pi )${\ into three components}, two of which only depend on the variational approximation of the global parameters $\theta $, and {a third component, }$\Upsilon _{n}(q)$ , {that} yang2020alpha refer to {as the average (with respect to }$ q_{\theta }(\theta )$){\ \textquotedblleft Jensen's gap\textquotedblright , which} encapsulates the error introduced by approximating the latent states using a given variational class. While the first and last term in the decomposition can easily be controlled by choosing an appropriate class for $ q_{\theta }(\theta )$, it is the average Jensen's gap that ultimately determines the behavior of the variational approximation.
The decomposition in (ref) has specific implications for variational inference in SSMs, which can be most readily seen by first considering the case where we only employ a variational approximation for the states, and consider point estimation of the parameters $\theta $. In this case, we can think of the variational family as $\mathcal{Q} :=\{q:q(\theta ,x_{1}^{n})=\delta _{\theta }\times q_{x}(x_{1}^{n}|\theta )\} $, where $\delta _{\theta }$ is the Dirac delta function at $\theta $, and we can then write
where we abuse notation and represent functions with arguments $\delta _{\theta }$ only by the parameter value $\theta \in \Theta $, and also make use of the short-hand notation $q_{x}$ for $q_{x}(x_{1}^{n}|\theta ).$ Define the variational point estimator as
At a minimum, we would hope that the variational estimator $\widehat{\theta } _{n}$ converges to the same point as the MLE. To deduce the behavior of $ \widehat{\theta }_{n}$, we employ the following high-level regularity conditions.
Low level regularity conditions that imply Assumption (ref) are given in douc2011consistency. Since the main thrust{\ of this paper is to deduce the accuracy of }variational methods{\ in SSMs, and not }to focus on{\ the technical details of the SSMs in particular, }we make use of high-level conditions to simplify the exposition{\ and reduce necessary technicalities that may otherwise obfuscate the main point.}
The following result shows that consistency of $\widehat{\theta }_{n}$ (for $ \theta _{0}$) is guaranteed if the variational family for the states is `good enough'.
The above result demonstrates that for the variational point estimator $ \widehat{\theta }_{n}$ to be consistent, the (infeasible) average Jensen's gap must converge to zero. Intuitively, this requires that the error introduced by approximating the states grows more slowly than the rate at which information accumulates in our observed sample, i.e., $n$. The condition $\kappa _{n}=o_{p}(1)$ is stated at the true value, $\theta _{0}$, rather than at the estimated value, as it will often be easier to deduce satisfaction of the condition, or otherwise, at convenient points in the parameter space.
As the following example illustrates, even in the simplest SSMs, the scaled (average) Jensen's gap need not vanish in the limit, and can ultimately pollute the resulting inference on $\theta _{0}$.
Lemma (ref) demonstrates that even in {this} simplest of SSMs, variational inference is {inconsistent in anything other than }the most vacuous cases. {In short}, so long as there is weak dependence in states the estimator of $\alpha $ is inconsistent; alternatively, if there is no relationship between $Y_{t}$ {and }$X_{t}$, i.e., $\alpha _{0}=0$, then the only way in general to obtain consistent inference for $\rho _{0}$ is if $ \rho _{0}=0$!
While the above results pertain to variational point estimators of $\theta _{0}$, a similar result can be stated in terms of the so-called `idealized'\ variational posterior. To state this result, we approximate the state posterior using the class of variational approximations,
where $\lambda \in \Lambda $ denotes the vector of {so-called `variational} parameters' that characterize the elements in $\mathcal{Q}$. With reference to ((ref)), making the dependence of $q_{\lambda }(x_{1}^{n})$ on the variational parameter $\lambda $ explicit leads to the criterion $ L_{n}(\theta ,\lambda )$, where $q_{x}(x_{1}^{n})$ in ((ref)) is replaced by $q_{\lambda }(x_{1}^{n} )$. Optimizing over $\lambda $ for fixed $\theta $ yields the profiled criterion,
and the `idealized'\ variational posterior for $\theta $,
We remark that, unlike with the frequentist optimization problem, the idealized VB posterior incorporates {a component of }Jensen's gap directly into the definition of that posterior. A sufficient condition for the `VB ideal' to concentrate onto $\theta _{0}$ is that $\theta _{0}$ is the {maximum } of a well-defined limit counterpart to $\widehat{L}_{n}(\theta )$. However, there is no reason to suspect this is the case a priori.
The `idealized' variational posterior $\widehat{q}(\theta |y_{1}^{n})$ is a generalized posterior, in the sense of bissiri2016general, based on the profiled criterion function $\widehat{L}_{n}(\theta )$. Given that $ \widehat{q}(\theta |y_{1}^{n})$ is constructed from a profiled criterion, the `idealized' variational posterior is then related to the {frequentist profiled variational inference }approach described in westling2019beyond. In their analysis, the authors view variational point estimators of the global parameters $\theta $ as $M$-estimators based on the profiled variational criterion function in (ref). They then explore conditions and examples under which the variational point estimator, based on maximizing $\widehat{L}_{n}(\theta )$, do, or do not, deliver consistent estimates of $\theta _{0}$.
While westling2019beyond focus on consistency of variational point estimators, we study concentration of the `idealized' posterior distribution $\widehat{q}(\theta |y_{1}^{n})$. The following result shows that, under regularity conditions similar to those maintained in westling2019beyond, the `idealized' variational posterior $\widehat{q} (\theta |y_{1}^{n})$ is Bayes consistent for some value that may or may not coincide with $\theta _{0}$.
Assumption (ref)(2.b) implies that $L_{n}(\theta ,\lambda )/n$ converges to $\mathcal{L}(\theta ,\lambda) $ (uniformly in $\theta $ and $ \lambda $); while part (2.a) is an identification condition and states that $ \mathcal{L}(\theta ,\lambda)${\ }is maximized at some $\theta _{\star }$, which may differ from $\theta _{0}$. This identification {condition }{makes clear} that if $\theta _{\star }\neq \theta _{0}$, then $\mathcal{L}[\theta _{\star },\lambda (\theta _{\star })]>\mathcal{L}[\theta _{0},\lambda (\theta _{0})]$ and the idealized posterior for $\theta $ will not concentrate onto $\theta _{0}$. This can be interpreted {explicitly }in terms of {Jensen's gap as defined in ((ref))} by recalling that under Assumption (ref), $ \ell _{n}(\theta _{0})\rightarrow _{p}H(\theta _{0})$, and by considering the limit of (the scaled) Jensen's gap evaluated at $\theta _{0}$,
for some $\delta \geq 0$. If Assumption (ref)(2.a) is satisfied at $ \theta _{\star }\neq \theta _{0}$, then $\delta >0$, and $\kappa _{n}:=\Upsilon _{n}(\theta _{0},q_{\widehat{\lambda }_{n}(\theta _{0})})/n\rightarrow _{p}C>0$.
Taken together, Lemmas (ref) and (ref) show that, regardless of whether one conducts variational frequentist or Bayesian inference in SSMs, consistent inference for $\theta _{0}$ will require that a version of Jensen's gap converges to zero. Moreover, as Example (ref) has demonstrated, this is not likely to occur even in simple SSMs. The point is further exemplified in the follow example, where we explore the Bayesian consistency of the idealized VB posterior in the same linear Gaussian SSM.
The above results suggest that VB methods can lead to inaccurate inference in the case of SSMs. In this section, we discuss further the implications of {these results for inference on the global parameters, plus their implications for predictive accuracy.}
{When conducting {VB in SSMs}, the need to approximate the posterior of $x_{1}^{n}$ introduces a discrepancy between the exact posterior, and {that which results from the VB approach}}. In this way, we can view the latent states $x_{1}^{n}$ as incidental or nuisance parameters { (see lancaster2000incidental, for a review)}, which are needed to make feasible the overall optimization problem, but which, in and of themselves, are not the object of interest. {A similar point is made by westling2019beyond in the case of independent states, and frequentist variational inference, where the authors demonstrate that inconsistency can occur, even in the case of independent observations, if delicate care is not taken with the {choice of variational class for} $x_{1}^{n}.$}
{However, the incidental parameter problem} has not stopped researchers from using VB methods to conduct inference on $\theta $ in SSMs. While the general conclusions elucidated above apply, in principle, to all such methods, we next discuss two specific categories of VB methods in greater detail, and comment on {their} ability to deliver consistent {inference} for $\theta _{0}$.
{A possible }VB approach is to first `integrate out' the latent states so that there is no need to perform joint inference on $(\theta ,x_{1}^{n})$. Such an approach can be motivated by the fact that if we take $ q_{x}(x_{1}^{n}|\theta )=\pi (x_{1}^{n}|y_{1}^{n},\theta ),$ (i.e. take the variational approximation for $x_{1}^{n}$\ to be equivalent\ to the exact posterior for $x_{1}^{n}$ conditional on $\theta $),\ then we can rewrite $\text{KL}(q||\pi )$ as
with the final line exploiting the fact that $\pi (x_{1}^{n}|y_{1}^{n},\theta )$ integrates to one for all\ $\theta .$\ Thus, if we are able to use as our variational approximation for the states the actual (conditional) posterior, {we can transform} a variational problem for $(\theta ,x_{1}^{n})$ into a variational problem for $\theta $ alone.
The above approach is adopted by loaiza2020fast, and is applicable in any case where draws from $p(x_{1}^{n}|y_{1}^{n},\theta )$ can be reliably and cheaply obtained, with the resulting draws then used to `integrate out'\ the states via the above KL divergence representation. While the approach of loaiza2020fast results in the above simplification, the real key to their approach is that it can be used to unbiasedly estimate the gradient of $\text{ELBO}[q_{\theta }||\pi (\theta |y_{1}^{n})]$\ (equivalent, in turn, to the gradient of the joint ELBO in ((ref)), by the above argument). This, in turn, allows optimization over $q_{\theta }$ to produce an approximation to the posterior $\pi (\theta |y_{1}^{n})$. Indeed, such an approach can be applied in many SSMs, such as unobserved component models like the LGSSM, in which draws from $\pi (x_{1}^{n}|y_{1}^{n},\theta )$ can be generated {exactly via, for example, forward (Kalman) filtering and backward sampling (carter1994gibbs,fruhwirth1994data); or various nonlinear models ({e.g. those featuring }stochastic volatility), in which { efficient }}Metropolis- Hastings-{{within-Gibbs algorithms are available ( kim1998stochastic,jacquier2002bayesian,primiceri2005time,huber2020inducing )}. }
In cases where we are not able to sample {readily }from $\pi (x_{1}^{n}|y_{1}^{n},\theta )$ it may still be possible to integrate out the states {using particle filtering methods.} To this end, assume that we can obtain an unbiased estimate of the observed data likelihood $p_{\theta }(y_{1}^{n})$ using a particle filter, which we denote by $\widehat{p} _{\theta }(y_{1}^{n})$. {We follow tran2017variational} and write $ \widehat{p}_{\theta }(y_{1}^{n})$ as $\widehat{p}(y_{1}^{n}|\theta ,z)$ to make the estimator's dependence on the random filtering explicit through the dependence on a random variable $z$, with $z$ subsequently defined by the condition $z=\log \widehat{p}(y_{1}^{n}|\theta ,z)-\log p_{\theta }(y_{1}^{n})$. For $g(z|\theta )$ denoting the density of $z|\theta $, tran2017variational consider VB for the augmented posterior
which, marginal of $z$, has the correct target posterior $\pi (\theta |y_{1}^{n})$ {due to the unbiasedness of the estimator }$\widehat{p} (y_{1}^{n}|\theta ,z)$. The {authors} refer to the resulting method as variational Bayes with an intractable likelihood function (VBIL). The VBIL posteriors can be obtained by considering a variational approximation to $ \pi (\theta ,z|y_{1}^{n})$ that minimizes the KL divergence between $ q(\theta ,z)=q_{\theta }(\theta )g(z|\theta )$ and $\pi (\theta ,z|y_{1}^{n}) $:
where in this case
For fixed $\theta $, $\mathbb{E}_{z}\left[ \widehat{p}(y_{1}^{n}|\theta \,z) \right] =p_{\theta }(y_{1}^{n})$, but in general $\log \widehat{p} (y_{1}^{n}|\theta ,z)$ is a biased estimator of $\log p_{\theta }(y_{1}^{n})$ , {from which it follows that} $\Upsilon _{n}[q(\theta ,z)]\geq 0.$ However, in contrast to the general approximation of the states discussed {in Section (ref)}, which intimately relies on the choice of the approximating density $q_{x}(x_{1}^{n}|\theta )$, VBIL can achieve consistent inference on $\theta _{0}$ by choosing {an appropriate} number of particles $N$ {in the production of }$\widehat{p}(y_{1}^{n}|\theta ,z)$ .
To see this, we recall that a maintained assumption in the literature on { PMCMC} methods is that, for all $n$ and $N$, the conditional mean and variance of the density $g(z|\theta )$ satisfy $\mathbb{E}[z|\theta ]=-\gamma (\theta )^{2}/2N$, and $\text{Var}\left[ z|\theta \right] =\gamma (\theta )^{2}/N$, where $\gamma (\theta )^{2}$ is bounded uniformly over $ \Theta $; see, e.g., Assumption 1 in doucet2015efficient and Assumption 1 in tran2017variational. However, in general, $N$ is assumed to be chosen so that $\mathbb{E}[z|\theta ]=-\sigma ^{2}/2$ and $ \text{Var}\left[ z|\theta \right] =\sigma ^{2}>0$, $0<\sigma <\infty $. Note that, under this choice for $N$, for any $\varepsilon >0$
assuming $q_{\theta }(\theta _{0}),\sigma ^{2}<\infty $.
From this condition, we see that the VBIL inference problem is asymptotically the same as the VB inference problem for $\theta $ alone. Consequently, existing results on the posterior concentration of VB methods for $\theta $ alone can be used to deduce posterior concentration of the VBIL posterior for $\theta $.
Yet {another approach} for dealing with variational inference {in the presence of} states is to consider a structured approximation that allows for a dynamic updating of the approximation for the posterior of the states. Such an approximation can be achieved by embedding in the class of variational densities an analytical filter, like the Kalman filter. koop2018variational propose the use of the Kalman filter within VB {(VBKF) }as a means of approximating the posterior density of the states using Kalman recursions. In particular, the authors approximate the posterior $\pi (x_{1}^{n}|y_{1}^{n},\theta )$ by approximating the relationship between $X_{t}$ and $X_{t-1}$, which may in truth be non-linear in $\theta $, by the random walk model $X_{t}=X_{t-1}+\epsilon _{t}$, with $ \epsilon _{t}\sim i.i.d.N(0,\sigma _{0}^{2})$, and then use Kalman filtering to update the states in conjunction with a linear approximation to the measurement equation. Using this formulation, the variational approximation is of the form $q(x_{1}^{n},\theta )=q_{\theta }(\theta )q_{x}(x_{1}^{n})$, where $q_{x}(x_{1}^{n})\propto \prod_{k\geq 1}\exp (-\left\{ x_{k}-\widehat{x }_{k|k}\right\} ^{2}(1-\mathcal{K}_{k})P_{k|k-1}/2)$ and where the terms $ \mathcal{K}_{k},P_{k|k-1},\widehat{x}_{k|k}$ are explicitly calculated using the Kalman recursion: $\widehat{x}_{k|k}=\widehat{x}_{k|k-1}+\mathcal{K} _{k}(y_{t}^{\star }-\widehat{x}_{k|k-1}),$ and where $\mathcal{K}_{k}$ is the Kalman gain, $P_{k|k-1}$ is the predicted variance of the state, and in the application of koop2018variational, $y_{t}^{\star }=\log (y_{t}^{2})$.
While the solution proposed by the VBKF is likely to lead to better inference on the states, especially when $x_{1}^{n}$ behaves like a random walk, ultimately we are still `conducting inference'\ on $x_{1}^{n}$, and thus we still encounter the incidental parameter problem {as a consequence}. Indeed, taking as the variational family for $x_{1}^{n}$ the Kalman filter approximation yields, at time $k\geq 1$, a conditionally normal density with mean $\hat{x}_{k|k}$ and variance $(1-\mathcal{K}_{k})P_{k|k-1}$. Hence, we have a variational density that has the same structure as in Lemma (ref), but which allows for a time varying mean and variance. Given this similarity, there is no reason to suspect that such an approach will yield inferences that are consistent. Indeed, further intuition can be obtained by noting that, in the VBKF formulation, the simplification of the state equation means that we disregard any dependence between the states and the values of $\theta $ that drive their dynamics.
The variational approach of chan2020fast can be viewed similarly: the suggested algorithm assumes and exploits a particular dynamic structure for the states that allows for analytical (posterior) updates and thus leads to computationally simple estimates for the variational densities of $ q_{x}(x_{1}^{n}).$ As with the VBKF approach, the assumed nature of the state process used by chan2020fast to estimate $q_{x}(x_{1}^{n})$ implies that, in general, it is unlikely that Bayesian consistency can be achieved. {Due to space restrictions, further discussion on the specifics of this approach are relegated to Sections (ref) and (ref) of the Supplementary Appendix}.
{VB provides, at best, an approximation to the posterior and, as a result, may well yield less accurate inferences than those produced by the exact posterior (see, e.g. koop2018variational; gunawan2020variational). However, VB} can perform admirably in predictive settings{, see, e.g., {quiroz2018gaussian\ and frazierloss, amongst others}, in the sense of replicating the out-of-sample accuracy achieved by exact predictives, when such comparators are available. }(See frazier2019approximate {for a comparable finding in the context of predictions based on }ABC.) Therefore, even though the VB posterior may not necessarily converge to the true value $\theta _{0}$, so long as the value onto which it is concentrating is not too far away from $\theta_0$, it may be that VB-based predictions perform well in practice.
Recall the conditional density of $Y_{n+1}$ given $x_{n+1}$ and $\theta $ is $g_{\theta }(Y_{n+1}|x_{n+1})$, so that the predictive pdf for $Y_{n+1}$ can be expressed as
where the last line follows from the Markovianity of the state transition equation (see equation (ref)). In many large SSMs, using MCMC methods to estimate (ref) is infeasible or {prohibitive computationally}, {due to the difficulty of sampling from }$\pi (x_{1}^{n+1},\theta |y_{1}^{n}).$ {Instead}, VB methods can produce an estimate of $p(Y_{n+1}|y_{1}^{n})$ by approximating, in various ways, the two pieces in equation (ref) underlined as (1) and (2). All such methods replace the second underlined term by some approximate posterior for $\theta $, but differ in how they access the first underlined term.
In all the cases {of which }we are aware, we can separate VB methods for prediction in SSMs into two {classes}: a class which makes explicit use of a variational approximation to the states, $\widehat{q}_{x}$ {to replace } $p(x_{1}^{n}|y_{1}^{n},\theta )$; and a class that uses an accurate simulation-based estimate of $p(x_{1}^{n}|y_{1}^{n},\theta )$. Due to space restrictions, we do not give a detailed discussion of how these VB predictives are produced, and instead refer the interested reader to Section (ref) in the Supplementary Appendix.
Any Bayesian method that replaces $\pi(\theta|y_1^n)$ in part (2) by an approximation, e.g., $\widehat{q}_\theta$ in the case of VB, will lead to some inaccuracy, however, as shown by frazier2019approximate in the case of ABC, this loss in accuracy is often minimal. Therefore, what really matters in terms of accurate prediction in SSMs using VB is the replacement of (1) in (ref). Replacing (1) with an accurate simulation-based estimate is likely to deliver more accurate estimators, at the cost of additional computation. However, it is not necessarily clear that the resulting predictions will perform much better than those approaches based on the approximation $\widehat{q}_x$. In the following section, we demonstrate that even through the inference that results from using $\widehat{q}_x$ instead of (1) can be poor, the resulting predictive performance is often quite reasonable, at least for sample sizes that are not too large.
In this section, we shed {further }light on the phenomenon of {the } predictive accuracy of VB methods, and connect the performance of these methods to the inconsistency for $\theta _{0}$ that can result as the sample size diverges. {The results suggest that, in terms of predictive accuracy, there is little difference between methods in small sample sizes or with a small number of out-of-sample observations. However, we document a clear hierarchy across methods as the sample} size becomes larger and as the out-of-sample evaluation increases.
{We now compare the inferential and predictive accuracy of the variational methods of quiroz2018gaussian and loaiza2020fast against an exact MCMC-based estimate of $\pi (\theta ,x_{1}^{n}|y_{1}^{n})$, referred to as `exact Bayes' hereafter, in a simulation exercise}. {In Section (ref) of the Supplementary Appendix, we provide complete details on the implementation of each of these methods under this particular simulation design.} However, we remark here that quiroz2018gaussian is an example of a VB method in which the states are approximated via a particular choice of variational family, whilst loaiza2020fast (as noted in Section (ref)) adopt a variational approximation for the posterior of the global parameters only, with the conditional posterior of the states accessed via simulation.
The assumed DGP is specified as an unobserved component model with stochastic volatility (UCSV):
where $(\varepsilon _{t},\eta _{t},u_{t})^{^{\prime }}\overset{i.i.d.}{\sim } N(0,I_{3})$. The unobserved component term $\mu _{t}$ is a latent variable that captures the persistence in the conditional mean of $Y_{t}$, while the stochastic volatility term $h_{t}$ captures the persistence in the conditional variance. We consider the following three set of values for the true parameters:
The specifications for DGP 1 produce a time series process that has substantial persistence in the conditional mean, and a constant variance; DGP 2 generates a process that has substantial persistence in the conditional variance, and a fixed marginal mean of zero; whilst DGP 3 corresponds to a process that exhibits persistence in both the conditional mean and variance. The true parameter vector in each case is defined as $\theta _{0}=\left( \bar{\mu}_{0},{\rho _{\mu }}_{0},{\sigma _{\mu }}_{0},\bar{h} _{0},{\rho _{h}}_{0},{\sigma _{h}}_{0},\right) ^{\prime }$.
For the predictive assessment we compare exact Bayes with the two variational methods cited above plus the method of chan2020fast. As discussed in Section (ref) of the Supplementary Appendix, the method of chan2020fast exploits a very specific structure in the construction of the variational algorithm, which in this case corresponds to DGP 2 under the parameter restrictions ${\rho _{h}}_{0}=1.0$, ${\sigma_{\mu }}_0=0.0$, and $\bar{h}_{0}=0.0$. Thus, application of this approach under any of the above true DGPs constitutes misspecified inference; hence, we do not include this technique in the inferential assessment. Due to space constraints, certain tables and figures are included in Section (ref) of the Supplementary Appendix.
We assess inferential accuracy through lens of state estimation. To this end, we generate a times series of length $T=11000${\ from each of the three true DGP specifications.} {The full sample} is used {to produce the} exact posterior as well as {the two} {approximate posteriors corresponding to the QNK and LSND methods};{\ hence, we are able to shed some light on the theoretical consistency results provided above. }We assess the inferential accuracy of each method {(exact and approximate) }by calculating the root mean squared error ({RMSE)} and mean absolute error (MAE) of {each sequence of} {marginal }posterior means, {for }$t=1,2,....,T$, for the unobserved component, $\mu _{t},$ and the stochastic {{standard deviation}, $\exp (h_{t}/2)$}, relative to the{\ marginal posterior means that }results when we condition on the true parameters, denoted respectively by $\mathbb{E}[\mu _{t}|\theta _{0},y_{1}^{T}]$ and $\mathbb{E}[\exp (h_{t}/2)|\theta _{0},y_{1}^{T}]$, $t=1,2,....,T.$ The results are presented in Table (ref).
As expected, exact Bayes produces the most accurate point estimates for the two sets of latent variables (both $\mu _t$ and $\exp (h_t/2)$), as tallies with the theoretical guarantees of this method. In terms of the VB methods, the LSND results closely match those of exact Bayes as this method does not suffer from the incidental parameter problem. In contrast, the QNK method does not deal directly with this problem and, as a consequence, exhibits - across all of the designs recorded in Table (ref) - inaccuracy that is between two and ten times greater than that of both exact Bayes and the LSND method. From the results recorded in Table (ref) in Supplementary Appendix (ref), we also note that the time taken to estimate the UCSV model, under all three DGPs, is approximately the same for exact Bayes and the LSND method, with the QNK approach taking roughly twice as long as both.
We further highlight the results in Table (ref) by plotting, in Figure (ref), the marginal posterior means for both $ \mu _{t}${\ and }$\exp (h_{t}/2)${, for each point in time across a given sample period, for all three methods; with the sequence of `true' posterior means (that condition on }$\theta _{0}$){\ included for comparison.} For the sake of brevity, we only present results for DGP 2, with the corresponding results for DGPs 1 and 3 placed in Supplementary Appendix (ref). Consistent with the summary results in Table (ref), the posterior means for exact Bayes and LSND are both very similar, for each $t$ , and visually very close to the corresponding true time posterior means, across the entire sample. In comparison, the QNK method consistently produces {point estimates of }the states that are {very different }from the { values that condition on the true parameters}, {as accords with the dependence of the method on a variational approximation for the states, and the consequent loss of Bayesian consistency for }$\theta _{0}$. We note that the additional figures in Supplementary Appendix (ref) demonstrate that, at least visually, the QNK method seems to produce more accurate estimates of $\mu _{t}$ under DGPs 1 and 3 than it does under DGP 2; however, it remains inaccurate in terms of estimating $\exp (h_{t}/2)$ under these alternative DGPs.
To assess the predictive accuracy of each method we conduct an expanding window prediction exercise using the same generated data as in the previous subsection. The exercise consists of constructing the Bayesian predictive density for $Y_{n+1}$, conditional on the sample $y_{1}^{n}$, for each of the competing approaches and for $n\in \{1000,\dots ,T-1\}$. For each method and each out-of-sample time point we evaluate eight measures of predictive accuracy: the logarithmic score, four censored scores, the continuously ranked probability score, the tail weighted continuously ranked probability score and the interval score. Details of all scoring rules, including appropriate references, are provided in Section (ref) of the Supplementary Appendix. We document results using 100, 1000 and 10000 out-of-sample evaluations respectively, remembering that the CY method is now included in the comparison, but only for the case of DGP 2. For reasons of space, we only present results for the largest number of out-of-sample evaluations (10000) in the main text, in Table (ref), while the results for the other evaluation periods are given in Section (ref) of the Supplementary Appendix, in Tables (ref) and (ref) respectively.
Focussing first on the results in Table (ref), based on the very large number of out-of-sample evaluations, we observe an interesting ranking. Across all designs, and according to all measures of accuracy, exact Bayes is the most accurate method. As accords with the inferential results discussed above, the LSND method has a predictive accuracy that often matches, or is extremely similar to, that of exact Bayes, followed, in order, by CY and QNK. A similar ranking holds for the results recorded in {Tables (ref) and (ref)} in Section (ref) of the Supplementary Appendix. However, the differences {between methods }are somewhat less stark over the smaller out-of-sample evaluation periods, which highlights the fact that it is ultimately the consistency properties of the different VB methods (in evidence for the largest evaluation period, given the large size of the expanding estimation windows) that is driving the discrepancies between the predictive accuracy of the competing methods.
Whilst a ranking is {in evidence} {in Table (ref)}, it can be argued that across certain DGP and scoring rule combinations, the predictive results across the different methods are {still quite} similar, both between the exact and (all) VB methods, and between the different VB methods. That is, for certain combinations of DGPs and scoring rules, all methods are {seen to perform well (relative to the benchmark of the true predictive)}, and the more substantial inferential discrepancies observed between certain of the methods are not reflected at the predictive level. This finding corroborates the point made earlier, and which has been supported by other findings in the literature, namely that computing a posterior {via an approximate method does not necessarily reduce predictive accuracy (relative to exact Bayes) by a }substantial{\ amount. }
{However, despite there being} certain DGP and scoring rule combinations where the methods perform similarly, this is not true across all DGPs and loss measures, {in particular for the larger out-of-sample evaluation period}. For example, {and with specific reference to Table (ref),} there is a clear trend that as model complexity increases {(i.e. moving from DGP 1 through to DGP 3)}, variational methods that work harder to correctly approximate the states have greater predictive accuracy. This finding is particularly marked for the log score and the interval score, which directly measure the dispersion of the posterior predictive. {In the case of DGP 3}, the all-purpose variational method of quiroz2018gaussian performs the worst across all the {methods} under analysis, {and most notably for the log score and the interval score}. This feature is most likely due to the fact that the posteriors associated with the method of quiroz2018gaussian have overly thin tails. Consequently, parameter uncertainty is not adequately accounted for when constructing the posterior predictive, {which} results in a predictive with thin tails, and ultimately translates into poor performance in scores that measure both location and/or dispersion.
We have systematically documented the behavior of variational methods, in terms of inference and prediction, within the class of state space models (SSMs). Sufficient conditions for (both frequentist and Bayesian) consistency of variational inference (VI) in SSMs have been presented in terms of the so-called Jensen's gap, which measures the discrepancy introduced within VI due to the approximation of the states. {Focusing on } variational {Bayes (VB) methods specifically, we show that only} methods that are capable of closing Jensen's gap yield Bayesian consistent inference for the global parameters and, in turn, deliver more accurate inferences for the states.
In the context of empirically relevant SSMs, we find numerical evidence of a clear hierarchy in terms of the accuracy of state inference across different variational methods: methods that can close Jensen's gap produce qualitatively more accurate inferences than those that do not. However, whilst this same hierarchy also holds for VB-based prediction, we find that the extent to which different variational approaches vary in terms for predictive accuracy depends on the data generating process (DGP), the loss in which the different methods are evaluated, {and - most importantly - the size of the out-of-sample evaluation period.} Indeed, we document that there are certain circumstances, i.e., {sample size, }DGP and loss combinations, where there is little to separate the various approaches. However, in large samples, methods that attain Bayesian consistent inference on the global parameters produce more accurate predictions.
To keep the length of this paper manageable, we have deliberately analysed and compared only a select few of the variational methods used to conduct inference and prediction in SSMs. Our findings, however, suggest that certain classes of approximations for the state posterior employed in the machine learning literature, e.g., classes based on normalising or autoregressive flows, may be flexible enough to deliver accurate inferences and predictions; we refer to, e.g., ryder2018black, and the references therein, for a discussion of such methods in SSMs. We leave a comparison between the approaches discussed herein and those commonly used in machine learning for future research.