Recursive variational Gaussian approximation with the Whittle likelihood for linear non-Gaussian state space models

An SSM consists of two sub-models: a data model that describes the relationship between the observations and the latent state variables, and a state evolution model that describes the (temporal) evolution of the latent states.

In SSMs, the $d$-dimensional latent state variables $\{\mathbf{x}_{t}\}_{t=1}^{T}$ are characterised by an initial density ${\mathbf{X}_{1}\sim p(\mathbf{x}_{1}\mid{\bm{\theta}})}$ and state transition densities

where $\bm{\theta}$ is a vector of model parameters. The observations $\{\mathbf{y}_{t}\}_{t=1}^{T}$ are related to the latent states via

where each $\mathbf{Y}_{t}$ is assumed to be conditionally independent of the other observations, when conditioned on the latent state $\mathbf{x}_{t}$. A linear SSM is one where $\mathbf{X}_{t}$ and $\mathbf{Y}_{t}$ in (1) and (2) can each be expressed as the sum of two components: a linear function of their conditioning variables and an exogenous random component.

Let $\mathbf{y}_{1:T}\equiv(\mathbf{y}_{1}^{\top},\dots,\mathbf{y}_{T}^{\top})^{\top}$ and $\mathbf{x}_{1:T}\equiv(\mathbf{x}_{1}^{\top},\dots,\mathbf{x}_{T}^{\top})^{\top}$ be vectors of observations and latent states, respectively. The posterior density of the model parameters is $p(\bm{\theta}\mid\mathbf{y}_{1:T})\propto p(\mathbf{y}_{1:T}\mid\bm{\theta})p(\bm{\theta})$, where the likelihood function is

This likelihood function is intractable, except for special cases such as the linear Gaussian SSM (e.g., Gibson and Ninness,, 2005). In non-Gaussian SSMs, the likelihood is typically estimated using a particle filter (Gordon et al.,, 1993), such as in the pseudo-marginal Metropolis-Hastings (PMMH) method of Andrieu and Roberts, (2009), which can be computationally expensive for large $T$.

The Whittle likelihood (Whittle,, 1953) is a frequency-domain approximation of the exact likelihood, and offers a computationally efficient alternative for parameter inference in SSMs. Transforming the data to the frequency domain “bypasses” the need for integrating out the states when computing the likelihood, as the Whittle likelihood only depends on the model parameters and not the latent variables. In the following sections, we review the Whittle likelihood and its properties, and then propose a sequential VB approach that uses the Whittle likelihood for parameter inference.

The Whittle likelihood has been widely used to estimate the parameters of complex time series or spatial models. It is more computationally efficient to evaluate than the exact likelihood, although it is approximate (e.g., Sykulski et al.,, 2019). Some applications of the Whittle likelihood include parameter estimation with long memory processes (e.g., Robinson,, 1995; Hurvich and Chen,, 2000; Giraitis and Robinson,, 2001; Abadir et al.,, 2007; Shao and Wu, 2007b, ; Giraitis et al.,, 2012), and with spatial models with irregularly spaced data (e.g., Fuentes,, 2007; Matsuda and Yajima,, 2009). The Whittle likelihood is also used for estimating parameters of vector autoregressive moving average (VARMA) models with heavy-tailed errors (She et al.,, 2022), and of vector autoregressive tempered fractionally integrated moving average (VARTFIMA) models (Villani et al.,, 2022).

The original R-VGA algorithm of Lambert et al., (2022) is a sequential VB algorithm for parameter inference in statistical models where the observations are independent, such as linear or logistic regression models. Given independent and exchangeable observations $\{\mathbf{y}_{i}:i=1,\dots,N\}$, this algorithm sequentially approximates the “pseudo-posterior” distribution of the parameters given data up to the $i$th observation, $p(\bm{\theta}\mid\mathbf{y}_{1:i})$, with a Gaussian distribution $q_{i}(\bm{\theta})=\mathrm{Gau}(\bm{\mu}_{i},\mathbf{\Sigma}_{i})$. The updates for $\bm{\mu}_{i}$ and $\mathbf{\Sigma}_{i}$ are available in closed form, as seen in Algorithm 1. We describe this algorithm in more detail in Section S1 of the online supplement.

