Ancillarity-Sufficiency Interweaving Strategy (ASIS) for Boosting MCMC Estimation of Stochastic Volatility Models

TL;DR

This paper introduces an interweaving strategy that significantly enhances MCMC sampling efficiency for stochastic volatility models, enabling more robust Bayesian inference across diverse parameter settings.

Contribution

It applies the ancillarity-sufficiency interweaving strategy to stochastic volatility models, improving sampling efficiency and allowing inference in previously infeasible parameter regimes.

Findings

Enhanced MCMC sampling efficiency across all parameters.

Effective inference without pre-selecting parameterizations.

Overcomes limitations of centered and non-centered models.

Abstract

Bayesian inference for stochastic volatility models using MCMC methods highly depends on actual parameter values in terms of sampling efficiency. While draws from the posterior utilizing the standard centered parameterization break down when the volatility of volatility parameter in the latent state equation is small, non-centered versions of the model show deficiencies for highly persistent latent variable series. The novel approach of ancillarity-sufficiency interweaving has recently been shown to aid in overcoming these issues for a broad class of multilevel models. In this paper, we demonstrate how such an interweaving strategy can be applied to stochastic volatility models in order to greatly improve sampling efficiency for all parameters and throughout the entire parameter range. Moreover, this method of "combining best of different worlds" allows for inference for parameter constellations that have previously been infeasible to estimate without the need to select a particular parameterization beforehand.

Figures

Tables

Table 1. Table 1: Inefficiency factors for 100 000 100000 100\,000 draws from p ​ ( μ | 𝐲 ) 𝑝 conditional 𝜇 𝐲 p(\mu|\mathbf{y}) in various parameterizations using different blocking strategies. Time series length T = 5000 𝑇 5000 T=5000 , the values reported are medians of 500 500 500 repetitions and T CPU subscript 𝑇 CPU T_{\text{CPU}} denotes the median time to complete 1000 1000 1000 iterations. Shading: 0 50 100 150 200 250 300 350 400 450 500+

0

50

100

150

200

250

300

350

400

450

500+

Table 2. Table 2: Percentage gains in effective sample size for the 2-block and the 3-block sampler. First and third table: ESS GIS-C subscript ESS GIS-C \text{ESS}_{\text{GIS-C}} vs.  max ⁡ ( ESS C , ESS NC ) subscript ESS C subscript ESS NC \max(\text{ESS}_{\text{C}},\text{ESS}_{\text{NC}}) . Second and fourth table: ESS GIS-C subscript ESS GIS-C \text{ESS}_{\text{GIS-C}} vs.  min ⁡ ( ESS C , ESS NC ) subscript ESS C subscript ESS NC \min(\text{ESS}_{\text{C}},\text{ESS}_{\text{NC}}) . Shading: 0 10 20 30 40 50 60 70 80 90 100+ (1 and 3) Shading: 0 100 200 300 400 500 600 700 800 900 1000+ (2 and 4)

0

10

20

30

40

50

60

70

80

90

100+

(1 and 3)

Table 3. Table 3: Inefficiency factors for 100 000 100000 100\,000 draws from p ​ ( ϕ | 𝐲 ) 𝑝 conditional italic-ϕ 𝐲 p(\phi|\mathbf{y}) in various parameterizations using different blocking strategies. Time series length T = 5000 𝑇 5000 T=5000 , the values reported are medians of 500 500 500 repetitions and T CPU subscript 𝑇 CPU T_{\text{CPU}} denotes the median time to complete 1000 1000 1000 iterations. Shading: 0 50 100 150 200 250 300 350 400 450 500+

0

50

100

150

200

250

300

350

400

450

500+

