The Stochastic Stationary Root Model

Abstract

We propose and study the stochastic stationary root model. The model resembles the cointegrated VAR model but is novel in that: (i) the stationary relations follow a random coefficient autoregressive process, i.e., exhibhits heavy-tailed dynamics, and (ii) the system is observed with measurement error. Unlike the cointegrated VAR model, estimation and inference for the SSR model is complicated by a lack of closed-form expressions for the likelihood function and its derivatives. To overcome this, we introduce particle filter-based approximations of the log-likelihood function, sample score, and observed Information matrix. These enable us to approximate the ML estimator via stochastic approximation and to conduct inference via the approximated observed Information matrix. We conjecture the asymptotic properties of the ML estimator and conduct a simulation study to investigate the validity of the conjecture. Model diagnostics to assess model fit are considered. Finally, we present an empirical application to the 10-year government bond rates in Germany and Greece during the period from January 1999 to February 2018.

1. Introduction

In this paper, we introduce the multivariate stochastic stationary root (SSR) model. The SSR model is a nonlinear state space model, which resembles the Granger-Johansen representation of the cointegrated vector autoregressive (CVAR) model, see inter alia Johansen (1996) and Juselius (2007). The SSR model decomposes a p–dimensional observation vector into r stationary components and $p-r$ nonstationary components, which is similar to the CVAR model. However, the roots of the stationary components are allowed to be stochastic; hence the name ‘stochastic stationary root’. The stationary and nonstationary dynamics of the model are observed with measurement error, which in this model prohibits close-form expressions for e.g., the log-likelihood, sample score and observed Information matrix. Likelihood-based estimation and inference therefore calls for non-standard methods.

Although the SSR model resembles the CVAR model, it is differentiated by its ability to characterize heavy-tailed dynamics in the stationary component. Heavy-tailed dynamics, and other types of nonlinear dependencies, are not amenable to analysis with the CVAR model, which has prompted work into nonlinear alternatives, see inter alia Bohn Nielsen and Rahbek (2014), Kristensen and Rahbek (2013), Kristensen and Rahbek (2010), and Bec et al. (2008). Similarly, cointegration in the state space setting has been considered in term of the common stochastic trend (CST) model by Chang et al. (2009) as well as the CVAR model with measurement errors by Bohn Nielsen (2016). Additionally, the SSR model is also related to the stochastic unit root literature, see inter alia Granger and Swanson (1997), Leybourne and McCabe (1996), Lieberman and Phillips (2014), Lieberman and Phillips (2017), McCabe and Tremayne (1995), and McCabe and Smith (1998). Relevant empirical applications where the SSR model could potentially provide a better fit than the CVAR model include, but are not limited to, (i) log-prices of assets that exhibit random walk behavior in the levels and heavy-tailed error-correcting dynamics in the no-arbitrage relations, and (ii) interest rates for which the riskless rate exhibits random walk-type dynamics and the risk premia undergo periods of high levels and high volatility.

The stationary and nonstationary components of the SSR model are treated as unobserved processes, and consequently need to be integrated out in order to compute the log-likelihood function and its derivatives. Due to the nonlinearity of the model, this cannot be accomplished analytically. We appeal to the incomplete data framework and the simulation-based approach known as particle filtering to approximate the log-likelihood function, sample score and observed Information matrix. See inter alia Gordon et al. (1993), Doucet et al. (2001), Cappé et al. (2005), and Creal (2012) for an overview of the particle filtering literature. Moreover, we rely on stochastic approximation methods to obtain the maximum likelihood (ML) estimator, see Poyiadjis et al. (2011). Summarizing, the main contributions of this paper are to

i

introduce and study the SSR model, and

ii

propose a method for approximate frequentist estimation and inference.

It is beyond the scope of this paper to provide a complete proof of the asymptotic properties of the ML estimator. The study of the asymptotic properties of the ML estimator in general state space models, such as the SSR model, is an emerging area of research. Most existing results rely on compactness of the state space, which excludes the SSR model and is generally restrictive. For results in this direction, see e.g., Olsson and Rydén (2008) who derive consistency and asymptotic normality for the ML estimator by discretizing the parameter space. Douc et al. (2011) have shown consistency of the ML estimator without assuming compactness, but the regularity conditions are nonetheless too restrictive to encompass the SSR model. Instead of providing a complete proof of the asymptotic properties of the ML estimator, we conjecture the asymptotic properties of the derivatives of the log-likelihood function. We base the conjecture on known properties of models that are closely related to the SSR model, and corroborate it by a simulation study. Given the conjecture holds, it allows us to establish the asymptotic properties of the ML estimator. We leave proving the conjecture for future work, and focus in this paper on developing methods for approximate frequentist estimation and inference.

The rest of the paper is organized as follows. We introduce the SSR model in Section 2, and study some properties of the process in Section 3. In Section 4 we introduce likelihood-based estimation and inference for the unknown model parameter. In Section 5 we introduce the incomplete data framework. In Section 6 we introduce the particle filter-based approximations to the log-likelihood function, sample score and Information matrix. In Section 7 we propose how to approximate the ML estimator and classic standard errors. In Section 8 we consider model diagnostics. In Section 9 we conduct a simulation study of the asymptotic distribution of the ML estimator. In Section 10 we apply the SSR model to monthly observations of 10-year government bond rates in Germany and Greece from January 1999 to February 2018. We conclude in Section 11. All proofs have been relegated to Appendix B, while Appendix A contains various auxiliary results.

Notation-wise, we adopt the convention that the ‘blackboard bold’ typeface, e.g., $\mathbb{E}$, denotes operators, and the ‘calligraphy’ typeface, e.g., $\mathcal{X}$, denotes sets. We thus let $\mathcal{R}$ and $\mathcal{N}$ denote the real and natural numbers, respectively. For any matrix A, we denote by $\left|A\right|$ the determinant, by $\parallel A\parallel =\sqrt{\mathrm{tr}\left(A^{\prime }A\right)}$ the Euclidean norm, and by $\rho \left(A\right)$ the spectral radius. For some positive definite matrix A, we let $A^{1/2}$ denote the lower triangular Cholesky decomposition. For some function $f:{\mathcal{R}}^{d_z}\mapsto {\mathcal{R}}^{d_f}$, let $\partial f\left(z\right)/\partial z$ denote the derivative of $f\left(z\right)$ with respect to z. For some stochastic variable $z\in {\mathcal{R}}^{d_z}$ with Gaussian distribution with mean $\mu$ and covariance $\Sigma$, let $N(z;\mu ,\Sigma )$ denote the Gaussian probability density function evaluated at z. We let $p\left(z\right)$ denote the probability density of stochastic variable $z\in {\mathcal{R}}^{d_z}$ with respect to the $d_z$–dimensional Lebesgue measure m, while $p\left(\mathrm{d}z\right)=p\left(z\right)\mathrm{d}m$ denotes the corresponding probability measure. Additionally, the letter ‘p’ is generic notation for probability density functions and measures induced by the model defined in (1)–(3) below. The ‘bold’ typeface, e.g., $\mathit{p}$, is generic notation for analytically intractable quantities, in the sense of having no closed-form expression. Finally, we denote a sequence of $n\in {\mathcal{N}}_{+}$ real $d_z$–dimensional vectors by $z_{1:n}:={\left[\begin{array}{ccc}{z}_1^{\prime }& \dots & z_n^{\prime }\end{array}\right]}^{\prime }\in {\mathcal{R}}^{n\times d_z}$.

2. The Model

The structure of the SSR model is similar to the Granger-Johansen representation of the CVAR model, cf. Johansen (1996, chp. 4), but departs from it in two respects. First, the stationary component is a random coefficient autoregressive process, cf. e.g., Feigin and Tweedie (1985), rather than an autoregressive process. Second, the stationary and nonstationary components are observed with measurement error. This makes the SSR model is a state space model, whereas the CVAR model is observation-driven. In addition to resembling the CVAR model, the SSR model constitutes an extension of the CST model, cf. Chang et al. (2009). However, while the CST model is a linear Gaussian state space model, the SSR model is a nonlinear Gaussian state space model as it allows the stationary component to be a random coefficient autoregressive process.

Formally, we consider the observable p-dimensional discrete time vector process $y_t$, for $t=1,2,\dots ,T$ given by,

$$\begin{array}{c}{y}_t=C\left(y_0\right)+B\sum _{i=1}^t{\eta }_i+A{\xi }_t+u_t\end{array}$$

(1)

$$\begin{array}{c}{\xi }_t=\mu +{\mathsf{\Phi }}_t{\xi }_{t-1}+{\nu }_t,\end{array}$$

(2)

for fixed initial values $y_0$ and ${\xi }_0$, and with $u_t$, ${\mathsf{\Phi }}_t$ and ${[{\eta }_t^{\prime },{\nu }_t^{\prime }]}^{\prime }$ mutually independent. We define ${\epsilon }_t:={\sum }_{i=1}^t{\eta }_i$ with ${\epsilon }_0=0_{p-r}$. The sequences ${\epsilon }_{1:T}$ and ${\xi }_{1:T}$ are unobserved and take values ${\epsilon }_t\in {\mathcal{R}}^{p-r}$ and ${\xi }_t\in {\mathcal{R}}^r$ for $0<r<p$. Additionally, the matrices are of dimensions $A\in {\mathcal{R}}^{p\times r}$ and $B\in {\mathcal{R}}^{p\times p-r}$, with $\left[AB\right]\in {\mathcal{R}}^{p\times p}$ and invertible. Let the random coefficient, ${\mathsf{\Phi }}_t$, be i.i.d. Gaussian,

$$\begin{array}{c}\hfill \mathrm{vec}\left({\mathsf{\Phi }}_t\right)\sim N(\mathrm{vec}\left(\mathsf{\Phi }\right),{\mathsf{\Omega }}_{\mathsf{\Phi }}).\end{array}$$

(3)

with ${\mathsf{\Omega }}_{\mathsf{\Phi }}$ a positive definite covariance matrix. Let the observation error be i.i.d. Gaussian, such that $u_t\sim N(0,{\mathsf{\Omega }}_u)$ with ${\mathsf{\Omega }}_u$ a positive definite matrix, and let the innovations ${\eta }_t$ and ${\nu }_t$ be jointly Gaussian such that ${\eta }_t\sim N(0,{\mathsf{\Omega }}_{\eta })$ and ${\nu }_t\sim N(0,{\mathsf{\Omega }}_{\nu })$ with cross-covariance $\mathbb{C}\mathrm{ov}\left[{\eta }_t,{\nu }_t\right]={\mathsf{\Omega }}_{\eta ,\nu }$, such that the joint covariance matrix,

$$\begin{array}{cc}\hfill \mathsf{\Lambda }& :=\left[\begin{array}{cc}{\mathsf{\Omega }}_{\eta }& {\mathsf{\Omega }}_{\eta ,\nu }\\ {\mathsf{\Omega }}_{\eta ,\nu }^{\prime }& {\mathsf{\Omega }}_{\nu }\end{array}\right],\hfill \end{array}$$

(4)

is positive definite. Let all the introduced matrices be of appropriate dimensions and full rank. Furthermore, we introduce the orthogonal complements to A and B, which we denote $b\in {\mathcal{R}}^{p\times r}$ and $a\in {\mathcal{R}}^{p\times p-r}$, such that $b^{\prime }B=0$ and $a^{\prime }A=0$ with b and a of full column rank. Finally, we let $C\left(y_0\right):=B{\left(a^{\prime }B\right)}^{-1}{a}^{\prime }{y}_0$.

Define the parameter vectors,

$$\begin{array}{c}\omega :={\left[\begin{array}{ccc}\mathrm{vec}{\left(B\right)}^{\prime }& \mathrm{vec}{\left(A\right)}^{\prime }& \mathrm{vech}{\left({\mathsf{\Omega }}_u\right)}^{\prime }\end{array}\right]}^{\prime }\end{array}$$

(5)

$$\begin{array}{c}\lambda :={\left[\begin{array}{cccc}{\mu }^{\prime }& \mathrm{vec}{\left(\mathsf{\Phi }\right)}^{\prime }& \mathrm{vech}{\left({\mathsf{\Omega }}_{\mathsf{\Phi }}\right)}^{\prime }& \mathrm{vech}{\left(\mathsf{\Lambda }\right)}^{\prime }\end{array}\right]}^{\prime },\end{array}$$