At each iteration ${i=1,\dots,N}$, the R-VGA algorithm requires the gradient and Hessian of the log-likelihood contribution from each observation, $\log p(\mathbf{y}_{i}\mid\bm{\theta})$. In a setting with independent observations, these likelihood contributions are easy to obtain; however, in SSMs where the latent states follow a Markov process and are hence temporally dependent, the log-likelihood function is intractable, as shown in Section 2.1. The gradient and Hessian of $\log p(\mathbf{y}_{t}\mid\bm{\theta})$, ${t=1,\dots,T}$, are hence also intractable. While these quantities can be estimated via particle filtering (e.g., Doucet and Shephard,, 2012), one would need to run a particle filter from the first time period to the current time period for every step in the sequential algorithm. This quickly becomes infeasible for long time series.

We now extend the R-VGA algorithm of Lambert et al., (2022) to state space models by replacing the time domain likelihood with a frequency domain alternative: the Whittle likelihood. From (6) and (10) we see that the log of the Whittle likelihood is a sum over individual log-likelihood components, a consequence of the asymptotic independence of the periodogram frequency components (Shao and Wu, 2007a, ). We replace the quantities ${\nabla_{\bm{\theta}}\log p(\mathbf{y}_{i}\mid\bm{\theta})}$ and ${\nabla_{\bm{\theta}}^{2}\log p(\mathbf{y}_{i}\mid\bm{\theta})}$ in Algorithm 1 with the gradient and Hessian of the Whittle log-likelihood evaluated at each frequency. In the univariate case, the individual log-likelihood contribution at frequency $\omega_{k}$ is

while in the multivariate case,

for $\omega_{k}=\frac{2\pi k}{T}$. The log-likelihood contributions $l_{W_{k}}(\bm{\theta})$ are available for models where the power spectral density is available analytically, such as models where the latent states follow a VARMA (She et al.,, 2022) or a VARTFIMA (Villani et al.,, 2022) process. The gradient ${\nabla_{\bm{\theta}}l_{W_{k}}(\bm{\theta})}$ and Hessian ${\nabla_{\bm{\theta}}^{2}l_{W_{k}}(\bm{\theta})}$ are therefore also available in closed form for these models. At each VB update, ${\nabla_{\bm{\theta}}l_{W_{k}}(\bm{\theta})}$ and ${\nabla_{\bm{\theta}}^{2}l_{W_{k}}(\bm{\theta})}$ are used to update the variational parameters for $k=1,\dots,\left\lfloor\frac{T-1}{2}\right\rfloor$. Their expectations with respect to $q_{k-1}(\bm{\theta})$ are approximated using Monte Carlo samples ${\bm{\theta}^{(s)}\sim q_{k-1}(\bm{\theta})}$, ${s=1,\dots,S}$:

for $k=1,\dots,\left\lfloor\frac{T-1}{2}\right\rfloor$. We set $S$ following an experiment on the sensitivity of the algorithm to the choice of $S$; see Section S4 of the online supplement for more details.

We also implement the damping approach of Vu et al., (2024) to mitigate instability that may arise when traversing the first few frequencies. Damping is done by splitting an update into $D$ steps, where $D$ is user-specific. In each step, the gradient and the Hessian of the log of the Whittle likelihood are multiplied by a “step size” $a=\frac{1}{D}$, and the variational parameters are updated $D$ times during one VB iteration. This has the effect of splitting a log-likelihood contribution into $D$ “parts”, and using them to update variational parameters one part at a time. In practice, we set $D$ through experimentation, and choose it so that initial algorithmic instability is reduced with only a small increase in computational cost. We summarise R-VGA-Whittle in Algorithm 2.