Table 4. Table 4: Percentage gains in effective sample size for the 2-block and the 3-block sampler. First and third table: ESS GIS-C subscript ESS GIS-C \text{ESS}_{\text{GIS-C}} vs.  max ⁡ ( ESS C , ESS NC ) subscript ESS C subscript ESS NC \max(\text{ESS}_{\text{C}},\text{ESS}_{\text{NC}}) . Second and fourth table: ESS GIS-C subscript ESS GIS-C \text{ESS}_{\text{GIS-C}} vs.  min ⁡ ( ESS C , ESS NC ) subscript ESS C subscript ESS NC \min(\text{ESS}_{\text{C}},\text{ESS}_{\text{NC}}) . Shading: 0 10 20 30 40 50 60 70 80 90 100+ (1 and 3) Shading: 0 100 200 300 400 500 600 700 800 900 1000+ (2 and 4)

0

10

20

30

40

50

60

70

80

90

100+

(1 and 3)

Table 5. Table 5: Inefficiency factors for 100 000 100000 100\,000 draws from p ​ ( σ | 𝐲 ) 𝑝 conditional 𝜎 𝐲 p(\sigma|\mathbf{y}) in various parameterizations using different blocking strategies. Time series length T = 5000 𝑇 5000 T=5000 , the values reported are medians of 500 500 500 repetitions and T CPU subscript 𝑇 CPU T_{\text{CPU}} denotes the median time to complete 1000 1000 1000 iterations. Shading: 0 50 100 150 200 250 300 350 400 450 500+

0

50

100

150

200

250

300

350

400

450

500+

Table 6. Table 6: Percentage gains in effective sample size for the 2-block and the 3-block sampler. First and third table: ESS GIS-C subscript ESS GIS-C \text{ESS}_{\text{GIS-C}} vs.  max ⁡ ( ESS C , ESS NC ) subscript ESS C subscript ESS NC \max(\text{ESS}_{\text{C}},\text{ESS}_{\text{NC}}) . Second and fourth table: ESS GIS-C subscript ESS GIS-C \text{ESS}_{\text{GIS-C}} vs.  min ⁡ ( ESS C , ESS NC ) subscript ESS C subscript ESS NC \min(\text{ESS}_{\text{C}},\text{ESS}_{\text{NC}}) . Shading: 0 10 20 30 40 50 60 70 80 90 100+ (1 and 3) Shading: 0 100 200 300 400 500 600 700 800 900 1000+ (2 and 4)

0

10

20

30

40

50

60

70

80

90

100+

(1 and 3)

Table 7. Table 7: Posterior means and inefficiency factors for various estimation methods of the SV model, applied to EUR exchange rate data.

	Posterior means			${\text{IF}}_{\text{C}}$			${\text{IF}}_{\text{NC}}$			${\text{IF}}_{\text{GIS-C}}$
	$\mu$	$\varphi$	$\sigma$	$\mu$	$\varphi$	$\sigma$	$\mu$	$\varphi$	$\sigma$	$\mu$	$\varphi$	$\sigma$
Australian dollar	-10.3	0.976	0.17	3	167	256	216	120	149	2	68	97
Canadian dollar	-10.1	0.987	0.09	3	350	539	211	128	159	3	89	120
Swiss franc	-12.0	0.985	0.21	3	55	129	807	99	162	3	33	73
Czech koruna	-11.5	0.953	0.28	3	140	200	151	121	154	3	72	96
Danish krone	-18.0	0.916	0.38	5	102	145	85	90	115	4	57	72
UK pound sterling	-10.8	0.992	0.10	2	95	233	802	120	150	2	39	87
Hong Kong dollar	-10.2	0.993	0.07	2	128	309	526	89	96	2	36	75
Indonesian rupiah	-9.9	0.966	0.23	4	234	313	216	201	242	3	114	142
Japanese yen	-10.0	0.989	0.12	3	110	227	513	98	140	3	47	91
Korean won	-10.0	0.987	0.14	2	92	192	501	96	134	2	40	79
Mexican peso	-9.8	0.977	0.16	2	137	221	220	101	128	2	60	86
Malaysian ringgit	-10.3	0.990	0.08	2	220	390	404	104	124	2	59	90
Norwegian krone	-11.1	0.970	0.18	2	139	217	150	83	107	2	53	76
New Zealand dollar	-10.0	0.963	0.17	4	344	432	105	153	177	3	114	135
Philippine peso	-10.1	0.981	0.12	3	376	529	209	189	220	2	123	160
Polish zloty	-10.4	0.975	0.19	2	96	171	261	85	117	2	43	69
Romanian leu	-11.1	0.970	0.31	2	51	100	446	87	137	2	32	60
Russian rouble	-10.6	0.988	0.15	4	75	172	891	121	151	3	38	82
Swedish krona	-11.3	0.991	0.11	1	52	156	752	75	99	1	23	60
Singapore dollar	-10.6	0.995	0.07	4	138	348	998	126	132	4	47	100
Thai bhat	-10.2	0.980	0.13	3	202	314	207	102	125	3	64	90
Turkish lira	-9.8	0.966	0.27	2	72	127	259	84	127	2	42	69
US dollar	-10.1	0.993	0.07	2	126	308	504	87	99	2	37	74

