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.
78,128 characters · 26 sections · 56 citation commands
Efficient variational approximations for state space models
{\bf Keywords}: State space models, Variational Bayes, Stochastic volatility, Multivariate Skellam model, Time-varying parameter vector autoregression \\ {\bf JEL Classification:} C11, C22, C32, C58
\thispagestyle{empty} \setcounter{page}{1}
Estimation of many state space models with nonlinear and/or non-Gaussian measurement equations is computationally challenging gribisch2022modeling,chan2022large,cross2021macroeconomic. The likelihood function of these models involves a high-dimensional integral with respect to the state variables which cannot be solved analytically, and hence renders maximum likelihood estimation to be infeasible. As an alternative, exact Bayesian estimation methods allow for the computation of the posterior distribution of the model parameters. These methods either use particle filtering chopin2020introduction, or sample from the augmented posterior of the model parameters and the states using analytical filtering carter1994gibbs. Both approaches can become computationally costly, especially with high-dimensional state vectors or when dependence between the states and the parameters is strong quiroz2018gaussian.
Variational Bayes (VB) methods can provide a scalable alternative to exact Bayesian methods. Instead of sampling exactly from the posterior, VB calibrates an approximation to the posterior via the minimization of a divergence function. However, off-the-shelf variational methods for state space models, such as mean-field variational approximations, are known to be poor wang2004lack. Gaussian VB methods as proposed by tan2018gaussian and quiroz2018gaussian use a variational family for the states that conditions on the model parameters and not on the data. The inaccuracy of these existing methods is due to the quality of the approximation to the conditional posterior distribution of the states frazier2022variational. More accurate VB methods are computationally infeasible for many state space models. For instance, tran2017variational exactly integrate out the states in the variational approximation using particle filtering. The method of loaiza2022fast is designed for the specific class of state space models where generation from the conditional posterior of the states is computationally feasible.
This paper proposes a novel VB method that is accurate and fast, and can be applied to a wide range of state space models for which estimation with existing methods is either inaccurate or computationally infeasible. Our method uses a variational approximation to the states that directly conditions on the observed data, and as such produces an accurate approximation to the exact posterior distribution. The approach is faster than existing VB methods for state space models due to the computationally efficient calibration steps it entails. The implementation only requires a measurement equation with a closed-form density representation, and a state transition distribution that belongs to the class of exponential distributions. This allows for a wide range of state space models, including ones with nonlinear measurement equations, certain types of nonlinear transition equations, and high-dimensional state vectors.
Our approximation to the states is the importance density proposed by richard2007efficient in the context of efficient importance sampling. Hence, we refer to our method as Efficient VB. scharth2016particle employ this importance distribution within a particle Markov chain Monte Carlo (PMCMC) sampler to reduce the variance of the estimate of the likelihood function. The use of this importance density inside PMCMC does not result in substantial computational gains, as it must be recalibrated at each iteration. Since VB poses an optimization problem, we can use stochastic gradient ascent (SGA) instead of a sampling algorithm. Our SGA algorithm requires draws from the approximation to the states, which is used to construct an estimate of the gradient of the objective function with respect to the parameters. Since the importance density is easy to generate from, and it does not have to be recalibrated at each SGA step, the optimization routine is fast and hence scalable to state space models with high-dimensional state vectors and a large number of observations.
Numerical experiments show that the proposed Efficient VB method provides accurate posterior densities, while it only takes a fraction of the computational cost of MCMC. The experiments employ the stochastic volatility model as the true data generating process. Since the exact posterior can be computed by MCMC methods, the accuracy of our method can be assessed for this model. We find that Efficient VB produces variational approximations to the states that are close to the exact posteriors, which result in accurate variational approximations to the parameters of the model. Efficient VB is faster than all benchmark methods with all sample sizes under consideration.
To illustrate the contributions of our method, we apply it in two empirical applications. The first application fits a multivariate Skellam stochastic volatility model with high-frequency tick-by-tick discrete price changes of four stocks. With the recent availability of high frequency trading data, the modelling of tick-by-tick price changes has become increasingly popular shephard2017continuous,koopman2018dynamic,catania2022dynamic. The model we are estimating is in the spirit of the univariate Skellam stochastic volatility model of koopman2017intraday but extended to the multivariate setting. This state space model with a non-linear measurement equation and a multivariate state vector cannot be accurately estimated with existing methods in a reasonable amount of time. Efficient VB produces posterior distributions close to the exact posteriors.
The second empirical application fits a state space model with Efficient VB that can also be estimated with existing, computationally costly, VB methods. This application fits a time-varying parameter vector autoregression with a stochastic volatility model to eight macroeconomic variables. This high-dimensional time series model is proposed by huber2021inducing, and related models are used by for instance clark2015macroeconomic and carriero2019large. This is a state space model with a nonlinear measurement equation and a high-dimensional state vector. In this complex model, our approach is accurate while the computation time is a fraction of the computation time of the benchmark methods.
The proposed VB method has the potential to produce fast and accurate estimation for a wide range of models. Computationally challenging state space models are currently estimated by VB methods that are limited to specific state space formulations. For instance, chan2022fast and gefang2019variational propose a VB method for a specific class of vector autoregression models. koop2018variational propose a VB method for a class of time-varying parameter models.
Existing variational inference methods for state space models that construct point estimates for model parameters, instead of posterior distributions, are computationally expensive. For instance, naesseth2018variational construct an approximation to the conditional posterior of the state using particle filtering, which is computationally costly and hence hinders the scalability to problems with high-dimensional state vectors. archer2015black use neural networks to construct an approximation to the posterior of the states. The parameters of these neural networks are calibrated jointly with the parameters of the model, which means that a high-dimensional gradient has to be computed in each iteration of the optimization algorithm.
The outline of the remainder of this paper is as follows. Section (ref) discusses specification and exact estimation of state space models, and Section (ref) develops our VB method. Section (ref) conducts numerical experiments to evaluate its accuracy and computational costs, and Section (ref) and (ref) apply our method to real data. Section (ref) concludes.
Let $\mathbf{y}=(\mathbf{y}_1^\top,\dots,\mathbf{y}_T^\top)^\top$ be an observed time series assumed to have been generated by a state space model with measurement and state densities
respectively, and where the prior density for $\mathbf{X}_1$ is $p(\mathbf{x}_1|\boldsymbol\theta)$, $\boldsymbol{\theta }\in\Theta$ is a $d$-dimensional parameter vector and $\mathbf{y}_t$ is an $N$-dimensional observation vector with $t=1,\dots,T$. The likelihood function for this model is given by
where $\mathbf{x} = (\mathbf{x}_1^\top,\dots,\mathbf{x}_T^\top)^\top$ and $p(\mathbf{y},\mathbf{x}|\boldsymbol\theta) = \prod_{t=1}^{T}p(\mathbf{y}_{t}|\mathbf{x}_{t},{\boldsymbol\theta })p(\mathbf{x}_{t}|\mathbf{x}_{t-1},{\boldsymbol\theta })$. Typically, the integral that characterises the likelihood is intractable, as it does not have an analytical solution. This is the case for state space models that assume non-linear or non-Gaussian measurement equations. These types of models are pervasive in econometrics and include, for instance, stochastic volatility models, some time-varying parameter models, and states space models for discrete data. Hence, maximum likelihood estimation is infeasible for a large class of econometric problems.
Bayesian analysis is concerned with computing the posterior density $p(\boldsymbol{\theta }|\mathbf{y})\propto p(\mathbf{y}|\boldsymbol\theta)p(\boldsymbol{\theta })$, where $p(\boldsymbol{\theta })$ is a given choice of prior density. The intractability in the likelihood function is tackled via two different avenues. First, for certain state space models it is feasible to use Markov chain Monte Carlo (MCMC) to generate from the augmented density
where analytical filtering methods are used to obtain draws from $p(\mathbf{x}|\boldsymbol\theta,\mathbf{y})$. MCMC effectively samples from $p(\boldsymbol\theta|\mathbf{y})$ which is a marginal density of $p(\boldsymbol\theta,\mathbf{x}|\mathbf{y})$. This approach is limited to certain classes of state space models, and the filtering techniques used can become computationally costly for large sample sizes or high-dimensional state vectors. The second Bayesian avenue generates samples from the posterior by replacing the likelihood function by its unbiased estimate $\widehat{p}_S(\mathbf{y}|\boldsymbol{\theta})$. This unbiased estimate, evaluated via particle methods, is then used inside a Metropolis-Hastings scheme. This approach, known as particle MCMC, trades off accuracy in estimation of $p(\boldsymbol\theta|\mathbf{y})$ by computational speed, via the choice in the number of particles $S$ andrieu2010particle,doucet2015efficient. While it can be applied to a broad class of state space models, this approach is known to be computationally costly and highly noisy for an inadequately low number of particles. This issue is exacerbated when the state vector is high-dimensional and a larger number of particles is required.
Variational Bayes may overcome the computational challenges in estimating state space models. The general idea behind VB is to approximate the exact posterior $p(\boldsymbol{\theta}|\mathbf{y})$ with an approximating density $q_{\hat{\lambda}}(\boldsymbol\theta)\in\mathcal{Q}$, where $\mathcal{Q} = \{q_{\lambda}(\boldsymbol\theta): \boldsymbol\lambda\in\Lambda\}$ is a class of tractable approximating densities indexed by the variational parameter $\boldsymbol{\lambda}\in\Lambda$. The most popular choice for $\mathcal{Q}$ is the Gaussian distribution class. The optimal variational parameter $\hat{\boldsymbol{\lambda}}$ is then calibrated by finding the element in $\mathcal{Q}$ that minimizes the Kullback-Leibler (KL) divergence - or any other divergence - to the exact posterior. Implementation of VB requires evaluation of the likelihood function $p(\mathbf{y}|\boldsymbol{\theta})$, which is infeasible for most state space models. tran2017variational circumvent this issue by replacing $p(\mathbf{y}|\boldsymbol{\theta})$ by the unbiased estimate $\widehat{p}_S(\mathbf{y}|\boldsymbol{\theta})$. While this approach is faster than PMCMC, it remains computationally costly due to its use of particle filtering.
VB can circumvent the computational challenges of exactly integrating out $\mathbf{x}$, by instead constructing an approximation to the augmented posterior in (ref). In this case, the approximating density is $q_{\hat{\lambda}}(\boldsymbol\theta,\mathbf{x})$ and $\mathcal{Q} = \{q_{\lambda}(\boldsymbol\theta,\mathbf{x}): \boldsymbol\lambda\in\Lambda\}$. Then, $\widehat{\boldsymbol{\lambda}}$ is obtained by minimising the KL divergence from $q_{\lambda}(\boldsymbol\theta,\mathbf{x})$ to $p(\boldsymbol\theta,\mathbf{x}|\mathbf{y})$, which is equivalent to maximising the evidence lower bound (ELBO) function $\mathcal{L}(\boldsymbol{\lambda}) = E_{q_\lambda}\left[\log p(\mathbf{y},\mathbf{x}|\boldsymbol\theta)p(\boldsymbol\theta)-\log q_{\lambda}(\boldsymbol\theta,\mathbf{x})\right]$:
VB methods that target the augmented posterior are much faster to implement relative to methods that approximate $p(\boldsymbol\theta|\mathbf{y})$ directly.
This paper proposes a variational approximation of the form
For the choice of $q_\lambda(\boldsymbol{\theta})$, we follow ong2018gaussian and employ a $d$-dimensional Gaussian density with mean $\boldsymbol{\mu}$ and a covariance matrix with a factor structure representation $\Omega = BB^\top+\text{diag}(\boldsymbol{d}^2)$, where $B$ is a $d\times p$ matrix and $\boldsymbol{d}$ is a $d$-dimensional vector. The variational parameter vector is $\boldsymbol{\lambda} = (\boldsymbol{\mu}^\top,\boldsymbol{d}^\top,\text{vech}(B)^\top)^\top$, where the vech denotes the half vectorization of a rectangular matrix.
loaiza2022fast show that the optimal choice of approximation for the latent states is $q(\mathbf{x}|\mathbf{y},\boldsymbol\theta) = p(\mathbf{x}|\mathbf{y},\boldsymbol{\theta})$, which guarantees exact integration of $\mathbf{x}$. However, the implementation of this approximation requires one to generate from $p(\mathbf{x}|\mathbf{y},\boldsymbol{\theta})$. This is computationally challenging or even infeasible for many state space model specifications, including nonlinear and high-dimensional state space models.
A faster approach which can be applied to a wider range of state space models, is to take $q(\mathbf{x}|\mathbf{y},\boldsymbol\theta) = q(\mathbf{x}|\boldsymbol{\theta})$. For instance, tan2018gaussian and quiroz2018gaussian take a multivariate Gaussian for $q(\mathbf{x}|\boldsymbol\theta)$ that does not condition on the data $\mathbf{y}$. frazier2022variational show that this type variational approximations to latent states may lead to inferential and predictive inaccuracies.
This paper develops an accurate and fast variational Bayes method for the state space model in (ref)--(ref), by proposing an approximation that can be expressed as $q(\mathbf{x}|\mathbf{y},\boldsymbol\theta)=q(\mathbf{x}|\mathbf{y})$. The proposed approximation is accurate due to the conditioning on $\mathbf{y}$. In addition, it is fast to implement as it does not directly condition on the parameter vector $\boldsymbol{\theta}$. The method is developed for models that have a closed-form measurement density $p(\mathbf{y}_t|\mathbf{x}_t,\boldsymbol\theta)$, and state transition density $p(\mathbf{x}_t|\mathbf{x}_{t-1},\boldsymbol\theta)$ that belongs to the exponential family of distributions. Generation from $p(\mathbf{x}|\mathbf{y},\boldsymbol{\theta})$ is not required. This makes our approach applicable to a wide range of different state space models. As will be discussed next, our approach is inspired by the literature on efficient importance sampling; hence we refer to it as Efficient VB.
We propose variational approximations to the states of the form
where $\boldsymbol{\varphi}$ is an auxiliary parameter vector, which works as a proxy for $\boldsymbol{\theta}$. The conditional densities $q(\mathbf{x}_t|\mathbf{x}_{t-1},\mathbf{y},\boldsymbol\varphi)$ are written in terms of a transition kernel $k(\mathbf{x}_t,\mathbf{x}_{t-1}|\boldsymbol{a}_t,\boldsymbol\varphi)$ and an integration constant $\chi(\mathbf{x}_{t-1}|\boldsymbol{a}_t,\boldsymbol\varphi)=\int k(\mathbf{x}_t,\mathbf{x}_{t-1}|\boldsymbol{a}_t,\boldsymbol{\varphi}) d \mathbf{x}_t$:
where $\boldsymbol{a}_t$ is a vector of parameters dependent on $\mathbf{y}$.
Denote $D[P||F]$ to be a divergence function between two distributions $P$ and $F$. The parameters $\boldsymbol{a} = (\boldsymbol{a}_1^\top,\dots,\boldsymbol{a}_T^\top)^\top\in A$ are calibrated so that $q(\mathbf{x}|\mathbf{y})$ accurately approximates $p(\mathbf{x}|\mathbf{y},\boldsymbol\varphi)$ as measured by $D$, that is
The optimization problem in (ref) is similar to the one considered in efficient importance sampling richard2007efficient,koopman2015numerically. Instead of sampling from the distribution of interest $p(\mathbf{x}|\mathbf{y},\boldsymbol\varphi)$, importance sampling replaces that distribution with an auxiliary distribution $q(\mathbf{x}|\mathbf{y})$. The parameters $\boldsymbol a$ are calibrated to minimize the variance of the ratio $p(\mathbf{x}|\mathbf{y},\boldsymbol\varphi)/q(\mathbf{x}|\mathbf{y})$. We follow this approach, and solve (ref) according to the steps proposed by richard2007efficient. Algorithm (ref) in scharth2016particle summarises this method for calibrating $\boldsymbol a$.
The calibration of $\boldsymbol a$ in (ref) is the most computationally expensive step in the VB optimization routine. However, if $q(\mathbf{x}|\mathbf{y})$ does not depend on $\boldsymbol\theta$, $\boldsymbol a$ does not have to be calibrated at each iteration of VB. Hence, we set $\boldsymbol\varphi$ to be a parameter vector close to $\boldsymbol\theta$ rather than $\boldsymbol\theta$ itself, and only update $\boldsymbol\varphi$ and $\boldsymbol a$ every 200 steps. The numerical experiments in Section (ref) and the empirical applications in Section (ref) and (ref) demonstrate that this does not hinder the accuracy of the approach.
Since the transition density belongs to the exponential family of distributions, it can be written as
where $\boldsymbol{T}(\mathbf{x}_t)$ denotes a vector of sufficient statistics. We select the transition kernel $k(\mathbf{x}_t,\mathbf{x}_{t-1}|\boldsymbol a_t,\boldsymbol\varphi)$ to be
The properties of the exponential family in combination with this choice of kernel guarantees the practical applicability of the proposed method because calibration of $q(\mathbf{x}|\mathbf{y})$ can be implemented via an algorithm that involves a fast recursive sequence of linear regressions. This algorithm is feasible because $ q(\mathbf{x}|\mathbf{y})$ is easy to generate from and the integration constant $\chi(\mathbf{x}_{t-1}|\boldsymbol{a}_t,\boldsymbol\varphi)$ can be evaluated, as shown in Appendix (ref).
The suggested approximation $q(\mathbf{x}|\mathbf{y})$ differs from existing Gaussian approximations in three ways. First, it allows for choices of exponential kernels that do not lead to a Gaussian approximation on the states. Second, note that while the approximation does not directly condition $\boldsymbol{\theta}$, it can condition on the proxy parameter $\boldsymbol{\varphi}$ in a non-linear fashion. That is, even if a Gaussian kernel was considered, the Efficient VB method has the ability to capture non-linear relationships between $\mathbf{x}$ and $\boldsymbol{\theta}$ via the proxy parameter, which the Gaussian approximation cannot. Third, the recursive nature of $q(\mathbf{x}|\mathbf{y})$ accurately captures the dependence structure of the exact conditional posterior $p(\mathbf{x}|\mathbf{y},\boldsymbol{\theta})$.
We solve the optimization problem in (ref) using SGA methods. SGA calibrates the variational parameter by iterating over
until convergence is achieved. The vector $\boldsymbol{\rho}^{[j]}$ contains the so called “learning parameters”, which we set according to the ADADELTA approach in zeiler2012adadelta. The vector $\widehat{\nabla_\lambda \mathcal{L}\left(\boldsymbol{\lambda}^{[j]}\right)}$ is an unbiased estimate of the gradient of the ELBO evaluated at $\boldsymbol{\lambda}^{[j]}$.
Any draw from $q_\lambda(\boldsymbol{\theta})$ can be expressed as $\boldsymbol{\theta} = \boldsymbol{\theta}(\boldsymbol{\varepsilon},\boldsymbol{\lambda})$, where $\boldsymbol{\varepsilon}\sim f_{\varepsilon}$ and $f_{\varepsilon}$ is a distribution that does not depend on $\boldsymbol{\lambda}$. Using the re-parametrization trick in kingma2013auto the ELBO gradient is written as
where the expectation is taken with respect to $q(\mathbf{x}|\mathbf{y})$ and $f_\varepsilon$. The gradient $\nabla_\theta\log p(\mathbf{y},\mathbf{x}|\boldsymbol\theta)p(\boldsymbol{\theta})$ is model specific. The expressions $\partial\boldsymbol{\theta} /\partial\boldsymbol{\lambda}$ and $\nabla_\theta \log q_\lambda(\boldsymbol{\theta} )$ are provided in ong2018gaussian. At each SGA iteration $[j]$, we calculate a sample estimate of (ref) based on only one draw for both $\mathbf{x}$ from $q(\mathbf{x}|\mathbf{y})$ and $\boldsymbol\varepsilon$ from $f_\varepsilon$. Note that $q(\mathbf{x}|\mathbf{y})$ only depends on the parameters $\boldsymbol\varphi$ and $\boldsymbol a$, which are updated by setting $\boldsymbol\varphi = \boldsymbol\mu^{[j]}$ and re-calibrating $\boldsymbol a$ as in (ref). The VB estimation routine is summarized in Algorithm (ref).
Since line 6 in Algorithm (ref) takes the most computation time, $\boldsymbol\varphi$ and $\boldsymbol a$ are only updated every 200 steps, as discussed in Section (ref). Appendix (ref) provides a detailed algorithm together with additional details on how line 6 is implemented.
This section presents numerical experiments to assess the accuracy and the computational costs of the proposed VB approach in a stochastic volatility model. This is a state space model that is widely used in economics, and that can also be estimated by MCMC and existing VB methods. This allows for the investigation of the properties of our method relative to benchmark methods for varying sample sizes.
The stochastic volatility model is defined as
where $\phi_1(x;\mu,s^2)$ denotes the univariate Gaussian density function with mean $\mu$ and variance $s^2$, $\boldsymbol{\theta} = (\bar{x},\rho,\sigma)'$ are the parameters of the model, and $x_t$ denotes the latent log-variances of the time series process at time $t$. We generate $T = 4000$ observations from the model in (ref) with the true parameter values set as $\bar{x}_0 = -1.3$, $\rho_0 = 0.95$, and $\sigma_0 = 0.3$.
The objective of the experiments is to assess the accuracy and computational costs of different VB methods in approximating the augmented posterior
for varying sample sizes, where $p(x_1|\boldsymbol{\theta})=\phi_1(x_1;\frac{\bar{x}}{1-\rho},\frac{\sigma^2}{1-\rho^2})$. Here, $\mathbf{x}=(x_1,\ldots,x_T)^\top$ and $p\left(\boldsymbol{\theta}\right) = p(\bar{x})p(\rho)p(\sigma)$ is the prior density for $\theta$, with $p(\bar{x})=N(0,1000)$, $p(\rho)=\text{Uniform}(0,0.995)$, and $p(\sigma)=\text{Inverse-Gamma}(1.001,1.001)$. VB is implemented by transforming all the parameters to the real line so that $\rho = 0.995/(1+\exp(-\kappa))$, and $\sigma = \exp(c/2)$.
The MCMC sampler generates from the exact posterior (ref). Appendix (ref) outlines the steps of the MCMC algorithm. Because of the low computational costs of implementing MCMC, we can estimate the model for multiple sample sizes. Specifically, we estimate the model using $T = 1,10,20,\dots,4000$ observations. These results provide insights into the accuracy and computational costs of the different VB methods in small and large samples.
The stochastic volatility model allows for the construction of the efficient variational approximation in (ref). Since the state transition density is Gaussian, set $\boldsymbol T( x)=(x,x^2)^\top$ and $\boldsymbol a_t = (b_t,c_t)^\top$ with $b_t$ and $c_t$ both scalars. Denote $\boldsymbol\varphi = (\bar{x}_\varphi,\rho_\varphi,\sigma_\varphi)$. The approximation to the states $q(\mathbf{x}|\mathbf{y}) = \prod_{t=1}^Tq(x_t|x_{t-1},\mathbf{y},\boldsymbol\varphi)$ is a product of Gaussian densities such that $q(x_t|x_{t-1},\mathbf{y},\boldsymbol\varphi) = \phi_1(x_t;\mu_t,\sigma_t^2)$ with $\sigma_{t} = (\sigma_\varphi^{-2}-2c_{t})^{-1/2}$, $\mu_{t} = \sigma_{t}^2\left[b_{t}+\frac{\bar{x}_\varphi+\rho_\varphi(x_{t-1}-\bar{x}_\varphi)}{\sigma_\varphi^2}\right]$, and normalising constant
Because generation of a random draw from $ q(\mathbf{x}|\mathbf{y}) $ is fast, so is line 10 in Algorithm 1.
In addition, we also implement the Hybrid VB approach proposed by loaiza2022fast and a Gaussian variational approximation. Hybrid VB requires draws from the conditional density $p(\mathbf{x}|\mathbf{y},\boldsymbol\theta)$, which can be generated using the filtering steps discussed in kim1998stochastic. Gaussian VB takes $q(\mathbf{x}|\boldsymbol\theta)=\phi_T(\mathbf{x},\boldsymbol\mu_x,C_xC_x^\top)$ to be a $T$-dimensional multivariate Gaussian density, where the Cholesky factor $C_x$ is a lower triangular matrix with three non-negative bands. In all three methods we use a Gaussian approximation with a factor covariance matrix for $\boldsymbol{\theta}$ and set the number of factors to one. The gradient expressions required for the implementation of all three VB methods are provided in Appendix (ref). We run the VB algorithms for a total of 10,000 iterations. MCMC is implemented using a burn-in sample of size 10,000 an inference sample of size 10,000.
First, we assess the accuracy of the posterior distribution for the states in $\mathbf{x}$ in a small sample with $T =500$ observations. Figure (ref) shows the Efficient VB, Gaussian VB, and MCMC posterior means of the states. Since Hybrid VB uses the exact conditional density of the states in its variational approximation, its posterior mean for the states is very similar to that of MCMC and is not included in the figure. Although Efficient VB uses an approximation to the conditional state density, its posterior mean is almost identical to the posterior mean of MCMC over the whole sample period. This is not the case for the posterior mean of Gaussian VB. The Gaussian VB posterior means overestimate the states compared to the posterior means from MCMC in almost each time period.
Additionally, we analyse the posterior dependence structure of the states. Panels (a) to (c) in Figure (ref) show the posterior correlations between $x_{100+i}$ and $x_{100+j}$ for $i,j=1,\dots,10$, for Gaussian VB, Efficient VB and MCMC, respectively. Gaussian VB underestimates most pairwise posterior correlations. Panel (a) in Figure (ref) shows that only the first and second order posterior correlations are nonzero, while the MCMC posterior correlations are positive up to at least the tenth order. The posterior correlations of Efficient VB and MCMC do not show any differences. We find the same patterns across different time periods and across longer time samples.
The differences in accuracy of the posterior distribution for the states can be explained by the choice of variational approximation for the states. Gaussian VB calibrates $q(\mathbf{x}|\boldsymbol\theta)$, which is far from the ideal distribution $p(\mathbf{x}|\mathbf{y},\boldsymbol\theta)$, as it does not condition on $\mathbf{y}$. On the other hand, Efficient VB uses $q(\mathbf{x}|\mathbf{y})$, which directly conditions on $\mathbf{y}$. Figures (ref) and (ref) suggest that when it comes to accurately representing the posterior distribution of the states, it is important to directly condition on the data. While estimation of the posterior of the states is not always the target of Bayesian analysis, accurate estimation of this posterior is critical to obtaining accurate variational approximations to the target posterior $p(\boldsymbol{\theta}|\mathbf{y})$. We demonstrate this in the next section.
Second, we assess the accuracy of the posterior distribution for the parameters. Figure (ref) shows the posterior parameter distributions for Efficient VB (solid yellow), Gaussian VB (dotted blue), Hybrid VB (dashed purple), and MCMC (solid black) in a small sample with 500 observations. The posterior of Efficient VB is close to the exact posterior of MCMC for all three parameters. The small differences between these posteriors can be summarized as a slight change in location for $\bar{x}$ and underestimation of the variance of the posteriors of $\rho$ and $\sigma$, which is a well-known property of VB. Hence, Hybrid VB also underestimates the posterior variance, but is slightly more accurate in the posterior locations. However, Gaussian VB only produces an accurate approximation for $\bar{x}$.
Additionally, we assess how the accuracy of the approximations changes with the sample size. The red lines in Panels (a.1) to (a.3) in Figure (ref) show the 99% posterior intervals for all the three parameters using Efficient VB. The shaded areas correspond to the MCMC posterior intervals. The x-axis indicates the sample size used for estimation. We find that the posterior intervals of Efficient VB are similar to those of MCMC, and the accuracy does not seem to be affected by the sample size. The posterior interval of Efficient VB concentrates to a location close to the true values with increasing sample size. Panels (c.1) to (c.3) in the figure show that the posterior intervals of Hybrid VB are more accurate for $\bar{x}$, but are not that different from Efficient VB for $\rho$ and $\sigma$. Panels (b.1) to (b.3) show that Gaussian VB is less accurate for all parameters and all sample sizes under consideration. Moreover, as the sample size increases its posteriors do not concentrate close to the true parameter values. This behaviour is likely to be related to the approximation errors to the posterior of the states that are exhibited by Gaussian VB.
Third, we compare the computational costs of the different estimation methods. Figure (ref) shows the estimation time with increasing sample sizes for Efficient VB, Gaussian VB, Hybrid VB, and MCMC. Efficient VB is substantially faster than the other methods: more than six times faster than MCMC, but also three times faster than Hybrid VB and more than twice as fast as Gaussian VB with a sample size of 4000 observations. Hence, Efficient VB is both more accurate and faster relative to Gaussian VB. Although the differences between the methods are smaller for smaller sample sizes, the ordering in the computational costs remains the same.
The substantial reduction in computational costs by Efficient VB can be explained by two properties of its variational approximation to the states. First, Efficient VB does not require computationally costly filtering steps to calibrate its approximation to the states. This explains the computational gains relative to Hybrid VB, which samples from the exact conditional density of the states by forward filtering and backward smoothing. Second, the number of variational parameters in Efficient VB depends only on the dimension of $\boldsymbol\theta$. On the other hand, the number of variational parameters in Gaussian VB also depends on the dimension of $\mathbf{x}$, making the required gradient and matrix operations computationally costly, especially as the sample size $T$ increases.
The computation time of Efficient VB is approximately equally distributed between three steps of Algorithm (ref): the calibration of the approximation to the states in line 6, the generation of the states in line 10 and the estimation of the gradient in line 11. The generation of the states is efficient since it does not require filtering steps, and so is evaluation of the gradient as its dimension is not affected by $T$. On the other hand, calibration of $q(\mathbf{x}|\mathbf{y})$ can be costly, since it involves a recursive sequence of linear regressions that increases in the number of observations. As a result, we only run step 6 of the algorithm every 200 iterations. The results in Figure (ref) indicate that this choice of update frequency does not hinder the accuracy of the model. This is also corroborated by the sensitivity analysis presented in Appendix (ref), which demonstrates that more frequent updates of $q(\mathbf{x}|\mathbf{y})$ do not produce higher accuracy of the variational approximation for either small or large sample sizes. We find similar results for the state space models with high-dimensional state vectors discussed in Section (ref) and (ref).
To illustrate Efficient VB with real data, we first fit a multivariate Skellam stochastic volatility model for high-frequency tick-by-tick discrete price changes. The Skellam stochastic volatility model is known to be a challenging model to estimate. The model is introduced by koopman2017intraday in a univariate setting. Since their inferential approach is not feasible for a multivariate setting, multivariate Skellam models have been proposed outside the state space framework koopman2018dynamic,catania2022dynamic. We extend the model to multiple assets by modelling each asset with independent Skellam distributions conditional on a multivariate transition equation that allows for volatility clustering and dependence between assets.
This application of Efficient VB to a model with a non-linear measurement equation and a multivariate state vector demonstrates that our approach is also fast and accurate in complex state space models. In contrast, the VB method of tran2017variational, which uses particle filtering, is computationally impractical with a multivariate state vector: a large number of particles would be required to produce a low variance estimate of the gradient required for its implementation. Hybrid VB is also infeasible in this model as it requires generation from the full-conditional $p(\mathbf{x}|\mathbf{y},\boldsymbol{\theta})$, for which no computationally efficient method is currently available.
The data is obtained from the New York Stock Exchange (NYSE) trades and quotes database, from which we consider trades at the NYSE of the four U.S. companies studied in koopman2017intraday: Walmart Stores Inc. (WMT), The Coca-Cola Company (KO), JPMorgan Chase & Co. (JPM), and Caterpillar Inc. (CAT). We consider price changes at a 15 second frequency in one of the most volatile weeks in the last decade: June 8th to June 12th in 2020. After applying the data cleaning steps P1, T1, T2 and T3 of barndorff2009realized, prices are constructed with the previous-tick method. Since trading hours run from 9:30:00am to 4:00:00pm, there are 1560 price observations at the 15 second frequency, which after differencing result in 1559 price changes. The final sample includes 5$\times$1559 = 7795 price changes for the week under consideration.
Let $y_{i,t}\in \mathbb{Z}$ denote the price change of asset $i=1,\dots,N$ at time $t$. We assume that these price changes follow a zero-inflated Skellam distribution as proposed by catania2022dynamic:
where $\delta_{0}(y_{i,t})$ is a Dirac point mass function at $0$ and $\exp\left(-\sigma^2_{i,t}\right)\mathcal{I}_{|y_{i,t}|}(\sigma^2_{i,t})$ denotes the Skellam probability mass function, with $\mathcal{I}_{|y_{i,t}|}(\cdot)$ the modified Bessel function of order $|y_{i,t}|$. To account for different amounts of zero returns, $\kappa_i$ governs the zero-inflation. The stochastic volatility is defined as $\sigma^2_{i,t} = \exp(s_{i,t}+x_{i,t})$, where $s_{i,t}$ captures the seasonal variation in intraday volatility, and $x_{i,t}$ is a state variable that captures the non-seasonal variation.
We assume that the price changes are conditionally independent and define the measurement equation
with $\mathbf{y}_t=(y_{1,t},y_{2,t},\ldots,y_{N,t})^\top$. The state vector of the model is $\mathbf{x}_t=(x_{1,t},x_{2,t},\ldots,x_{N,t})^\top$ for which we assume the transition equation
where $\bar{\boldsymbol x}=(\bar{x}_1,\dots,\bar{x}_N)^\top$, $\Omega$ is a diagonal matrix with elements $\boldsymbol\omega=(\boldsymbol\omega_{1},\dots,\boldsymbol\omega_N)^\top$, and $\Sigma =(LL^\top)^{-1}$ a covariance matrix. Persistence in volatility is captured by $\boldsymbol\omega $ and cross-sectional dependence between the volatilities of the stocks by $\Sigma$. The full parameter vector of this model is $\boldsymbol\theta=(\boldsymbol\kappa^\top,\bar{\boldsymbol x}^\top,\boldsymbol\omega^\top,\text{vech}(L)^\top)^\top$, where $\boldsymbol\kappa=(\kappa_1,\dots,\kappa_N)^\top$.
The seasonal component is constructed as a cubic spline function $s_{i,t}=\tilde{W}_t^\top \boldsymbol{\beta}_i$ with basis $\tilde{W}_t$. The basis functions are constructed in two steps. First, we use the steps in section 9.1.1 of greenberg2012introduction with four knots located at times 9:30:30am, 10:00:30am, 12:30:00pm and 4:00:00pm to construct the basis $W_t = (W_{1,t},\dots,W_{4,t})^\top$. Second, we apply the transformation in harvey1993forecasting to form the zero sum basis vector $\tilde{W}_t = (\tilde{W}_{1,t},\tilde{W}_{2,t},\tilde{W}_{3,t})^\top$, where $\tilde{W}_{j,t}={W}_{j,t}-{W}_{4,t}\frac{\bar{W}_j}{\bar{W}_4}$ and $\bar{W}_j=\sum_{t=1}^{1559}{W}_{j,t}$. This second step allows us to identify $\bar{x}_i$ by ensuring that $\sum_t \tilde{W}_t^\top \boldsymbol{\beta}_i=0$.
The augmented posterior of $\boldsymbol\theta$ and $\mathbf{x}$ conditional on $\mathbf{y}$ is
where $\mathbf{x}_0 = (\log\text{var}(y_1),\dots,\log \text{var}(y_N)))^\top$. The prior density for $\text{\boldmath$\theta$}$ is specified as $p(\text{\boldmath$\theta$})= p(L)\prod_{i=1}^Np(\kappa_i)$ $p(\bar{x}_i)p(\omega_i)p(\beta_i)$ with $\kappa_i\sim \text{Uniform}(0,1)$, $\bar{x}_i\sim N(0,100)$, $\omega_i\sim \text{Uniform}(0,1)$, $\boldsymbol{\beta}_i\sim N(\boldsymbol{0}_3,100I_3)$ and $p(L) = |LL^\top|^{-\frac{N+1}{2}}\prod_{i=1}^NL_{i,i}(N+1-i)$, which implies a Jeffrey's prior on $\Sigma$.
The augmented posterior of the model admits an approximation as proposed in (ref). Since the state transition density is a multivariate Gaussian, $\boldsymbol T( \mathbf{x}_t)=(\mathbf{x}_{t}^\top,\text{vec}(\mathbf{x}_{t}\mathbf{x}_{t}^\top))^\top$. We define the vector of kernel parameters to be $\boldsymbol a_{t} = (\mathbf{b}_{t},\text{vec}(C_{t}))^\top$ with $\boldsymbol{b}_{t}$ an $N$-dimensional vector and $C_{t} = \text{diag}(\boldsymbol{c}_{t})$ specified as a diagonal matrix for computational efficiency, where $\boldsymbol{c}_{t}$ is an $N$-dimensional vector.
The approximation to the states $q(\mathbf{x}|\mathbf{y}) = \prod_{t=1}^Tq(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{y},\boldsymbol\varphi)$ is a product of Gaussian densities such that $q(\mathbf{x}_{t}|\mathbf{x}_{t-1},\mathbf{y},\boldsymbol\varphi) = \phi_{N}(\mathbf{x}_{t};\boldsymbol\mu_{t},V_{t})$ with $\boldsymbol{\mu}_{t} = V_{t} \left(\boldsymbol{b}_{t}+\Sigma^{-1}\left(\bar{\mathbf{x}}+\Omega\mathbf{x}_{t-1}\right)\right)$ and $V_{t} = \left(\Sigma^{-1}-2C_{t}\right)^{-1}$. The integration constant of the transition kernel equals
While not made explicit in the notation, the parameters $\bar{\mathbf{x}}$, $\Omega$ and $\Sigma$ in the approximation are determined by the proxy parameter vector $\boldsymbol\varphi$. Since we assume $C_{t}$ to be diagonal, only the kernel parameters $\boldsymbol{b}_{t}$ and $\boldsymbol{c}_{t}$ have to be calibrated. Hence, $\boldsymbol{\gamma}_{t}^\top\boldsymbol{T}(\mathbf{x}_{t}^{[s]})$ boils down to $\boldsymbol{\tilde{\gamma}}_{t}^\top(\mathbf{x}_{t}^\top,(\mathbf{x}_{t}^2)^\top )^\top$ in Algorithm (ref) in Appendix (ref).
In addition to our method, we also consider the Gaussian VB approach. Gaussian VB takes $q(\mathbf{x}|\boldsymbol\theta)=\phi_T(\mathbf{x},\boldsymbol\mu_x,C_xC_x^\top)$ to be a $T$-dimensional multivariate Gaussian density, where the Cholesky factor $C_x$ is a lower triangular matrix with three non-negative bands. Both VB methods use a Gaussian approximation with a factor covariance matrix for $\boldsymbol{\theta}$ and set the number of factors to two. To judge the accuracy of the approximations we also implement PMCMC. Appendix (ref) presents further implementation details of all methods in this example. PMCMC is implemented using a burn-in sample of size 15,000, an inference sample of size 15,000, and a total of 1,000 particles. We run both VB algorithms for a total of 15,000 iterations. The ELBO figures in Appendix (ref) shows that the number of VB iterations suffices to achieve convergence.
We assess the accuracy of Efficient VB by comparing its posterior parameter distributions to the posteriors produced by PMCMC. Because the PMCMC algorithm requires the use of a particle filter, which is computationally infeasible for the multivariate Skellam SV model, this analysis is conducted in a univariate setting. First, consider the price changes for WMT. The computation time of PMCMC is substantial: 41 hours compared to 145 seconds of Gaussian VB and 87 seconds of Efficient VB.
Figure (ref) shows the posterior parameter distributions for $\kappa_1$, $\bar{x}_1$, $\omega_{1}$, and $\Sigma_{11}$, and the posterior mean for $\{s_{1t}\}_{t=1}^T$ in the Skellam stochastic volatility model. The posterior of Efficient VB is close to the location of the posterior of PMCMC for all four parameters, but underestimates the posterior variances. The posterior mean of Efficient VB is close to the posterior mean of PMCMC for the seasonality component across all intra-day periods. The location of the posterior distribution of Gaussian VB is less accurate for the four parameters and the seasonality.
We find similar results for the univariate Skellam stochastic volatility models for the other stocks. One way to summarize the improvement in accuracy of Efficient VB over Gaussian VB, is by comparing the ELBO values for the augmented posterior of both methods. The ELBO values averaged over the final 100 VB iterations divided by one thousand for the stocks WMT, KO, JPM, CAT equal -21.033, -14.852, -24.754 and -25.002 for Gaussian VB and -20.737, -14.508, -24.477 and -24.745 for Efficient VB, respectively. These numbers summarize the accuracy of the variational approximations to the augmented posterior, with larger numbers being preferred. Efficient VB produces larger ELBO’s for all stocks.
Next, we show that our method can be implemented efficiently to the multivariate Skellam stochastic volatility model with four stocks. The computation time is 18.2 minutes. This is substantially faster than MCMC with only one stock, which takes 41 hours as discussed in the previous section. Here we discuss the posterior parameter distributions and the posterior time-varying conditional volatilities of Efficient VB in the multivariate case.
First, Table (ref) shows the posterior means for the multivariate Skellam stochastic volatility model estimated by Efficient VB. The $\kappa_i$ can be interpreted as the percentage of price changes that are set to zero as a result of a temporary market freezing. CAT has a relatively large amount of zero-inflation compared to almost no zero-inflation in JPM. CAT and JPM are also the stocks with the highest volatility level according to the posterior means of $\bar{x}_i$. The persistence ($\omega_i$) and variance ($\Sigma_{ii}$) in the volatilities are similar across assets. The posterior means of the correlations show the strong dependence between the volatility of the assets, which cannot be captured with univariate models. The posterior mean of the correlations vary between 0.579 (WMT and JPM) and 0.712 (CAT and JPM).
Second, Figure (ref) shows the posterior mean of the time-varying conditional volatilities, together with the absolute value of the price changes. The conditional volatilities exhibit the correct volatility level and follow the daily seasonality in the data, with a spike in the volatility at the start of each trading day. Additionally, conditional volatilities peak at times of outliers in the observed price changes and show persistence over time. We conclude that the intradaily patterns in the magnitude of the price variations are captured in the conditional volatilities.
This second empirical application of Efficient VB fits a time-varying parameter vector autoregression with a stochastic volatility model (TVP-VAR-SV) to eight macroeconomic variables. This application shows that our approach is also fast and accurate in a model with a nonlinear measurement equation and a high-dimensional state vector. In contrast to the empirical application in Section (ref), Hybrid VB is also feasible in this model, which allows us to show that estimation by Efficient VB is faster with a negligible loss in accuracy.
The data contains 150 quarterly observations from 1980:Q3 to 2017:Q4 on eight macroeconomic variables. The FRED mnemonics for these variables are GDPC1, PCECC96, FPIx, CE16OV, CES0600000007, GDPCTPI, CES0600000008, and FEDFUNDS. We fit a TVP-VAR-SV with a lag length of 2. The data set is described in detail by huber2021inducing.
The VAR representation of the model is
where $\mathbf{y}_t=(y_{1,t},y_{2,t},\ldots,y_{N,t})^\top$ represents the $N$ macroeconomic variables at time $t$, $L_t^{-1}$ is a lower triangular matrix with unit-valued diagonal elements and lower-diagonal elements denoted as $\boldsymbol{l}_t$, $\text{\boldmath$\beta$}_{0,t}$ is the intercept vector, $B_{1,t},\ldots,B_{p,t}$ are $(N\times N)$ autoregressive coefficient matrices, and $H_t=\mbox{diag}(e^{h_{1,t}},\ldots,e^{h_{N,t}})$ is a diagonal matrix. The $K=(pN^2+N+N(N-1)/2)$ time-varying coefficients are collected in the $K$-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,\boldsymbol{l}_t^\top)$ and $V=\mbox{diag}(v_1,\ldots,v_K)$ is a diagonal matrix. The logarithms of the volatilities $h_{i,1},\ldots,h_{i,T}$ follow a stationary first order autoregression with mean $\bar{h}_i$ and autoregressive parameter $|\rho_i|<1$. We use a horseshoe prior to regularize the time-varying parameters, as proposed by huber2021inducing.
Estimation of the joint model in (ref) is difficult, and therefore it is common to transform the VAR model to $N$ unrelated regressions carriero2019large,kastner2020sparse. Moreover, horseshoe priors are known to result in posterior densities that are difficult to approximate ghosh2019model, which can be solved by adopting the re-parametrization proposed by ingraham2017variational. Appendix (ref) shows that after these two transformations, (ref) can be represented by $i=1,\dots,N$ state space models:
where $\mathbf{x}_{i,t} = (\tilde{\boldsymbol{\eta}}_{i,t}^\top,h_{i,t})^\top$ is the $(Np+i+1)-$dimensional state vector, with $\tilde{\boldsymbol{\eta}}_{i,t}^\top$ a function of the coefficient vector $\boldsymbol\beta_t$. The parameter vector for equation $i$ is defined as $\boldsymbol{\theta}_i=\left(\boldsymbol{\tau}_i^\top,\boldsymbol{\chi}_i^\top,\xi_i,\bar{h}_i,\rho_i,\sigma_i^2\right)^\top$, with $\boldsymbol{\alpha}_i = \sqrt{\xi_i}(\boldsymbol{\tau}_i\circ\sqrt{\boldsymbol{\chi}_i})$ and $J_i$-dimensional parameter vectors $\text{\boldmath$\tau$}_i = (\tau_{i,1},\dots,\tau_{i,J_i})^\top$ and $\sqrt{\boldsymbol{\chi}_i} = \left(\chi_{i,1}^{1/2},\dots,\chi_{i,J_i}^{1/2}\right)^\top$ with $J_i=2(pN+i)$, and scalar parameter $\xi_i$. The ($p+N-1$)-dimensional vector ${\boldsymbol{z}}_{i,t}=\left(\boldsymbol{y}_{t-1}^\top,\dots,\boldsymbol{y}_{t-p}^\top,1,-\boldsymbol{y}_{1:i-1,t}^\top\right)^\top$ with $\boldsymbol{y}_{1:i-1,t} = \left(y_{1,t},\dots,y_{i-1,t}\right)^\top$, represents the covariates in the measurement density. The parameters in the state density $\bar{\mathbf{x}}_i$, $A_{1,i}$ and $A_{2,i}$ are functions of $\rho_i$, $\sigma_i$ and $\bar{h}_i$. Here, $A_{2,i}^2$ denotes the operation of squaring each of the elements in $A_{2,i}.$
Since the state vector $\mathbf{x}_{i,t}$ is high-dimensional and enters the measurement equation non-linearly via $h_{i,t}$, the states cannot be analytically integrated out of the likelihood function. Hence, we consider the augmented posterior distribution. Let $\mathbf{y}_{(i)} \equiv (y_{i,1},\ldots,y_{i,T})^\top$ be the observations on the $i$th macroeconomic variable, $\mathbf{y}_{(\backslash i)}$ be the observations on the other $N-1$ macroeconomic variables, and $\mathbf{x}_{(i)}\equiv (\mathbf{x}_{i,1}^\top,\ldots,\mathbf{x}_{i,T}^\top)^\top$ the latent states in the $i$th equation. The augmented posterior is
where $\bar{\mathbf{x}}_{i,1} = (\boldsymbol{0}_{Np+i}^\top,\bar{h}_i)^\top$ and $V_{i,1} = \text{diag}((\boldsymbol{1}_{Np+i}^\top,\frac{\sigma_i^2}{1-\rho_i^2}))$. The prior density for $\boldsymbol\theta_i$ is specified as $p(\text{\boldmath$\theta$}_i)=p(\xi_{i}|\kappa_{i})p(\kappa_i)p(\bar{h}_i)p(\rho_i)p(\sigma_i^2)\prod_{j=1}^{J_i}p(\tau_{i,j})p(\chi_{i,j}|\nu_{i,j})p(\nu_{i,j})$, with $p(\xi_{i}|\kappa_{i})=$Inverse-Gamma(0.5,$\kappa_i^{-1})$, $p(\kappa_i)=$Inverse-Gamma(0.5,1), $p(\bar{h}_i)=N(0,100)$, $p((\rho_i+1)/2)=$Beta(25,5), $p(\sigma_i^2)$=Gamma(0.5,0.5), $p(\tau_{i,j})=N(0,1)$, $p(\chi_{i,j}|\nu_{i,j})=$Inverse-Gamma(0.5,$\nu_{i,j}^{-1})$, and $p(\nu_{i,j})=$Inverse-Gamma(0.5,1). The MCMC sampler from the augmented posterior in (ref) is discussed by huber2021inducing.
Each separate augmented posterior of the model admits an approximation as proposed in (ref). Since the state transition is a multivariate Gaussian, the sufficient summary vector is $\boldsymbol T( \mathbf{x}_t)=(\mathbf{x}_{i,t}^\top,\text{vec}(\mathbf{x}_{i,t}\mathbf{x}_{i,t}^\top))^\top$. We define the vector of kernel parameters to be $\boldsymbol a_{i,t} = (\mathbf{b}_{i,t},\text{vec}(C_{i,t}))^\top$ with $\boldsymbol{b}_{i,t}$ an $(Np+i+1)$-dimensional vector and $C_{i,t} = \text{diag}(\boldsymbol{c}_{i,t})$ specified as a diagonal matrix for computational efficiency, where $\boldsymbol{c}_{i,t}$ is an $(Np+i+1)$-dimensional vector.
The approximation to the states $q(\mathbf{x}_{(i)}|\mathbf{y}) = \prod_{t=1}^Tq(\mathbf{x}_{i,t}|\mathbf{x}_{i,t-1},\mathbf{y},\boldsymbol\varphi_i)$ is a product of Gaussian densities such that $q(\mathbf{x}_{i,t}|\mathbf{x}_{i,t-1},\mathbf{y},\boldsymbol\varphi_i) = \phi_{(Np+i+1)}(\mathbf{x}_{i,t};\boldsymbol\mu_{i,t},\Sigma_{i,t})$ with $\Sigma_{i,t} = \left(A_{2,i}^{-2}-2C_{i,t}\right)^{-1}$, and $\boldsymbol{\mu}_{i,t} = \Sigma_{i,t} \left(\boldsymbol{b}_{i,t}+A_{2,i}^{-2}\left(\bar{\mathbf{x}}_i+A_{i,1}\mathbf{x}_{i,t-1}\right)\right)$. The integration constant of the transition kernel equals
While not made explicit in the notation, the parameters $\bar{\mathbf{x}}_i$, $A_{1,i}$ and $A_{2,i}$ are determined by the proxy parameter vector $\boldsymbol\varphi_i$. Since we assume $C_{i,t}$ to be diagonal, only the kernel parameters $\boldsymbol{b}_{i,t}$ and $\boldsymbol{c}_{i,t}$ have to be calibrated. Hence, $\boldsymbol{\gamma}_{i,t}^\top\boldsymbol{T}(\mathbf{x}_{i,t}^{[s]})$ boils down to $\boldsymbol{\tilde{\gamma}}_{i,t}^\top(\mathbf{x}_{i,t}^\top,(\mathbf{x}_{i,t}^2)^\top )^\top$ in Algorithm (ref) in Appendix (ref).
In addition to our method, we also consider Gaussian and Hybrid variational approximations. All three VB methods use a Gaussian approximation with a factor covariance matrix for $\boldsymbol{\theta}_i$ and set the number of factors to one. The gradients for our variational approximation together with implementation details of the benchmark methods are provided in loaiza2022fast. MCMC is implemented using a burn-in sample of size 15,000 and inference sample of size 15,000. We run all three VB algorithms for a total of 10,000 iterations. Appendix (ref) shows that this number of iterations is suffices to achieve convergence.
The computation time of Efficient VB is 2.170 minutes. This is faster than Gaussian VB and Hybrid VB, which takes 2.452 and 8.791 minutes, respectively. Since MCMC takes 26.390 minutes, Efficient VB uses less than 9% of the time required for MCMC.
The TVP-VAR-SV model in (ref) contains a total of 172 states at each of the 150 time periods. To illustrate the posterior estimates for the state vectors, we consider the posterior mean of one of the time-varying VAR coefficients and of one of the time-varying volatilities.
First, Figure (ref) shows the posterior mean of one of the time-varying VAR coefficients in (ref), $B_{2,t}(1,3)$, across time $t$, for Efficient VB, Gaussian VB, and MCMC. Remember that Hybrid VB uses the exact conditional density of the states in its variational approximation, and hence its posterior mean for the states is very similar to MCMC and not included in the figure. Many of the time-varying VAR coefficients are regularized to zero with all methods. Therefore we illustrate the posterior state distributions by the posterior mean of $B_{2,t}(1,3)$, a coefficient that actually has time dynamics and differences in these dynamics across the different methods.
The posterior mean for Efficient VB is similar to the posterior mean for MCMC over the whole sample period. Both show substantial variation over time, with a posterior mean close to zero in 1980, increasingly positive between 1980 and 1995, decreasing between 1995 and 2010, and close to zero again between 2010 and 2017. The posterior mean for Gaussian VB follows a different time path with a small amount of variation and close to zero across the whole sample.
Second, Figure (ref) shows the posterior mean of one of the time-varying volatilities in (ref), $\exp(h_{4,t}/2)$, across time $t$. The results are similar to the posterior means for the time-varying coefficients. Efficient VB more accurately approximates the posterior mean for MCMC than Gaussian VB. Similar to the findings for Figure (ref) in the numerical experiment, Gaussian VB seems to overestimate states compared to the posterior means for MCMC. Figure (ref) in Appendix (ref) shows that a similar conclusion can be drawn when looking at the posterior mean for the time-varying volatilities of the remaining equations.
To assess the accuracy of the posterior distribution for the parameters, we focus on the parameters in the stochastic volatility model for Real Gross Domestic Product (GDPC1), which is the first variable in (ref). Figure (ref) shows the posterior parameter distributions for $\bar{h}_1$, $\rho_1$, and $\sigma_1^2$. The posterior of Efficient VB is close to the location of the posterior of MCMC for all three parameters, but underestimates the posterior variances. Although slightly more accurate, we find the same for Hybrid VB. The location of the posterior distribution of Gaussian VB is less accurate for $\rho_1$ and $\sigma_1^2$.
Table (ref) shows the ELBO values for the augmented posterior for both Efficient VB and Gaussian VB, for each equation in (ref) averaged over the final 100 VB iterations divided by one thousand. These numbers summarize the accuracy of the variational approximations to the augmented posterior, with larger numbers being preferred. We find that Efficient VB produces larger ELBO's for all equations.
The results are similar to the findings for Figure (ref) in the numerical experiment, with a smaller sample and a more complex model. If the empirical results are consistent with the numerical analyses in Figure (ref), the bias in Gaussian VB is expected to increase as the sample grows, while Efficient VB is not expected to lose accuracy.
This paper proposes a variational Bayes method for state space models, that uses a new variational approximation to the states. This approximation conditions on the observed data, which results in an accurate approximation to the posterior distribution of both the states and the model parameters. Since the approximation is calibrated in a computationally efficient way, the method is fast and scalable to a large number of states and a large number of observations. The combination of accuracy and speed of the variational approximation is illustrated in numerical experiments with a simple stochastic volatility model and empirical applications to a novel multivariate Skellam stochastic volatility model for high-frequency tick-by-tick discrete price changes and a modern macroeconomic time-varying parameter vector autoregression with stochastic volatility.
The proposed efficient variational Bayes method is applicable to a wide range of state space models, including models for which accurate estimation is computationally infeasible using existing methods. First, our method can be applied to many models with nonlinear or non-Gaussian measurement equations and/or high-dimensional state vectors. For instance, two potential additional applications of the approach are dynamic stochastic copula models hafner2012dynamic, and multivariate stochastic volatility models with realised volatility yamauchi2020multivariate.
Second, the method can be applied to models with any transition density that is a member of the exponential family of distributions. This opens up the possibility of applying state space modelling with a non-Gaussian transition density to, for instance, realised covariance matrices of asset returns. The uptake of state space models in this literature has been limited, since high-dimensional state vectors with nonlinear restrictions are generally required gribisch2022modeling. One future extension of our method is the accommodation of transition equations that are not a member of the exponential family, as is the case with the Heston model eraker2004stock.