The spectral density of a process captures the distribution of signal power (variance) across frequency (Percival and Walden,, 1993). The periodogram is an estimate of the spectral density. In Figure 1(a), we show the periodogram for data generated from the linear Gaussian SSM that we consider in Section 3.1. The variance (power) of the signal gradually decreases as the frequency increases, which means that lower frequencies contribute more to the total power of the process. A similar pattern is seen in Figure 1(b), which shows the periodogram for data generated from a univariate stochastic volatility model that we consider in Section 3.2. In this case there is an even larger concentration of power at the lower frequencies.

Frequencies that carry less power have less effect on R-VGA-Whittle updates. To speed up the algorithm and reduce the number of gradient/Hessian calculations, we therefore process angular frequency contributions with low power, namely those beyond the 3dB cutoff point, together in “blocks”. We define the 3dB cutoff point as the angular frequency at which the power drops to half the maximum power in the (smoothed) periodogram, which we obtain using Welch’s method (Welch,, 1967). Section S3 studies several cutoff points and empirically justifies our choice of a 3dB cutoff.

Assume the 3dB cutoff point occurs at $\omega_{\tilde{n}}$; we split the remaining $\lfloor\frac{T-1}{2}\rfloor-\tilde{n}$ angular frequencies into blocks of approximately equal length $B$. Suppose that there are $n_{b}$ blocks, with $\Omega_{b},b=1,\dots,n_{b}$, denoting the set of angular frequencies in the $b$th block. The variational mean and precision matrix are then updated based on each block’s Whittle log-likelihood contribution. For $b=1,\dots,n_{b}$, this is given by

that is, by the sum of the Whittle log-likelihood contributions of the frequencies in the block.

For a multivariate time series, we compute the periodogram for each individual time series along with its corresponding 3dB cutoff point, and then use the highest cutoff point in R-VGA-Whittle. An example of bivariate time series periodograms along with their 3dB cutoffs points is shown in Figure 2.

We summarise R-VGA-Whittle with block updates in Algorithm 3.

We now discuss the use of HMC with the Whittle likelihood for estimating parameters in SSMs. Suppose we want to sample from the posterior distribution $p(\bm{\theta}\mid\mathbf{y}_{1:T})$. We introduce an auxiliary momentum variable $\mathbf{r}$ of the same length as $\bm{\theta}$ and assume $\mathbf{r}$ comes from a density $\mathrm{Gau}(\mathbf{0},\mathbf{M})$, where $\mathbf{M}$ is a mass matrix that is often a multiple of the identity matrix. The joint density of $\bm{\theta}$ and $\mathbf{r}$ is $p(\bm{\theta},\mathbf{r}\mid\mathbf{y})\propto\exp(-H(\bm{\theta},\mathbf{r}))$, where

is called the Hamiltonian, and $\mathcal{L}(\bm{\theta})=\log p(\bm{\theta}\mid\mathbf{y}_{1:T})$. HMC proceeds by moving $\bm{\theta}$ and $\mathbf{r}$ according to the Hamiltonian equations

to obtain new proposed values $\bm{\theta}^{*}$ and $\mathbf{r}^{*}$. By Bayes’ theorem, $p(\bm{\theta}\mid\mathbf{y}_{1:T})\propto p(\mathbf{y}_{1:T}\mid\bm{\theta})p(\bm{\theta})$, so the gradient $\nabla_{\bm{\theta}}\mathcal{L}(\bm{\theta})$ can be written as