(6)

which contain the parameters governing the observations $y_t$, and unobserved components ${\epsilon }_t$ and ${\xi }_t$, respectively. The parameter vectors take values in $\omega \in {\mathsf{\Theta }}_{\omega }$ and $\lambda \in {\mathsf{\Theta }}_{\lambda }$, respectively. Additionally, we define the full parameter vector as

$$\begin{array}{c}\theta :={\left[\begin{array}{cc}{\omega }^{\prime }& {\lambda }^{\prime }\end{array}\right]}^{\prime }\in {\mathsf{\Theta }}_{\omega }\times {\mathsf{\Theta }}_{\lambda }=:\mathsf{\Theta }.\end{array}$$

(7)

which indexes the model, and we refer to $\mathsf{\Theta }$ as the parameter space. Note that $\omega$ and $\lambda$ in $\theta$ are variation free in the sense of Engle et al. (1983). The parameter space is a subset of the $d_{\theta }$-dimensional Euclidean space $\mathsf{\Theta }\subseteq {\mathcal{R}}^{d_{\theta }}$, where $d_{\theta }$ denotes the number of elements in $\theta$. In the case where no restrictions are imposed on $\theta$, the dimension $d_{\theta }$ increases rapidly in r due to the $\frac{1}{2}(r^2+1)r^2$ parameters in ${\mathsf{\Omega }}_{\mathsf{\Phi }}$. We suggest restricting the off-diagonal elements of ${\mathsf{\Omega }}_{\mathsf{\Phi }}$ to zero to avoid over-parameterization. The number of parameters is then $d_{\theta }=2p^2+p+2r^2+r$ when the model is otherwise unrestricted.

The log-likelihood function for any parameter vector $\theta \in \mathsf{\Theta }$, fixed initial values $y_0\in {\mathcal{R}}^p$, ${\epsilon }_0=0_{p-r}$ and ${\xi }_0\in {\mathcal{R}}^r$, and observation sequence $y_{1:T}\in {\mathcal{R}}^{p\times T}$ is given by,

$$\begin{array}{cc}\hfill {\ell }_T\left(\theta \right)& :=log{\mathit{p}}_{\theta }({\epsilon }_0,{\xi }_0,y_{0:T}).\hfill \end{array}$$

(8)

The sample score is given by the first derivative of (8),

$$\begin{array}{cc}\hfill {\mathit{S}}_T\left(\theta \right)& :=\frac{\partial }{\partial \theta }{\ell }_T\left(\theta \right),\hfill \end{array}$$

(9)

and the observed Information matrix is given by minus the second derivative of (8),

$$\begin{array}{cc}\hfill {\mathit{I}}_T\left(\theta \right)& :=-\frac{{\partial }^2}{\partial \theta \partial {\theta }^{\prime }}{\ell }_T\left(\theta \right).\hfill \end{array}$$

(10)

Due to the nonlinear dynamics of the unobserved process (2), the log-likelihood function (8) and its derivatives (9)–(10) do not have closed-form solutions. In the following, we suppress the dependence on the initial values ${\epsilon }_0$, ${\xi }_0$ and $y_0$, but note they remain fixed.

3. Properties of the Process

In this section we consider some properties of the process defined by Equations (1)–(3) for a given parameter value $\theta \in \mathsf{\Theta }$. Specifically, we study the nonstationary and stationary components, including conditions on the parameter $\theta$ that ensure strict stationarity of the stationary component. Additionally, we decompose the observation $y_t$ into nonstationary and stationary directions.

3.1. The Unobserved Components

The first component of the model, ${\epsilon }_t$, is a random walk (RW) in $p-r$ dimensions, equivalently expressed as an autoregressive process with a unit root. That is, for $t=1,\dots ,T$,

$$\begin{array}{c}\hfill {\epsilon }_t={\epsilon }_{t-1}+{\eta }_t,\end{array}$$

(11)

with ${\epsilon }_0=0_{p-r}$. The process (11) admits the transition density $p_{\lambda }({\epsilon }_t\mid {\epsilon }_{t-1})$ with respect to the $p-r$–dimensional Lebesgue measure; however, it does not have a stationary distribution. This type of process has been studied extensively, see e.g., Dickey and Fuller (1979). In summary, the RW process is linear and Gaussian, but nonstationary.

The second unobserved component of the model, ${\xi }_t$, is a random coefficient autoregressive (RCAR) process of lag order one in r dimensions. The RCAR process (2)–(3) is observationally equivalent to a double autoregressive (DAR) process with one lag, cf. Ling (2007), which we formalize in Lemma 1.

Lemma 1.

For $\theta \in \mathsf{\Theta }$, the random coefficient autoregressive process (2)–(3) with $k=1$ has the following double autoregressive process representation, $t=1,2,\dots ,T$

$$\begin{array}{cc}\hfill {\xi }_t& =\mu +\mathsf{\Phi }{\xi }_{t-1}+{\mathsf{\Omega }}_{\nu ,t}^{1/2}{z}_t\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_{\nu ,t}& ={\mathsf{\Omega }}_{\nu }+\left({\xi }_{t-1}^{\prime }\otimes I_r\right){\mathsf{\Omega }}_{\mathsf{\Phi }}{\left({\xi }_{t-1}^{\prime }\otimes I_r\right)}^{\prime },\hfill \end{array}$$

(13)

for ${\xi }_0$ fixed, $z_t\sim N(0,I_r)$, cross-covariance $\mathbb{C}\mathrm{ov}\left[{\eta }_t,z_t\right]={\mathsf{\Omega }}_{\eta \nu }$, and with the joint innovation process ${\left[{\eta }_t^{\prime },z_t^{\prime }\right]}^{\prime }$ independent and identically distributed.

The DAR representation in Lemma 1 of the RCAR process in (2)–(3) characterizes the process dynamics in terms of the conditional mean and variance. The conditional mean ${\mathbb{E}}_{\lambda }\left[{\xi }_t\mid {\xi }_{t-1}\right]$ is autoregressive. However, the conditional variance $\mathbb{V}{\mathrm{ar}}_{\lambda }\left[{\xi }_t\mid {\xi }_{t-1}\right]$ depends positively on the lagged level ‘squared’. The conditional variance is heteroskedastic, but not in the well-known ARCH sense of e.g., Engle (1982); rather, the lagged level of the process ${\xi }_{t-1}$ enters the variance, not the lagged innovation ${\nu }_{t-1}$. To illustrate the point, we consider for a moment the conditional variance in the univariate case $r=1$, which is given by ${\omega }_{\nu ,t}^2={\omega }_{\nu }^2+{\omega }_{\varphi }^2{\xi }_{t-1}^2$. Here we see that a relatively large (in absolute terms) lagged level $|{\xi }_{t-1}|$ will result in a relatively large volatility ${\omega }_{\nu ,t}$ in the present period, and vice versa.

We make the following assumption on the random coefficients (3) in order to ensure strict stationarity of the RCAR process (2)–(3).

Assumption 1.

Assume that the top Lyapunov exponent is strictly negative,

$$\begin{array}{c}\mathit{\gamma }:=\underset{n\to \infty }{lim}\frac{1}{n}{\mathbb{E}}_{\lambda }\left[log∥\prod _{t=1}^n{\mathsf{\Phi }}_t∥\right]<0.\end{array}$$

(14)

Remark 1.

The top Lyapunov exponent (14) is intractable but can be approximated to arbitrary precision via simulation, cf. inter alia Ling (2007) and Francq and Zakoian (2010). The following approximation converges almost surely

$$\begin{array}{cc}\hfill {\widehat{\gamma }}_n& :=\frac{1}{n}log∥\prod _{t=1}^n{\mathsf{\Phi }}_t∥\stackrel{a.s.}{\to }\mathit{\gamma },\hfill \end{array}$$

(15)

as $n\to \infty$. In turn, ${\widehat{\gamma }}_n$ can be computed efficiently via the QR-decomposition, cf. Dieci and Van Vleck (1995).

Assumption 1 ensures that the RCAR process can be characterized as a geometrically ergodic Markov chain, cf. Meyn and Tweedie (2005). This is formalized in the following theorem.

Theorem 1 

(Feigin and Tweedie (1985), Theorem 3). Under Assumption 1, the process ${\left\{{\xi }_t\right\}}_{t=0,1,\dots }$ is geometrically ergodic. In particular, the initial value ${\xi }_0$ can be given an initial distribution ${\mathit{p}}_{\theta }\left({\xi }_0\right)$ such that ${\left\{{\xi }_t\right\}}_{t=0,1,\dots }$ is stationary and geometrically ergodic with some fractional moment.

Remark 2.

The stationary component, ${\xi }_t$, exhibits heavy-tailed behavior since it satisfies a stochastic recurrence equation. Pedersen and Wintenberger (2018) have recently considered the tail properties of processes of the form (2) for a more general specification of the random coefficient, ${\mathsf{\Phi }}_t$, that includes BEKK-ARCH and DAR-type processes as special cases. It should be possible to show that the stationary distribution of ${\xi }_t$ as defined in (2)–(3) also has power-law tails under suitable conditions.

The RCAR process (2)–(3) admits the transition density $p_{\lambda }({\xi }_t\mid {\xi }_{t-1})$ with respect to the r-dimensional Lebesgue measure. Moreover, the process has the stationary distribution ${\mathit{p}}_{\theta }\left({\xi }_t\right)$ under Assumption 1. In summary, the RCAR process is Gaussian and strictly stationary, but nonlinear.

3.2. The Observed Process

The observations ${\left\{y_t\right\}}_{t=1,2,\dots }$ are conditionally independent given the sequence of unobserved components ${\{{\epsilon }_t,{\xi }_t\}}_{t=1,2,\dots }$. Thus, the dynamics of the observed process are determined by the dynamics of the unobserved components.

We use the orthogonal complements $b^{\prime }$ and $a^{\prime }$ of the loading matrices B and A, respectively, and the skew-projection identity of Johansen (1996) to decompose the observation vector $y_t$ as follows,

$$\begin{array}{cc}\hfill y_t& =B_a{a}^{\prime }{y}_t+A_b{b}^{\prime }{y}_t,\hfill \end{array}$$

(16)

where we define $B_a:=B{\left(a^{\prime }B\right)}^{-1}$ and $A_b:=A{\left(b^{\prime }A\right)}^{-1}$. Here $a^{\prime }B$ and $b^{\prime }A$ are invertible thanks to our assumption that $\left[AB\right]$ is square and invertible. By premultiplying $y_t$ by $a^{\prime }$ we eliminate the stationary directions, while leaving the nonstationary directions,

$$\begin{array}{cc}\hfill a^{\prime }{y}_t& =a^{\prime }C\left(y_0\right)+a^{\prime }B{\epsilon }_t+a^{\prime }{u}_t.\hfill \end{array}$$

(17)

What is left after the linear transformation (17) is a random walk with Gaussian measurement error. Similarly, premultiplying $y_t$ by $b^{\prime }$ eliminates the nonstationary directions while the stationary directions remain,

$$\begin{array}{cc}\hfill b^{\prime }{y}_t& =b^{\prime }A{\xi }_t+b^{\prime }{u}_t.\hfill \end{array}$$

(18)

The process given by (18) is a stationary random coefficient autoregressive process with Gaussian measurement error.

The decomposition of the observation process (16) allows for a cointegration interpretation of the SSR model. The p observed variables in $y_t$ share $p-r$ common stochastic trends (17) with loading matrix $B_a$, while the r linear combinations (18) are stationary and load into the levels with the matrix $A_b$. The observed process admits the conditional density ${\mathit{p}}_{\theta }(y_t\mid y_{1:t-1})$ with respect to the p–dimensional Lebesgue measure; however, this density does not have a closed-form expression. Moreover, the observed process does not have a stationary distribution.

4. Likelihood-Based Estimation and Inference

