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.
126,863 characters · 16 sections · 61 citation commands
Ancillarity-Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Estimation of Stochastic Volatility Models
\hrule
\hrule
{\bf Keywords:} Markov Chain Monte Carlo, Non-Centering, Auxiliary Mixture Sampling, Massively Parallel Computing, State Space Model, Exchange Rate Data
Returns of financial and economic time series often exhibit time-varying volatilities. To account for this behavior, tay:fin suggests in his pioneering paper to model the logarithm of the squared volatilities by latent autoregressive processes of order one. This specification, commonly referred to as the stochastic volatility (SV) model, presents itself as a competitive alternative to GARCH-type designs by modeling the volatilities non-deterministically. Also, it arises naturally as a discretization of continuous-time models frequently appearing in the mathematical finance literature hul-whi:pri.
Following e.g.\ jac-etal:bayJBES or kim-etal:sto, observed log-returns are denoted $\mathbf y = (y_1, y_2, \dots, y_T)'$ and the SV model is specified as
where it is assumed that the iid standard normal innovations $\epsilon_t$ and $\eta_s$ are independent for $t, s \in \{1,\dots,T\}$. The unobserved process ${\mathbf h}=(h_{0}, h_{1}, \ldots,h_{T})$ appearing in state equation ((ref)) is usually interpreted as the latent time-varying volatility process with initial state distributed according to the stationary distribution, i.e. $ h_0|\mu,\phi,\sigma \sim \mathcal{N}\!\left(\mu,\sigma^2/(1-\phi^2)\right). $ From now on, we will refer to equations ((ref)) and ((ref)) as the SV model in its centered parameterization (C).
Simulation efficiency in state-space models can often be improved through model reparameterization. Papers related to this matter include gel-etal:eff, pit-she:ana, fru:eff, rob-etal:bay04, fru-soe:bayhes, and str-etal:par. The pioneering paper by tay:fin as well as several other papers like kim-etal:sto or lie-ric:cla consider a partially non-centered parameterization, where the level $\mu$ of $h_{t}$ -- which defines the scale of $y_{t}$ -- is shifted from the state equation ((ref)) to the observation equation ((ref)) by setting $\bar h_t = h_t - \mu$. kim-etal:sto compare both parameterizations within a Bayesian inference. They show that the partially non-centered parameterization leads to very high inefficiency when sampling $\mu$ and recommend choosing the centered parameterization in any case. Nevertheless, the centered parameterization has several disadvantages. Firstly, inefficiency when drawing $\sigma$ is still high, see e.g. Table 1 in kim-etal:sto. Secondly, the conclusions are only valid if $\phi$ is close to one, which is commonly the case when the SV model is applied to capture conditional heteroskedasticity of observed financial times series. However, this is not necessarily true when the SV model is applied in more general contexts such as capturing conditional heteroskedasticity in latent variables or regression residuals.
For the purpose of this paper, the (fully) non-centered parameterization (NC), given through
where $\omega=\operatorname{e}^\mu$, is of particular importance. The initial value of $\tilde h_{0}|\phi$ is once again drawn from the stationary distribution of the latent process, i.e.\ $\tilde h_{0}|\phi \sim \mathcal{N}\!\left(0, 1/(1- \phi^2)\right)$. Note that $\tilde h_t = (h_t-\mu)/\sigma$. For a moderate parameter range where $\phi_\text{true} \in \{0.8, 0.9, 0.95\}$ and $\sigma_\text{true} \in \{0.2, 0.3, 0.4\}$, str-etal:par illustrate that this type of non-centering typically yields lowest inefficiency factors when estimating stochastic volatility and stochastic conditional duration models with randomly sized block updating. Also, in similar contexts, there are several papers showing that MCMC sampling improves a lot by considering a non-centered version of a state space model, see e.g. fru:eff and fru-wag:sto. These authors show that non-centering is especially useful if the error variance in the state equation is considerably smaller than the error variance in the observation equation. pit-she:ana show for linear Gaussian state space models that the speed of convergence in the centered parameterization decreases as $|\phi|$ increases when the signal to noise ratio is fixed.
No matter which parameterization is chosen, the likelihood in the SV model has an intractable form. Thus, Bayesian estimation commonly relies on sampling the latent states $\mathbf{h}$ and treat these as known for updating the parameters $\mu$, $\phi$, and $\sigma$. In their seminal paper, jac-etal:bayJBES propose a single-move Metropolis-Hastings (MH) algorithm. Each individual $h_t$ is sampled conditional on past and future, i.e. drawn from $ p(h_t|\mathbf{h}_{[-t]},\sigma,\phi,\mu,{\mathbf y}), $ where $\mathbf{h}_{[-t]}$ denotes all elements of $\mathbf{h}$ except $h_t$. Due to the commonly high persistence of the latent process, she-kim:bay note that the draws obtained from this sampler are also highly correlated and thus only slowly converging to the stationary distribution. Alternatively, she-pit:lik propose a multi-move sampler, where volatility blocks of random length are updated at a time, while she:par and omo-etal:sto propose a method to draw directly from $ p(\mathbf{h}|\sigma,\phi,\mu,{\mathbf y}). $ This becomes possible through a normal mixture approximation of $\log(\epsilon_t^2)$ and requires forward filtering backward sampling (FFBS) methods car-koh:ong,fru:dat,dur-koo:sim. Within a more general Gaussian state-space framework, rue:fas and mcc-etal:sim propose sampling the latent volatilities through Cholesky-factorization of the precision matrix by exploiting its band-diagonal structure. We adopt this method to sample the latent volatilities “all without a loop” (AWOL). For a more extensive review of both Bayesian and non-Bayesian SV estimation methods, see bos:rel.
The contribution of this paper is threefold. Firstly, we explore the impact of alternative parameterizations for a wide parameter range including empirically plausible values and more extreme ones that can be relevant for applications of SV models within more general frameworks such as SV factor models or regression analysis. It turns out that simulation efficiency heavily depends on the true parameter values of the data generating process, thus no single “best” parameterization exists. Secondly, we provide a strategy to overcome this deficiency by interweaving C and NC utilizing an ancillarity-sufficiency interweaving strategy (ASIS) introduced by yu-men:cen. This results in a robustly efficient sampler that always outperforms the more efficient parameterization with respect to all parameters at little extra cost in terms of design and computation. Thirdly, we provide evidence that empirical sampling efficiency heavily depends on the realization of the data generating process by massively parallel simulation experiments.
The paper is structured as follows: Section (ref) gives insight into the estimation procedure for each of the two selected parameterizations. Section (ref) explains how ASIS can be applied in order to interweave these parameterizations. Extensive simulation results presented in Section (ref) compare sampling efficiency for all parameters amongst the different parameterizations, Section (ref) provides real-data results for several daily exchange rates, and Section (ref) concludes.
To perform Bayesian inference, a prior distribution $p(\mu,\phi,\sigma)$ needs to be specified. For both parameterizations, we choose the same independent components for each parameter. The level $\mu \in \mathbb{R}$ is equipped with the usual normal prior $\mu \sim \mathcal{N}\!\left(b_{\mu}, B_{\mu}\right)$. For the persistence parameter $\phi \in (-1,1)$, we choose $(\phi+1)/2 \sim \mathcal{B}\left(a_0, b_0\right)$ as in kim-etal:sto, implying
Clearly, the support of this distribution is the unit ball and thus guarantees stationarity of the autoregressive volatility process. For the volatility of volatility $\sigma \in \mathbb{R}^+$, we choose $\sigma^2 \sim B_\sigma\cdot \chi^2_1= \mathcal{G}\left(\frac{1}{2},\frac{1}{2B_\sigma}\right)$. Note that this specification differs from the commonly employed conjugate Inverse-Gamma prior $\sigma^2 \sim \mathcal{G}^{-1} \left(c_0,C_0\right)$ and is motivated by fru-wag:sto, who equivalently stipulate the prior for $\pm \sqrt { \sigma^2}$ to follow a centered normal distribution, i.e.\ $\pm \sqrt { \sigma^2} \sim \mathcal{N}\!\left(0,B_\sigma\right)$. It turns out that this choice is less influential when the true volatility of volatility is small because $\sigma$ is not bound away from zero a priori.
Observation equation ((ref)) can easily be rewritten as
where $\tilde y_t$ denotes $\log y_t^2$. Alternatively, $\tilde y_t$ can be interpreted as the transformed de-meaned returns $\log (y_t-\bar y)^2$, but might just as well be taken $\log ( (y_t-\bar y)^2 + c)$ with a fixed offset constant $c=10^{-3}$ as in kim-etal:sto or $\log ( y_t^2 + c)$ with $c= 10^{-4}$ as in omo-etal:sto in order to avoid values equal to zero. Equation ((ref)) now takes the form of a linear but non-Gaussian state space model. Moreover, one can approximate the distribution of $\log(\epsilon_t^2)$ by a mixture of normal distributions, i.e.\ $ \log(\epsilon_t^2)|r_t \sim \mathcal{N}\!\left(m_{r_t},s_{r_t}^2\right). $ Here, $r_t\in\{1,\ldots,10\}$ defines the mixture component indicator at time $t$, while $m_{r_t}$ and $s^2_{r_t}$ denote mean and variance of the $r_t$th mixture component as tabulated in omo-etal:sto. This representation allows rewriting ((ref)) as a linear and conditionally Gaussian state space model,
where speedy MCMC sampling becomes possible in three steps.
Conditional on all other variables, the joint density for ${\mathbf h}$ (and $\tilde {\mathbf h}$) is multivariate normal. Due to the order-one autoregressive nature of the latent volatility process, this distribution can be written in terms of the tridiagonal precision matrix $\mathbf{\Omega}$, giving rise to sampling all without a loop (AWOL). This method is employed in rue:fas and mcc-etal:sim and does not require the “end-user” to implement any loops -- hence the name. Thus, it is very convenient in terms of implementation and fast in terms of computation. No FFBS methods are needed, there is no need to invert the tridiagonal precision matrix $\mathbf{\Omega}$ and it is fast due to the availability of band back-substitution already implemented in practically all widely used programming libraries.
In the centered parameterization, we draw from ${\mathbf h}_{[-0]}|\mu, \sigma, \phi, \mathbf{r}, {\mathbf y} \sim N _{T}\left(\mathbf{\Omega}^{-1}\mathbf{c}, \mathbf{\Omega}^{-1}\right)$ with \[ \mathbf{\Omega}=
, \] and \[ \mathbf{c}=
. \]
Analogously, in the noncentered case, we draw from $\tilde {\mathbf h}_{[-0]}|\mu, \sigma, \phi, \mathbf{r}, {\mathbf y} \sim N _{T}\left(\mathbf{\Omega}^{-1}\mathbf{c}, \mathbf{\Omega}^{-1}\right)$ with \[ \mathbf{\Omega}=
, \] and \[ \mathbf{c}=
. \]
For both parameterizations, this is accomplished by first computing the Cholesky decomposition $\mathbf{\Omega}=\mathbf{LL}'$. Due to the band structure of $\mathbf{\Omega}$, this is computationally inexpensive and can either be implemented directly or via the LAPACK-routine \verb.dpbtrf. lapack, to name but one of the many widely available (and thoroughly tested) linear algebra routines designed for this task. Note that only main diagonal and lower first off-diagonal elements of $\mathbf{L}$ will be nonzero. Next, we draw $\boldsymbol{\epsilon} \sim N _{T}\left(\boldsymbol{0},\mathbf{I}_T\right)$ and then efficiently solve $\mathbf{L}\mathbf{a}=\mathbf{c}$ for $\mathbf{a}$ and $\mathbf{L}'{\mathbf h}=\mathbf{a}+\boldsymbol{\epsilon}$ for ${\mathbf h}$ by using band back-substitution instead of actually calculating $\mathbf{L}^{-1}$. Finally, the initial value can be sampled from $h_0|h_1,\mu,\phi,\sigma \sim \mathcal{N}\!\left(\mu+\phi(h_{1}-\mu), \sigma^2\right)$ in C and from $\tilde h_0|\tilde h_1,\phi \sim N(\tilde h_1\phi,1)$ in NC.
For sampling $\boldsymbol{\theta} = (\mu, \phi, \sigma^2)$, it is helpful to rewrite the conditional AR(1) model as a conditional regression model with the lagged latent variables as regressors, \[ h_t=\gamma + \phi h_{t-1} + \eta_t, \quad \eta_t \sim \mathcal{N}\!\left(0,\sigma^2\right), \] via $\gamma=(1-\phi)\mu$. Note that the implied conditional prior $p(\gamma|\phi)$ follows a normal distribution with mean $b_{\mu}(1-\phi)$ and variance $B_\mu(1-\phi)^2$. In this Subsection, we will discuss three common blocking strategies for sampling $\boldsymbol{\theta}$.
For a one block update of $\boldsymbol{\theta}$, we use a single MH step. The posterior arising from an auxiliary regression model with conjugate priors is used as the proposal density: \[ p_\text{aux}( \boldsymbol{\theta}_\text{new} |{\mathbf h}) = p_\text{aux}(\gamma_\text{new}, \phi_\text{new} |{\mathbf h}, \sigma^2_\text{new}) p_\text{aux}(\sigma^2_\text{new}|{\mathbf h}). \] We choose $p_\text{aux}(\sigma^2) \propto\sigma^{-1} $ to denote the density of an auxiliary improper conjugate $\mathcal{G}^{-1} \left(-\frac{1}{2}, 0\right)$ prior, and $p_\text{aux}(\gamma,\phi|\sigma^2) $ to denote the density of an auxiliary conjugate $N _{2}\left(\mathbf{0}, \sigma^2\mathbf{B}_0\right)$ prior with $\mathbf{B}_0=\text{diag}(B_0^{11}, B_0^{22})$. More specifically, $p_\text{aux}(\gamma|\sigma) \sim \mathcal{N}\!\left(0, \sigma^2B_0^{11}\right)$ and $p_\text{aux}(\phi|\sigma) \sim \mathcal{N}\!\left(0, \sigma^2B_0^{22}\right)$. In order to avoid collinearity problems when $\sigma^2$ is close to zero (and thus $h_t$ almost constant for all $t$), we pick slightly informative variances, i.e. $B_0^{11} = 10^{12}$ and $B_0^{22} = 10^{8}$. This yields
with $\mathbf{B}_T = (\mathbf{X'X} + \mathbf{B}_0^{-1})^{-1}$ and $ \mathbf{b}_T = \mathbf{B}_T\mathbf{X}'\mathbf{h}_{[-0]}$, where $\mathbf{X}$ is the $T \times 2$ design matrix with ones in the first column and $\mathbf{h}_{[-T]}$ in the second. The marginalized auxiliary posterior distribution for $\sigma^2$ is given through $ \sigma^2 | {\mathbf h} \sim \mathcal{G}^{-1} \left(c_T, C_T\right), $ with $c_T = (T-1)/2$ and $C_T = \frac{1}{2}\left( \sum_{i=1}^T h_i^2 - \mathbf{b}_T' \mathbf{X}' {\mathbf h}_{[-0]} \right)$. The acceptance probability is given through $\min(1,R)$, with
In the two-block sampler, we draw the first block from the full conditional distribution $\gamma,\phi|{\mathbf h}, \sigma^2$ given in ((ref)) and accept with probability $\min(1,R)$, where \[ R= \frac {p(h_0|\gamma_\text{new},\phi_\text{new})p(\gamma_\text{new}|\phi_\text{new})p(\phi_\text{new})} {p(h_0|\gamma_\text{old},\phi_\text{old})p(\gamma_\text{old}|\phi_\text{old})p(\phi_\text{old})} \times \frac {p_\text{aux}(\gamma_{\text{old}},\phi_{\text{old}})} {p_\text{aux}(\gamma_{\text{new}},\phi_{\text{new}})}. \]
In order to construct a suitable proposal for the -- now full conditional -- density $p(\sigma^2|{\mathbf h},\mu,\phi)$, we again use the auxiliary conjugate prior $p_\text{aux}(\sigma^2) \propto\sigma^{-1} $, under which we straightforwardly obtain
where $c_T=T/2$ and $ C_T=\frac{1}{2}\left(\sum_{t=1}^T ((h_t-\mu)-\phi(h_{t-1}-\mu))^2 +(h_0-\mu)^2(1-\phi^2)\right). $ The acceptance probability simplifies to $\min(1,R)$ with
In the three-block sampler, each individual parameter is drawn from the full conditional distribution $\mu|\cdot$, $\phi|\cdot$, and $\sigma^2|\cdot$, respectively. Thus, $\sigma^2$ is drawn from ((ref)). For sampling $\phi$, we obtain a proposal from \[ \phi|{\mathbf h}, \gamma,\sigma^2 \sim \mathcal{N}\!\left( \frac{\left[\sum_{t=1}^Th_{t-1}h_t\right]-\gamma\sum_{t=0}^{T-1}h_t}{\sum_{t=0}^{T-1}h_t^2 + 1/B_0^{22}}, \frac{\sigma^2}{\sum_{t=0}^{T-1}h_t^2 + 1/B_0^{22}} \right). \] The acceptance probability is equal to $\min(1,R)$ with \[ R=\frac{p(h_0|\phi_\text{new}, \mu, \sigma^2)p(\phi_\text{new})}{p(h_0|\phi_\text{old}, \mu, \sigma^2)p(\phi_\text{old})}\times \frac{p_\text{aux}(\phi_\text{old}|\sigma^2)}{p_\text{aux}(\phi_\text{new}|\sigma^2)}. \] For sampling $\gamma$ from the full conditional posterior distribution, we obtain a proposal from \[ \gamma|{\mathbf h}, \phi,\sigma^2 \sim \mathcal{N}\!\left( \frac{\sum_{t=1}^T h_t-\phi\sum_{t=0}^{T-1}h_t}{T + 1/B_0^{11}}, \frac{\sigma^2}{T + 1/B_0^{11}}\right) \] and an acceptance probability equaling $\min(1,R)$ with \[ R=\frac{p(h_0|\gamma_\text{new}, \phi, \sigma^2)p(\gamma_\text{new}|\phi)}{p(h_0|\gamma_\text{old}, \phi, \sigma^2)p(\gamma_\text{old}|\phi)}\times \frac{p_\text{aux}(\gamma_\text{old}|\sigma^2)}{p_\text{aux}(\gamma_\text{new}|\sigma^2)}. \]
In the noncentered parameterization, only $\phi$ is left in the state equation. To sample this parameter, we employ a flat auxiliary prior $p_\text{aux}(\phi)\propto c$, yielding the proposal \[ \phi|\tilde {\mathbf h} \sim \mathcal{N}\!\left( \frac{\sum_{t=0}^{T-1}\tilde h_t\tilde h_{t+1}}{\sum_{t=0}^{T-1}\tilde h_t^2}, \frac{1}{\sum_{t=0}^{T-1}\tilde h_t^2} \right), \] and an acceptance probability of $\min(1,R)$, where $ R= {p(\tilde h_0|\phi^\text{new})p(\phi^\text{new})}/ { p(\tilde h_0|\phi^\text{old})p(\phi^\text{old})}. $
For sampling $\mu$ and $\sigma$, one can straightforwardly rewrite the conditional observation equation ((ref)) as a regression model with homoskedastic errors, i.e.\
where $\boldsymbol{\epsilon} \sim N _{K}\left(\mathbf{0},\mathbf{I}_K\right)$, and \[ \mathbf{\breve y}=
, \quad \mathbf{X}=
. \]
The joint posterior distribution is again bivariate Gaussian with variance-covariance matrix $\mathbf{B}_T=(\mathbf{B}_0^{-1}+\mathbf{X}'\mathbf{X})^{-1}$ and mean $\mathbf{b}_T = \mathbf{B}_T(\mathbf{B}_0^{-1}\mathbf{b}_0+\mathbf{X}'\mathbf{\breve y})$, where $\mathbf{b}_0=(b_\mu,0)'$ and $\mathbf{B}_0=\text{diag}(B_\mu, B_\sigma)$ denote mean and variance of the joint prior density $p(\mu,\sigma)$, respectively.
Alternatively, one could sample both parameters from the full conditional posteriors (three-block sampling), yielding $ \mu|{\mathbf y}, \tilde {\mathbf h},\mathbf{r},\sigma \sim \mathcal{N}\!\left(b_{T,\mu}, B_{T,\mu}\right) $ with \[ b_{T,\mu} = B_{T,\mu}\left( \sum_{t=1}^T \frac{\tilde y_t -m_{r_t}-\sigma \tilde h_t}{s^2_{r_t}} + \frac{b_\mu}{B_\mu} \right), \qquad B_{T,\mu} = 1/\left({\displaystyle \sum_{t=1}^T 1/s_{r_t}^2 +\frac{1}{B_\mu}}\right), \] and $ \sigma|{\mathbf y}, \tilde {\mathbf h},\mathbf{r},\mu \sim \mathcal{N}\!\left(b_{T,\sigma}, B_{T,\sigma}\right) $ with \[ b_{T,\sigma} = B_{T,\sigma} \sum_{t=1}^T \frac{\tilde h_t(\tilde y_t-m_{r_t}-\mu)}{s^2_{r_t}}, \qquad B_{T,\sigma} = 1/\left({\sum_{t=1}^T\frac{\tilde h_t^2}{s_{r_t}^2} +\frac{1}{B_\sigma}}\right). \]
We proceed exactly as omo-etal:sto. Observing that $\tilde y_t-h_t=\epsilon_t^*$ with $\epsilon_t^* \sim \mathcal{N}\!\left(m_{r_t},s^2_{r_t}\right)$, one easily obtains the posterior probabilities $\mathbb{P}(r_t=k|\cdot)$ for $k \in \{1,\dots, 10\}$ and $t \in \{1,\dots, T\}$ according to \[ \mathbb{P}(r_t=k|\cdot) \propto \mathbb{P}(r_t=k)\frac{1}{s_k}\exp \left\{-\frac{(\epsilon_t^*-m_k)^2}{2s_k^2} \right\}, \] where $\mathbb{P}(r_t=k)$ denotes the mixture weights of the $k$th component. In our implementation, we do the calculations on a $\log$-scale and normalize with respect to the maximum as required. The actual drawing is then conducted via inverse transform sampling. Note that due to $T\times10$ exponential function calls, this step is computationally rather expensive but can easily be parallelized.
To provide some intuition about the sampling efficiency in C, let $\phi=0$ for a moment. This implies that the state equation ((ref)) reduces to $h_t \sim \mathcal{N}\!\left(\mu,\sigma^2\right)$ iid for all $t \in \{1,\dots,T\}$. In this setting, ${\mathbf h}$ becomes more informative about $\mu$ when the conditional variance $\sigma^2$ gets smaller. Thus, more information is missing when treating ${\mathbf h}$ as latent data. Consequently, when sampling $\mu$ under the assumption that $\phi = 0$ and $\sigma^2$ is small, we expect C to be inefficient. On the other hand, if $\phi$ approaches $1$, the latent process converges towards a random walk and ${\mathbf h}$ will be very uninformative about $\mu$. Thus, only little information is lost when treating ${\mathbf h}$ as latent data and C has better chances to work fine. In NC, no major troubles are to be expected if $\phi=0$, since the state equation ((ref)) reduces to $\tilde h_t \sim \mathcal{N}\!\left(0,1\right)$ iid for all $t \in \{1,\dots,T\}$, which is obviously independent of the value of $\sigma$. Thus, sampling $\mu$ and $\sigma$ in the linearized equation ((ref)) reduces to simple linear regression with independent regressors. If, however, $\phi$ goes towards $1$, we are prone to running into spurious regression problems. Certainly these arguments rely on massive oversimplification (e.g.\ by not taking into account the impact of the mixture approximation or spillover effects by inefficient proposal densities and different blocking strategies) and can only provide a faint idea of what is going on in the general case.
Nevertheless, due to the fact that in the context of the model at hand, the latent variables ${\mathbf h}$ in C form a sufficient statistic for $\mu$ and $\sigma$, while the transformed volatilities $\tilde {\mathbf h}$ in NC form an ancillary statistic for these parameters, there is hope that interweaving C and NC helps to increase sampling efficiency. yu-men:cen propose an ancillary-sufficiency interweaving strategy (ASIS) which, in certain situations, converges geometrically even when C and/or NC fail to do so. They explain this “seemingly magical property” by relating to Basu's theorem bas:ons on the independence of complete sufficient and ancillary statistics and show in a quite general context that the geometric convergence rate of the sampler interweaving ${\mathbf h}$ and $\tilde {\mathbf h}$ is always bound by $R\sqrt{r_\text{C}r_\text{NC}}$, where $R$ is the maximal correlation between ${\mathbf h}$ and $\tilde {\mathbf h}$ in their joint posterior distribution $p({\mathbf h}, \tilde {\mathbf h}|\mathbf{y})$ and $r_\text{C},r_\text{NC}$ denote the geometric rate of convergence of C and NC, respectively. This means that the rate of convergence of the interwoven sampler is mainly governed by the individual convergence rates and the posterior correlation $R$, implying that ancillary-sufficiency pairs of latent variables are likely to be good candidates for reducing sampling inefficiency. It is worth noting that the original ASIS notation $Y_{obs}$ for the observed data directly transforms to ${\mathbf y}$ for the model at hand, while $Y_{mis}$ -- denoting the “missing” part of the data -- equals $\mathbf{h}$.
The idea of interweaving is surprisingly simple. It is based on sampling the parameters in question -- in our case $\mu$, $\sigma$ (and $\phi$) -- twice: once utilizing C and again utilizing NC. Ad hoc, it is not clear whether one should start C and redraw NC (“baseline C”) or vice versa (“baseline NC”). We will discuss both strategies and assess their performance individually. Algorithm (ref) below describes the former, i.e.\ both the latent volatilities and the indicators are sampled once with baseline C, while the parameters $\mu, \phi, \sigma$ are sampled once in each parameterization within each iteration of the sampler. It is termed “GIS-C”, where the first three letters are borrowed from yu-men:cen and stand for global interweaving strategy. C simply denotes the fact that we use the centered baseline.
The individual sampling steps are implemented exactly as described in subsections (ref) to (ref). Note that since $\phi$ is not involved in the reparameterization, in step (b**), one might as well redraw $\mu$ and $\sigma$ only; the difference concerning sampling efficiency is however negligible. Also note that although additional sampling steps are introduced as (b*) to (b***), overall sampling time is only affected minimally because these steps are very cheap in terms of computation cost.
The sampler with noncentered baseline is of course very similar. As before, for each iteration the parameters $\mu$, $\phi$, and $\sigma$ are sampled twice (once in C and once in NC), while the latent volatilities and the indicators are sampled in NC only.
To conclude, note that the strategy of interweaving is intrinsically different to alternating the parameterizations, for instance by (randomly) choosing one parameterization and running a complete MCMC cycle within that parameterization. Also, it is distinct from compromising between two parameterizations, e.g.\ by partial noncentering.
In order to assess simulation efficiency of our algorithms, we simulate data from the model specified in equations ((ref)) and ((ref)). For the sake of simplicity and readability, $\mu_\text{true}$ is set to $-10$ for all runs. Results not reported here show that this choice is of minor influence. The parameters $\phi_\text{true}$ and $\sigma_\text{true}$ vary on a $\{0,0.5,0.8,0.9,0.95,0.96,0.97,0.98,0.99\}\times\{0.5,0.4,0.3,0.2,0.1\}$ grid, resulting in $45$ distinct parameter settings. This choice includes previously investigated and empirically plausible values, see e.g.\ jac-etal:bayJBES, kim-etal:sto, lie-ric:cla, and str-etal:par. Moreover, the range is chosen to also include more extreme values that frequently arise when univariate SV is applied to capture conditional heteroskedasticity in latent variables such as factors or residuals of regression-type problems. We repeat this exercise for $500$ data sets and apply four sampling schemes (C, NC, GIS-C, GIS-NC) by using $M=100\,000$ MCMC draws after a burn-in of $10\,000$ for each data set. Time series length is fixed to $T=5000$, which corresponds to just above 20 years of daily data. Overall, this results in $90\ 000$ chains of length $110\ 000$, or a total of around $50$ trillion latent instantaneous volatility draws. Nevertheless, due to parallel implementation of native C code on our local computer cluster using 500 cores, sampling can easily be done overnight. Throughout all simulations we use priors with means equaling the true values, more specifically $b_\mu=\mu_{\text{true}}$, $B_\mu = 10$, $a_0=40$, $b_0=80/(1+\phi_\text{true})-40$, $B_\sigma = \sigma_\text{true}^2$, and starting values are set to true values to avoid values outside the stationary distribution after the burn-in period.
Computation of parallel MCMC chains for each parameter constellation was conducted on a cluster of workstations consisting of 44 IBM dx360M3 nodes with a total of 544 cores running R\ 2.15.1 r:r and OpenMPI 1.4.3 gab-etal:ope. For high-level-parallelization and parallel random number generation according to lec-etal:obj, the R packages \verb.parallel. (part of R) and \verb.snow. r:sno were used. Ex-post analysis and timing was done on a Laptop with a 2.67GHz Intel i5 M560 CPU running the same R version. For the actual sampling, the R package \verb.stochvol. r:sto, available on CRAN, was created. The core implementation is written in C, interfaced to R via \verb.Rcpp. edd-fra:rcp. Inefficiency factors and effective sample sizes were computed with the R package \verb.coda. plu-etal:cod.
The mean time for running $1000$ simulation draws varies between $2.3$ seconds for C and $2.4$ seconds for GIS-NC on a Laptop with a 2.67GHz Intel i5 M560 CPU using one core. Note that these numbers are fairly constant for all true parameter values and grow linearly with $T$. As an example, the time to run $1000$ simulations for $T=500$ varies between $0.23$ and $0.24$ seconds.
Simulation efficiency of the two raw parameterizations mainly depends on the values of the parameters $\phi$ (persistence) and $\sigma$ (volatility of volatility). To illustrate the latter, Figure (ref) shows autocorrelations of an exemplary parameter setup with small volatility of volatility $\sigma_\text{true}=0.1$ for a single time series that has been randomly selected from the pool of all $500$ time series. Here, C “fails” in the sense that the draws from $p(\mu|\mathbf{y})$ and $p(\sigma|\mathbf{y})$ exhibit large autocorrelation, while NC performs substantially better. This observation is in line with findings of pit-she:ana and fru:eff, who observe that simulation efficiency in the centered parameterization decreases with decreasing $\sigma_\text{true}$ for linear Gaussian state space models.
On the other hand, Figure (ref) portraits a parameter setup with larger volatility of volatility $\sigma_\text{true}=0.5$, while persistence $\phi_\text{true}$ and level $\mu_\text{true}$ are the same as before. Here, we see that draws from C show little autocorrelation, while MCMC chains obtained from NC do not mix well.
For assessing simulation efficiency, the inefficiency factor (IF) is employed as a benchmark. It is an estimator for the integrated autocorrelation time $\tau$ of a stochastic process given through $ \tau = 1+2\sum_{s=1}^\infty\rho(s), $ where $\rho(s)$ is the autocorrelation function for lag $s$. We estimate $\tau$ through the spectral density of the Markov chain, i.e.\ $\text{IF}=\gamma_0/s^2$, where $\gamma_0$ denotes the estimated spectral density evaluated at zero and $s^2$ denotes the sample variance of the MCMC draws. The inefficiency factor is directly proportional to the squared Monte Carlo standard error MCSE$^2$ through the relationship $\text{MCSE}^2 = \frac{s^2}{M}\times \text{IF}$. In other words, $100\,000$ draws from a Markov chain with an IF of $100$ have roughly the same MCSE as $1000$ draws from an independent sample. Consequently, the effective sample size ESS is given by $M/\text{IF}$. Clearly, the aim is to provide samplers with small IFs, thus large ESSs, at smallest possible computational cost.
Even for artificially created datasets of length $T=5000$ or larger, estimation results may depend substantially on the actual realization of the underlying process. Also, other factors -- most importantly the initial seed for drawing pseudo random variables in the individual MCMC steps -- can influence both sample statistics from the posterior distribution as well as sample statistics for evaluating simulation efficiency. To compensate for this fact, we repeat each simulation with $500$ independently generated artificial data sets. The boxplots provided in Figure (ref) and Figure (ref) illustrate the variation of IFs for the same parameter constellations as above. While the overall picture about non-centering remains the same, we can now observe some substantial deviation from the median for certain realizations. Furthermore, the plots show that both interweaving strategies GIS-C and GIS-NC help to avoid woe by working well no matter which raw parameterization fails.
In order to gain insight into the entire parameter range of interest, Tables (ref), (ref) and (ref) provide a summary of median inefficiency factors across all $45$ parameter constellations.
Median IFs obtained from draws from $p(\mu|\mathbf{y})$ in Table (ref) confirm clearly that the centered parameterization is quite capable of efficiently estimating the level $\mu$ of the latent process throughout a wide parameter range, no matter which blocking strategy is used. Only a combination of both small $\sigma_\text{true}$ and small $\phi_\text{true}$ leads to large inefficiency. As was to be expected, this is exactly the area where the non-centered parameterization performs comparably well; median IFs are small to moderately large. On the other end of the scale -- where we find both highly persistent and highly varying latent variables -- NC becomes close to useless with very large IFs of $1000$ and above for both blocking strategies.
The lower three panels of Table (ref) show the performance of the interwoven samplers with different baselines, GIS-C and GIS-NC. It stands out that in terms of simulation efficiency, both variants are always better than or en par with the ideal parameterization, while there are practically no differences between the sampler with baseline C and the one with baseline NC. Comparing CPU time of the raw samplers with their interwoven counterparts reveals the computational cost of interweaving, which amounts to merely around $2\%$ in our setup. Thus, even when taking into account the extra cost, interweaving is hardly ever a bad choice. Also note that GIS-C is practically as fast as NC.
Table (ref) shows a direct comparison of the interwoven sampler with the raw parameterizations in terms of increase in effective sample size. All numbers are positive, showing that interweaving is more efficient than the ideal parameterization, but sometimes only slightly. Note that in comparison to the suboptimal parameterization, GIS is always at least twice as effective.
Next, we turn to assessing simulation efficiency for the persistence parameter $\phi$, summarized in Table (ref). Because $\phi$ is not involved in the reparameterization, the differences between C and NC (and consequently also between the raw and the interwoven samplers) are much less pronounced. For a summary of efficiency gains, see Table (ref). It stands out that one- and two-block samplers show very similar IFs, whereas the three-block sampler deteriorates due to massive overconditioning for moderate and small $\phi_\text{true}$ or $\sigma_\text{true}$. Note however that again the interwoven sampler is exempt from these defects due to the fact that NC performs solidly. Results not reported here show that for shorter time series with $T=500$, sampling inefficiency is uniformly smaller for all parameterizations. The interwoven 2-block samplers for instance show IFs of 30 or below for all underlying true parameter values.
Finally, we investigate sampling efficiency for $\sigma$. Table (ref) summarizes median IFs for draws from $p(\sigma|\mathbf{y})$. We observe a similar overall picture to the one presented in Table (ref): C performs poorly when $\sigma_\text{true}$ and $\phi_\text{true}$ are small, and NC performs poorly when $\sigma_\text{true}$ and $\phi_\text{true}$ are large, while interweaving strategies perform well for all underlying parameter values. This result partially contrasts the conclusions of str-etal:par, who associate better mixing with larger $|\phi|$ for all parameterizations and recommend the non-centered parameterization in any setup. It should be noted, however, that these authors use a different sampling algorithm that does not rely on Gaussian mixture approximation. Moreover, the parameter range investigated in their paper does not span the range of parameters examined in our paper. For a summary of percentage gains in terms of effective sample size, see Table (ref).
We apply our estimation methodology to daily Euro exchange rates. The data stems from the European Central Bank's Statistical Data Warehouse and comprises 3140 observations of 23 currencies ranging from January 3, 2000 to April 4, 2012. In choosing the prior for $\phi$ we follow kim-etal:sto, i.e. $(\phi+1)/2 \sim \mathcal{B}\left(20, 1.5\right)$, and for the other parameters we pick rather vague priors: $\mu \sim \mathcal{N}\!\left(-10, 100\right)$ and $\sigma^2 \sim \mathcal{G}\left(\frac{1}{2},\frac{1}{2}\right)$. After a burn-in of $10\ 000$, we use $1\ 000\ 000$ draws from the respective distributions in each parameterization for posterior inference.
To exemplify, Figure (ref) shows exchange rates of EUR/US\$ along with absolute de-meaned log-returns, which are then used to estimate the time-varying volatilities displayed below. The transformed latent process and the absolute log-returns exhibit a similar overall pattern. Nevertheless, the volatility path is much smoother, which is due to the highly persistent autoregressive process (the posterior mean of $\phi$ is $0.993$, the posterior mean of $\sigma$ is $0.07$). Marginal posterior density estimates and two-way scatterplots can be found in Figure (ref). Note that only $p(\mu|\mathbf{y})$ is symmetric, while both $p(\phi|\mathbf{y})$ and $p(\sigma|\mathbf{y})$ are skewed. Moreover, the parameter draws are (sometimes nonlinearly) correlated.
Results for all 23 examined exchange rates are displayed in Table (ref). It stands out that for currencies which are closely tied to the Euro, posterior parameter means differ substantially to those found above. Most notably, the Danish krone exhibits very low overall level of volatility ($\mu_\text{mean} = -18$), paired with moderate persistence ($\phi_\text{mean}=0.916$) and moderately high volatility of volatility ($\sigma_\text{mean}=0.38$). Looking at the inefficiency factors for the raw parameterizations, one observes striking superiority of C in terms of sampling efficiency of $\mu$, while NC usually performs better in terms of sampling efficiency of $\sigma$. Again, interweaving overcomes these problems by showing lowest IFs uniformly for all parameters and all time series. Even though not reported here in detail due to space constraints, the choice of the baseline (GIS-C vs. GIS-NC) is negligible.
Previous studies have shown that simple reparameterizations often turn out to have substantial impact on MCMC simulation efficiency in state-space specifications. This paper contributes to the literature by exploring the influence of choosing between two selected parameterizations for Bayesian estimation of SV models. Moreover, it provides evidence that inefficiency factors obtained from simulation experiments can heavily depend on the realization of the data generating process. Through the findings of this paper it becomes clear that employing an ancillarity-sufficiency interweaving strategy (ASIS) introduced by yu-men:cen helps to overcome shortcomings of either the centered or the non-centered parameterization by outperforming those in terms of sampling efficiency with respect to all parameters at very little extra computational cost, whereas the baseline of the interweaving strategy is of minor influence.
The concept of interweaving different parameterizations of state-space models is clearly very general, and there is good reason to hope for similar magic when applying ASIS to extension of the basic SV model such as more general innovation distributions lie-jun:sto, del-gri:bay, asymmetry yu:on, omo-etal:sto or both chi-etal:mar, wan-etal:sto, tsi-on, ish-omo:eff, nak-omo:sto. Preliminary results for an SV model with leverage, where a centered parameterization from yu:on is compared with a non-centered version based on transforming $h_t$ into $\tilde h_t = (h_t-\mu)/\sigma$ as in the present paper, show that this hope is in fact an actual possibility. A thorough investigation of this issue is however beyond the scope of this article.
The authors would like to thank the editor and two referees for their perspicacious comments on an earlier draft of this paper, and Stefan Theu\ss l for helpful advice concerning coding and implementation.