The gradient $\nabla_{\bm{\theta}}\log p(\mathbf{y}_{1:T}\mid\bm{\theta})$ on the right hand side of (16) is typically intractable in SSMs and needs to be estimated using a particle filter (e.g., Nemeth et al.,, 2016). However, if the exact likelihood is replaced with the Whittle likelihood, an approximation to ${\nabla_{\bm{\theta}}\mathcal{L}(\bm{\theta})}$ can be obtained without the need for running an expensive particle filter repeatedly. Rather than targeting the true posterior distribution $p(\bm{\theta}\mid\mathbf{y}_{1:T})$, HMC-Whittle proceeds by targeting an approximation $\tilde{p}(\bm{\theta}\mid\mathbf{y}_{1:T})\propto p_{W}(\mathbf{y}_{1:T}\mid\bm{\theta})p(\bm{\theta})$, where $p_{W}(\mathbf{y}_{1:T}\mid\bm{\theta})\equiv\exp(l_{W}(\bm{\theta}))$ is the Whittle likelihood as defined in (6) for a univariate time series, or (10) for multivariate time series. Then $\mathcal{L}(\bm{\theta})$ in (14) is replaced by the quantity $\mathcal{L}_{W}(\bm{\theta})=\log\tilde{p}(\bm{\theta}\mid\mathbf{y}_{1:T})$, and the term $\nabla_{\bm{\theta}}\mathcal{L}(\bm{\theta})$ in (15) is replaced with

For models where the spectral density is available in closed form, the Whittle likelihood is also available in closed form, and the gradient $\nabla_{\bm{\theta}}\mathcal{L}_{W}(\bm{\theta})$ can be obtained analytically.

In practice, the continuous-time Hamiltonian dynamics in (15) are simulated over $L$ discrete time steps of size $\epsilon$ using a leapfrog integrator (see Algorithm 4). The resulting $\bm{\theta}^{*}$ and $\mathbf{r}^{*}$ from this integrator are then used as proposals and are accepted with probability

This process of proposing and accepting/rejecting is repeated until the desired number of posterior samples has been obtained. We summarise HMC-Whittle in Algorithm 5.

We generate $T=10000$ observations from the linear Gaussian SSM,

We set the true auto-correlation parameter to $\phi=0.9$, the true state disturbance standard deviation to $\sigma_{\eta}=0.7$ and the true measurement error standard deviation to $\sigma_{\epsilon}=0.5$.

The parameters that need to be estimated, $\phi,\sigma_{\eta}$ and $\sigma_{\epsilon}$, take values in subsets of $\mathbb{R}$. We therefore instead estimate the transformed parameters ${\bm{\theta}=(\theta_{\phi},\theta_{\eta},\theta_{\epsilon})^{\top}}$, where $\theta_{\phi}=\tanh^{-1}(\phi)$, $\theta_{\eta}=\log(\sigma_{\eta}^{2})$, and $\theta_{\epsilon}=\log(\sigma_{\epsilon}^{2})$, which are unconstrained.

The spectral density of the observed time series $\{y_{t}\}_{t=1}^{T}$ is

where

is the spectral density of the AR(1) process $\{x_{t}\}_{t=1}^{T}$, and

is the spectral density of the white noise process $\{\epsilon_{t}\}_{t=1}^{T}$. The Whittle log-likelihood for this model is given by (6), where the periodogram $\{I(\omega_{k}):k=-\left\lceil\frac{T}{2}\right\rceil+1,\dots,\left\lfloor\frac{T}{2}\right\rfloor\}$ is the discrete Fourier transform (DFT) of the observations $\{y_{t}\}_{t=1}^{T}$, which we compute using the fft function in R. We evaluate the gradient and Hessian of the Whittle log-likelihood at each frequency, $\nabla_{\theta}l_{W_{k}}(\bm{\theta})$ and $\nabla_{\theta}^{2}l_{W_{k}}(\bm{\theta})$, using automatic differentiation via the tensorflow library (Abadi et al.,, 2015) in R. The initial/prior distribution we use is

This leads to 95% probability intervals of $(-0.96,0.96)$ for $\phi$, and $(0.23,1.59)$ for $\sigma_{\eta}$ and $\sigma_{\epsilon}$.

Figure 3 shows the trajectory of the R-VGA-Whittle variational mean for each parameter. There are only 226 R-VGA-Whittle updates, as the 3dB cutoff point is the 176th angular frequency, and there are 50 blocks thereafter. The figure shows that the trajectory of $\phi$ moves towards the true value within the first 10 updates. This is expected, since most of the information on $\phi$, which is reasonably large in our example, is in the low-frequency components of the signal. The trajectories of $\sigma_{\eta}$ and $\sigma_{\epsilon}$ are more sporadic before settling close to the true value around the last 25 (block) updates. This later convergence is not unexpected since high-frequency log-likelihood components contain non-negligible information on these parameters that describe fine components of variation (the spectral density of white noise is a constant function of frequency). This behaviour is especially apparent for the measurement error standard deviation $\sigma_{\epsilon}$.