Equations

$$y_t$$

$$h_t$$

$$y_t\sim \mathcal{N}\left(0,\omega e^{\sigma {\tilde{h}}_t}\right),$$

$${\tilde{h}}_t=\varphi {\tilde{h}}_{t-1}+{\eta }_t,{\eta }_t\sim \mathcal{N}\left(0,1\right),$$

$$p(\varphi )=\frac{1}{2B(a_0,b_0)}{\left(\frac{1+\varphi }{2}\right)}^{a_0-1}{\left(\frac{1-\varphi }{2}\right)}^{b_0-1}.$$

$${\tilde{y}}_t$$

$$\mathbf{\Omega }=\left[\begin{array}{ccccc}\frac{1}{s_{r_1}^2}+\frac{1}{{\sigma }^2}& \frac{-\varphi }{{\sigma }^2}& 0& \dots & 0\\ \frac{-\varphi }{{\sigma }^2}& \frac{1}{s_{r_2}^2}+\frac{1+{\varphi }^2}{{\sigma }^2}& \frac{-\varphi }{{\sigma }^2}& \ddots & ⋮\\ 0& \frac{-\varphi }{{\sigma }^2}& \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \frac{1}{s_{r_{T-1}}^2}+\frac{1+{\varphi }^2}{{\sigma }^2}& \frac{-\varphi }{{\sigma }^2}\\ 0& \dots & 0& \frac{-\varphi }{{\sigma }^2}& \frac{1}{s_{r_T}^2}+\frac{1}{{\sigma }^2}\end{array}\right],$$

$$\mathbf{c}=\left[\begin{array}{c}\frac{1}{s_{r_1}^2}({\tilde{y}}_1-m_{r_1})+\frac{\mu (1-\varphi )}{{\sigma }^2}\\ \frac{1}{s_{r_2}^2}({\tilde{y}}_2-m_{r_2})+\frac{\mu (1-\varphi {)}^2}{{\sigma }^2}\\ ⋮\\ \frac{1}{s_{r_{T-1}}^2}({\tilde{y}}_{T-1}-m_{r_{T-1}})+\frac{\mu (1-\varphi {)}^2}{{\sigma }^2}\\ \frac{1}{s_{r_T}^2}({\tilde{y}}_T-m_{r_T})+\frac{\mu (1-\varphi )}{{\sigma }^2}\end{array}\right].$$

$$\mathbf{\Omega }=\left[\begin{array}{ccccc}\frac{{\sigma }^2}{s_{r_1}^2}+1& -\varphi & 0& \dots & 0\\ -\varphi & \frac{{\sigma }^2}{s_{r_2}^2}+1+{\varphi }^2& -\varphi & \ddots & ⋮\\ 0& -\varphi & \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \frac{{\sigma }^2}{s_{r_{T-1}}^2}+1+{\varphi }^2& -\varphi \\ 0& \dots & 0& -\varphi & \frac{{\sigma }^2}{s_{r_T}^2}+1\end{array}\right],$$