In this section, we introduce the ML estimator and consider its asymptotic properties. We wish to conduct estimation and inference based on the true, but intractable, model likelihood. Due to the intractability of the likelihood, we can neither compute the ML estimator via numerical optimization of (8), nor compute classic standard errors via the observed Information matrix (10). We refer to the ML estimator as being ‘doubly intractable’, with reference to the concept from the literature in Bayesian statistics on models with intractable likelihoods, see e.g., Murray et al. (2006). It is beyond the scope of this paper to derive a full asymptotic theory for the SSR model. Instead, we conjecture the limiting properties of the likelihood function (8) and its derivatives (9)–(10). We obtain the asymptotic properties for the ML estimator based on the conjecture.

We recall preliminarily that the ML estimator is defined as the parameter vector $\theta \in \mathsf{\Theta }$ that maximizes the log-likelihood function (8),

$$\begin{array}{cc}\hfill {\widehat{\mathit{\theta }}}_T& :=\underset{\theta \in \mathsf{\Theta }}{argsup}{\ell }_T\left(\theta \right),\hfill \end{array}$$

(19)

noting that the ML estimator (19) is a function of the observation sequence $y_{1:T}$. We denote by ${\theta }^{*}\in \mathsf{\Theta }$ the true parameter value for the data generating process (1)–(3). In the following, we make the below conjecture on the asymptotic properties of (8)–(10). Note that, having assumed that $B^{*}$ is known, the score, information, and likelihood in the conjecture refer to the unknown parameters only; that is, all elements in $\theta$ excluding $\mathrm{vec}\left(B\right)$.

Conjecture 1.

If Assumption 1 holds, $B^{*}$ is known, and ${\theta }^{*}\in \mathsf{\Theta }\subseteq {\mathcal{R}}^{d_{\theta }}$, then the log-likelihood function ${\ell }_T(·):{\mathcal{R}}^{d_{\theta }}\mapsto \mathcal{R}$ is three times continuously differentiable in θ, and

1. 

$\frac{1}{\sqrt{T}}{\mathit{S}}_T\left({\theta }^{*}\right)\stackrel{D}{\to }N(0,{\mathsf{\Omega }}_S)$ as $T\to \infty$, with ${\mathsf{\Omega }}_S>0$,

2. 

$\frac{1}{T}{\mathit{I}}_T\left({\theta }^{*}\right)\stackrel{P}{\to }{\mathsf{\Omega }}_I$ as $T\to \infty$, with ${\mathsf{\Omega }}_I>0$, and

3. 

${max}_{h,i,j=1,\dots ,d_{\theta }}{sup}_{\theta \in \mathcal{N}\left({\theta }^{*}\right)}|{\partial }^3{\ell }_T\left(\theta \right)/\partial {\theta }_h\partial {\theta }_i\partial {\theta }_j|\le c_T$,

where $\mathcal{N}\left({\theta }^{*}\right)$ is a neighborhood of ${\theta }^{*}$ and $0\le c_T\stackrel{P}{\to }c$, $0<c<\infty$, as $T\to \infty$.

Remark 3.

Theorem 3 in Bohn Nielsen and Rahbek (2014) shows that Conjecture 1 holds in the case of the strictly stationary bivariate double autoregressive model with BEKK-type time-varying covariance. With $B^{*}$ known, the SSR model corresponds closely to this model plus Gaussian measurement errors.

It should be noted that we propose Conjecture 1 despite lack of finite moments of the RCAR process, cf. Theorem 1. This is in line with the results of inter alia Bohn Nielsen and Rahbek (2014) for the bivariate DAR model, and Ling (2004, 2007) for the univariate DAR model.

The result in Theorem 2 below states that if Conjecture 1 holds true, then the ML estimator (19) is unique, $\sqrt{T}$–consistent and asymptotically Gaussian. The result follows from applying Lemma 1 in Jensen and Rahbek (2004), the conditions of which correspond to (1.)–(3.) of Conjecture 1.

Theorem 2 

(Jensen and Rahbek (2004), Lemma 1). If Conjecture 1 holds, then there exists a fixed open neighborhood $\mathcal{U}\left({\theta }^{*}\right)\subseteq \mathcal{N}\left({\theta }^{*}\right)$ of the true parameter ${\theta }^{*}$, which is an interior point of Θ, such that with probability tending to one as $T\to \infty$, there exists a minimum point ${\widehat{\mathit{\theta }}}_T$ in $\mathcal{U}\left({\theta }^{*}\right)$ and ${\ell }_T\left(\theta \right)$ is convex in $\mathcal{U}\left({\theta }^{*}\right)$. In particular, ${\widehat{\mathit{\theta }}}_T$ is unique and satisfies the score equation

$$\begin{array}{c}\hfill {\mathit{S}}_T\left({\widehat{\mathit{\theta }}}_T\right)=0.\end{array}$$

(20)

Additionally, the ML estimator is consistent ${\widehat{\mathit{\theta }}}_T\to {\theta }^{*}$, and asymptotically Gaussian,

$$\begin{array}{cc}\hfill \sqrt{T}({\widehat{\mathit{\theta }}}_T-{\theta }^{*})& \stackrel{D}{\to }N(0,{\mathsf{\Omega }}_I^{-1}{\mathsf{\Omega }}_S{\mathsf{\Omega }}_I^{-1}),T\to \infty .\hfill \end{array}$$

(21)

Proof. 

Conjecture 1 satisfies the Cramer-type conditions of Lemma 1 in Jensen and Rahbek (2004), which provides the result. ☐

We assume that the true value of B is known, because Chang et al. (2009) showed that the ML estimator of the loading matrix B exhibits T-convergence and is asymptotically mixed Gaussian in the CST model. The CST model corresponds to the SSR model with $p-r=1$, but without the stationary components, i.e., $A=0_{p\times r}$ for any p. We find it reasonable to believe that this result carries over to the SSR model. Moreover, fixing B is conceptually similar to classic cointegration analysis with known cointegrating vectors, which is an accepted starting point for new methodological developments, see e.g., Bec and Rahbek (2004). In applications we often have a predefined set of cointegrating vectors that we are interested in. In the context of the SSR model, the cointegrating vectors correspond to the rows of the orthogonal complement $b^{\prime }$. As an example, for the empirical illustration in Section 10 we consider an interest rate spread in a bivariate system with one common stochastic trend, i.e., $p=2$ and $p-r=1$. The spread implies $b^{\prime }=\left[\begin{array}{cc}1& -1\end{array}\right]$, which in turn corresponds to the loading matrix $B={\left[\begin{array}{cc}1& 1\end{array}\right]}^{\prime }$ when normalizing on the first element.

The Fisher Information matrix, ${\mathsf{\Omega }}_I$, is consistently estimated by the (scaled) observed Information matrix evaluated at ${\widehat{\mathit{\theta }}}_T$, cf. Conjecture 1.(3.). Moreover, the asymptotic variance of the score, ${\mathsf{\Omega }}_S$, is equal to the Fisher Information matrix when the model is well-specified; the information matrix equality holds, cf. e.g., Hamilton (1994, sct. 14.4). In this case, the asymptotic variance of the ML estimator (19) is simply the inverse Fisher Information matrix. Thus, we can use classic standard errors, that are based on the observed Information matrix (10), to conduct inference on the ML estimates.

5. The Incomplete Data Framework

In this section, we appeal to the incomplete data framework of Dempster et al. (1977) to deal with the unobserved components of the SSR model. We first formulate the state space representation of the model in (1)–(3) and its associated optimal filtering problem. Secondly, we formulate the intractable sample score (9) and observed information matrix (10) in terms of the optimal filtering problem. In Section 6 we introduce a particle filter algorithm with which we can approximate the optimal filtering problem. This enables approximation of the intractable sample score and observed information matrix via the particle filter algorithm.

5.1. The State Space Form and the Optimal Filtering Problem

Preliminarily, we collect the unobserved components in the vector $x_t:={\left[\begin{array}{cc}{\epsilon }_t^{\prime }& {\xi }_t^{\prime }\end{array}\right]}^{\prime }$, which we refer to as the state vector. The unobserved components are Markov, see (11)–(13), and the observation depends only on the contemporary values of the unobserved components. Thus, the SSR model in (1)–(3) has the dependency structure of a state space model. Formally, for $t=1,\dots ,T$, the SSR model in (1)–(3) has the following state space representation,

$$\begin{array}{c}{y}_t=C\left(y_0\right)+\mathsf{\Pi }{x}_t+{\mathsf{\Omega }}_u^{1/2}{u}_t\end{array}$$

(22)

$$\begin{array}{c}{x}_t=\alpha +\mathsf{\Gamma }{x}_{t-1}+{\mathsf{\Lambda }}_t^{1/2}{v}_t,\end{array}$$

(23)

with $y_0$ and $x_0$ fixed, $u_t\sim N(0,I_p)$ and $v_t\sim N(0,I_p)$, and $u_t$ and $v_t$ mutually independent. We define accordingly,

$$\begin{array}{cc}\hfill \mathsf{\Pi }& :={\left[\begin{array}{c}{B}^{\prime }\\ A^{\prime }\end{array}\right]}^{\prime },\alpha :=\left[\begin{array}{c}0\\ \mu \end{array}\right],\mathsf{\Gamma }:=\left[\begin{array}{cc}{I}_{p-r}& 0\\ 0& \mathsf{\Phi }\end{array}\right]\mathrm{and}{\mathsf{\Lambda }}_t:=\left[\begin{array}{cc}{\mathsf{\Omega }}_{\eta }& {\mathsf{\Omega }}_{\eta ,\nu }\\ {\mathsf{\Omega }}_{\eta ,\nu }^{\prime }& {\mathsf{\Omega }}_{\nu ,t}\end{array}\right],\hfill \end{array}$$

(24)

and recall that ${\mathsf{\Omega }}_{\nu ,t}$ is defined in Lemma (1). We refer to (22) as the observation equation, and to (23) as the transition equation. It is easy to verify that the state space representation in (22) and (23) is observationally equivalent to the SSR model as presented in (1)–(3). The observation and transition equations admit the densities with respect to the p-dimensional Lebesgue measure,

$$\begin{array}{c}{p}_{\omega }\left(y_t\mid x_t\right)=N(y_t;C\left(y_0\right)+\mathsf{\Pi }{x}_t,{\mathsf{\Omega }}_u)\end{array}$$

(25)

$$\begin{array}{c}{p}_{\lambda }\left(x_t\mid x_{t-1}\right)=N(x_t;\alpha +\mathsf{\Gamma }{x}_{t-1},{\mathsf{\Lambda }}_t),\end{array}$$

(26)

respectively. We refer to (25) as the observation density and to (26) as the transition density. As mentioned previously, we suppress the dependence on the initial observation $y_0$.

One approach to conducting inference on the unobserved components, i.e., the state vector $x_t$, is the optimal filtering problem, cf. Anderson and Moore (1979). The optimal filtering problem refers to the general problem of computing the conditional expectation of some sequence of unobserved states given some sequence of observations. In the following, we consider the specific instance of the optimal filtering problem known as the smoothing problem. Formally, the smoothing problem is a conditional expectation of the form,

$$\begin{array}{cc}\hfill {\mathbb{E}}_{\theta }\left[{\gamma }_t\left(x_{1:t}\right)\mid y_{1:t}\right]& =\int {\gamma }_t\left(x_{1:t}\right){\mathit{p}}_{\theta }\left(x_{1:t}\mid y_{1:t}\right)\mathrm{d}{x}_{1:t},\hfill \end{array}$$

(27)

for any function ${\gamma }_t\left(x_{1:t}\right)\in L^1\left[{\mathcal{R}}^{tp},{\mathit{p}}_{\theta }\left(x_{1:t}\mid y_{1:t}\right)\right]$ and point in time $t\in \{1,\dots ,T\}$. We refer to the function ${\gamma }_t\left(x_{1:t}\right)$ as the test function and to the density ${\mathit{p}}_{\theta }\left(x_{1:t}\mid y_{1:t}\right)$ as the smoothing density. The test function may be time-varying, but of known form for a fixed observation sequence $y_{1:T}$. The smoothing density in (27) can be expressed as the recursion of the lagged smoothing density,

$$\begin{array}{cc}\hfill {\mathit{p}}_{\theta }\left(x_{1:t}\mid y_{1:t}\right)& =\frac{p_{\omega }\left(y_t\mid x_t\right)p_{\lambda }\left(x_t\mid x_{t-1}\right)}{{\mathit{p}}_{\theta }\left(y_t\mid y_{1:t-1}\right)}{\mathit{p}}_{\theta }\left(x_{1:t-1}\mid y_{1:t-1}\right),\hfill \end{array}$$