Figures 8(a) and 8(b) in Section S5 of the online supplement show the trace plots for each parameter obtained from HMC-exact and HMC-Whittle. Trace plots from both methods show good mixing for all three parameters but, for $\sigma_{\eta}$ and $\sigma_{\epsilon}$ in particular, HMC-Whittle has slightly better mixing than HMC-exact.

Figure 4 shows the marginal posterior distributions of the parameters, along with bivariate posterior distributions. The three methods yield very similar posterior distributions for all parameters. We also show in Section S3 of the online supplement that the posterior densities obtained from R-VGA-Whittle using block size $B=100$ are highly similar to those obtained without block updates. This suggests that block updates are useful for speeding up R-VGA-Whittle without compromising the accuracy of the approximate posterior densities.

A comparison of the computing time between R-VGA-Whittle, HMC-Whittle and HMC-exact for this model can be found in Table 1. The table shows that R-VGA-Whittle is nearly 10 times faster than HMC-Whittle and more than 30 times faster than HMC-exact.

SV models are a popular class of non-Gaussian financial time series models used for modelling the volatility of stock returns. An overview of SV models and traditional likelihood-based methods for parameter inference for use with these models can be found in Kim et al., (1998); an overview of MCMC-based inference for a univariate SV model can be found in Fearnhead, (2011).

In this section, we follow Fearnhead, (2011) and simulate $T=10000$ observations from the following univariate SV model:

We set the true parameters to $\phi=0.99$, $\sigma_{\eta}=0.1$, and $\kappa=2$.

Following Ruiz, (1994), we apply a log-squared transformation,

where $\mu\equiv\log(\kappa^{2})+\mathbb{E}(\log(\epsilon_{t}^{2}))$ and $\xi_{t}\equiv\log(\epsilon_{t}^{2})-\mathbb{E}(\log(\epsilon_{t}^{2}))$, for $t=1,\dots,T$. The terms $\{\xi_{t}\}_{t=1}^{T}$ are independent and identically distributed (iid) non-Gaussian terms with mean zero and variance $\frac{\pi^{2}}{2}$ (see Ruiz, (1994) for a derivation, which is reproduced in Section S2 of the online supplement for completeness).

Rearranging (27) yields the de-meaned series

Since $\mathbb{E}(\log(y_{t}^{2}))=\mu$, the parameter $\mu$ can be estimated using the sample mean:

As $\mu=\log(\kappa^{2})+\mathbb{E}(\log(\epsilon_{t}^{2}))$, once an estimate $\hat{\mu}$ is obtained, $\kappa$ can be estimated as

where $\mathbb{E}(\log(\epsilon_{t}^{2})))=\psi(1/2)+\log(2)$, and $\psi(\cdot)$ is the digamma function (Abramowitz and Stegun,, 1968). As we show below, the estimate $\hat{\kappa}$ is not required for R-VGA-Whittle and HMC-Whittle as the Whittle likelihood does not depend on $\kappa$. However, $\kappa$ appears in the likelihood used by HMC-exact. Since $\kappa$ is not estimated by the Whittle-based methods, for comparison purposes we fix it to its estimate $\hat{\kappa}$ when using HMC-exact.

The remaining parameters to estimate are $\phi$ and $\sigma_{\eta}$. As with the linear Gaussian state space model, we work with the transformed parameters $\bm{\theta}=(\theta_{\phi},\theta_{\eta})^{\top}$, where $\theta_{\phi}=\tanh^{-1}(\phi)$ and $\theta_{\eta}=\log(\sigma_{\eta}^{2})$, which are unconstrained. We fit the spectral density

to the log-squared, demeaned series $\{z_{t}\}_{t=1}^{T}$, where