$$\mathbf{c}=\left[\begin{array}{c}\frac{\sigma }{s_{r_1}^2}({\tilde{y}}_1-m_{r_1}-\mu )\\ \frac{\sigma }{s_{r_2}^2}({\tilde{y}}_2-m_{r_2}-\mu )\\ ⋮\\ \frac{\sigma }{s_{r_{T-1}}^2}({\tilde{y}}_{T-1}-m_{r_{T-1}}-\mu )\\ \frac{\sigma }{s_{r_T}^2}({\tilde{y}}_T-m_{r_T}-\mu )\end{array}\right].$$

$$h_t=\gamma +\varphi h_{t-1}+{\eta }_t,{\eta }_t\sim \mathcal{N}\left(0,{\sigma }^2\right),$$

$$p_{\text{aux}}({\mathbf{\theta }}_{\text{new}}|\mathbf{h})=p_{\text{aux}}({\gamma }_{\text{new}},{\varphi }_{\text{new}}|\mathbf{h},{\sigma }_{\text{new}}^2)p_{\text{aux}}({\sigma }_{\text{new}}^2|\mathbf{h}).$$

$$\gamma ,\varphi |\mathbf{h},{\sigma }^2\sim N_2\left({\mathbf{b}}_T,{\sigma }^2{\mathbf{B}}_T\right),$$

$$R=\frac{p(h_0|{\mathbf{\theta }}_{\text{new}})p({\gamma }_{\text{new}}|{\varphi }_{\text{new}})p({\varphi }_{\text{new}})p({\sigma }_{\text{new}}^2)}{p(h_0|{\mathbf{\theta }}_{\text{old}})p({\gamma }_{\text{old}}|{\varphi }_{\text{old}})p({\varphi }_{\text{old}})p({\sigma }_{\text{old}}^2)}\times \frac{p_{\text{aux}}({\varphi }_{\text{old}},{\gamma }_{\text{old}}|{\sigma }_{\text{old}}^2)p_{\text{aux}}({\sigma }_{\text{old}}^2)}{p_{\text{aux}}({\varphi }_{\text{new}},{\gamma }_{\text{new}}|{\sigma }_{\text{new}}^2)p_{\text{aux}}({\sigma }_{\text{new}}^2)}.$$

$$R=\frac{p(h_0|{\gamma }_{\text{new}},{\varphi }_{\text{new}})p({\gamma }_{\text{new}}|{\varphi }_{\text{new}})p({\varphi }_{\text{new}})}{p(h_0|{\gamma }_{\text{old}},{\varphi }_{\text{old}})p({\gamma }_{\text{old}}|{\varphi }_{\text{old}})p({\varphi }_{\text{old}})}\times \frac{p_{\text{aux}}({\gamma }_{\text{old}},{\varphi }_{\text{old}})}{p_{\text{aux}}({\gamma }_{\text{new}},{\varphi }_{\text{new}})}.$$

$${\sigma }^2|\mathbf{h},\mu ,\varphi \sim {\mathcal{G}}^{-1}\left(c_T,C_T\right),$$

$$R=\frac{p({\sigma }_{\text{new}}^2)}{p({\sigma }_{\text{old}}^2)}\times \frac{p_{\text{aux}}({\sigma }_{\text{old}}^2)}{p_{\text{aux}}({\sigma }_{\text{new}}^2)}=\exp \left\{\frac{{\sigma }_{\text{old}}^2-{\sigma }_{\text{new}}^2}{2B_{\sigma }}\right\}.$$