(28)

initialized with ${\mathit{p}}_{\theta }\left(x_1\mid x_0,y_0,y_1\right)$. The normalizing constant in (28) is the likelihood contribution, which is given by the integral,

$$\begin{array}{cc}\hfill {\mathit{p}}_{\theta }\left(y_t\mid y_{1:t-1}\right)& =\int p_{\omega }\left(y_t\mid x_t\right)p_{\lambda }\left(x_t\mid x_{t-1}\right){\mathit{p}}_{\theta }\left(x_{1:t-1}\mid y_{1:t-1}\right)\mathrm{d}{x}_{1:t}.\hfill \end{array}$$

(29)

We note the smoothing density recursion (28) is intractable due to the intractability of the likelihood contribution (29). In the following, we will use the smoothing problem (27) to address computation of the sample score (9) and observed Information matrix (10).

5.2. The Sample Score and Observed Information as Smoothing Problems

The incomplete data framework is closely associated with the classic expectation maximization (EM) algorithm, introduced in Dempster et al. (1977). The EM algorithm is a common approach to maximizing the log-likelihood function (8) to obtain the ML estimator (19) for models with unobserved variables. When the EM algorithm is applicable, it is also possible to evaluate the sample score (9) and observed Information matrix (10). For the SSR model, however, the EM algorithm does not apply directly, yet we may use the incomplete data framework to reformulate the sample score and observed Information in terms of intractable smoothing problems of the form (27).

A central concept of the EM algorithm is the auxiliary function called the intermediate quantity, which is defined as,

$$\begin{array}{cc}\hfill {\mathit{Q}}_T\left(\theta \mid \vartheta \right)& :=\int log{p}_{\theta }\left(y_{1:T},x_{1:T}\right){\mathit{p}}_{\vartheta }\left(x_{1:T}\mid y_{1:T}\right)\mathrm{d}{x}_{1:T}\hfill \\ & ={\mathit{\ell }}_T\left(\theta \right)-{\mathit{H}}_T(\theta \mid \vartheta ),\hfill \end{array}$$

(30)

where

$$\begin{array}{cc}\hfill {\mathit{H}}_T(\theta \mid \vartheta )& :=-\int log{p}_{\theta }(x_{1:T}\mid y_{1:T}){\mathit{p}}_{\vartheta }(x_{1:T}\mid y_{1:T})\mathrm{d}{x}_{1:T},\hfill \end{array}$$

(31)

for any parameter values $\theta ,\vartheta \in \mathsf{\Theta }$. We refer to $log{p}_{\theta }\left(y_{1:T},x_{1:T}\right)$ as the complete data log-likelihood. By the state space model structure (22)–(23) and variation freeness of $\theta$ defined in (7), we have that the complete data log-likelihood is given by,

$$\begin{array}{cc}\hfill log{p}_{\theta }\left(y_{1:T},x_{1:T}\right)& =\sum _{t=1}^T\left[log{p}_{\omega }\left(y_t\mid x_t\right)+log{p}_{\lambda }\left(x_t\mid x_{t-1}\right)\right].\hfill \end{array}$$

(32)

The intermediate quantity (30) is sometimes also called the expected log-likelihood, since it is interpretable as the conditional expectation of the complete data log-likelihood (32) given the observations $y_{1:T}$. We note the term separating the log-likelihood (8) and the intermediate quantity (30) is the entropy of the smoothing density (28) with parameters $\vartheta$ and $\theta$, defined in (31).

We are interested in the intermediate quantity (30) because it provides a convenient way to derive the sample score and observed Information matrix in terms of the derivatives of the complete data log-likelihood (32). The first and second derivatives of the complete data log-likelihood function in (32) are the sum of the first and second order derivatives of the observation and transition log-densities with respect to $\omega$ and $\lambda$, respectively. These can be computed by either analytical or numerical differentiation of (32). For $\vartheta \in \mathsf{\Theta }$, we define the derivatives of (32) in terms of the functions,

$$\begin{array}{cc}\hfill U_T\left(x_{1:T};\vartheta \right)& :=\frac{\partial }{\partial \theta }log{p}_{\theta }\left(y_{1:T},x_{1:T}\right){|}_{\theta =\vartheta }=\sum _{t=1}^T{u}_t\left(x_t,x_{t-1};\vartheta \right)\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill V_T\left(x_{1:T};\vartheta \right)& :=\frac{{\partial }^2}{\partial \theta \partial {\theta }^{\prime }}log{p}_{\theta }\left(y_{1:T},x_{1:T}\right){|}_{\theta =\vartheta }=\sum _{t=1}^T{v}_t\left(x_t,x_{t-1};\vartheta \right),\hfill \end{array}$$

(34)

where, taking advantage of the variation freeness of the model parameter, $\theta$, we define the summands of (33) and (34), respectively, as

$$\begin{array}{cc}\hfill u_t\left(x_t,x_{t-1};\vartheta \right)& :=\left[\begin{array}{c}\frac{\partial }{\partial \omega }log{p}_{\omega }\left(y_t\mid x_t\right)\\ \frac{\partial }{\partial \lambda }log{p}_{\lambda }\left(x_t\mid x_{t-1}\right)\end{array}\right]{|}_{\theta =\vartheta },\hfill \end{array}$$

(35)

and

$$\begin{array}{cc}\hfill v_t\left(x_t,x_{t-1};\vartheta \right)& :=\left[\begin{array}{cc}\frac{{\partial }^2}{\partial \omega \partial {\omega }^{\prime }}log{p}_{\omega }\left(y_t\mid x_t\right)& 0_{d_{\omega }\times d_{\lambda }}\\ 0_{d_{\lambda }\times d_{\omega }}& \frac{{\partial }^2}{\partial \lambda \partial {\lambda }^{\prime }}log{p}_{\lambda }\left(x_t\mid x_{t-1}\right)\end{array}\right]{|}_{\theta =\vartheta },\hfill \end{array}$$

(36)

We note that the functions (35) and (36) should not be confused with the measurement error in (22) and innovations in (23), respectively.

If the first and second order derivatives of the complete data log-likelihood in (33) and (34), respectively, are integrable with respect to the smoothing density (28), then we may appeal to Fisher’s and Louis’ identities (defined below) to express the sample score (9) and observed Information matrix (10) in terms of smoothing problems of the form (27).

Conjecture 2.

For any $\theta \in \mathsf{\Theta }$ and observation sequence $y_{1:T}\in {\mathcal{R}}^{p\times T}$, it holds that $U_T(x_{1:T};\theta )\in L^2[{\mathcal{R}}^{p\times T},{\mathit{p}}_{\theta }(x_{1:T}\mid y_{1:T})]$ and $V_T(x_{1:T};\theta )\in L^1[{\mathcal{R}}^{p\times T},{\mathit{p}}_{\theta }(x_{1:T}\mid y_{1:T})]$.

For the same reasons we conjectured the asymptotic properties of the true log-likelihood function, sample score, observed information matrix, we conjecture integrability of the derivatives of the complete data log-likelihood (33) and (34).

Fisher’s identity, cf. Dempster et al. (1977), states the first derivative of the intermediate quantity (30) is equivalent to the sample score (9). Similarly, Louis’ identity of Louis (1982) establishes a relation between the first and second derivatives of the intermediate quantity (30) and the observed Information matrix (10).

Lemma 2 

(Fisher’s and Louis’ identities, cf. Cappé et al. (2005), Proposition 10.1.6). If Conjecture 2 holds and $\theta \in \mathsf{\Theta }$, then the sample score (9) is equivalently given by

$$\begin{array}{cc}\hfill {\mathit{S}}_T\left(\theta \right)& =\int U_T(x_{1:T};\theta ){\mathit{p}}_{\theta }\left(x_{1:T}\mid y_{1:T}\right)\mathrm{d}{x}_{1:T},\hfill \end{array}$$

(37)

and the observed Information (10) is equivalently given by

$$\begin{array}{cc}\hfill {\mathit{I}}_T\left(\theta \right)& ={\mathit{S}}_T\left(\theta \right){\mathit{S}}_T{\left(\theta \right)}^{\prime }-{\mathit{G}}_T\left(\theta \right)-{\mathit{K}}_T\left(\theta \right),\hfill \end{array}$$

(38)

where

$$\begin{array}{cc}\hfill {\mathit{G}}_T\left(\theta \right)& :=\int V_T(x_{1:T};\theta ){\mathit{p}}_{\theta }\left(x_{1:T}\mid y_{1:T}\right)\mathrm{d}{x}_{1:T}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill {\mathit{K}}_T\left(\theta \right)& :=\int U_T(x_{1:T};\theta )U_T{(x_{1:T};\theta )}^{\prime }{\mathit{p}}_{\theta }\left(x_{1:T}\mid y_{1:T}\right)\mathrm{d}{x}_{1:T},\hfill \end{array}$$

(40)

and the functions $U_T(x_{1:T};\theta )$ and $V_T(x_{1:T};\theta )$ are defined in (33) and (34), respectively.

Although Lemma 2 shows the sample score (9) and observed Information (10) can be restated as smoothing problems of the form (27), we still cannot obtain closed-form expressions due to the intractability of the optimal filtering problem, cf. Section 5.1. In the next section, we introduce a particle filter algorithm that can approximate smoothing problems for appropriately chosen test functions, such as the functions $U_T(x_{1:T};\theta )$ and $V_T(x_{1:T};\theta )$ under Conjecture 2.

6. Particle Filter-Based Approximations

In this section, we introduce a particle filter algorithm that produces pointwise approximations to the true but intractable log-likelihood function (8), sample score (9), and observed Information matrix (10) for any parameter $\theta \in \mathsf{\Theta }$ and fixed observation sequence $y_{1:T}\in {\mathcal{R}}^{p\times T}$. In Section 7, we show how to apply the particle filter-based approximations introduced in this section to approximate the true, intractable ML estimator and classic standard errors, which we introduced in Section 4.

6.1. Particle Filtering

A particle filter is a simulation-based algorithm that produces approximations to smoothing problems of the form (27) for state space models. We introduce here a standard particle filter, which produces empirical measures that recursively approximate the smoothing density (28) for each time point in the observed sample $t\in \{1,\dots ,T\}$. The empirical measures consist of point masses, which we refer to as particles, and we use these for Monte Carlo integration in order to approximate the smoothing problem (27). Additionally, the particle filter produces a point-wise approximation of the log-likelihood function as a by-product. For an introduction to particle filtering in the context of economics and finance see Creal (2012).

The particle filter algorithm relies on an importance density, denoted $q_{\theta }\left(x_{1:t}\mid y_{1:t}\right)$, that has the same support and recursive structure as the smoothing density (28). Formally, for $t=1,\dots ,T$, we define the importance density as,

$$\begin{array}{cc}\hfill q_{\theta }\left(x_{1:t}\mid y_{1:t}\right)& :=q_{\theta }\left(x_t\mid x_{t-1},y_t\right)q_{\theta }\left(x_{1:t-1}\mid y_{1:t-1}\right),\hfill \end{array}$$

(41)

initialized by $q_{\theta }\left(x_1\mid x_0,y_0,y_1\right)$. We note the importance density (41) is defined recursively by $q_{\theta }\left(x_t\mid x_{t-1},y_t\right)$, which we refer to as the importance transition density.

Assuming the smoothing density (28) is absolutely continuous with respect to the importance density (41), we can write the former as a the product of the importance density and a weight function,

$$\begin{array}{cc}\hfill {\mathit{p}}_{\theta }\left(x_{1:t}\mid y_{1:t}\right)& ={\overline{\mathit{w}}}_t\left(x_{1:t}\right)q_{\theta }\left(x_{1:t}\mid y_{1:t}\right),{\overline{\mathit{w}}}_t\left(x_{1:t}\right):=\frac{{\mathit{p}}_{\theta }\left(x_{1:t}\mid y_{1:t}\right)}{q_{\theta }\left(x_{1:t}\mid y_{1:t}\right)}.\hfill \end{array}$$

(42)