and $f_{\xi}(\omega_{k};\bm{\theta})$ is the Fourier transform of the covariance function $C_{\xi}(h;\bm{\theta})$ of $\{\xi_{t}\}_{t=1}^{T}$:

and is thus given by

The prior/initial variational distribution we use is

which leads to 95% probability intervals of $(0.57,0.99)$ for $\phi$, and $(0.11,0.44)$ for $\sigma_{\eta}$. The interval for $\phi$ aligns closely with those used in similar studies of the stochastic volatility (SV) model. For example, Kim et al., (1998) employed the transformation $\phi=2\phi^{*}-1$ with $\phi^{*}\sim\textrm{Beta}(20,1.5)$, resulting in a 95% probability interval for $\phi$ of $(0.59,0.99)$. Additionally, Kim et al., (1998) utilized an inverse gamma prior for $\sigma_{\eta}^{2}$ with shape and rate parameters of $2.5$ and $0.025$, respectively, corresponding to a 95% probability interval $\sigma_{\eta}$ for $(0.06,0.24)$, which is comparable to the interval implied by our chosen prior.

Figure 5 shows the trajectories of the R-VGA-Whittle variational means for all parameters. The 3dB cutoff point for this example is the 75th frequency, after which there are 50 blocks, for a total of 125 updates. Similar to the case of the LGSS model in Section 3.1, the variational mean of $\phi$ is close to the true value from the early stages of the algorithm, while that of $\sigma_{\eta}$ requires more iterations to achieve convergence.

Figures 9(a) – 9(c) in Section S5 of the online supplement show the HMC-exact and HMC-Whittle trace plots. For this example, HMC-exact produces posterior samples with a much higher correlation than HMC-Whittle, particularly for the variance parameter $\sigma_{\eta}$. This is because HMC-Whittle only targets the posterior density of the parameters, while HMC-exact targets the joint posterior density of the states and parameters, which are high dimensional. We thin the posterior samples from HMC-exact by a factor of 100 for the remaining samples to be treated as uncorrelated draws from the posterior distribution.

We compare the posterior densities from R-VGA-Whittle, HMC-Whittle, and HMC-exact in Figure 6. The posterior densities from R-VGA-Whittle and HMC-Whittle are very similar, while the posterior densities from HMC-exact are slightly narrower. Table 1 shows that R-VGA-Whittle is 8 times faster than HMC-Whittle, and more than 1000 times faster than HMC-exact in this example.

We now consider the following bivariate extension of the SV model in Section 3.2 and simulate $T=5000$ observations as follows:

where $\mathbf{x}_{t}\equiv(x_{1,t},x_{2,t})^{\top}$, $\mathbf{y}_{t}\equiv(y_{1,t},y_{2,t})^{\top}$, ${\mathbf{V}_{t}\equiv\operatorname{diag}(\exp(\tfrac{x_{1,t}}{2}),\exp(\tfrac{x_{2,t}}{2}))}$, and $\mathbf{I}_{m}$ is the $m\times m$ identity matrix. The covariance matrix $\mathbf{\Sigma}_{\eta_{1}}$ is the solution to the equation $\operatorname{vec}(\mathbf{\Sigma}_{\eta_{1}})=(\mathbf{I}_{4}-\mathbf{\Phi}\otimes\mathbf{\Phi})^{-1}\operatorname{vec}(\mathbf{\Sigma}_{\eta})$, where $\textrm{vec}(\cdot)$ is the vectorisation operator and $\otimes$ denotes the Kronecker product of two matrices. We let $\mathbf{\Phi}$ be a diagonal matrix but let $\mathbf{\Sigma}_{\eta}$ have off-diagonal entries so that there is dependence between the two series. We set these matrices to

As with the univariate model, we take a log-squared transformation of (32) to find

where

The terms $\xi_{1,t}$ and $\xi_{2,t}$ are iid non-Gaussian uncorrelated terms, each with mean zero and variance $\sigma_{\xi}^{2}=\frac{\pi^{2}}{2}$. We express (36) as