$$\varphi |\mathbf{h},\gamma ,{\sigma }^2\sim \mathcal{N}\left(\frac{\left[{\sum }_{t=1}^T{h}_{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).$$

$$R=\frac{p(h_0|{\varphi }_{\text{new}},\mu ,{\sigma }^2)p({\varphi }_{\text{new}})}{p(h_0|{\varphi }_{\text{old}},\mu ,{\sigma }^2)p({\varphi }_{\text{old}})}\times \frac{p_{\text{aux}}({\varphi }_{\text{old}}|{\sigma }^2)}{p_{\text{aux}}({\varphi }_{\text{new}}|{\sigma }^2)}.$$

$$\gamma |\mathbf{h},\varphi ,{\sigma }^2\sim \mathcal{N}\left(\frac{{\sum }_{t=1}^T{h}_t-\varphi {\sum }_{t=0}^{T-1}{h}_t}{T+1/B_0^{11}},\frac{{\sigma }^2}{T+1/B_0^{11}}\right)$$

$$R=\frac{p(h_0|{\gamma }_{\text{new}},\varphi ,{\sigma }^2)p({\gamma }_{\text{new}}|\varphi )}{p(h_0|{\gamma }_{\text{old}},\varphi ,{\sigma }^2)p({\gamma }_{\text{old}}|\varphi )}\times \frac{p_{\text{aux}}({\gamma }_{\text{old}}|{\sigma }^2)}{p_{\text{aux}}({\gamma }_{\text{new}}|{\sigma }^2)}.$$

$$\varphi |\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),$$

$$\breve{\mathbf{y}}=\mathbf{X}\left[\begin{array}{c}\mu \\ \sigma \end{array}\right]+\mathbf{ϵ},$$

$$\breve{\mathbf{y}}=\left[\begin{array}{c}({\tilde{y}}_1-m_{r_1})/s_{r_1}\\ ⋮\\ ({\tilde{y}}_T-m_{r_T})/s_{r_T}\end{array}\right],\mathbf{X}=\left[\begin{array}{cc}{\tilde{h}}_1/s_{r_1}& 1/s_{r_1}\\ ⋮& ⋮\\ {\tilde{h}}_T/s_{r_T}& 1/s_{r_T}\end{array}\right].$$

$$b_{T,\mu }=B_{T,\mu }\left(\sum _{t=1}^T\frac{{\tilde{y}}_t-m_{r_t}-\sigma {\tilde{h}}_t}{s_{r_t}^2}+\frac{b_{\mu }}{B_{\mu }}\right),B_{T,\mu }=1/\left(\sum _{t=1}^{T}1/s_{r_t}^2+\frac{1}{B_{\mu }}\right),$$

$$b_{T,\sigma }=B_{T,\sigma }\sum _{t=1}^T\frac{{\tilde{h}}_t({\tilde{y}}_t-m_{r_t}-\mu )}{s_{r_t}^2},B_{T,\sigma }=1/\left(\sum _{t=1}^T\frac{{\tilde{h}}_t^2}{s_{r_t}^2}+\frac{1}{B_{\sigma }}\right).$$

$$\mathbb{P}(r_t=k|\cdot )\propto \mathbb{P}(r_t=k)\frac{1}{s_k}\exp \left\{-\frac{({ϵ}_t^{\ast }-m_k{)}^2}{2s_k^2}\right\},$$

Full text

Keywords: Markov Chain Monte Carlo, Non-Centering, Auxiliary Mixture Sampling, Massively Parallel Computing, State Space Model, Exchange Rate Data

1 Introduction

Returns of financial and economic time series often exhibit time-varying volatilities. To account for this behavior, Taylor (1982) 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 (see e.g. Hull & White, 1987).

Following e.g. Jacquier et al. (1994) or Kim et al. (1998), observed log-returns are denoted $\mathbf{y}=(y_{1},y_{2},\dots,y_{T})^{\prime}$ and the SV model is specified as

[TABLE]

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 (2) 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 (1) and (2) 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 Gelfand et al. (1995), Pitt & Shephard (1999), Frühwirth-Schnatter (2004), Roberts et al. (2004), Frühwirth-Schnatter & Sögner (2008), and Strickland et al. (2008). The pioneering paper by Taylor (1982) as well as several other papers like Kim et al. (1998) or Liesenfeld & Richard (2006) 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 (2) to the observation equation (1) by setting $\bar{h}_{t}=h_{t}-\mu$. Kim et al. (1998) 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 et al. (1998). 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

[TABLE]

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\}$, Strickland et al. (2008) 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. Frühwirth-Schnatter (2004) and Frühwirth-Schnatter & Wagner (2010). 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. Pitt & Shephard (1999) 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, Jacquier et al. (1994) 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, Shephard & Kim (1994) note that the draws obtained from this sampler are also highly correlated and thus only slowly converging to the stationary distribution. Alternatively, Shephard & Pitt (1997) propose a multi-move sampler, where volatility blocks of random length are updated at a time, while Shephard (1994) and Omori et al. (2007) 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 (Carter & Kohn, 1994; Frühwirth-Schnatter, 1994; Durbin & Koopman, 2002). Within a more general Gaussian state-space framework, Rue (2001) and McCausland et al. (2011) 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 (2012).

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 & Meng (2011). 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 2 gives insight into the estimation procedure for each of the two selected parameterizations. Section 3 explains how ASIS can be applied in order to interweave these parameterizations. Extensive simulation results presented in Section 4 compare sampling efficiency for all parameters amongst the different parameterizations, Section 5 provides real-data results for several daily exchange rates, and Section 6 concludes.

2 Bayesian Inference in the SV Model

2.1 Prior Distributions

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 et al. (1998), implying

[TABLE]

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 Frühwirth-Schnatter & Wagner (2010), 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.

2.2 MCMC Methodology

Observation equation (1) can easily be rewritten as

[TABLE]

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 et al. (1998) or $\log(y_{t}^{2}+c)$ with $c=10^{-4}$ as in Omori et al. (2007) in order to avoid values equal to zero. Equation (6) 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 Omori et al. (2007). This representation allows rewriting (6) as a linear and conditionally Gaussian state space model,

[TABLE]

where speedy MCMC sampling becomes possible in three steps.

Algorithm 1** (AWOL Sampler).**

Choose appropriate starting values for the parameters $\mu$, $\phi$, $\sigma$ and the indicators $\mathbf{r}=(r_{1},r_{2},\dots,r_{T})^{\prime}$ – e.g. start with components with high weights – and repeat the following steps:

1. (a)

Sample the latent volatilities AWOL by drawing from ${\mathbf{h}}_{[-0]}|{\mathbf{y}},\mathbf{r},\mu,\phi,\sigma^{2}$ or $\tilde{\mathbf{h}}_{[-0]}|{\mathbf{y}},\mathbf{r},\mu,\phi,\sigma^{2}$, respectively. The initial value is drawn from $h_{0}|h_{1},\mu,\phi,\sigma^{2}$ or from $\tilde{h}_{0}|\tilde{h}_{1},\phi$. 2. (b)

Sample $\mu,\phi,\sigma^{2}$ via Bayesian regression.

For C, we investigate a 1-block sampler, drawing from $\mu,\phi,\sigma^{2}|{\mathbf{h}}$, a 2-block sampler, where $\sigma^{2}$ is drawn from $\sigma^{2}|{\mathbf{h}},\mu,\phi$, while $\mu$ and $\phi$ are sampled jointly from $\mu,\phi|{\mathbf{h}},\sigma^{2}$, and a 3-block sampler, where all parameters are individually drawn from the full conditionals. Due to non-conjugacy of the chosen priors, MH updates are used in all variants.

• -

In NC, MH is needed only for updating $\phi$ by drawing from $\phi|\tilde{\mathbf{h}}$, while $\mu$ and $\sigma^{2}$ can be Gibbs-updated jointly from $\mu,\sigma^{2}|\mathbf{y},\tilde{\mathbf{h}},\mathbf{r}$ (2-block) or individually from $\mu|\mathbf{y},\tilde{\mathbf{h}},\mathbf{r},\sigma^{2}$ and $\sigma^{2}|\mathbf{y},\tilde{\mathbf{h}},\mathbf{r},\mu$ (3-block). 3. (c)

Update the indicators $\mathbf{r}$ from $\mathbf{r}|\mathbf{y},{\mathbf{h}}$ in C, or $\mathbf{r}|\mathbf{y},\tilde{\mathbf{h}},\mu,\sigma^{2}$ in NC, via inverse transform sampling.

2.3 Step (a): Sampling the Latent Volatilities AWOL

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 (2001) and McCausland et al. (2011) 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

[TABLE]

and

[TABLE]

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

[TABLE]

and

[TABLE]

For both parameterizations, this is accomplished by first computing the Cholesky decomposition $\mathbf{\Omega}=\mathbf{LL}^{\prime}$. Due to the band structure of $\mathbf{\Omega}$, this is computationally inexpensive and can either be implemented directly or via the LAPACK-routine dpbtrf (Anderson et al., 1999), 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}^{\prime}{\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.

2.4 Step (b)-C: Sampling of $\mu$, $\phi$ and $\sigma$ in C

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,

[TABLE]

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:

[TABLE]

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^{2}B_{0}^{11}\right)$ and $p_{\text{aux}}(\phi|\sigma)\sim\mathcal{N}\!\left(0,\sigma^{2}B_{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

[TABLE]

with $\mathbf{B}_{T}=(\mathbf{X^{\prime}X}+\mathbf{B}_{0}^{-1})^{-1}$ and $\mathbf{b}_{T}=\mathbf{B}_{T}\mathbf{X}^{\prime}\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}^{\prime}\mathbf{X}^{\prime}{\mathbf{h}}_{[-0]}\right)$. The acceptance probability is given through $\min(1,R)$, with