We refer to the weight function ${\overline{\mathit{w}}}_t\left(x_{1:t}\right)$ as the normalized importance weight. We note that (42) constitutes a change of measure from the smoothing density to the importance density, and the normalized importance weight is a Radon-Nikodym derivative between the two densities.

Substituting the recursive expressions for the smoothing density (28) and importance density (41) into the expression for the normalized importance weight in (42), we obtain a recursive expression for the normalized importance weight,

$$\begin{array}{cc}\hfill {\overline{\mathit{w}}}_t\left(x_{1:t}\right)& =\frac{{\tilde{w}}_t\left(x_{t-1:t}\right)}{{\mathit{p}}_{\theta }\left(y_t\mid y_{1:t-1}\right)}{\overline{\mathit{w}}}_{t-1}\left(x_{1:t-1}\right),\hfill \end{array}$$

(43)

where we define

$$\begin{array}{cc}\hfill {\tilde{w}}_t\left(x_{t-1:t}\right)& :=\frac{p_{\omega }\left(y_t\mid x_t\right)p_{\lambda }\left(x_t\mid x_{t-1}\right)}{q_{\theta }\left(x_t\mid x_{t-1},y_t\right)}.\hfill \end{array}$$

(44)

We refer to (44) as the incremental importance weights. The recursion for the normalized importance weight (43) is normalized by the likelihood contribution (29) and is therefore also intractable.

For particle filtering in general, the importance transition density is subject to choice under mild regularity conditions, cf. e.g., Assumption 9.4.1 in Cappé et al. (2005). We let the importance transition density be the corresponding model density; formally,

$$\begin{array}{cc}\hfill q_{\theta }\left(x_t\mid x_{t-1},y_t\right)& :=p_{\theta }\left(x_t\mid x_{t-1},y_t\right).\hfill \end{array}$$

(45)

We refer to (45) as the locally optimal transition density. This choice of importance transition density is optimal in the sense that it is conditional on the the contemporary observation $y_t$, cf. Doucet et al. (2000). This is sometimes also referred to as ‘fully adapted’, cf. e.g., Pitt and Shephard (1999b). If we instead let the importance transition density be the model transition density (26), we omit the information about $x_t$ that is contained in $y_t$. The locally optimal transition density is not necessarily available in closed-form for nonlinear state space models. It is, however, available for the SSR model and we present it in Lemma 3.

Lemma 3.

For $\theta \in \mathsf{\Theta }$, the locally optimal transition density has the closed-form expression

$$\begin{array}{c}\hfill p_{\theta }\left(x_t\mid x_{t-1},y_t\right)=N(x_t;{\mu }_{t|t}^x,{\Sigma }_{t|t}^x),\end{array}$$

(46)

where the conditional mean and variance are given by,

$$\begin{array}{c}{\mu }_{t|t}^x={\mu }_{t|t-1}^x+{\Sigma }_{t|t-1}^x{\mathsf{\Pi }}^{\prime }{\left[{\Sigma }_{t|t-1}^y\right]}^{-1}\left(y_t-{\mu }_{t|t-1}^y\right)\end{array}$$

(47)

$$\begin{array}{c}{\Sigma }_{t|t}^x={\Sigma }_{t|t-1}^x-{\Sigma }_{t|t-1}^x{\mathsf{\Pi }}^{\prime }{\left[{\Sigma }_{t|t-1}^y\right]}^{-1}\mathsf{\Pi }{\Sigma }_{t|t-1}^x,\end{array}$$

(48)

with

$$\begin{array}{c}{\mu }_{t|t-1}^y=C\left(y_0\right)+\mathsf{\Pi }{\mu }_{t|t-1}^x\end{array}$$

(49)

$$\begin{array}{c}{\Sigma }_{t|t-1}^y=\mathsf{\Pi }{\Sigma }_{t|t-1}^x{\mathsf{\Pi }}^{\prime }+{\mathsf{\Omega }}_u\end{array}$$

(50)

$$\begin{array}{c}{\mu }_{t|t-1}^x=\alpha +\mathsf{\Gamma }{x}_{t-1}\end{array}$$

(51)

$$\begin{array}{c}{\Sigma }_{t|t-1}^x={\mathsf{\Lambda }}_t,\end{array}$$

(52)

and the state space form definitions given in (24).

Remark 4.

The locally optimal transition density (46) is related to the Kalman (1960) filter, which solves the optimal filtering problem analytically for linear and Gausian models. Equations (49)–(52) correspond the Kalman filter for a known value of $x_{t-1}$. Related methods for efficient particle filtering include the mixture Kalman filter and Rao-Blackwellisation, cf. Chen and Liu (2000) and Andrieu and Doucet (2002).

It is straightforward to use the general expression for the incremental importance weight in (44) to show that letting the importance transition density be the locally optimal transition density, i.e., (45), results in the following specific expression for incremental importance weights,

$$\begin{array}{cc}\hfill {\tilde{w}}_t\left(x_{t-1}\right)& =p_{\theta }\left(y_t\mid x_{t-1}\right).\hfill \end{array}$$

(53)

We refer to the density in (53) as the predictive observation density. It has a closed-form expression that follows from the closed-form expression of the locally optimal transition density in Lemma 3.

Corollary 1.

For $\theta \in \mathsf{\Theta }$, the predictive observation density has the closed-form expression

$$\begin{array}{cc}\hfill p_{\theta }\left(y_t\mid x_{t-1}\right)& =N(y_t;{\mu }_{t|t-1}^y,{\Sigma }_{t|t-1}^y),\hfill \end{array}$$

(54)

recalling the definitions in (49)–(52).

Proof. 

Contained in the proof of Lemma 3. ☐

Remark 5.

The choice of importance transition density (45) is locally optimal in the sense that the conditional variance of the incremental importance weights (53) given $x_{t-1}$ is zero, cf. Doucet et al. (2000).

The particle filter, presented in Algorithm 1 below, produces weighted particle samples approximately distributed as the smoothing density (28) at each point in time $t=1,\dots ,T$. The algorithm consists of iterating over three steps. At point t in time, the first step is to sample N particles, denoted ${\{{\tilde{x}}_{1:t}^{\left(i\right)}\}}_{i=1}^N$, from the importance density (41) given the particle sample from $t-1$. This is called the propagation step. Step two consists of computing self-normalized importance weights, denoted ${\{{\overline{w}}_t^{\left(i\right)}\}}_{i=1}^N$, that approximate the normalized importance weights (43). This is the weighting step. The third step is to sample N particle indices, denoted ${\{I^{\left(i\right)}\}}_{i=1}^N$, with replacement. We sample index j with probability ${\overline{w}}_t^{\left(j\right)}$ for $j\in \{1,\dots ,N\}$. We retain the number of particles indicated by the resulting sample of particle indices, denoted ${\{x_{1:t}^{\left(i\right)}\}}_{i=1}^N$, and let the importance weights be uniform. This is the resampling step. After resampling, we store the particle samples and proceed to $t+1$.

For a fixed parameter value $\theta \in \mathsf{\Theta }$ and observation sequence $y_{1:T}\in {\mathcal{R}}^{p\times T}$, we run the locally optimal particle filter for the SSR model as specified in Algorithm 1 below.

Algorithm 1: Locally Optimal Particle Filter.

Given a parameter $\theta \in \mathsf{\Theta }$, initialize by setting $x_0^{\left(i\right)}:=x_0$ and ${\overline{w}}_0^{\left(i\right)}:=1/N$ for $i=1,\dots ,N$. For $t=0,1,\dots ,T$: Sample particles ${\{{\tilde{x}}_t^{\left(i\right)}\}}_{i=1}^N$ with distribution $$\begin{array}{c}\hfill {\tilde{x}}_t^{\left(i\right)}\sim p_{\theta }(x_t\mid x_{t-1}^{\left(i\right)},y_t),\end{array}$$ (55) and set ${\tilde{x}}_{1:t}^{\left(i\right)}:=\left[\begin{array}{cc}{x}_{1:t-1}^{\left(i\right)}& {\tilde{x}}_t^{\left(i\right)}\end{array}\right]$ for $i=1,2,\dots ,N$. Calculate the unnormalized importance weights, ${\{w_t^{\left(i\right)}\}}_{i=1}^N$, $$\begin{array}{cc}\hfill w_t^{\left(i\right)}& =p_{\theta }\left(y_t\right|{\tilde{x}}_{t-1,t}^{\left(i\right)}){\overline{w}}_{t-1}^{\left(i\right)},\hfill \end{array}$$ (56) for $i=1,\dots ,N$. Then compute the normalized importance weights $$\begin{array}{cc}\hfill {\overline{w}}_t^{\left(i\right)}& =\frac{w_t^{\left(i\right)}}{W_t^N},W_t^N:=\sum _{i=1}^N{w}_t^{\left(i\right)},\hfill \end{array}$$ (57) for $i=1,\dots ,N$. Sample N particle indices ${\{I^{\left(i\right)}\}}_{i=1}^N$, $I^{\left(i\right)}\in \{1,\dots ,N\}$, with probabilities $$\begin{array}{cc}\hfill Pr(I^{\left(i\right)}=j\mid {\tilde{\mathcal{F}}}_t,y_{1:t})& ={\overline{w}}_t^{\left(j\right)},j\in \{1,\dots ,N\}\hfill \end{array}$$ (58) for $i=1,\dots ,N$. Set the resampled particles $x_{1:t}^{\left(i\right)}:={\tilde{x}}_{1:t}^{\left(I^{\left(i\right)}\right)}$, and the normalized importance weights ${\overline{w}}_t^{\left(i\right)}:=1/N$ for $i=1,\dots ,N$.

Remark 6.

The resampling method applied in step (3.) of Algorithm 1 is known as multinomial resampling. Alternative methods that are guaranteed to produce lower Monte Carlo variance exists, cf. Douc et al. (2005). We consider multinomial resampling for its analytical tractability, and recommend applying one of the more efficient alternatives in practice.

Remark 7.

The notation $x_{1:t}^{\left(i\right)}$ is ambiguous due to the resampling step of Algorithm 1, since the elements of the ith particle path at time $t-1$, denoted $x_{1:t-1}^{\left(i\right)}$, are not necessarily the same as the first $t-1$ elements of the ith particle path at time t, denoted $x_{1:t}^{\left(i\right)}$. By convention, $x_{1:t}^{\left(i\right)}$ always refers to the particle chain after resampling at time t (similarly ${\tilde{x}}_{1:t}^{\left(i\right)}$ refers to the chain before resampling). We refer to elements k to l of the ith particle chain after resampling at time t as $x_{l:k,t}^{\left(i\right)}$.

The particle filter in Algorithm 1 produces two particle samples at each point in time, t. The first set, ${\{{\tilde{x}}_{1:t}^{\left(i\right)}\}}_{i=1}^N$, is produced at the propagation step (1.) and is associated with importance weights in the weighting step (2.), ${\{{\overline{w}}_t^{\left(i\right)}\}}_{i=1}^N$. The second set, ${\{x_{1:t}^{\left(i\right)}\}}_{i=1}^N$, is produced at the resampling step (3.). Both sets are approximately drawn from the smoothing density (28). We note the resampling step introduces additional sampling error, cf. Chopin (2004), so we calculate approximations using the weighted sample unless otherwise specified.

The particle filter iterates over over the propagation, weighting and resampling steps throughout the sequence, $t=1,\dots ,T$, after which the algorithm terminates. We note the two sets of particles produced during each iteration are themselves random variables measurable with respect to the sub-$\sigma$-algebras ${\tilde{\mathcal{F}}}_t$ and ${\mathcal{F}}_t$, defined next.

Definition 1.

Define the sub-σ-algebras ${\tilde{\mathcal{F}}}_t:={\mathcal{F}}_{t-1}\cup \sigma ({\tilde{x}}_t^{\left(1\right)},\dots ,{\tilde{x}}_t^{\left(N\right)})$, ${\mathcal{F}}_t:={\tilde{\mathcal{F}}}_t\cup \sigma (x_t^{\left(1\right)},\dots ,x_t^{\left(N\right)})$ for $t=1,\dots ,T$, initialized by ${\mathcal{F}}_0:=\varnothing$.

At each point in time, we associate an empirical measure with the weighted particle sample generated by the propagation (1.) and reweighting (2.) steps in Algorithm 1. Formally, for $t=1,2,\dots ,T$, we define the empirical measure,

$$\begin{array}{cc}\hfill {\tilde{p}}_{\theta }^N(\mathrm{d}{x}_{1:t}\mid y_{1:t})& :=\sum _{i=1}^N{\overline{w}}_t^{\left(i\right)}{\delta }_{{\tilde{x}}_{1:t}^{\left(i\right)}}\left(\mathrm{d}{x}_{1:t}\right),\hfill \end{array}$$

(59)

where ${\delta }_{x^{\prime }}\left(\mathrm{d}x\right)$ denotes the point measure at $x^{\prime }\in {\mathcal{R}}^p$ with respect to $\mathrm{d}x$. The weighted particles that constitute the empirical measure (59) are approximately distributed according to the smoothing density (28). We emphasize the weighted particles are not independent draws from (28), because the resampling step introduces dependence between the particles at each iteration of the algorithm. We use the empirical measure (59) to define a particle filter-based approximation of the intractable smoothing problem in (27),

$$\begin{array}{cc}\hfill {\tilde{\mathbb{E}}}_{\theta }^N[{\gamma }_t\left(x_{1:t}\right)\mid y_{1:t}]& :=\int {\gamma }_t\left(x_{1:t}\right){\tilde{p}}_{\theta }^N(\mathrm{d}{x}_{1:t}\mid y_{1:t})=\sum _{i=1}^N{\overline{w}}_t^{\left(i\right)}{\gamma }_t\left({\tilde{x}}_{1:t}^{\left(i\right)}\right),\hfill \end{array}$$

(60)

for any point in time $t\in \{1,\dots ,T\}$. Due to dependence between the weighted particles, we cannot establish the asymptotic properties of the approximation (60) based on the law of large numbers and central limit theorem for independent random variables. For appropriately chosen test functions ${\gamma }_t\left(x_{1:t}\right)$, the approximation (60) is both consistent and asymptotically Gaussian as the number of particles tends to infinity, $N\to \infty$, cf. Theorem 9.4.5 in Cappé et al. (2005).

The particle filter in Algorithm 1 also produces an approximation of the log-likelihood function (8) evaluated at the parameter value $\theta$ and the observation sequence $y_{1:T}$,

$$\begin{array}{cc}\hfill {\tilde{\ell }}_T^N\left(\theta \right)& :=\sum _{t=1}^{T}log{W}_t^N.\hfill \end{array}$$

(61)

We note that the approximate log-likelihood function (61) consists of the logarithm of the product of normalizing constants produced by Algorithm 1. The approximate log-likelihood (61) is consistent in the sense that it converges in probability to the true log-likelihood function, as the number of particles tends to infinity, see Lemma 4.

Lemma 4.

For the model (1)–(3) and $\theta \in \mathsf{\Theta }$, the approximate log-likelihood function (61) produced by Algorithm 1 is a consistent estimator of the true log-likelihood (8),

$$\begin{array}{c}\hfill {\tilde{\ell }}_T^N\left(\theta \right)\stackrel{P}{\to }{\ell }_T\left(\theta \right),\end{array}$$

(62)

as $N\to \infty$.

In addition to producing an approximation of the intractable log-likelihood function (8), we apply the approximation (60) of the intractable smoothing problem in (27) to produce approximations of the sample score and observed Information matrix via Fisher’s and Louis’ identities in Lemma 2.

6.2. The Approximate Sample Score and Observed Information Matrix

We showed in Section 5 that the sample score and observed Information matrix can be expressed in terms of smoothing problems of the form (27). Appealing to Fisher’s identity (37) in Lemma 2, and to the approximation of the smoothing problem (60), we define the particle filter-based approximate sample score as,

$$\begin{array}{cc}\hfill {\tilde{S}}_T^N\left(\theta \right)& :=\sum _{i=1}^N{U}_T({\tilde{x}}_{1:T}^{\left(i\right)};\theta ){\overline{w}}_T^{\left(i\right)},\hfill \end{array}$$

(63)

for any parameter $\theta \in \mathsf{\Theta }$, with the function $U_T(x_{1:T};\theta )$ as defined in (33). If Conjecture 2 holds, then the approximate sample score in (63) is both consistent and asymptotically normal.

Lemma 5.

If Conjecture 2 holds and $\theta \in \mathsf{\Theta }$, then the approximate sample score (63) is asymptotically normal,

$$\begin{array}{c}\hfill \sqrt{N}\left\{{\tilde{S}}_T^N\left(\theta \right)-{\mathit{S}}_T\left(\theta \right)\right\}\stackrel{D}{\to }N(0,{\tilde{\mathbb{S}}}_T\left[U_T(x_{1:T};\theta )\right]),\end{array}$$

(64)

as $N\to \infty$. An intractable expression for the asymptotic covariance matrix ${\tilde{\mathbb{S}}}_T\left[U_T(x_{1:T};\theta )\right]$ is given in Lemma A.5 by setting $t=T$ and ${\gamma }_T\left(x_{1:T}\right)=U_T(x_{1:T};\theta )$.

Similarly, by appealing to Louis’ identity (38) in Lemma 2, and to the approximation of the smoothing problem (60), we define the particle filter-based approximate observed Information matrix as,

$$\begin{array}{cc}\hfill {\tilde{I}}_T^N\left(\theta \right)& :={\tilde{S}}_T^N\left(\theta \right){\tilde{S}}_T^N{\left(\theta \right)}^{\prime }-{\tilde{G}}_T^N\left(\theta \right)-{\tilde{K}}_T^N\left(\theta \right),\hfill \end{array}$$

(65)

for any parameter $\theta \in \mathsf{\Theta }$, where we define the approximations to (39) and (40) as

$$\begin{array}{cc}\hfill {\tilde{G}}_T^N\left(\theta \right)& :=\sum _{i=1}^N{V}_T({\tilde{x}}_{1:T}^{\left(i\right)};\theta ){\overline{w}}_T^{\left(i\right)}\hfill \end{array}$$

(66)

$$\begin{array}{cc}\hfill {\tilde{K}}_T^N\left(\theta \right)& :=\sum _{i=1}^N{U}_T({\tilde{x}}_{1:T}^{\left(i\right)};\theta )U_T{({\tilde{x}}_{1:T}^{\left(i\right)};\theta )}^{\prime }{\overline{w}}_T^{\left(i\right)},\hfill \end{array}$$

(67)

and the functions $U_T(x_{1:T};\theta )$ and $V_T(x_{1:T};\theta )$ are defined in (33) and (34), respectively. If Conjecture 2 holds, then the approximate observed Information in (65) is consistent, stated in the following lemma.

Lemma 6.

If Conjecture 2 holds and $\theta \in \mathsf{\Theta }$, then the approximate observed Information matrix (65) is consistent,

$$\begin{array}{c}\hfill {\tilde{I}}_T^N\left(\theta \right)\stackrel{P}{\to }{\mathit{I}}_T\left(\theta \right)\end{array}$$

(68)

as $N\to \infty$.

Both the approximate sample score (63) and observed Information matrix (65) are biased for finite N. This is a general issue related to the particle filter-based approximation of the smoothing problem (60). At each iteration, the particle filter in Algorithm 1 relies on an approximation of the normalized constant, i.e., likelihood contribution. This induces a finite-sample bias in (60) that gradually disappears as the number of particles N tends to infinity and is negligible for large enough N, cf. e.g., Robert and Casella (2010, sct. 3.3.2).

The particle filter-based approximation of the sample score (63) and observed Information matrix (65) correspond to a batch version of Algorithm A in Poyiadjis et al. (2011), which is of computational cost $O\left(N\right)$, but exhibits quadratically increasing variance of the approximate sample score as a function of the sample size T. We note that Poyiadjis et al. (2011) also suggest an alternative algorithm, that exhibits linearly increasing variance as a function of T, but at the computational cost $O\left(N^2\right)$. For smaller sample sizes, such as monthly observations as usually encountered in economics, we have found that the $O\left(N\right)$ algorithm is adequate.

7. Particle Filter-Based Estimation and Inference

In this section, we show how the approximate sample score (63) and observed Information matrix (65) can be used to perform parameter estimation and inference. We apply a stochastic approximation method based on the approximate sample score to approximate the ML estimator (19). This has recently been suggested in Poyiadjis et al. (2011). We then use the approximate observed Information matrix to obtain approximate standard errors for the approximate ML estimates. Although these quantities are ‘approximate’, we note that they can be made arbitrarily precise by increasing the number of particles, N, at the expense of increased computational effort.

Recall from Section 4 that the ML estimator (19) is doubly intractable. Consequently, we cannot apply gradient-based optimization algorithms to maximize the log-likelihood function (8). Originally proposed in Robbins and Monro (1951), stochastic approximation methods are conceptually similar to gradient-based optimization methods, but rely on noisy rather than exact evaluations of the sample score to optimize the objective function. The basic idea is that appropriately decreasing the step sizes provides an averaging of the random errors induced by the noisy evaluations of the sample score. For a book-length treatment of stochastic approximation, we refer to Kushner and Yin (2003).

The stochastic approximation algorithm proposed in Poyiadjis et al. (2011, sct. 3.1) consists of a recursion that is conceptually similar to the steepest descent method, cf. e.g., Nocedal and Wright (2006, chp. 3). Prior to executing the algorithm, we choose a fixed initial parameter value ${\theta }_0\in \mathsf{\Theta }$, a sequence of particle counts ${\{N_j\}}_{j=1}^{\infty }$, a sequence of step sizes ${\{{\gamma }_j\}}_{j=1}^{\infty }$, and a sequence of weight matrices ${\{B_j\}}_{j=1}^{\infty }$. The particle counts must be monotonically increasing positive integers, the step sizes must be strictly positive, non-summable but square summable,

$$\begin{array}{c}\hfill \sum _{j=1}^{\infty }{\gamma }_j=\infty \mathrm{and}\sum _{j=1}^{\infty }{\gamma }_i^2<\infty ,\end{array}$$

(69)

and the weight matrices must be positive definite. Having chosen the initial parameter, particle counts, step sizes, and weight matrices, we run the recursion,

$$\begin{array}{c}\hfill {\theta }_{j+1}={\theta }_j+{\gamma }_j{B}_j{\tilde{S}}_T^{N_j}\left({\theta }_j\right),\end{array}$$

(70)

for $j=0,1,\dots ,K$. Here K has to be sufficiently large in the sense that the sequence of parameter values generated by the recursion (70) has stabilized in a neighborhood of the true ML estimate. Additionally, if the particle count $N_j$ is large enough, the approximation error affecting the stochastic approximation recursion (70) will be approximately normal, cf. Lemma 5. In this case large disturbances will be rare, such that the parameter sequence ${\{{\theta }_j\}}_{j=1}^K$ is likely to stabilize without exhibiting large jumps.

We denote by ${\{{\tilde{x}}_{1:T,j}^{\left(i\right)},{\overline{w}}_{t,j}^{\left(i\right)}\}}_{i=1}^{N_j}$ the particle paths produced by the particle filter in Algorithm 1 at iteration j of the stochastic approximation recursion (70). The iteration index j is notationally identical to time index of the particle path, cf. Remark 7. Although this is abuse of notation, it is clear from the context whether we refer to the parameter iteration or particle path time index. The parameter ${\theta }_{j+1}$ produced by iteration j of (70) is a random variable that is measurable with respect to the sub-$\sigma$-algebra ${\mathcal{G}}_j$, defined next.

Definition 2.

Let ${\mathcal{F}}_{T,j}:=\sigma (x_{1:T,j}^{\left(1\right)},\dots ,x_{1:T,j}^{\left(N_j\right)})$ denote the sub-σ-algebra in Definition 1 generated with the parameter value ${\theta }_j$, and define the sub-σ-algebras ${\mathcal{G}}_j:={\mathcal{G}}_{j-1}\cup {\mathcal{F}}_{T,j}$ for $j=1,\dots$, initialized by ${\mathcal{G}}_0:={\mathcal{F}}_{T,0}$.