[TABLE]

In the two-block sampler, we draw the first block from the full conditional distribution $\gamma,\phi|{\mathbf{h}},\sigma^{2}$ given in (8) and accept with probability $\min(1,R)$, where

[TABLE]

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

[TABLE]

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

[TABLE]

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 (9). For sampling $\phi$, we obtain a proposal from

[TABLE]

The acceptance probability is equal to $\min(1,R)$ with

[TABLE]

For sampling $\gamma$ from the full conditional posterior distribution, we obtain a proposal from

[TABLE]

and an acceptance probability equaling $\min(1,R)$ with

[TABLE]

2.5 Step (b)-NC: Sampling of $\mu$, $\phi$ and $\sigma$ in NC

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

[TABLE]

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 (7) as a regression model with homoskedastic errors, i.e.

[TABLE]

where $\boldsymbol{\epsilon}\sim N_{K}\left(\mathbf{0},\mathbf{I}_{K}\right)$, and

[TABLE]

The joint posterior distribution is again bivariate Gaussian with variance-covariance matrix $\mathbf{B}_{T}=(\mathbf{B}_{0}^{-1}+\mathbf{X}^{\prime}\mathbf{X})^{-1}$ and mean $\mathbf{b}_{T}=\mathbf{B}_{T}(\mathbf{B}_{0}^{-1}\mathbf{b}_{0}+\mathbf{X}^{\prime}\mathbf{\breve{y}})$, where $\mathbf{b}_{0}=(b_{\mu},0)^{\prime}$ 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