One of the main benefits of the stochastic approximation method is that the method is known to stabilize for a wide variety of initial values, sample counts, step sizes, and weight matrices. In practice, all of these choices affect the number of iterations needed to bring the parameter sequence into the neighborhood of the true ML estimator. The choice of step sizes is particularly important, since large step sizes generally speed up the convergence, but fail to dampen the approximation-induced noise. Small step sizes reduce the noise, but cause slow convergence. The particle count has a similar effect, since a low number of particles will result in a computationally cheap but noisy approximation of the sample score, while a large number of particles reduces the noise but increases the computational cost. Heuristically, it is appropriate to use a combination of large step sizes and small particle counts until the parameter sequence has reached a neighborhood of the ML estimator, and then switch to a combination of smaller step sizes and larger particle counts to reduce the noise. The intuition is that, while far away from the ML estimator, a relatively noisy approximation of the sample score will still on average lead the algorithm in the right direction.

The presence of noise in the sample score is not an impediment when applying stochastic approximation, since the use of decreasing step sizes provides an averaging of the errors. However, the finite sample bias of the particle filter-based approximate sample score, cf. Section 6.2, poses a problem since its effect is not mitigated by decreasing the step sizes. Bias reduction is possible by increasing the particle count $N_j$ together with the iteration number j.

The stochastic approximation method is presented in Algorithm 2 below.1

Algorithm 2: Stochastic Approximation.

Choose the initial parameter ${\theta }_0\in \mathsf{\Theta }$, the particle counts ${\left\{N_j\right\}}_{j=1}^{\infty }$, the step sizes ${\left\{{\gamma }_j\right\}}_{j=1}^{\infty }$ and weighting matrices ${\left\{B_j\right\}}_{j=1}^{\infty }$. For $j=0,1,\dots ,K$: Run Algorithm 1 for ${\theta }_j$ to generate $N_j$ weighted particle paths, denoted ${\{x_{1:T,j}^{\left(i\right)},{\overline{w}}_{t,j}^{\left(i\right)}\}}_{i=1}^{N_j}$. Compute the approximate sample score (63), denoted $$\begin{array}{cc}\hfill {\tilde{S}}_T^{N_j}\left({\theta }_j\right)& =\sum _{i=1}^{N_j}{U}_T(x_{1:T,j}^{\left(i\right)};{\theta }_j){\overline{w}}_{t,j}^{\left(i\right)}.\hfill \end{array}$$ (71) With step size ${\gamma }_j$, ascend along the direction $B_j$, $$\begin{array}{cc}\hfill {\theta }_{j+1}& ={\theta }_j+{\gamma }_j{B}_j{\tilde{S}}_T^{N_j}\left({\theta }_j\right).\hfill \end{array}$$ (72)

Polyak (1990) and Polyak and Juditsky (1992) showed that if the step sizes ${\{{\gamma }_j\}}_{j=1}^{\infty }$ satisfy the summability conditions (69) and tend to zero slower than $j^{-1}$, then the average of the last $j-K_0$ iterations converges at an optimal rate. Here $K_0<K$ denotes the iteration number at which the averaging begins; implicitly, we discard the initial $K_0$ iterations. We define the approximate ML estimator as,

$$\begin{array}{cc}\hfill {\tilde{\theta }}_T& :=\frac{1}{K-K_0}\sum _{j=K_0}^K{\theta }_j,\hfill \end{array}$$

(73)

suppressing the dependence on the particle count. Establishing convergence of the approximate ML estimator (73) to the true ML estimator (19) is outside the scope of this paper. However, if (73) converges in probability to (19) for any fixed T, then (73) inherits the consistency property, cf. Theorem 2, of the true ML estimator.

Convergence of the particle filter-based stochastic approximation method proposed in Poyiadjis et al. (2011) has, to the author’s knowledge, not been studied yet. The finite-sample bias of the approximate sample score (63) presents the primary obstacle to establishing convergence results. Intuition suggests that increasing the number of particles $N_j$ with the iteration number j solves the problem. However, convergence of such schemes has not been carefully established, cf. Douc et al. (2014, sct. 12.1.2). Poyiadjis et al. (2011) report stabilization of the particle filter-based stochastic approximation method with constant particle count. In Section 10, we report similar stabilization with increasing particle counts.

If the model is correctly specified, we would conduct inference on the ML estimator via the observed Information matrix, cf. Section 4. Analogously, since the approximate observed Information matrix (65) converges in probability to the true observed Information matrix (10), we can conduct inference for the approximate ML estimator (73) via the approximate observed Information matrix (65), the same way we would conduct inference given the true observed Information matrix (10).

8. Model Diagnostics

In this section, we introduce a method to conduct model diagnostics, such that we may assess whether the SSR model is well-specified for a given parameter $\theta$ and observation sequence $y_{1:T}$. Recall that the disturbances $u_t$, ${\eta }_t$ and ${\nu }_t$ are normally distributed and serially independent with mean zero and unit variances. Because the components ${\epsilon }_t$ and ${\xi }_t$ are hidden to us, we cannot directly compute the residuals corresponding to the disturbances. Instead, we introduce the normalized one-step prediction errors, cf. Durbin and Koopman (2012, sct. 2.12), that can be approximated via particle filtering. This approach to model diagnostics for state space models has also previously been considered in Pitt and Shephard (1999a).

We define the normalized one-step prediction errors as,

$$\begin{array}{cc}\hfill {\mathit{e}}_t& :=\mathbb{V}{\mathrm{ar}}_{\theta }{\left[y_t\mid y_{1:t-1}\right]}^{-1/2}\left(y_t-{\mathbb{E}}_{\theta }\left[y_t\mid y_{1:t-1}\right]\right),\hfill \end{array}$$

(74)

for $t=1,\dots ,T$. For a well-specified model, the sequence of normalized one-step prediction errors should be serially independent with mean zero with unit variance. Any deviation from these characteristics are indicative of model misspecification.

The conditional mean and variance in (74) can be stated in terms of smoothing problems, where the test functions are the conditional mean and variance of the predictive observation density,

$$\begin{array}{c}{\mathbb{E}}_{\theta }\left[y_t\mid y_{1:t-1}\right]={\mathbb{E}}_{\theta }\left[{\mathbb{E}}_{\theta }\left[y_t\mid x_{t-1}\right]\mid y_{1:t-1}\right],\end{array}$$

(75)

$$\begin{array}{c}\mathbb{V}{\mathrm{ar}}_{\theta }\left[y_t\mid y_{1:t-1}\right]={\mathbb{E}}_{\theta }\left[\mathbb{V}{\mathrm{ar}}_{\theta }\left[y_t\mid x_{t-1}\right]\mid y_{1:t-1}\right]+\mathbb{V}{\mathrm{ar}}_{\theta }\left[{\mathbb{E}}_{\theta }\left[y_t\mid x_{t-1}\right]\mid y_{1:t-1}\right].\end{array}$$

(76)

We note that the conditional mean and variance of the predictive observation density are given in Lemma 3. Using the locally optimal particle filter in Algorithm 1, we define approximations to (75) and (76) as

$$\begin{array}{c}{\tilde{\mathbb{E}}}_{\theta }^N\left[y_t\mid y_{1:t-1}\right]:=\sum _{i=1}^N{\tilde{\mu }}_{t|t-1}^{y,\left(i\right)}{\overline{w}}_{t-1}^{\left(i\right)}\end{array}$$

(77)

$$\begin{array}{c}{\widetilde{\mathbb{V}\mathrm{ar}}}_{\theta }^N\left[y_t\mid y_{1:t-1}\right]:=\sum _{i=1}^N{\tilde{\Sigma }}_{t|t-1}^{y,\left(i\right)}{\overline{w}}_{t-1}^{\left(i\right)}+\sum _{i=1}^N\left({\tilde{\mu }}_{t|t-1}^{y,\left(i\right)}\right){\left({\tilde{\mu }}_{t|t-1}^{y,\left(i\right)}\right)}^{\prime }{\overline{w}}_{t-1}^{\left(i\right)}-{\tilde{\mathbb{E}}}_{\theta }^N\left[y_t\mid y_{1:t-1}\right]{\tilde{\mathbb{E}}}_{\theta }^N{\left[y_t\mid y_{1:t-1}\right]}^{\prime },\end{array}$$

(78)

respectively, where we have defined the conditional moments given each individual particle as,

$$\begin{array}{cc}\hfill {\tilde{\mu }}_{t|t-1}^{y,\left(i\right)}& :={\mathbb{E}}_{\theta }[y_t\mid {\tilde{x}}_{t-1}^{\left(i\right)}]\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill {\tilde{\Sigma }}_{t|t-1}^{y,\left(i\right)}& :=\mathbb{V}{\mathrm{ar}}_{\theta }[y_t\mid {\tilde{x}}_{t-1}^{\left(i\right)}],\hfill \end{array}$$

(80)

for $i=1,\dots ,N$. Finally, we use the approximations (77) and (78) to define the approximate normalized likelihood contributions as follows,

$$\begin{array}{cc}\hfill {\tilde{e}}_t^N& :={\widetilde{\mathbb{V}\mathrm{ar}}}_{\theta }^N{\left[y_t\mid y_{1:t-1}\right]}^{-1/2}\left(y_t-{\tilde{\mathbb{E}}}_{\theta }^N\left[y_t\mid y_{1:t-1}\right]\right),\hfill \end{array}$$

(81)

for $t=1,\dots ,T$. Thus, by applying the particle filter in Algorithm 1, we obtain the sequence of approximate normalized one-step prediction errors ${\tilde{e}}_{1:T}$ via (77)–(81). For N sufficiently large, we can use the sequence ${\tilde{e}}_{1:T}$ to test whether the true sequence of normalized one-step prediction errors ${\mathit{e}}_{1:T}$ is serially independent with mean zero and unit variance. For common tests for serial dependence and ARCH effects see e.g., Doornik and Hendry (2013, sct. 11.9.2–3).

9. Simulation Study

In this section, we conduct a simulation study of the asymptotic properties of the ML estimator, stated in Theorem 2. We limit our treatment to B, A, $\mathsf{\Phi }$ and ${\mathsf{\Omega }}_{\mathsf{\Phi }}$, leaving aside the remaining parameters ${\mathsf{\Omega }}_u$, $\mu$, ${\mathsf{\Omega }}_{\eta }$, ${\mathsf{\Omega }}_{\nu }$ and ${\mathsf{\Omega }}_{\eta ,\nu }$. Recall, the loading matrix for the stationary components A is conjectured to be asymptotically normal, while the loading matrix of the nonstationary components B is kept fixed. Due to the results of Chang et al. (2009), we expect the asymptotic distribution of B to be mixed normal, and we tentatively investigate this. Moreover, we consider the case where ${\mathsf{\Phi }}_t$ is a stochastic unit root. A deterministic unit root is associated with the Dickey-Fuller distribution, cf. Dickey and Fuller (1979), while a stochastic unit root has been shown to be asymptotically normal, see e.g., Ling (2007) and Bohn Nielsen and Rahbek (2014).

Recall, Theorem 2 is based on the conjectured properties of the true, intractable log-likelihood function and its derivatives, cf. Conjecture 1. The aim is to substantiate this conjecture by obtaining the distribution of the approximate ML estimator based on simulated data sets. Usually, the number of realizations in a simulation study of this type is in excess of 1000 and the sample length in excess of 2500 observations. Due to the computational intensity of the particle filter-based stochastic approximation method in Algorithm 2, we limit ourselves to 250 realizations and 500 observations.

We let each of the simulated data sets be a bivariate $p=2$ series of length $T=500$ observations with $r=1$ stationary component and $p-r=1$ nonstationary component. We use the parameter

$$\begin{array}{c}B=\left[\begin{array}{c}1\\ 1\end{array}\right],A=\left[\begin{array}{c}0\\ 1\end{array}\right],{\mathsf{\Omega }}_u=\left[\begin{array}{cc}{2.5}^2& 0\\ 0& {2.5}^2\end{array}\right],\end{array}$$

(82)

$$\begin{array}{c}\mu =0,\varphi =1,{\omega }_{\varphi }^2={0.25}^2,{\omega }_{\eta }^2={15}^2,{\omega }_{\eta ,\nu }=0,\mathrm{and}{\omega }_{\nu }^2={2.5}^2,\end{array}$$

(83)

to generate the simulated data sets. We note the parameter values (83) result in a top Lyapunov coefficient of ${\gamma }_n=-0.035$, computed via (15) with $n={10}^6$, such that the RCAR process ${\left\{{\xi }_t\right\}}_{t=0,1,\dots }$ is strictly stationary.

Having simulated 250 series with the data generating process given by (1)–(3) and (83), we apply Algorithm 2 with $K=600$ iterations to obtain the approximate ML estimate for the parameter in question, e.g., $\varphi$, keeping all other parameters fixed at the true values in (83). We initialize the algorithm at the true parameter value, and initiate Polyak averaging at iteration $K_0=100$.2 Moreover, we let the particle count increase as

$$\begin{array}{c}\hfill N_j=50+\lfloor 1/20j\rfloor ,\end{array}$$

(84)

where $\lfloor ·\rfloor$ denotes the largest integer that is smaller than the argument. We let the step size sequence to decrease as

$$\begin{array}{c}\hfill {\gamma }_j=100{(j+500)}^{-2/3},\end{array}$$

(85)

and set the weight matrix to

$$\begin{array}{c}\hfill B_j=T^{-1}\mathrm{diag}\left(\left[\begin{array}{cccccccccc}{10}^{-5}& 1& 1& 1& 1& {10}^{-2}& 1& 1& 1& {10}^{-3}\end{array}\right]\right),\end{array}$$

(86)

for $j=1,2,\dots ,K$. Note the particle count (84) tends to infinity as $j\to \infty$, eliminating the finite-sample bias of (63)–(65), the step sizes satisfy (69), and the weight matrix is constant.3

The results from the simulation experiment are presented in Figure 1. Despite the relatively low number of realizations and observations, Figure 1 is instructive of the asymptotic distributions of $A_1$, $\varphi$ and ${\omega }_{\varphi }^2$, cf. Panels (a), (c) and (d). These all appear to be normal. Recall, Theorem 2 does not state the asymptotic distribution of the ML estimator for $B_2$, and from Panel (b) it does not appear to be normal. Rather, the realizations in Panel (b) are consistent with mixed normality, as we would expect from the closely-related CST model, cf. Chang et al. (2009). To investigate further, one could to simulate the t-ratios of $B_2$, which should be standard normal. This involves the approximation of the observed Information matrix for each realization, which further increases the computational cost. For this reason, and because we consider B fixed, we do not pursue this further here.

In summary, the findings of the simulation study tentatively support the conjecture made in Section 4. Namely, the ML estimator for A, $\mathsf{\Phi }$ and ${\mathsf{\Omega }}_{\mathsf{\Phi }}$ is asymptotically normal. The ML estimator for $B_2$ appears to be consistent with mixed normality. We have not investigated the remaining parameters.

10. An Illustration

In this section, we illustrate the use of the SSR model by applying it to the monthly 10-year government bond rates for Germany and Greece from January 1999 to February 2018.4 We denote the German and Greek bond rates $y^{GE}$ and $y^{GR}$, respectively, and measure these in basis points per year. The sample begins at the introduction of the euro area and ends at present day. During this period, the rates initially exhibit convergence towards a common ‘euro area rate’, until interrupted by the euro area crisis beginning in 2009 and culminating in 2011. The rates, the spread and the changes in the spread are illustrated in Figure 2 below. Because the spread is up to 75 times larger during the second half of the sample than during the first half, we split the display of the sample into the first and second half, respectively.

Panels (a) and (b) in Figure 2 show the bond rates, Panels (c) and (d) show the spread, and Panels (e) and (f) show the changes in the spread in the two periods. We note two features of the observations. First, Panel (a) suggests the rates can be characterized by a shared common stochastic trend, since these tend to move in tandem. Second, Panels (d) and (f) suggest the spread can be characterized by a RCAR process, since the changes in the spread, cf. Panel (f), are clearly positively associated with the level of the spread itself, cf. Panel (d).

We define the observation vector as $y_t:={\left[\begin{array}{cc}{y}_t^{GE}& y_t^{GR}\end{array}\right]}^{\prime }$. We condition on the observation for January 1999, which we denote $y_0$, such that the effective sample spans $t=1,\dots ,229$. From visual inspection of Figure 2, our working assumption is that the spread $y_t^{GR}-y_t^{GE}$ is strictly stationary, while the rates $y_t$ share a common stochastic trend. With a $p=2$ dimensional system, we thus have $r=1$ stationary component and $p-r=1$ nonstationary component. Moreover, we fix $B={\left[\begin{array}{cc}1& 1\end{array}\right]}^{\prime }$, such that the orthogonal complement $b={\left[\begin{array}{cc}-1& 1\end{array}\right]}^{\prime }$ produces the spread. To ensure the model is just-identified, we normalize on the second element of A, such that $A_2=1$.

We apply the particle filter-based stochastic approximation method in Algorithm 2 to obtain the approximate ML estimate of the model parameter $\theta$. For this illustration, we run the algorithm for $K=\mathrm{10,000}$ iterations. We let the particle count increase as (84), the step size sequence decrease as (85), and the weighting matrix as (86). We initiate Polyak averaging at iteration $K_0=5000$.

Figure 3 shows the results of running the particle filter-based stochastic approximation method. Panel (a) displays the iterations for the parameters in the observation Equation (22), Panel (b) displays the iterations for the parameters in the transition Equation (23), and Panel (c) displays the sequence of realized approximate log-likelihoods together with a moving average of lag order 500. The algorithm has been implemented in the Ox 7 programming language, cf. Doornik (2012), using analytical derivatives of the complete data log-likelihood (32) for the evaluation of the function (33). The elements of the parameter sequence shown in Panels (a) and (b) have stabilized after the initial 7500 iterations. At the $\mathrm{10,000}$th iteration, the particle count has increased to 550, the step size decreased to $0.2085$, and the sequences have stabilized. By inspection of the sequence of the approximate log-likelihood in Panel (c), we see that the value has also stabilized after approximately 7500 iterations.

The estimation results are presented in

Table 1. Approximate ML estimate, ${\tilde{\theta }}_T$.

Parameter	Estimate	Std.err.	Parameter	Estimate	Std.err.
$B_1$	$1.0000$	-	$\mu$	$0.3449$	$0.5526$
$B_2$	$1.0000$	-	$\varphi$	$1.0085$	$0.0152$
$A_1$	$-0.0154$	$3.4\times {10}^{-5}$	${\omega }_{\varphi }^2$	$0.0306$	$0.0031$
$A_2$	$1.0000$	-	${\omega }_{\eta }^2$	$360.1600$	$33.6250$
${\omega }_{u11}^2$	$2.1063$	$0.1974$	${\omega }_{\eta ,\nu }$	$-22.2400$	$0.8119$
${\omega }_{u12}$	$3.5924$	$0.3435$	${\omega }_{\nu }^2$	$1.7880$	$2.0728$
${\omega }_{u22}^2$	$6.6327$	$0.6214$

Note: The approximate log-likelihood is ${\tilde{\ell }}_T=-2094.1$. The approximate ML estimate has been obtained by running Algorithm 2 for $K=\mathrm{10,000}$ iterations with the particle count increasing to $N=550$ particles, as described in the main text. The standard errors are based on the inverse of the approximate observed Information matrix computed with $N=1000$ particles.

, together with approximate classic standard errors.5 Before considering inference, we assess the model fit. We compute the normalized one-step prediction errors ${\tilde{e}}_{1:T}^N$ via (81) using $N=1000$ particles.

Table 2. Model diagnostics.

Univariate tests for AR 1-2:	${\tilde{\mathit{e}}}_{\mathit{t},1}$	$\mathit{F}\left(2227\right)=1.4523$	$\mathit{p}=0.2362$
	${\tilde{\mathit{e}}}_{\mathit{t},2}$	$\mathit{F}\left(2227\right)=1.2086$	$\mathit{p}=0.3005$
Multivariate test for AR 1-2:		$F\left(8448\right)=1.6084$	$p=0.1200$
Univariate tests for ARCH:	${\tilde{e}}_{t,1}$	$F\left(1227\right)=4.7008$	$p=0.0312$
	${\tilde{e}}_{t,2}$	$F\left(1227\right)=0.58861$	$p=0.4438$

Note: The approximate normalized one-step prediction errors ${\tilde{e}}_{1:T}$ have been computed with $N=1000$ particles for the approximate ML estimate ${\tilde{\theta }}_T$, cf. Table 1.

presents univariate tests for autocorrelation (AR) of order one and two, autoregressive conditional heteroskedasticity (ARCH) of order one, and a multivariate test for AR of order one and two, cf. Doornik and Hendry (2013, sct. 11.9.2–3). We cannot reject the null hypothesis of no-AR of order one and two in the univariate as well as multivariate tests at a $5\%$ critical level. Nor can we reject the null hypothesis of no-ARCH for the residuals at a $5\%$ critical level. However, we note the test for the German rate is close to, but below, our chosen critical level. This could suggest unmodeled heteroskedasticity in the German bond rate. In conclusion, the overall specification of the model is acceptable. Moreover, computing the top Lyapunov coefficient via (15) with $n={10}^5$ produces a coefficient of ${\widehat{\gamma }}_n=-0.007$, which indicates the stationary direction is strictly stationary for ${\tilde{\theta }}_T$.

The model is reasonably well-specified, and we therefore proceed to use the approximate classic standard errors to conduct inference on the approximate ML estimates. First, we note the standard error of the estimate of $A_1$ is extremely small. Since the test for no-ARCH for the residuals associated with the German rate is rejected at the $5\%$ critical level, this could affect the approximate classic standard errors.6 Nevertheless, it is economically plausible that the stationary component also loads into the German rate, given that a large increase in the Greek rate would in this case coincide with a small drop in the German rate, which is consistent with risk-averse investors seeking safer assets in times of uncertainty, such as the euro area crisis. Second, we cannot reject the null hypothesis that $H_0:\varphi =1$ at a $5\%$ critical level with $p=0.577$. Third, the estimate of ${\omega }_{\varphi }^2$ is significantly different from zero at any commonly used critical level. However, the constant term $\mu$ is not significantly different from zero with $p=0.533$. Fourth, the measurement errors are highly positively correlated with coefficient $0.961$, and the innovations of the unobserved components are highly negatively correlated with coefficient $-0.876$. The results in Table 1 suggest the level of the stationary direction is a stochastic unit root process without a constant term. An approximate likelihood ratio test for the joint null hypothesis $H_0:\varphi =1,\mu =0$ fails to reject the null at a $5\%$ critical level with $p=0.374$.

Based on the estimates in Table 1, we use the orthogonal complements b and a to compute the changes of the nonstationary and stationary components, given by $b^{\prime }\Delta y_t$ and $a^{\prime }\Delta y_t$, respectively. These are illustrated in Figure 4. First, we note the magnitude of the changes in Panels (a) and (b) of Figure 4 are slightly larger during the second half of the sample than during the first (standard deviations $18.01$ and $20.16$, respectively). Otherwise, the series in Panels (a) and (b) in Figure 4 are consistent with a homoskedastic random walk plus measurement error, cf. (17). The magnitude of the changes in Panels (c) and (d) of Figure 4 is positively associated with the level, just as observed in Figure 2. This is consistent with a random coefficient autoregressive process plus measurement error, cf. (18).

Summarizing, the empirical illustration suggests that the SSR model successfully characterizes the 10-year government bond rates for Germany and Greece during the period from January 1999 to February 2018. During this sample, the spread exhibits bubble-like behavior, which is captured by the random coefficient autoregressive dynamics of the stationary component. Additionally, the levels exhibit a shared common stochastic trend, which is captured by the random walk dynamics of the nonstationary component.

11. Conclusions

In this paper, we have proposed and studied the stochastic stationary root model, which is a multivariate nonlinear state space model. We introduced particle filter-based approximations of the intractable log-likelihood function, sample score and observed Information matrix. In turn, we used these to approximate the ML estimator via stochastic approximation, and showed how to perform inference via the approximate observed Information matrix. We considered model diagnostics to assess the model fit. Additionally, we conducted a simulation study to investigate the asymptotic properties of the ML estimator. Finally, we presented an empirical application to the 10-year government bond rates in Germany and Greece in the period from January 1999 to February 2018 to illustrate the usefulness of the SSR model.