[TABLE]

and $\sigma|{\mathbf{y}},\tilde{\mathbf{h}},\mathbf{r},\mu\sim\mathcal{N}\!\left(b_{T,\sigma},B_{T,\sigma}\right)$ with

[TABLE]

2.6 Step (c): Sampling the Indicators $\mathbf{r}$

We proceed exactly as Omori et al. (2007). 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

[TABLE]

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\times 10$ exponential function calls, this step is computationally rather expensive but can easily be parallelized.

3 Interweaving C and NC by ASIS

To provide some intuition about the sampling efficiency in C, let $\phi=0$ for a moment. This implies that the state equation (2) 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 (4) 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 (10) 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 & Meng (2011) 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 (Basu, 1955) 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 2 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 & Meng (2011) and stand for global interweaving strategy. C simply denotes the fact that we use the centered baseline.

Algorithm 2** (GIS-C).**

Choose appropriate starting values and repeat the following steps:

1. (a)

Draw $\mathbf{h}$ (C). 2. (b)

Draw $\mu,\phi,\sigma$ (C). 3. (b*)

Move to NC by the simple deterministic transformation $\tilde{h}_{t}=\frac{h_{t}-\mu}{\sigma}$ for all $t.$ 4. (b**)

Redraw $\mu,\phi,\sigma$ (NC). 5. (b***)

Move back to C by calculating $h_{t}=\mu+\sigma\tilde{h}_{t}$ for all $t$. 6. (c)

Draw the indicators $\mathbf{r}$ (C).

The individual sampling steps are implemented exactly as described in subsections 2.3 to 2.6. 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.

Algorithm 3** (GIS-NC).**

Choose appropriate starting values and repeat the following steps:

1. (a)

Draw $\mathbf{\tilde{h}}$ (NC). 2. (b)

Draw $\mu,\phi,\sigma$ (NC). 3. (b*)

Move to C by the simple deterministic transformation $h_{t}=\mu+\sigma\tilde{h}_{t}$ for all $t$. 4. (b**)

Redraw $\mu,\phi,\sigma$ (C). 5. (b***)

Move back to NC by transforming back: $\tilde{h}_{t}=\frac{h_{t}-\mu}{\sigma}$ for all $t$. 6. (c)

Draw the indicators $\mathbf{r}$ (NC).

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.

4 Simulation Results

In order to assess simulation efficiency of our algorithms, we simulate data from the model specified in equations (1) and (2). 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. Jacquier et al. (1994), Kim et al. (1998), Liesenfeld & Richard (2006), and Strickland et al. (2008). 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 Core Team, 2013) and OpenMPI 1.4.3 (Gabriel et al., 2004). For high-level-parallelization and parallel random number generation according to L’Ecuyer et al. (2002), the R packages parallel (part of R) and snow (Tierney et al., 2011) 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 stochvol (Kastner, 2013), available on CRAN, was created. The core implementation is written in C, interfaced to R via Rcpp (Eddelbuettel & François, 2011). Inefficiency factors and effective sample sizes were computed with the R package coda (Plummer et al., 2006).

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.

4.1 To Center Or Not to Center?

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 1 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 Pitt & Shephard (1999) and Frühwirth-Schnatter (2004), 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 2 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.

4.2 Sampling Efficiency

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 MCSE2 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 3 and Figure 4 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.

4.3 Efficiency Overview

In order to gain insight into the entire parameter range of interest, Tables 1, 3 and 5 provide a summary of median inefficiency factors across all $45$ parameter constellations.

Median IFs obtained from draws from $p(\mu|\mathbf{y})$ in Table 1 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 1 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 2 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 3. 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 4. 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 5 summarizes median IFs for draws from $p(\sigma|\mathbf{y})$. We observe a similar overall picture to the one presented in Table 1: 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 Strickland et al. (2008), 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 6.

5 Application to Exchange Rate Data

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 et al. (1998), 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 5 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 [6](#S5.F6). 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 7. 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.

6 Concluding Remarks

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 & Meng (2011) 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 (e.g. Liesenfeld & Jung, 2000; Delatola & Griffin, 2011), asymmetry (e.g. Yu, 2005; Omori et al., 2007) or both (e.g. Chib et al., 2002; Wang et al., 2011; Tsiotas, 2012; Ishihara & Omori, 2012; Nakajima & Omori, 2012). Preliminary results for an SV model with leverage, where a centered parameterization from Yu (2005) 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.

7 Acknowledgments

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ßl for helpful advice concerning coding and implementation.