2.1 The Model, Assumptions, and Notation

Consider the following MS state-space model:

(1) $$ \begin{align} \text{measurement equation: } & y_t= \beta_{s_t} x_t+\sigma_{s_t}u_t, \\[-36pt]\nonumber\end{align} $$

(2) $$ \begin{align} \text{transition equation: } & x_t = \alpha_{s_t} x_{t-1}+\varsigma_{s_t}v_t, \end{align} $$

where the transition equation describes the evolution of a latent variable $x_t$ , and the measurement equation describes how the observed variable $y_t$ responds to $x_t$ . The disturbances, $v_t$ and $u_t$ , are mutually and serially independent and follow standard normal distributions. The parameters $\beta _{s_t}$ , $\sigma _{s_t}$ , $\alpha _{s_t}$ , and $\varsigma _{s_t}$ depend on a latent regime variable $s_t$ . We assume that $s_t$ follows a first-order Markov chain on $\mathbb {S}=\{1,2,\ldots ,J\}$ with transition matrix $(q_{ij})_{1\leq i,j\leq J}$ and $q_{iJ}=1-\sum _{j=1}^{J-1}q_{ij}$ , $i=1,2,\ldots ,J$ . Moreover, we assume $\sigma _j,\varsigma _j>0$ for all $j=1,\ldots ,J$ . The unknown parameter is $\theta = \big ( (\beta _i,\sigma _i,\alpha _i,\varsigma _i)_{i\in \mathbb {S}}, (q_{ij})_{1\le i\le J,\,1\le j\le J-1} \big )$ , and $\theta _0$ denotes the true parameter. We write $\mathbb {E}_\theta $ , $\mathbb {P}_\theta $ , and $p_\theta $ when the parameter is $\theta $ , and $\mathbb {E}_0$ , $\mathbb {P}_0$ , and $p_0$ when the parameter is $\theta _0$ .

In the special case of $J=1$ , the model reduces to a state-space model. The Kalman filter can be applied to compute the likelihood, because the conditional distribution in (3) is Gaussian,

(3) $$ \begin{align} (y_t,x_t)' \big| (y_{t-1},y_{t-2},\ldots,y_1,x_0)=_d \text{normal distribution}, \end{align} $$

and the linear projections in the Kalman filter coincide with the mean and variance of (3). When $J\geq 2$ , Gaussianity holds only if the distribution is also conditioned on the entire path of the regime. The Kalman filter needs to be applied to all possible regime paths to deliver the likelihood function. Then the problem of path dependence arises, which is explained in detail in Section 2.2.

2.2 Path Dependence in the Application of the Kalman Filter

$\mathbf {y}_1^t$ and $x_t$ can be written as a linear combination of the Gaussian shocks $\mathbf {u}_1^t$ and $\mathbf {v}_1^t$ with coefficients that depend on the regimes:

$$ \begin{align*} \begin{pmatrix} y_1\\ \vdots\\y_t\\x_t \end{pmatrix} =& \begin{pmatrix} \beta_{s_1}\varsigma_{s_1}&0&\dots&0\\ & &\ddots& \\ \beta_{s_t}\prod_{k=2}^t\alpha_{s_k}\varsigma_{s_1}&\beta_{s_t}\prod_{k=3}^{t}\alpha_{s_k}\varsigma_{s_2} & \dots & \beta_{s_t}\varsigma_{s_t}\\ \prod_{k=2}^t\alpha_{s_k}\varsigma_{s_1}&\prod_{k=3}^{t}\alpha_{s_k}\varsigma_{s_2} & \dots & \varsigma_{s_t} \end{pmatrix} \begin{pmatrix} v_1\\ \vdots\\ v_t \end{pmatrix} + \begin{pmatrix} \sigma_{s_1}&\dots&0\\ &\ddots& \\ 0 &\dots& \sigma_{s_t}\\ 0&\dots & 0 \end{pmatrix} \begin{pmatrix} u_1\\ \vdots\\u_t \end{pmatrix}\\ & +\begin{pmatrix} \beta_{s_1}\alpha_{s_1}& \dots & \beta_{s_t}\prod_{k=1}^t\alpha_{s_k}& \prod_{k=1}^t\alpha_{s_k} \end{pmatrix}'x_0. \end{align*} $$

By the property of the Gaussian distribution,

(4) $$ \begin{align} (y_t, x_t)' \big| (\mathbf{y}_1^{t-1},\mathbf{s}_1^t,x_0)=_d \text{normal distribution}. \end{align} $$

If the Kalman filter is applied, it must be computed for each entire regime path $\mathbf {s}_1^t$ in order to derive the mean and variance in (4). If we truncate the regime path and only consider the most recent r periods, then the distribution of $(x_t,y_t)' \mid (\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t,x_0)$ becomes a Gaussian mixture, rather than Gaussian. In this case, the outputs of the Kalman filter no longer coincide with the mean and variance of the Gaussian mixture distribution, and the exact distribution cannot be recovered. Hence, to obtain the exact distribution, one must condition on the entire regime path.

Tracking all regime paths is computationally infeasible owing to the curse of dimensionality. More specifically, the Kalman filter given path $\mathbf {s}_1^t$ is

(5) $$ \begin{align} \text{Prediction: } \qquad x_{t|t-1}^\theta(\mathbf{s}_1^t) &= \alpha_{s_t} x^\theta_{t-1|t-1}(\mathbf{s}_1^{t-1}), \end{align} $$

(6) $$ \begin{align} y_{t|t-1}^\theta(\mathbf{s}_1^t) &= \beta_{s_t} x^\theta_{t|t-1}(\mathbf{s}_1^t), \end{align} $$

(7) $$ \begin{align} \Omega_{t|t-1}^\theta(\mathbf{s}_1^t) &= \alpha_{s_t}^2\Omega^\theta_{t-1|t-1}(\mathbf{s}_1^{t-1})+\varsigma_{s_t}^2, \end{align} $$

(8) $$ \begin{align} \Sigma_{t|t-1}^\theta(\mathbf{s}_1^t) &= \beta_{s_t}^2\Omega^\theta_{t|t-1}(\mathbf{s}_1^t)+\sigma_{s_t}^2, \end{align} $$

(9) $$ \begin{align} \text{Update:} \qquad \qquad x_{t|t}^\theta(\mathbf{s}_1^t) &= x^\theta_{t|t-1}(\mathbf{s}_1^t)+\beta_{s_t}\Omega^\theta_{t|t-1}(\mathbf{s}_1^t)(\Sigma^\theta_{t|t-1}(\mathbf{s}_1^t))^{-1}(y_t-y^\theta_{t|t-1}(\mathbf{s}_1^t)), \end{align} $$

(10) $$ \begin{align} \Omega_{t|t}^\theta(\mathbf{s}_1^t) &= \Omega^\theta_{t|t-1}(\mathbf{s}_1^t)-\beta_{s_t}^2 \Omega^\theta_{t|t-1}(\mathbf{s}_1^t)^2(\Sigma^\theta_{t|t-1}(\mathbf{s}_1^t))^{-1}, \end{align} $$

$t=1,2,\ldots $ , where $x_{t|k}^\theta (\mathbf {s}_1^t)$ (resp. $y_{t|k}^\theta (\mathbf {s}_1^t)$ ) denotes the linear projection of $x_t$ (resp. $y_t$ ) given the observation history $\mathbf {y}_1^k$ and the regime history $\mathbf {s}_1^t$ , and $\Omega _{t|k}^\theta (\mathbf {s}_1^t)$ (resp. $\Sigma _{t|k}^\theta (\mathbf {s}_1^t)$ ) denotes the mean squared error of $x_t$ (resp. $y_t$ ) given the observation history $\mathbf {y}_1^k$ and the regime history $\mathbf {s}_1^t$ . Then (5)–(10) need to be computed for all possible regime paths $\mathbf {s}_1^t$ , the number of which is $J^t$ . Because the number of paths increases geometrically, the computation quickly becomes practically infeasible.

It is worth pointing out that path dependence only arises when the Kalman filter is applied to compute the likelihood. It does not arise if we use simulation-based filters, such as particle filters. For simulation-based filters, the period likelihood $p_\theta (y_t|\mathbf {y}_1^{t-1})$ can be sequentially computed from

(11) $$ \begin{align} p_{\theta}(w_t|\mathbf{y}_1^{t-1}) &=\int p_{\theta}(w_{t}|w_{t-1})\mathbb{P}_{\theta}(dw_{t-1}|\mathbf{y}_1^{t-1}),\nonumber \\ p_{\theta}(y_t|\mathbf{y}_1^{t-1}) &=\int p_{\theta}(y_t|w_t)\mathbb{P}_\theta(dw_t|\mathbf{y}_1^{t-1}),\nonumber\\ p_{\theta}(w_t|\mathbf{y}_1^{t}) &=\frac{p_\theta(y_t|w_t)p_\theta(w_t|\mathbf{y}_1^{t-1})}{p_\theta(y_t|\mathbf{y}_1^{t-1})}, \end{align} $$

where $w_t=(x_t,s_t)'$ . The integral can be dealt with numerically in the simulation-based filters. At time t, information about previous unobservable variables, $w_{t-1}$ , is integrated out in (11), and thus we do not need to keep track of past regimes. The computation of simulation-based filters, though feasible, is still expensive because of the random sampling procedure at each sample and each iteration. In order to reduce computation time and still enjoy the simplicity of the Kalman filter, many empirical works apply an approximation strategy, which is explained in Section 2.3.

2.3 The Approximation Algorithm

This section explains the approximation algorithm in Kim (Reference Kim1994) and our collapsing method to average out the regimes. During the application of the Kalman filter, the approximation algorithm truncates the regime path and tracks only the most recent r periods of regimes, $\mathbf {s}_{t-r+1}^t$ , for $t\geq r$ . The likelihood function computed with this algorithm is approximate. We use $\overline {p}_r^\theta $ to denote the approximate likelihood function, with the subscript r emphasizing that it depends on the regime path length r.

The Kalman filter given the most recent r periods of regimes can be described as in (12)–(17).

(12) $$ \begin{align} \text{Prediction: } \qquad \overline{x}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) &= \alpha_{s_t} \overline{x}^\theta_{t-1|t-1}(\mathbf{s}_{t-r+1}^{t-1}), \end{align} $$

(13) $$ \begin{align} \overline{y}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) &= \beta_{s_t} \overline{x}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t), \end{align} $$

(14) $$ \begin{align} \overline{\Omega}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) &= \alpha_{s_t}^2\overline{\Omega}^\theta_{t-1|t-1}(\mathbf{s}_{t-r+1}^{t-1}) +\varsigma_{s_t}^2, \end{align} $$

(15) $$ \begin{align} \overline{\Sigma}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) &= \beta_{s_t}^2\overline{\Omega}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) +\sigma_{s_t}^2, \end{align} $$

(16) $$ \begin{align} \text{Update: } \qquad \qquad \overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+1}^t) &= {} \overline{x}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) \nonumber \\ & + \beta_{s_t} \overline{\Omega}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) \overline{\Sigma}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t)^{-1} \big(y_t-\overline{y}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t)\big), \end{align} $$

(17) $$ \begin{align} \overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+1}^t) &= \overline{\Omega}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t) -\beta_{s_t}^2 \overline{\Omega}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t)^2 \overline{\Sigma}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t)^{-1}. \end{align} $$

We put lines over the variables in (12)–(17) to distinguish them from the ones used in the exact algorithm (5)–(10). After the iteration at time t, $\overline {x}^\theta _{t|t}(\mathbf {s}_{t-r+1}^t)$ and $\overline {\Omega }^\theta _{t|t}(\mathbf {s}_{t-r+1}^t)$ depend on r periods of regimes. If we directly go to the filter at time $t+1$ , then we will include one more regime $s_{t+1}$ and track $r+1$ regimes. If we keep doing it, the number of regimes increases geometrically. To avoid the increase, we need to average out, or collapse, the most dated regime. Then the number of periods of regimes is $r-1$ at the end of iterations at time t and stays at r at time $t+1$ . Kim (Reference Kim1994) adopts the following collapsing procedures:

(18) $$ \begin{align} &\begin{aligned} \overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t}) =\sum_{s_{t-r+1}}\overline{p}^{\theta}_r(s_{t-r+1}|\mathbf{s}_{t-r+2}^{t},\mathbf{y}_1^{t})&\overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+1}^t), \end{aligned} \end{align} $$

(19) $$ \begin{align} &\begin{aligned} \overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t}) =\sum_{s_{t-r+1}}\overline{p}^{\theta}_r(s_{t-r+1}|\mathbf{s}_{t-r+2}^{t},\mathbf{y}_1^{t}) &\Big[\overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+1}^t)\nonumber \\ &\quad+\big(\overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-\overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+1}^t)\big)^2\Big]. \end{aligned} \end{align} $$

We propose the following new method to collapse $\overline {\Omega }^\theta _{t|t}(\mathbf {s}_{t-r+2}^{t})$ by discarding the nonlinear adjustment term in (19):

(20) $$ \begin{align} \overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t}) =\sum_{s_{t-r+1}}p_r^\theta(s_{t-r+1}|\mathbf{s}_{t-r+2}^{t},\mathbf{y}_1^{t})\overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+1}^t). \end{align} $$

To see more clearly how the two collapsing rules differ, note that their gap is

$$\begin{align*}\sum_{s_{t-r+1}}\overline{p}_r^\theta(s_{t-r+1}\mid\mathbf{s}_{t-r+2}^t,\mathbf{y}_{1}^t) \big(\overline{x}_{t|t}^\theta(\mathbf{s}_{t-r+2}^t)-\overline{x}_{t|t}^\theta(\mathbf{s}_{t-r+1}^t)\big)^2, \end{align*}$$

which measures the dispersion of $\overline {x}_{t|t}^\theta (\mathbf {s}_{t-r+1}^t)$ across different values of $s_{t-r+1}$ . The gap is small when $\max _{i,j\in \mathbb {S}}|\overline {x}_{t|t}^\theta (\mathbf {s}_{t-r+2}^t,s_{t-r+1}=i)-\overline {x}_{t|t}^\theta (\mathbf {s}_{t-r+2}^t,s_{t-r+1}=j)|$ is small. As shown in the proof of Theorem 11, in the simple collapsing method, this per-period gap is small when $\max _{j\in \mathbb {S}}\alpha _j$ is small and r is large. This suggests that the gap between the two collapsing rules may likewise depend on $\max _{j\in \mathbb {S}}\alpha _j$ and r. Moreover, because the gap occurs in every period, the cumulative difference over the sample of size n increases with n.

Kim (Reference Kim1994) includes the nonlinear term based on the argument that if $\overline {x}^\theta _{t|t}(\mathbf {s}_{t-r+1}^t)$ and $\overline {\Omega }^\theta _{t|t}(\mathbf {s}_{t-r+1}^t)$ represented the exact conditional mean and variance, respectively, then (19) would give the exact $\mathrm {Var}(x_t|\mathbf {y}_1^{t},\mathbf {s}_{t-r+2}^{t})$ . $\overline {x}^\theta _{t|t}(\mathbf {s}_{t-r+1}^t)$ and $\overline {\Omega }^\theta _{t|t}(\mathbf {s}_{t-r+1}^t)$ are, however, approximations instead of the exact conditional moments. There is no theoretical support that (19) must yield a smaller error than the simple collapsing (20). In this study, we focus on the simple collapsing (20) so that we can focus on the real theoretical challenges without the nonlinear part in (19) that complicates the proof.

The predictive density given r periods of regimes is imposed to be

(21) $$ \begin{align} \overline{p}_r^\theta(y_t|\mathbf{y}_1^{t-1},\mathbf{s}_{t-r+1}^t) =\frac{1}{(2\pi\overline{\Sigma}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t))^{\frac{1}{2}}} \exp\bigg\{\!\!\!-\frac{(y_t-\overline{y}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t))^2}{2\overline{\Sigma}^\theta_{t|t-1}(\mathbf{s}_{t-r+1}^t)}\bigg\}. \end{align} $$

$\overline {p}_r^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ approximates the exact predictive density in two aspects. First, the exact distribution of $y_t$ conditional on $(\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ is a Gaussian mixture, while we impose the Gaussian density function in (21). Moreover, owing to the non-Gaussianity of $x_t|(\mathbf {y}_1^{t},\mathbf {s}_{t-r+1}^t)$ , the outputs $\overline {y}^\theta _{t|t-1}(\mathbf {s}_{t-r+1}^t)$ and $\overline {\Omega }^\theta _{t|t-1}(\mathbf {s}_{t-r+1}^t)$ do not give the exact conditional mean and variance. We put $\overline {y}^\theta _{t|t-1}(\mathbf {s}_{t-r+1}^t)$ and $\overline {\Omega }^\theta _{t|t-1}(\mathbf {s}_{t-r+1}^t)$ in (21) in place of the true conditional mean and variance.

Now that $\overline {p}_r^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ and $p_\theta (s_t|s_{t-1})$ are ready, the algorithm in Hamilton (Reference Hamilton1989) can be applied to compute the likelihood function. Due to the potential difficulty in obtaining the closed-form expression of the stationary distribution, we use the conditional likelihood function given the initial latent variables $(x_0,s_0)=(\tilde {x},\tilde {s})$ . The choice of $\tilde {x}$ and $\tilde {s}$ can be arbitrary. In Section 4, we show the asymptotic properties of MLE with any choice of $\tilde {x}$ and $\tilde {s}$ , even if they are not the true underlying state and regime. The complete algorithm is stated as follows. For initialization, $\overline {x}_{0|0}^\theta =\tilde {x}$ and $\overline {\Omega }_{0|0}^\theta =0$ , $\overline {p}_r^\theta (s_0=\tilde {s})=1$ , and $\overline {p}_r^\theta (s_0=i)=0$ for $i\neq \tilde {s}$ . For $t=1,2,\ldots ,n$ ,

1. Step I. Calculate $\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t-1})=p_\theta (s_t|s_{t-1})\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^{t-1}|\mathbf {y}_1^{t-1})$ , where $p_\theta (s_t|s_{t-1})$ is given by the transition matrix $(q_{ij})$ , and $\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^{t-1}|\mathbf {y}_1^{t-1})$ is given by the initialization $\overline {p}^\theta _r(s_0)$ for $t=1$ and the previous Step V for $t\geq 2$ .

2. Step II. Apply (12)–(17). Use the outputs to compute $\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ according to (21).

3. Step III. Calculate $\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1})=\sum _{s_t,\dots ,s_{t-r+1}}\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t-1})$ , where $\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ is given by II, and $\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t-1})$ is given by Step I.

4. Step IV. Calculate $\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t})=\frac {\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1},\ \mathbf {s}_{t-r+1}^t)\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t-1})}{\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1})}$ , where $\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ is given by Step III, $\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t-1})$ is given by Step I, and $\overline {p}_{r}^\theta (y_t|\mathbf {y}_1^{t-1})$ is given by Step III.

5. Step V. Calculate $\overline {p}_{r}^\theta (\mathbf {s}_{t-r+2}^{t}|\mathbf {y}_1^{t})=\sum _{s_{t-r+1}}\overline {p}_{r}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t})$ and $\overline {p}_{r,\tilde {x}}^\theta (s_{t-r+1}|\mathbf {y}_1^{t},\mathbf {s}_{t-r+2}^{t})=\frac {\overline {p}_{r,\tilde {x}}^\theta (\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t})}{\overline {p}_{r,\tilde {x}}^\theta (\mathbf {s}_{t-r+2}^{t}|\mathbf {y}_1^{t})}$ .

6. Step VI. Apply (18) and (20) with the outputs in Steps II, IV, and V.

Let $\tilde {w}=(\tilde {x},\tilde {s})$ . We use the subscript $\tilde {w}$ to denote the initialization. The approximate log-likelihood function is

$$ \begin{align*} \overline{\ell}^{(r)}_{n,\tilde{w}}(\theta)=\log \overline{p}_{r}^\theta(y_n,\dots,y_{1}|x_0=\tilde{x},s_0=\tilde{s}) =\log \overline{p}_{r,\tilde{w}}^\theta(y_n,\dots,y_{1}) =\sum_{t=1}^{n} \log\overline{p}_{r,\tilde{w}}^\theta(y_t|\mathbf{y}_1^{t-1}), \end{align*} $$

where $\overline {p}_{r,\tilde {w}}^\theta (y_t|\mathbf {y}_1^{t-1})$ is obtained in Step III. The approximate MLE $\hat {\theta }^{(r)}_{n,\tilde {w}}=\mathrm{arg\,max}_{\theta \in \Theta }\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta )$ .

3.1 Expressions of Exact and Approximate Likelihood Functions

This section derives the expressions for the exact and approximate likelihood functions, from which we can analyze their difference explicitly. The exact likelihood function conditional on the initial latent variables is

(22) $$ \begin{align} p^\theta_{\tilde{w}}(y_n,\ldots,y_1) =&\sum_{s_n,\ldots,s_1}p_{\tilde{w}}^\theta(y_n,\ldots,y_1,s_n,\ldots,s_1)\notag\\ =&\sum_{s_n,\ldots,s_1}\prod_{t=1}^n (p_{\tilde{w}}^\theta(y_t|\mathbf{y}_1^{t-1},\mathbf{s}_1^{t})p_\theta(s_t|s_{t-1})), \end{align} $$

where

(23) $$ \begin{align} \begin{aligned} p_{\tilde{w}}^\theta(y_t|\mathbf{y}_1^{t-1},\mathbf{s}_1^t) =\frac{1}{(2\pi\Sigma^\theta_{t|t-1}(\mathbf{s}_1^t))^{\frac{1}{2}}} \exp\bigg\{\!\!\!-\frac{(y_t-y^\theta_{t|t-1}(\mathbf{s}_1^t))^2}{2\Sigma^\theta_{t|t-1}(\mathbf{s}_1^t)}\bigg\}. \end{aligned} \end{align} $$

The approximate likelihood function can be expressed as

(24) $$ \begin{align} \overline{p}_{r,\tilde{w}}^\theta(y_n,\dots,y_1) =\sum_{s_n,\dots,s_1} \prod_{t=1}^n\big[\overline{p}_{r,\tilde{w}}^\theta(y_t|\mathbf{y}_1^{t-1},\mathbf{s}_{t-r+1}^t)p_\theta(s_t|s_{t-1})\big], \end{align} $$

where $\overline {p}_{r,\tilde {w}}^\theta (y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ is given in (21). (24) can be derived by expanding the approximation algorithm step by step. More specifically,

(25) $$ \begin{align} &\overline{p}^{\theta}_{r,\tilde{w}}(y_n,y_{n-1},\ldots,y_1) =\overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1})\prod_{t=1}^{n-1} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1})\nonumber\\ &=\sum_{s_n,\dots,s_{n-r+1}}\overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1},\mathbf{s}_{n-r+1}^{n})\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{n-r+1}^{n}|\mathbf{y}_1^{n-1})\prod_{t=1}^{n-1} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1}) \end{align} $$

(26) $$ \begin{align} &=\sum_{s_n,\dots,s_{n-r+1}}\overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1},\mathbf{s}_{n-r+1}^{n})p_{\theta}(s_n|s_{n-1})\sum_{s_{n-r}}\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{n-r}^{n-1}|\mathbf{y}_1^{n-1})\prod_{t=1}^{n-1} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1}) \end{align} $$

(27) $$ \begin{align} &\begin{aligned} &=\sum_{s_n,\dots,s_{n-r}}\Big[\overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1},\mathbf{s}_{n-r+1}^{n})p_{\theta}(s_n|s_{n-1})\\ &\hspace{3.5cm}\times\frac{\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-1}|\mathbf{y}_1^{n-2},\mathbf{s}_{n-r}^{n-1})\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{n-r}^{n-1}|\mathbf{y}_1^{n-2})}{\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-1}|\mathbf{y}_1^{n-2})}\prod_{t=1}^{n-1} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1})\Big] \end{aligned} \end{align} $$

$$ \begin{align*} &\begin{aligned} &=\sum_{s_n,\dots,s_{n-r}}\Big[\overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1},\mathbf{s}_{n-r+1}^{n})p_\theta(s_n|s_{n-1})\\&\hspace{3.5cm}\times\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-1}|\mathbf{y}_1^{n-2},\mathbf{s}_{n-r}^{n-1})\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{n-r}^{n-1}|\mathbf{y}_1^{n-2})\prod_{t=1}^{n-2} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1})\Big], \end{aligned}\nonumber \end{align*} $$

where (25) follows from Step III, (26) follows by applying Steps I and V, and (27) follows from Step IV. If we continue to apply Steps I, V, and IV to $\overline {p}_{r,\tilde {w}}^\theta (\mathbf {s}_{n-r}^{n-1}|\mathbf {y}_1^{n-2})$ , then the likelihood function can be written as

$$ \begin{align*} &\overline{p}^{\theta}_{r,\tilde{w}}(y_n,y_{n-1},\ldots,y_1)\\ &=\sum_{s_n,\dots,s_{n-r-1}} \overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1},\mathbf{s}_{n-r+1}^{n})p_{\theta}(s_n|s_{n-1})\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-1}|\mathbf{y}_1^{n-2},\mathbf{s}_{n-r}^{n-1})p_{\theta}(s_{n-1}|s_{n-2})\\ &\hspace{3cm}\times\frac{\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-2}|\mathbf{y}_1^{n-3},\mathbf{s}_{n-r-1}^{n-2})\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{n-r-1}^{n-2}|\mathbf{y}_1^{n-3})}{\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-2}|\mathbf{y}_1^{n-3})}\prod_{t=1}^{n-2} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1})\\ &=\sum_{s_n,\dots,s_{n-r-1}} \overline{p}^{\theta}_{r,\tilde{w}}(y_n|\mathbf{y}_1^{n-1},\mathbf{s}_{n-r+1}^{n})p_{\theta}(s_n|s_{n-1})\overline{p}^{\theta}_{r,\tilde{w}}(y_{n-1}|\mathbf{y}_1^{n-2},\mathbf{s}_{n-r}^{n-1})p_\theta(s_{n-1}|s_{n-2})\\ &\hspace{3cm}\times \overline{p}^{\theta}_{r,\tilde{w}}(y_{n-2}|\mathbf{y}_1^{n-3},\mathbf{s}_{n-r-1}^{n-2})\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{n-r-1}^{n-2}|\mathbf{y}_1^{n-3})\prod_{t=1}^{n-3} \overline{p}^{\theta}_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1}). \end{align*} $$

If we repeat Steps I, V, and IV until $t=1$ , then (24) follows.

A close look at (22) and (24) shows that the difference in the two likelihood functions arises from the difference between $p^\theta _{\tilde {w}}(y_t|\mathbf {y}_1^{t-1},\mathbf {s}_1^t)$ and $\overline {p}^\theta _{r,\tilde {w}}(y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ for $t\geq r+1$ . Section 3.2 theoretically derives an upper bound of the difference.

3.2 Difference between Exact and Approximate Likelihood Functions

This section bounds the approximation error. Let $\Theta $ denote the parameter space. For a given $\theta \in \Theta $ , let $\alpha _j(\theta )$ denote the $\alpha _j$ -component of $\theta $ . For $i\in \mathbb S$ and $j=1,\ldots ,J-1$ , let $q_{ij}(\theta )$ denote the $q_{ij}$ -component of $\theta $ , and define $q_{iJ}(\theta )=1-\sum _{\ell =1}^{J-1}q_{i\ell }(\theta )$ . Define $\alpha _{\max }(\theta )\triangleq \max _{j\in \mathbb S}|\alpha _j(\theta )|$ as the maximum of the absolute values of the $\alpha _j$ -components across regimes under $\theta $ , and $\alpha _+\triangleq \textrm {sup}_{\theta \in \Theta }\alpha _{\max }(\theta )$ as the supremum of this quantity over all parameter values in $\Theta $ .

Throughout the article, we work on stationary and ergodic $(x_t,s_t,y_t)$ .

As stated at the end of Section 3.1, the difference between the approximate and exact likelihood functions arises from the difference in $p^{\theta }_{\tilde {w}}(y_t|\mathbf {y}_1^{t-1},\mathbf {s}_1^t)$ and $\overline {p}^{\theta }_{r,\tilde {w}}(y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ . Because both $p^{\theta }_{\tilde {w}}(y_t|\mathbf {y}_1^{t-1},\mathbf {s}_1^t)$ and $\overline {p}^{\theta }_{r,\tilde {w}}(y_t|\mathbf {y}_1^{t-1},\mathbf {s}_{t-r+1}^t)$ take the form of Gaussian densities, their difference comes from what we put in the places of means and variances or, more specifically, $\bar {y}^{\theta }_{t|t-1}(\mathbf {s}_{t-r+1}^t)-y^{\theta }_{t|t-1}(\mathbf {s}_1^t)$ and $\overline {\Sigma }^{\theta }_{t|t-1}(\mathbf {s}_{t-r+1}^t)-\Sigma ^{\theta }_{t|t-1}(\mathbf {s}_1^t)$ . Define

$$ \begin{align*} \Delta_{r,t}^{x}(\mathbf{s}_1^t;\theta)\triangleq\overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-x^\theta_{t|t}(\mathbf{s}_1^t),\qquad \Delta_{r,t}^\Omega(\mathbf{s}_1^t;\theta)\triangleq\overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-\Omega^\theta_{t|t}(\mathbf{s}_1^t). \end{align*} $$

From (6), (8), (13), and (15),

$$ \begin{align*} &\bar{y}^{\theta}_{t|t-1}(\mathbf{s}_{t-r+1}^t)-y^{\theta}_{t|t-1}(\mathbf{s}_1^t)=\beta_{s_t}\alpha_{s_t}\Delta_{r,t-1}^x(\mathbf{s}_1^{t-1};\theta),\\ &\overline{\Sigma}^{\theta}_{t|t-1}(\mathbf{s}_{t-r+1}^t)-\Sigma^{\theta}_{t|t-1}(\mathbf{s}_1^t)=\beta_{s_t}^2\alpha_{s_t}^2\Delta_{r,t-1}^\Omega(\mathbf{s}_1^{t-1};\theta). \end{align*} $$

Then it suffices to derive the upper bounds of $|\Delta _{r,t}^x(\mathbf {s}_1^t;\theta )|$ and $|\Delta _{r,t}^\Omega (\mathbf {s}_1^t;\theta )|$ . We first give the following lemma useful for understanding the derivation.

Lemma 2. Assume Assumptions 1–4. During the application of the Kalman filter (5)–(10), for one sequence of iterations $x^{\theta }_{t|t}(\mathbf {s}_{1}^{t})$ and $\Omega ^{\theta }_{t|t}(\mathbf {s}_{1}^{t})$ starting from $x_{0|0}^\theta $ and $\Omega _{0|0}^\theta $ and another sequence $\check {x}_{t|t}^\theta (\mathbf {s}_{1}^{t})$ and $\check {\Omega }_{t|t}^\theta (\mathbf {s}_{1}^{t})$ starting from $\check {x}_{0|0}^\theta $ and $\check {\Omega }_{0|0}^\theta $ , there exists a sequence of positive random variables $\xi _t$ and constants $c_\xi $ and c such that, for all $t=1,2,\ldots $ , $\mathbb {E}_{0}\xi _t^2\leq c_\xi $ , and for all $\theta \in \Theta $ ,

$$ \begin{align*} \max_{s_t,\ldots,s_1,s_0,\check{s}_0\in\mathbb{S}} |x_{t|t}^\theta(\mathbf{s}_{1}^{t})-\check{x}_{t|t}^\theta(\mathbf{s}_{1}^{t})|\leq& \alpha_{\max}(\theta)^{t} |x_{0|0}^\theta-\check{x}_{0|0}^\theta|+\alpha_{\max}(\theta)^{t}(\xi_t+c|\check{x}_{0|0}^\theta|)|\Omega_{0|0}^\theta-\check{\Omega}_{0|0}^\theta|,\\ \max_{s_t,\ldots,s_1,s_0,\check{s}_0\in\mathbb{S}}|\Omega_{t|t}^\theta(\mathbf{s}_{1}^{t})-\check{\Omega}_{t|t}^\theta(\mathbf{s}_{1}^{t})|\leq & \alpha_{\max}(\theta)^{2t}|\Omega_{0|0}^\theta-\check{\Omega}_{0|0}^\theta|. \end{align*} $$

Lemma 2 states that the Kalman filter is convergent. More specifically, for two sequences of Kalman filter iterations different only in their starting points, the difference in their outputs decays exponentially at the rate of $\alpha _{\max }(\theta )$ . Then we bound $|\Delta _{r,t}^x(\mathbf {s}_1^t;\theta )|$ and $|\Delta _{r,t}^\Omega (\mathbf {s}_1^t;\theta )|$ .

Theorem 1. Assume Assumptions 1–4. Constant $c_\Delta>0$ exists such that for any choice of initialization $w_0=\tilde {w}$ , there exist a constant $c_M>0$ and a sequence of positive random variables $M_t$ satisfying for all $t\geq r$ , $\mathbb {E}_0 \ \mathrm {sup}_{\theta \in \Theta } M_t^2<c_{M}$ , and for all $\theta \in \Theta $ ,

(28) $$ \begin{align} &\max_{\mathbf{s}_{1}^{t}\in\mathbb{S}^t}|\Delta_{r,t}^x(\mathbf{s}_1^t;\theta)|\leq c_\Delta\alpha_{\max}(\theta)^{r-1} M_{t-r+1}, \end{align} $$

(29) $$ \begin{align} &\max_{\mathbf{s}_{1}^{t}\in\mathbb{S}^t}|\Delta_{r,t}^\Omega(\mathbf{s}_1^t;\theta)|\leq c_\Delta\alpha_{\max}(\theta)^{2(r-1)}. \end{align} $$

We give the general idea of the proof of (29). The discrepancy at $t-1$ is $\Delta _{r,t-1}^\Omega (\mathbf {s}_1^{t-1};\theta )$ . Notice that, in spite of potentially different inputs, the prediction and update steps (12)–(17) in the approximation algorithm are the same as (5)–(10) in the Kalman filter. After the update step at time t, $|\Delta _{r,t}^\Omega (\mathbf {s}_1^t;\theta )|\leq \alpha _{\max }(\theta )^2|\Delta _{r,t-1}^\Omega (\mathbf {s}_1^{t-1};\theta )|$ , because the Kalman filter is convergent, as shown in Lemma 2. Then the collapsing step averages out the regime r periods ahead. The resulting difference is the discrepancy between the average and each path, which is bounded by the maximum difference between any two possible paths, differing only at the regime r periods ahead. It is reflected in (30):

(30) $$ \begin{align} |\Delta_{r,t}^\Omega(\mathbf{s}_1^t;\theta)|&\leq\alpha_{\max}(\theta)^2|\Delta_{r,t-1}^\Omega(\mathbf{s}_1^{t-1};\theta)| \nonumber\\ &+\max_{i,j\in\mathbb{S}}|\Omega^\theta_{t|t}(\mathbf{s}_{t-r+2}^t,s_{t-r+1}=i,\mathbf{s}_{1}^{t-r})-\Omega^\theta_{t|t}(\mathbf{s}_{t-r+2}^t,s_{t-r+1}=j,\mathbf{s}_{1}^{t-r})|. \end{align} $$

From Lemma 2, (30) is bounded by

$$ \begin{align*}\alpha_{\max}(\theta)^{2(r-1)}|\Omega^\theta_{t-r+1|t-r+1}(s_{t-r+1}=i,\mathbf{s}_{1}^{t-r})-\Omega^\theta_{t-r+1|t-r+1}(s_{t-r+1}=j,\mathbf{s}_{1}^{t-r})|.\end{align*} $$

Then (29) follows from upper bounds of $\Omega ^\theta _{t|t}(\mathbf {s}_1^t)$ and iteration of the inequality. The proof of (28) follows the same idea.

Equipped with Theorem 1, we can show an upper bound of the approximation error.

Proposition 1. Assume Assumptions 1–4. Constant $c>0$ exists such that for any initialization $w_0=\tilde {w}$ , there exist constants $c_M,c_N>0$ and sequences of positive random variables $M_t$ and $N_t$ satisfying for all $t\geq r$ , $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta }M_t^2<c_M$ , $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta }N_t<c_N$ , and for all $\theta \in \Theta $ ,

(31) $$ \begin{align} \begin{aligned} \max_{s_1,\ldots,s_t\in\mathbb{S}}|\log \overline{p}^\theta_{r,\tilde{w}}(y_t|\mathbf{y}_1^{t-1},\mathbf{s}_{t-r+1}^t)&-\log p_{\tilde{w}}^\theta(y_t|\mathbf{y}_1^{t-1},\mathbf{s}_1^t)|\leq c\alpha_{\max}(\theta)^{r-1}(N_t+M_{t-r+1}^2). \end{aligned} \end{align} $$

Moreover,

(32) $$ \begin{align} \big|\log \overline{p}^\theta_{r,\tilde{w}}(y_1,\dots,y_n)-\log p^\theta_{\tilde{w}}(y_1,\dots,y_n)\big| \leq c\alpha_{\max}(\theta)^{r-1} \sum_{t=r+1}^n(N_t+M_{t-r+1}^2). \end{align} $$

Kim and Kang (Reference Kim and Kang2019) used simulation to show that the approximation error is small when regimes are persistent. We provide an analytical explanation for this observation. When regimes are highly persistent and seldom switch, the likelihood is close to that of a single-regime state-space model, and the approximation error tends to be small. To see this analytically, note that the total approximation error can be decomposed across all possible regime paths. For persistent regimes, paths involving frequent regime switches carry low probability weights and thus contribute little to the overall error. The main contribution comes instead from paths where the regime remains invariant. For such a path with $s_t=s$ for all t, the approximation error for the filtered covariance matrix can be expressed as

(33) $$ \begin{align}& \big|\Delta_{r,t}^\Omega(\mathbf{s}_1^t;\theta)\big|\leq \alpha_{\max}(\theta)^2|\Delta_{r,t-1}^\Omega(\mathbf{s}_1^{t-1};\theta)|\nonumber\\ &\quad+\sum_{i\in\mathbb{S},i\neq s}\overline{p}_r^\theta(s_{t-r+1}=i|\mathbf{s}_{t-r+2}^t,\mathbf{y}_1^t)\big|\Omega_{t|t}^\theta(\mathbf{s}_{t-r+2}^t,s_{t-r+1}=i,\mathbf{s}_{1}^{t-r})-\Omega_{t|t}^\theta(\mathbf{s}_1^t)\big|. \end{align} $$

In (33), the term $\overline {p}_r^\theta (s_{t-r+1}=i|\mathbf {s}_{t-r+2}^t,\mathbf {y}_1^t)$ is small for $i\neq s$ , so the overall approximation error remains small under high regime persistence. The same reasoning applies to the approximation error for the filtered mean. Thus, the approximation error is small.

Regime persistence does not enter the theoretical bound in Proposition 1 explicitly. This is because, when deriving this bound, we condition on the full regime path and compute the maximum bound across all paths. As a result, Proposition 1 provides a uniform bound that holds for any regime sequence, regardless of its persistence.

Our simulation study in Section A.3 of the Supplementary Material adds one further observation. Although the approximation error is small when regimes are persistent, this relationship is not necessarily monotonic. Under certain specifications, the error may initially increase before decreasing as regimes become more persistent.

3.3 Difference between Exact and Approximate Inferred Latent Variables

In MS state-space models, $x_t$ and $s_t$ cannot be directly observed and thus require inference. This section examines the difference between the latent variable inference based on the exact likelihood function and that based on the approximation algorithm. We consider three types of inferences for the regime and the state at time t: the predictive inference using the information up to time $t-1$ , the updated inference using the information up to time t, and the smoothed inference using the full sample. The approximate inference is explained in Kim (Reference Kim1994). For readers’ convenience, we explain the computation method before we introduce our main findings.

For the regime process, its predictive inference $\{\overline {p}^{\theta }_{r,\tilde {w}}(s_t|\mathbf {y}_1^{t-1})\}_{1\leq t\leq n}$ and updated inference $\{\overline {p}^{\theta }_{r,\tilde {w}}(s_t|\mathbf {y}_1^{t})\}_{1\leq t\leq n}$ can be obtained in Steps I and IV of the algorithm to compute the likelihood function. For the smoothed probability, when $t=n$ , $\overline {p}^{\theta }_{r,\tilde {w}}(\mathbf {s}_{n-r+1}^{n}|\mathbf {y}_1^{n})$ is the last entry of $\{\overline {p}^{\theta }_{r,\tilde {w}}(\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t})\}_{1\leq t\leq n}$ obtained at Step IV. For $t<n$ , $\overline {p}^{\theta }_{r,\tilde {w}}(\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{n})$ can be computed backward given $\overline {p}^{\theta }_{r,\tilde {w}}(\mathbf {s}_{t-r+2}^{t+1}|\mathbf {y}_1^{n})$ followingFootnote ¹

(34) $$ \begin{align} \overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+1}^t|\mathbf{y}_1^n) =\sum_{s_{t+1}}\frac{\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+2}^{t+1}|\mathbf{y}_1^n)p_{\theta}(s_{t+1}|s_t)\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+1}^t|\mathbf{y}_1^{t})}{\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+2}^{t+1}|\mathbf{y}_1^{t})}. \end{align} $$

Theorem 2 gives the findings about the exact and approximate inferred probabilities.

Theorem 2. Assume Assumptions 2–4. Constant $c>0$ exists such that for any choice of initialization $w_0=\tilde {w}$ , there exist constants $c_M,c_N>0$ and sequences of positive random variables $M_t$ and $N_t$ satisfying $\mathbb {E}_0 \ \mathrm {sup}_{\theta \in \Theta } M_t^2<c_M$ and $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta } N_t<c_N$ for all $t\geq r$ , and for all $\theta \in \Theta $ ,

(35) $$ \begin{align} \max_{s_t\in\mathbb{S}}|\overline{p}^{\theta}_{r,\tilde{w}}(s_t|\mathbf{y}_1^{t-1}) -p^{\theta}_{\tilde{w}}(s_t|\mathbf{y}_1^{t-1})|&\leq c\alpha_{\max}(\theta)^{r-1}\sum_{k=r+1}^{t-1}(N_k+M_{k-r+1}^2), \end{align} $$

(36) $$ \begin{align} \max_{s_t\in\mathbb{S}}|\overline{p}^{\theta}_{r,\tilde{w}}(s_t|\mathbf{y}_1^{t}) -p^{\theta}_{\tilde{w}}(s_t|\mathbf{y}_1^{t})|&\leq c\alpha_{\max}(\theta)^{r-1}\sum_{k=r+1}^{t}(N_k+M_{k-r+1}^2), \end{align} $$

(37) $$ \begin{align} \max_{s_t\in\mathbb{S}}|\overline{p}^{\theta}_{r,\tilde{w}}(s_t|\mathbf{y}_1^{n}) -p^{\theta}_{\tilde{w}}(s_t|\mathbf{y}_1^{n})|&\leq c\alpha_{\max}(\theta)^{r-1}\sum_{k=r+1}^{n}(N_k+M_{k-r+1}^2). \end{align} $$

To show Theorem 2, we first derive the expressions of approximate inferred probabilities

(38) $$ \begin{align} \overline{p}^{\theta}_{r,\tilde{w}}(s_t|\mathbf{y}_1^{t-1}) &=\sum_{s_{t-1},\ldots,s_{1}}\frac{p_{\theta}(s_t|s_{t-1})\prod_{k=1}^{t-1}(\overline{p}^{\theta}_{r,\tilde{w}}(y_{k}|\mathbf{y}_1^{k-1},\mathbf{s}_{k-r+1}^{k})p_{\theta}(s_k|s_{k-1}))} {\prod_{k=1}^{t-1}\overline{p}^{\theta}_{r,\tilde{w}}(y_{k}|\mathbf{y}_1^{k-1})}, \end{align} $$

(39) $$ \begin{align} \overline{p}^{\theta}_{r,\tilde{w}}(s_t|\mathbf{y}_1^{t}) &=\sum_{s_{t-1},\ldots,s_{1}}\frac{\prod_{k=1}^{t}(\overline{p}^{\theta}_{r,\tilde{w}}(y_{k}|\mathbf{y}_1^{k-1},\mathbf{s}_{k-r+1}^{k})p_{\theta}(s_k|s_{k-1}))} {\prod_{k=1}^{t}\overline{p}^{\theta}_{r,\tilde{w}}(y_{k}|\mathbf{y}_1^{k-1})}, \end{align} $$

(40) $$ \begin{align} \overline{p}^{\theta}_{r,\tilde{w}}(s_t|\mathbf{y}_1^n) &=\sum_{s_n,\ldots,s_{t+1},s_{t-1},\ldots,s_1} \frac{\prod_{k=1}^n(\overline{p}^{\theta}_{r,\tilde{w}}(y_k|\mathbf{y}_1^{k-1},\mathbf{s}_{k-r+1}^{k})p_{\theta}(s_k|s_{k-1}))} {\prod_{k=1}^n \overline{p}^{\theta}_{r,\tilde{w}}(y_k|\mathbf{y}_1^{k-1})}. \end{align} $$

The derivation of the expressions essentially parallels the derivation of the approximate likelihood expression in Section 3.1 by the repeated application of Steps I, V, and IV. The details are deferred to the Supplementary Material. The expressions of the exact inferred probabilities can be obtained by replacing $\overline {p}^{\theta }_{r,\tilde {w}}(y_k|\mathbf {y}_1^{k-1},\mathbf {s}_{k-r+1}^{k})$ and $\overline {p}^{\theta }_{r,\tilde {w}}(y_k|\mathbf {y}_1^{k-1})$ with $p^{\theta }_{\tilde {w}}(y_k|\mathbf {y}_1^{k-1},\mathbf {s}_{1}^{k})$ and $p^{\theta }_{\tilde {w}}(y_k|\mathbf {y}_1^{k-1}),$ respectively. Then Theorem 2 follows if we apply Proposition 1.

Next, we explain the state inference algorithm. Let $x_{t|k}^{\theta }$ and $\overline {x}_{t|k}^{(r,\theta )}$ denote the exact and approximate inference about the state given $\{y_k,\ldots ,y_1,x_0,s_0\}$ , respectively. We might omit the superscript $(r)$ and simply write $\overline {x}^{\theta }_{t|k}$ when there is no confusion. The approximate predictive and updated inference can be computed from

$$ \begin{align*} \overline{x}^{\theta}_{t|t-1}&=\sum_{\mathbf{s}_{t-r+1}^t}\overline{x}^{\theta}_{t|t-1}(\mathbf{s}_{t-r+1}^t)\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+1}^t|\mathbf{y}_1^{t-1}),\\ \overline{x}^{\theta}_{t|t}&=\sum_{\mathbf{s}_{t-r+2}^{t}}\overline{x}^{\theta}_{t|t}(\mathbf{s}_{t-r+2}^{t})\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+2}^{t}|\mathbf{y}_1^{t}), \end{align*} $$

where $\overline {x}^{\theta }_{t|t-1}(\mathbf {s}_{t-r+1}^t)$ and $\overline {x}^{\theta }_{t|t}(\mathbf {s}_{t-r+2}^{t})$ can be obtained from (12) and (18), and $\overline {p}^{\theta }_{r,\tilde {w}}(\mathbf {s}_{t-r+1}^t|\mathbf {y}_1^{t-1})$ and $\overline {p}^{\theta }_{r,\tilde {w}}(\mathbf {s}_{t-r+2}^{t}|\mathbf {y}_1^{t})$ can be obtained in Steps I and IV. For the smoothed inference, when $t=n$ , $\overline {x}^{\theta }_{n|n}(\mathbf {s}_{n-r+2}^n)$ can be obtained from (18). For $t<n$ ,

$$ \begin{align*}\overline{x}^{\theta}_{t|n}=\sum_{\mathbf{s}_{t-r+2}^{t}}\overline{p}^{\theta}_{r,\tilde{w}}(\mathbf{s}_{t-r+2}^{t}|\mathbf{y}_1^n)\overline{x}^{\theta}_{t|n}(\mathbf{s}_{t-r+2}^{t}),\end{align*} $$

where $\overline {x}^{\theta }_{t|n}(\mathbf {s}_{t-r+2}^{t})$ can be computed backward following

$$ \begin{align*} &\overline{x}^{\theta}_{t|n}(\mathbf{s}_{t-r+2}^{t+1}) =\overline{x}^{\theta}_{t|t}(\mathbf{s}_{t-r+2}^{t})+\overline{J}^{\theta}_t(\mathbf{s}_{t-r+2}^{t+1})(\overline{x}^{\theta}_{t+1|n}(\mathbf{s}_{t-r+3}^{t+1})-\overline{x}^{\theta}_{t+1|t}(\mathbf{s}_{t-r+2}^{t+1})),\\ &\overline{x}^{\theta}_{t|n}(\mathbf{s}_{t-r+2}^{t})=\sum_{s_{t+1}}\overline{p}^{\theta}_{r,\tilde{w}}(s_{t+1}|\mathbf{s}_{t-r+2}^{t},\mathbf{y}_1^n)\overline{x}^{\theta}_{t|n}(\mathbf{s}_{t-r+2}^{t+1}), \end{align*} $$

where $\overline {J}^{\theta }_t(\mathbf {s}_{t-r+2}^{t+1})=\frac {\alpha _{s_{t+1}}\overline {\Omega }^{\theta }_{t|t}(\mathbf {s}_{t-r+2}^{t})}{\overline {\Omega }^{\theta }_{t+1|t}(\mathbf {s}_{t-r+2}^{t+1})}$ . The difference between exact and approximate inference of $x_t$ can be summarized as follows.

Theorem 3. Assume Assumptions 1–4. Constant $c_\Delta>0$ exists such that for any choice of initialization $w_0=\tilde {w}$ , there exist constants $c_{n},c_{N},c_{M},c_A>0$ and sequences of positive random variables $\{n_t\}_{t\geq 1}$ , $\{M_t\}_{t\geq 1}$ , $\{N_t\}_{t\geq 1}$ , and $\{A_{t,k}\}_{t\geq 1,k\geq 1}$ satisfying $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta }n_t^2<c_{n}$ , $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta }N_t^2<c_{N}$ , $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta }M_t^4<c_{M}$ , and $\mathbb {E}_0\ \mathrm {sup}_{\theta \in \Theta }\ A_{t,k}\leq c_{A}$ for all $t,k\geq 1$ , and for all $\theta \in \Theta $ ,

(41) $$ \begin{align} |\overline{x}_{t|t-1}^{\theta}-x_{t|t-1}^{\theta}|\leq &\alpha_{\max}(\theta)^{r-1}rn_{t-1}\sum_{k=r+1}^{t-1} (N_k+M_{k-r+1}^2)+c_\Delta\alpha_{\max}(\theta)^{r-1} M_{t-r+1}, \end{align} $$

(42) $$ \begin{align} |\overline{x}_{t|t}^{\theta}-x_{t|t}^{\theta}|\leq &\alpha_{\max}(\theta)^{r-1}rn_{t}\sum_{k=r+1}^{t} (N_k+M_{k-r+1}^2)+c_\Delta\alpha_{\max}(\theta)^{r-1} M_{t-r},\end{align} $$

(43) $$ \begin{align} |\overline{x}_{t|n}^{\theta}-x_{t|n}^{\theta}|\leq& \alpha_{\max}(\theta)^{r-1}r n_t\sum_{k=r+1}^n (N_k+M_{k-r+1}^2)\nonumber\\ &\qquad+2c_\Delta\alpha_{\max}(\theta)^{r-1} M_{t-r+1} +c_\Delta\alpha_{\max}(\theta)^{r-1} \sum_{k=t+1}^n A_{t+1,k}. \end{align} $$

4.1 Asymptotic Properties of the Exact MLE

The asymptotic properties of the exact MLE relate to the literature of hidden Markov models. The MS state-space model belongs to the class of hidden Markov models because the latent process $(x_t,s_t)$ is a Markov chain, and given $(x_t,s_t)$ , $(y_t)$ is a sequence of conditionally independent random variables. Consistency follows from Douc and Moulines (Reference Douc and Moulines2012), which shows the consistency of the MLE in a wide range of hidden Markov models.

Asymptotic normality of the exact MLE follows from a central limit theorem for the Fisher score function

(44) $$ \begin{align} \nabla_\theta \ell_{n,\tilde{w}}(\theta)=\sum_{t=1}^n \nabla_\theta\log p^{\theta}_{\tilde{w}}(y_t|\mathbf{y}_1^{t-1}) \end{align} $$

and a locally uniform law of large numbers for the observed Fisher information

(45) $$ \begin{align} \nabla_\theta^2 \ell_{n,\tilde{w}}(\theta)=\sum_{t=1}^n \nabla_\theta^2\log p^{\theta}_{\tilde{w}}(y_t|\mathbf{y}_1^{t-1}). \end{align} $$

The summands in (44) and (45) do not form stationary ergodic sequences, because $p_{\theta ,\tilde {w}}(y_t|\mathbf {y}_1^{t-1})$ depends on the accumulating $\mathbf {y}_1^{t-1}$ . Thus, we cannot directly apply the law of large numbers or the central limit theorem. The strategy, originated from Baum and Petrie (Reference Baum and Petrie1966), is to approximate the summands with stationary ergodic sequences. Bickel et al. (Reference Bickel, Ritov and Rydén1998), Jensen and Petersen (Reference Jensen and Petersen1999), Douc et al. (Reference Douc, Moulines and Rydén2004), Kasahara and Shimotsu (Reference Kasahara and Shimotsu2019), Pouzo et al. (Reference Pouzo, Psaradakis and Sola2022), and Li and Liu (Reference Li and Liu2023) then take this approach. A key step in their approximation is that the upper bounds of the difference between the summands and stationary ergodic sequences involve

$$ \begin{align*} \mathop{\mathrm{sup}}\limits_{w_{k-1},w_k} |\nabla_\theta^j\log p_\theta(y_k,w_k|w_{k-1})|\Big(1-\frac{\inf_{\theta}\inf_{w_{k-1},w_k}p_\theta(w_k|w_{k-1})}{\mathrm{sup}_\theta\mathrm{sup}_{w_{k-1},w_k}p_\theta(w_k|w_{k-1})}\Big)^k, \end{align*} $$

$j=1,2$ . Because they assumed that $\textrm {sup}_{w_{k-1},w_k} |\nabla _\theta ^j\log p_\theta (y_k,w_k|w_{k-1})|$ is bounded and that $\frac {\inf _{\theta }\inf _{w_{k-1},w_k}p_\theta (w_k|w_{k-1})}{\textrm {sup}_\theta \textrm {sup}_{w_{k-1},w_k}p_\theta (w_k|w_{k-1})}>0$ , the difference decays exponentially with k. For MS state-space models, however, $\textrm {sup}_{w_{k-1},w_k} |\nabla _\theta ^j\log p_\theta (y_k,w_k|w_{k-1})|=+\infty $ and $\inf _{\theta }\inf _{w_{k-1},w_k}p_\theta (w_k|w_{k-1})=0$ , and these two assumptions no longer hold. Thus, the previous results do not apply in this study.

We contribute to the literature by showing asymptotic normality of the exact MLE in MS state-space models. Two instruments play an important role in showing the approximation. The first one is the exponential forgetting of the filtering distribution established by Douc and Moulines (Reference Douc and Moulines2012, Prop. 5). It imposes a much less stringent assumption on the transition density and only requires a positive lower bound on a large enough set, which holds in MS state-space models. The second instrument is the relaxed upper bound of the score function and Hessian. We show that the score function and Hessian at each period are bounded above by random variables with finite second moments, which suffices to show the approximation in probability. With these two instruments combined, we can show the influence of the information far in the past vanishes quickly. Then the original sequences can be approximated with the stationary sequences starting from past infinity. Next, we introduce the notation and present the two instruments in Lemmas 3 and 4.

Because our subsequent analysis involves operations on transition kernels, we first recall the notation for these operations on a generic measurable space. Let $(\mathsf {X}, \mathcal {X})$ be an arbitrary measurable space, and let L be a transition kernel on $(\mathsf {X}, \mathcal {X})$ . That is, for any $x \in \mathsf {X}$ , $L(x, \cdot )$ is a finite measure on $(\mathsf {X}, \mathcal {X}),$ and for any $A \in \mathcal {X}$ , $x \mapsto L(x, A)$ is measurable function from $(\mathsf {X}, \mathcal {X})$ to $([0, 1], \mathcal {B}([0, 1]))$ . L acts on bounded functions f on $\mathsf {X}$ and on $\sigma $ -finite positive measures $\mu $ on $(\mathsf {X}, \mathcal {X})$ via

$$\begin{align*}Lf(x) = \delta_x Lf \triangleq \int L(x, \mathrm{d}y) f(y), \qquad \mu L(A) = \mu L\mathbb{1}_A \triangleq \int \mu(\mathrm{d}x) L(x, A). \end{align*}$$

If $L_1$ and $L_2$ are two transition kernels on $(\mathsf {X}, \mathcal {X})$ , then $L_1 L_2$ is the transition kernel on $(\mathsf {X}, \mathcal {X})$ , given, for any $x \in \mathsf {X}$ and $A \in \mathcal {X}$ by

$$ \begin{align*} L_1 L_2(x, A) = \int L_1(x, \mathrm{d}y) L_2(y, A). \end{align*} $$

For the MS state-space model, let $w_t=(x_t,s_t)'$ denote the latent process taking values on $\mathbb {W}=\mathbb {R}\times \mathbb {S}$ and with the associated Borel $\sigma $ -field $\mathcal {W}=\mathcal {B}(\mathbb {R}\times \mathbb {S})$ , and $\varpi $ denote the product measure of Lebesgue measure and counting measure on $(\mathbb {W},\mathcal {W})$ . We use $g_\theta (y_t|w_t)$ to denote the conditional density of $y_t$ given $w_t$ . We use $\mathcal {H}_\theta (w,A)$ with $w\in \mathbb {W}$ and $A\in \mathcal {W}$ to denote the transition kernel of the Markov process $w_t$ . Let $h_\theta (w,\cdot )$ denote the density of $\mathcal {H}_\theta (w,\cdot )$ with respect to $\varpi $ . We follow Douc and Moulines (Reference Douc and Moulines2012) and define the unnormalized kernel

$$ \begin{align*} \mathbf{L}^\theta\langle y_{-m}^{t}\rangle(w_{-m},A) =\int\dots\int\big[\prod_{k=-m}^{t}g_\theta(y_k|w_k)\mathcal{H}_\theta(w_k,dw_{k+1})\big]\mathbb{1}_A(w_{t+1}) \end{align*} $$

for $w_{-m}\in \mathbb {W}$ and $A\in \mathcal {W}$ . We use $\phi _{\theta ,k}$ as the shorthand notation of

(46) $$ \begin{align} \phi_{\theta}(y_k,(w_k,w_{k-1}))\triangleq \nabla_\theta\log p_\theta(y_{k},w_{k}|w_{k-1}) =\nabla_\theta \log(g_\theta(y_{k}|w_{k})h_\theta(w_{k-1},w_k)) \end{align} $$

and use $\dot {\phi }_{\theta ,k}$ as the shorthand notation of

(47) $$ \begin{align} \dot{\phi}_{\theta,k}(y_k,(w_k,w_{k-1})) \triangleq\nabla_\theta^2\log p_\theta(y_k,w_k|w_{k-1}) =\nabla_\theta^2 \log(g_\theta(y_k|w_k)h_\theta(w_{k-1},w_k)). \end{align} $$

Lemma 3 (Proposition 5 in Douc and Moulines, Reference Douc and Moulines2012).

Assume Assumptions 1–4. Let $0\leq \gamma ^-< \gamma ^+\leq 1$ . There exists a $\mathbb {P}_0$ -a.s. finite random variable K such that for any $\eta>0$ , there exists $\rho \in (0,1)$ satisfying

(48) $$ \begin{align} \begin{aligned} \Big|\frac{\xi \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{f}} {\xi \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{g}} -\frac{\xi' \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{f}} {\xi' \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{g}}\Big| \leq \rho^{\lfloor (k+m+1)(\gamma-\gamma^-)\rfloor} \Big[\frac{\xi \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle |\mathrm{f}|} {\xi \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{g}} +\frac{\xi' \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle|\mathrm{f}|} {\xi' \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{g}}\Big]&\\ +~2\eta^{\lfloor (k+m+1)(\gamma^+-\gamma)\rfloor/2} \frac{\big( \prod_{\ell={-m}}^k |\mathbf{L}^\theta\langle y_\ell\rangle(\cdot,\mathbb{W})|_\infty^2 \big)|\mathrm{f}|_\infty|\mathrm{g}|_\infty}{(\xi \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{g})( \xi' \mathbf{L}^\theta\langle \mathbf{y}_{-m}^{k} \rangle \mathrm{g})},& \end{aligned} \end{align} $$

for any $k+m\geq K$ , any $\gamma \in (\gamma ^-,\gamma ^+)$ , any bounded function $\mathrm {f}$ and any nonnegative bounded function $\mathrm {g}$ , any probability measure $\xi $ and $\xi '$ on $(\mathbb {W},\mathcal {W}),$ and any $\theta \in \Theta $ .

The left-hand side of (48) is the difference between two variables starting from distinct initial distributions $\xi $ and $\xi '$ . Because we can set $\rho ,\eta \in (0,1)$ , Lemma 3 states that the influence of the initial distributions decays geometrically if the random variables on the right-hand side of (48) can be bounded properly. Douc et al. (Reference Douc, Moulines and Rydén2004, Cor. 1), Kasahara and Shimotsu (Reference Kasahara and Shimotsu2019, Lem. 1), Pouzo et al. (Reference Pouzo, Psaradakis and Sola2022, Thm. 2), and Li and Liu (Reference Li and Liu2023, Lem. 3.1) serve the same function. They assumed $\frac {\inf _{\theta }\inf _{w_{k-1},w_k}h_\theta (w_{k-1},w_k)}{\textrm {sup}_\theta \textrm {sup}_{w_{k-1},w_k}h_\theta (w_{k-1},w_k)}>0$ and that the decaying rate is $1-\frac {\inf _{\theta }\inf _{w_{k-1},w_k}h_\theta (w_{k-1},w_k)}{\textrm {sup}_\theta \textrm {sup}_{w_{k-1},w_k}h_\theta (w_{k-1},w_k)}<1$ . Their rate does not apply here because $x_t$ takes values on the entire real line and $\inf _{\theta }\inf _{w_{k-1},w_k}p_\theta (w_k|w_{k-1})=0$ . Lemma 3 extends the previous results. It considers the difference on a compact set C and its complement separately. The decaying rate on C is $\rho \triangleq 1-\frac {\inf _{\theta }\inf _{w_{k-1},w_k\in C}h_\theta (w_{k-1},w_k)^2}{\textrm {sup}_\theta \textrm {sup}_{w_{k-1},w_k\in C}h_\theta (w_{k-1},w_k)^2}<1$ , which parallels the previous results. To control the difference on C’s complement, it is assumed $\textrm {sup}_{w\in C^c}g_\theta (y|w)<\eta \textrm {sup}_{w\in \mathbb {W}}g_\theta (y|w)$ for all $\theta $ and y. This assumption is satisfied if C is chosen to be large enough. The decaying rate on C’s complement is $\eta $ , which extends the previous results to noncompact sets.

Lemma 4 shows the relaxed upper bounds of the period score function and Hessian. We need the following additional assumptions.

We need the high-level Assumption 5 to show the approximation in $L^p$ . In the special cases of bounded state space, Douc et al. (Reference Douc, Moulines and Rydén2004), Pouzo et al. (Reference Pouzo, Psaradakis and Sola2022), and Li and Liu (Reference Li and Liu2023) could show Assumption 5 with verifiable assumptions. In our model with a non-compact state space, Lemmas 3 and 4 help to show the approximation only in probability.

4.1.1 A Central Limit Theorem for the Score Function

This section shows asymptotic normality of the score function. First, we show that the score function can be approximated by a sequence of integrable martingale increments. Using the Fisher identity,Footnote ² which states that the score function of a model with missing data equals the conditional expectation of the complete score given the observed data, the score function can be written as

$$ \begin{align*} \nabla_\theta\ell_{n,\tilde{w}}(\theta) =\sum_{t=1}^n \Phi_{t,0,\tilde{w}}(\theta), \end{align*} $$

where

(49) $$ \begin{align} \Phi_{t,m,\tilde{w}}(\theta) &\triangleq \nabla_\theta\log p_\theta(y_t|\mathbf{y}_{-m+1}^{t-1},w_{-m}=\tilde{w})\nonumber\\ &=\mathbb{E}_{\theta}\big[\sum_{k=-m+1}^{t}\phi_{\theta,k}|\mathbf{y}_{-m+1}^t,w_{-m}=\tilde{w}\big] -\mathbb{E}_{\theta}\big[\sum_{k=-m+1}^{t-1}\phi_{\theta,k}|\mathbf{y}_{-m+1}^{t-1},w_{-m}=\tilde{w}\big]. \end{align} $$

We can similarly define the period score function that starts from the stationary initial distribution

(50) $$ \begin{align} \Phi_{t,m}(\theta) \triangleq \nabla_\theta\log p_\theta(y_t|\mathbf{y}_{-m+1}^{t-1}) =\mathbb{E}_{\theta}\big[\sum_{k=-m+1}^{t}\phi_{\theta,k}|\mathbf{y}_{-m+1}^t\big] -\mathbb{E}_{\theta}\big[\sum_{k=-m+1}^{t-1}\phi_{\theta,k}|\mathbf{y}_{-m+1}^{t-1}\big]. \end{align} $$

Lemma 5 shows that the difference in the period score function with different starting points diminishes at a geometric rate. It parallels Douc et al. (Reference Douc, Moulines and Rydén2004, Lem. 8), Kasahara and Shimotsu (Reference Kasahara and Shimotsu2019, Lem. 4), Pouzo et al. (Reference Pouzo, Psaradakis and Sola2022, Lem. 7), and Li and Liu (Reference Li and Liu2023, Lem. 5.1). In their theories, the decaying rate $\tilde {\rho }=1-\frac {\inf _{\theta }\inf _{w_{k-1},w_k}h_\theta (w_{k-1},w_k)}{\textrm {sup}_\theta \textrm {sup}_{w_{k-1},w_k}h_\theta (w_{k-1},w_k)}$ is assumed to be strictly smaller than unity. For MS state-space models, $\inf _{\theta }\inf _{w_{k-1},w_k}h_\theta (w_{k-1},w_k)=0$ , and thus $\tilde {\rho }=1$ , the difference is non-decaying, and the previous theories no longer apply. We use Lemmas 3 and 4 to establish the geometric decaying rate. More specifically, we write (49) and (50) when $m=0,$ respectively, as

$$ \begin{align*} \mathbb{E}_\theta[\phi_{\theta,k}|\mathbf{y}_{1}^t,w_{0}=\tilde{w}] =\frac{\chi\mathbf{L}^\theta\langle \mathbf{y}_{1}^{k-2}\rangle \mathrm{f}_{k,t}^{\theta}} {\chi\mathbf{L}^\theta\langle \mathbf{y}_{1}^{k-2}\rangle \mathrm{g}_{k,t}^\theta},\quad \mathbb{E}_\theta[\phi_{\theta,k}|\mathbf{y}_{1}^t] =\frac{\overline{\chi}\mathbf{L}^\theta\langle \mathbf{y}_{1}^{k-2}\rangle \mathrm{f}_{k,t}^{\theta}} {\overline{\chi}\mathbf{L}^\theta\langle \mathbf{y}_{1}^{k-2}\rangle \mathrm{g}_{k,t}^{\theta}}, \end{align*} $$

for $3\leq k\leq t$ , where $\chi =\mathbb {P}_\theta (w_1\in A|w_0=\tilde {w})$ , $\overline {\chi }$ denotes the stationary distribution of $w_t$ ,

(51) $$ \begin{align} \mathrm{f}_{k,t}^{\theta}\langle \mathbf{y}_{k-1}^t\rangle (w_{k-1}) &=\int\dots\int \phi_{\theta,k}\prod_{\ell=k-1}^t\big[g_\theta(y_\ell|w_\ell)\mathcal{H}_\theta(w_\ell,dw_{\ell+1})\big]\mathbb{1}_{\mathbb{W}}(w_{t+1}), \end{align} $$

(52) $$ \begin{align} \mathrm{g}_{k,t}^\theta\langle \mathbf{y}_{k-1}^t\rangle (w_{k-1}) &=\int\dots\int \prod_{i=k-1}^t\big[g_\theta(y_\ell|w_\ell)\mathcal{H}_\theta(w_\ell,dw_{\ell+1})\big]\mathbb{1}_{\mathbb{W}}(w_{t+1}). \end{align} $$

Then Lemma 3 provides an upper bound of $\big |\mathbb {E}_\theta [\phi _{\theta ,k}|\mathbf {y}_{1}^t,w_{0}=\tilde {w}]-\mathbb {E}_\theta [\phi _{\theta ,k}|\mathbf {y}_{1}^t]\big |$ , which leads to the upper bound of $|\Phi _{t,0,\tilde {w}}(\theta _0)-\Phi _{t,0}(\theta _0)|$ in Lemma 5(a). Because $\rho ,\eta \in (0,1)$ in Lemma 3, the difference is decaying geometrically as k increases. Lemma 4 helps bound the other terms on the right-hand side of (48). Lemma 5(b) follows similarly.

Lemma 5. Assume Assumptions 1–5. Then there exist constants $\varrho \in (0,1)$ , $c>~0,$ and a random variable $B\in L^2(\mathbb {P}_0)$ such that for some $\mathbb {P}_0$ -a.s. finite random variable K,

1. (a) $|\Phi _{t,0,\tilde {w}}(\theta _0)-\Phi _{t,0}(\theta _0)|\leq \varrho ^t(t^2B+c\tilde {x})$ , $\mathbb {P}_0$ -a.s. for $t>K$ and any $\tilde {w}\in \mathbb {W}$ ;

2. (b) $|\Phi _{t,m}(\theta _0)-\Phi _{t,m'}(\theta _0)| \leq Bt^2(m^2+1)\varrho ^{t+m}$ , $\mathbb {P}_0$ -a.s. for $t+m>K$ and $m'\geq m $ .

Lemma 5 implies that $\{\Phi _{t,m}(\theta _0)\}_{m\geq 0}$ is a Cauchy sequence and converges to $\Phi _{t,\infty }(\theta _0)\triangleq \lim _{m\rightarrow \infty }\Phi _{t,m}(\theta _0)$ . Define the filtration $\mathcal {F}_t=\sigma (y_k,-\infty < k\leq t)$ . From Assumption 5, Durrett (Reference Durrett2019, Exer. 2.3.7 and Thm. 5.5.9), and $\mathbb {E}_0[\Phi _{t,m}(\theta _0)|\mathbf {y}_{-m+1}^{t-1}]=0$ , it follows that $\mathbb {E}_0[\Phi _{t,\infty }(\theta _0)|\mathbf {y}_{-\infty }^{t-1}]=0$ and $I(\theta _0)\triangleq \mathbb {E}_{0}[\Phi _{t,\infty }(\theta _0)\Phi _{t,\infty }(\theta _0)^T]<\infty $ . Therefore, $\{\Phi _{t,\infty }(\theta _0)\}_{t=-\infty }^\infty $ is an $(\mathcal {F},\mathbb {P}_0)$ -adapted stationary ergodic and square integrable martingale increment sequence. The central limit theorem implies $n^{-\frac {1}{2}}\sum _{t=1}^n\Phi _{t,\infty }(\theta _0)\Rightarrow \mathbb {N}(0,I(\theta _0))$ .

Because the influence of the initial distribution decays geometrically, as shown in Lemma 5, $\Phi _{t,0,\tilde {w}}(\theta )$ can be approximated with $\Phi _{t,\infty }(\theta )$ , which is summarized in Lemma 6.

Lemma 6. Assume Assumptions 1–5. Then, for any $\tilde {w}\in \mathbb {W}$ ,

$$ \begin{align*} \Big|\frac{1}{\sqrt{n}}\sum_{t=1}^n\big(\Phi_{t,0,\tilde{w}}(\theta_0)-\Phi_{t,\infty}(\theta_0)\big)\Big|\overset{a.s.}{\rightarrow}0,\qquad \text{ as }n\rightarrow \infty. \end{align*} $$

Then $n^{-\frac {1}{2}}\sum _{t=1}^n\Phi _{t,0,\tilde {w}}(\theta _0)$ has the same limiting distribution as $n^{-\frac {1}{2}}\sum _{t=1}^n\Phi _{t,\infty }(\theta _0)$ . The score function is asymptotically normally distributed.

Theorem 5. Assume Assumptions 1–5. Then, for any $\tilde {w}\in \mathbb {W}$ , $n^{-\frac {1}{2}}\nabla _\theta \ell _{n,\tilde {w}}(\theta _0) \Rightarrow \mathbb {N}(0,I(\theta _0))$ .

4.1.2 Law of Large Numbers for the Observed Fisher Information

This section presents the law of large numbers for the observed Fisher information. Using Louis missing information principle,Footnote ³ we write

$$ \begin{align*} \nabla_\theta^2\ell_{n,\tilde{w}}(\theta) =\mathbb{E}_{\theta}\Big[\sum_{t=1}^n\dot{\phi}_{\theta,t}\big|\mathbf{y}_{1}^n,w_0=\tilde{w}\Big] +\mathrm{Var}_\theta\Big[\sum_{t=1}^n\phi_{\theta,t}\big|\mathbf{y}_{1}^{n},w_0=\tilde{w}\Big] \end{align*} $$

with $\dot {\phi }_{\theta ,t}$ defined in (47). Similar to Section 4.1.1, we first approximate the Hessian with a stationary ergodic sequence. The proof parallels that in Section 4.1.1 and we defer the details to the Supplementary Material.

Theorem 6. Assume Assumptions 1–5. Let $\{\theta _n^*\}$ be any, possibly stochastic, sequence in $\Theta $ such that $\theta _n^*\overset {a.s.}{\rightarrow }\theta _0$ . Then, for any $\tilde {w}\in \mathbb {W}$ , $-n^{-1}\nabla _\theta ^2\ell _{n,\tilde {w}}(\theta _n^*)\overset {a.s.}{\rightarrow }I(\theta _0)$ .

Theorems 5 and 6 together provide the following theorem of asymptotic normality.

Theorem 7. Assume Assumptions 1–5 hold and that $I(\theta _0)$ is positive definite. Then, for any $\tilde {w}\in \mathbb {W}$ , $n^{1/2}(\tilde {\theta }_{n,\tilde {w}}-\theta _0)\Rightarrow \mathbb {N}(0,I(\theta _0)^{-1})$ .

4.2 Asymptotic Properties of the Approximate MLE

This section establishes the asymptotic properties of the approximate MLE and provides a consistent estimate of the asymptotic covariance matrix. We first show in Theorems 8 and 9 that the normalized likelihood function, score function, and Hessian matrix converge to the exact ones uniformly if $r_n$ increases at proper rates with n. Then Theorem 10 establishes the asymptotic normality of the approximate MLE and consistency of the asymptotic covariance matrix estimates.

Theorem 8 shows that the difference between the normalized approximate likelihood and the exact one diminishes to zero. Moreover, a larger $\alpha _+$ requires $r_n$ to grow at a faster rate.

Theorem 8. Assume Assumptions 1–4. If $r_n$ satisfies $\alpha _+^{r_n}\rightarrow 0$ as $n\rightarrow \infty $ , then for any initial $\tilde {w}$ ,

$$ \begin{align*} \mathop{\mathrm{sup}}\limits_{\theta\in\Theta}|\frac{1}{n}\log \overline{p}_{r_n,\theta}(y_1,\dots,y_n|w_0=\tilde{w})-\frac{1}{n} \log p_{\theta}(y_1,\dots,y_n|w_0=\tilde{w})|\overset{p}{\rightarrow} 0, \text{ as }n\rightarrow\infty. \end{align*} $$

As stated in Section 3.2, the discrepancy between the approximate and exact likelihood functions arises essentially from $\overline {x}_{t|t}^\theta (\mathbf {s}_{t-r+2}^{t})-x_{t|t}^\theta (\mathbf {s}_1^t)$ and $\overline {\Omega }_{t|t}^\theta (\mathbf {s}_{t-r+2}^{t})-\Omega _{t|t}^\theta (\mathbf {s}_1^t)$ . Then the difference between the approximate and exact score functions (resp. Hessian matrices) comes from the first (resp. second) order derivatives of $\overline {x}_{t|t}^\theta (\mathbf {s}_{t-r+2}^{t})-x_{t|t}^\theta (\mathbf {s}_1^t)$ and $\overline {\Omega }_{t|t}^\theta (\mathbf {s}_{t-r+2}^{t})-\Omega _{t|t}^\theta (\mathbf {s}_1^t)$ . Define

$$ \begin{align*} \dot{\Delta}_{r,t}^{x}(\mathbf{s}_1^t;\theta)\triangleq\nabla_\theta\overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-\nabla_\theta x^\theta_{t|t}(\mathbf{s}_1^t),\qquad \dot{\Delta}_{r,t}^\Omega(\mathbf{s}_1^t;\theta)\triangleq\nabla_\theta\overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-\nabla_\theta\Omega^\theta_{t|t}(\mathbf{s}_1^t),\\ \ddot{\Delta}_{r,t}^{x}(\mathbf{s}_1^t;\theta)\triangleq\nabla_\theta^2\overline{x}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-\nabla_\theta^2 x^\theta_{t|t}(\mathbf{s}_1^t),\qquad \ddot{\Delta}_{r,t}^\Omega(\mathbf{s}_1^t;\theta)\triangleq\nabla_\theta^2\overline{\Omega}^\theta_{t|t}(\mathbf{s}_{t-r+2}^{t})-\nabla_\theta^2\Omega^\theta_{t|t}(\mathbf{s}_1^t). \end{align*} $$

Their upper bounds are shown in Lemma 7.

Lemma 7. Assume Assumptions 1–5. Constant $c>0$ exists such that for any choice of initialization $w_0=\tilde {w}$ , there exist constants $c_{\tilde {M}},c_{\tilde {\tilde {M}}}$ and sequences of positive random variables $M_t$ , $\tilde {M}_t$ , and $\tilde {\tilde {M}}_t$ satisfying for all $t\geq r$ , $\mathbb {E}_0\mathrm {sup}_{\theta \in \Theta } M_t^4<c_{\tilde {M}}$ , $\mathbb {E}_0\mathrm {sup}_{\theta \in \Theta }\tilde {M}_t^4<c_{\tilde {M}}$ , $\mathbb {E}_0\mathrm {sup}_{\theta \in \Theta }\tilde {\tilde {M}}_t^4<c_{\tilde {\tilde {M}}}$ , and for all $\theta \in \Theta $ ,

(53) $$ \begin{align} \max_{\mathrm{s}_1^t\in\mathbb{S}^t}|\dot{\Delta}_{r,t}^x(\mathbf{s}_1^t;\theta)|&\leq c\alpha_{\max}(\theta)^{r-1}\tilde{M}_{t-r+1} +c(r-1)\alpha_{\max}(\theta)^{r-2}\mathbb{1}\{r\geq2\}M_{t-r+1}, \end{align} $$

(54) $$ \begin{align} \max_{\mathrm{s}_1^t\in\mathbb{S}^t}|\dot{\Delta}_{r,t}^\Omega(\mathbf{s}_1^t;\theta)|&\leq c\alpha_{\max}(\theta)^{2(r-1)} +c(r-1)\alpha_{\max}(\theta)^{2r-3}\mathbb{1}\{r\geq 2\}, \end{align} $$

(55) $$ \begin{align} \max_{\mathrm{s}_1^t\in\mathbb{S}^t}|\ddot{\Delta}_{r,t}^x(\mathbf{s}_1^t;\theta)|&\leq c\alpha_{\max}(\theta)^{r-1}\tilde{\tilde{M}}_{t-r+1} +c(r-1)\alpha_{\max}(\theta)^{r-2}\mathbb{1}\{r\geq2\}\tilde{M}_{t-r+1}\nonumber\\ &\qquad+c(r-2)^2\alpha_{\max}(\theta)^{r-3}\mathbb{1}\{r\geq 3\}M_{t-r+1}, \end{align} $$

(56) $$ \begin{align} \max_{\mathrm{s}_1^t\in\mathbb{S}^t}|\ddot{\Delta}_{r,t}^\Omega(\mathbf{s}_1^t;\theta)|&\leq c\alpha_{\max}(\theta)^{2(r-1)} +c(r-1)\alpha_{\max}(\theta)^{2r-3}\mathbb{1}\{r\geq 2\}\nonumber\\ &\qquad+c(r-2)^2\alpha_{\max}(\theta)^{2r-4}\mathbb{1}\{r\geq 3\}. \end{align} $$

Building on Lemma 7, we derive the convergence of the approximate normalized score functions and Hessian matrices to the exact ones.

Theorem 9. Assume Assumptions 1–5. If $r_n$ satisfies $n^{3/2}\alpha _+^{r_n}\rightarrow 0$ as $n\rightarrow \infty $ , then for any $\tilde {w}$ , as $n\rightarrow \infty $ ,

(57) $$ \begin{align} \mathop{\mathrm{sup}}\limits_{\theta\in G}|\frac{1}{\sqrt{n}}\nabla_\theta\log \overline{p}_{r_n,\theta}(y_1,\dots,y_n|w_0=\tilde{w})-\frac{1}{\sqrt{n}}\nabla_\theta\log p_{\theta}(y_1,\dots,y_n|w_0=\tilde{w})|\overset{p}{\rightarrow} 0. \end{align} $$

Moreover, if $r_n$ satisfies $n^3\alpha _+^{r_n}\rightarrow 0$ , then for any initial $\tilde {w}$ , as $n\rightarrow \infty $ ,

(58) $$ \begin{align} \mathop{\mathrm{sup}}\limits_{\theta\in G}|\frac{1}{n}\nabla_\theta^2\log \overline{p}_{r_n,\theta}(y_1,\dots,y_n|w_0=\tilde{w})-\frac{1}{n}\nabla_\theta^2\log p_{\theta}(y_1,\dots,y_n|w_0=\tilde{w})|\overset{p}{\rightarrow} 0. \end{align} $$

Theorem 10 establishes the asymptotic properties of the approximate MLE. It shows consistency in (i) based on Theorems 4 and 8, shows asymptotic normality based on Theorems 7 and 9 in (ii), and gives consistent estimates of the asymptotic covariance matrix based on the Hessian matrix and the outer product of scores in (iii) and (iv), respectively.

7.1 Difference between Collapsing Methods

Section 7.1 examines the difference in the log-likelihood functions with the collapsing steps in (19) and (20). We generate 1,000 datasets from the baseline data-generating process for each of the two sample sizes, $n=200$ and $n=2,000$ . For each dataset, we compute the approximate log-likelihood function based on the simple collapsing (20) $\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta )$ and that with Kim’s collapsing method (19) $\overline {\ell }^{(r)\text {Kim}}_{n,\tilde {w}}(\theta )$ with truncation number $r=2, 4$ .

For $n=2,000$ , Figure 1 plots the kernel density estimates of $\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta )-\overline {\ell }^{(r)\text {Kim}}_{n,\tilde {w}}(\theta )$ with parameters at the specified values and the other parameters at their true values. Figure 1 shows that $|\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta )-\overline {\ell }^{(r)\text {Kim}}_{n,\tilde {w}}(\theta )|$ is generally no larger than 0.02 when $r=2$ and no larger than 0.0001 when $r=4$ . For reference, the averages of $\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta _0)$ across 1,000 simulation experiments are $-5,622.033$ for both $r=2$ and $r=4$ . Thus, the difference in likelihood functions caused by different collapsing methods is small. The results of $n=200$ are presented in Figure B.1 in the Supplementary Material. The difference at $n=200$ is smaller than that for $n=2,000$ , indicating that the difference increases with n and decreases with r. The difference is also large when $(\alpha _1,\alpha _2)$ is large.

Figure 1

$\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta )-\overline {\ell }^{(r)\text {Kim}}_{n,\tilde {w}}(\theta )$ , $n=2,000$ . Note: Each panel plots the kernel densities of $\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta )-\overline {\ell }^{(r)\text {Kim}}_{n,\tilde {w}}(\theta )$ with one parameter at the specified value and the other parameters at their true values. For reference, the average of $\overline {\ell }^{(r)}_{n,\tilde {w}}(\theta _0)$ across 1,000 replications is $-5,622.03$ for both $r=2$ and $r=4$ .

Figure 1 Long description

The figure contains ten panels total, split into two groups of five.

Section a, r equals 2, contains five panels:

1. Top-left, alpha sub 1, alpha sub 2: Three curves centered at 0. The solid line is a sharp peak, while dashed and dotted lines are wider.

2. Top-center, beta sub 1, beta sub 2: Shows a broad distribution on the left and a very sharp, high peak near 0.002 on the right.

3. Top-right, sigma sub 1, sigma sub 2: Three distinct bell curves shifted along the x-axis, which is scaled by 10 super minus 3.

4. Bottom-left, zeta sub 2: Three overlapping bell curves centered near 0.

5. Bottom-center, q sub 1 1, q sub 2 2: Three overlapping bell curves centered near 0, with the x-axis scaled by 10 super minus 3.

Section b, r equals 4, contains five panels with similar parameter labels but different scales:

1. Top-left, alpha sub 1, alpha sub 2: A single extremely sharp vertical spike at 0, with the y-axis reaching 4 times 10 super 6 and x-axis scaled by 10 super minus 5.

2. Top-center, beta sub 1, beta sub 2: A sharp vertical spike at 0 reaching 5 times 10 super 6.

3. Top-right, sigma sub 1, sigma sub 2: A sharp vertical spike at 0 reaching 15 times 10 super 6.

4. Bottom-left, zeta sub 2: Three overlapping bell curves centered near 0, with the y-axis reaching 2.5 times 10 super 5.

5. Bottom-center, q sub 1 1, q sub 2 2: Three overlapping bell curves centered near 0, with the y-axis reaching 4 times 10 super 5 and x-axis scaled by 10 super minus 5.

7.2 Approximation Error, Computation Cost, and Error Bounding Method

This section examines the bias–complexity trade-off of the approximation algorithm and the performance of the error bounding method from Section 5.

We first examine how the approximation error and computation cost vary with r. Using a small sample size of $n=20$ , we compute the exact likelihood function, which becomes computationally infeasible for larger n. We generate 1,000 datasets from the baseline data-generating process with $n=20$ . For each dataset, we compute the approximation error for $r = 2,\ldots ,5$ at $(\alpha _1,\alpha _2)=(0.3,0.4)$ and $(0.95,0.99)$ , and record the computation time. Figures 2 and 3 plot the average approximation error and its logarithm across 1,000 replications. The logarithm of the error is linear in r, indicating that the approximation error decreases exponentially with r. Figure 4 plots the average computation time and its logarithm. The logarithm of computation time increases almost linearly with r, indicating that the computation time increases exponentially with r.

Figure 2

Approximation errors and error bounds, $n=20$ . Note: The figure displays the average approximation errors and the error bounds across 1,000 simulated datasets for $n=20$ . The left and right panels are evaluated at $(\alpha _1,\alpha _2)=(0.3,0.4)$ and $(0.95,0.99)$ , respectively. The other parameters are fixed at their true values.

Figure 2 Long description

Two panels show data for r on the x-axis ranging from 2 to 5.

Left panel: Titled (alpha sub 1, alpha sub 2) = (0.3, 0.4). The y-axis is scaled by 10 super minus 5, ranging from 0 to 3. Three lines (bound of error r-bar = r, bound of error r-bar = 5, and error) are nearly identical. They show a sharp linear decrease from approximately 2.6 at r = 2 to near zero at r = 3, remaining flat through r = 5.

Right panel: Titled (alpha sub 1, alpha sub 2) = (0.95, 0.99). The y-axis ranges from 0 to 0.02.

- The dashed line (bound of error r-bar = r) starts at 0.017 at r = 2, drops to 0.002 at r = 3, and reaches zero by r = 4.

- The dotted line (bound of error r-bar = 5) starts at 0.0115 at r = 2 and drops to 0.001 at r = 3.

- The solid line (error) starts lowest at 0.01 at r = 2 and converges with the dotted line at r = 3, both reaching zero by r = 4.

Figure 3

Log approximation error, $n=20.$ Note: The figure displays the average logarithm of approximation errors across 1,000 simulated datasets for $n=20$ . The left and right panels are evaluated at $(\alpha _1,\alpha _2)=(0.3,0.4)$ and $(0.95,0.99)$ , respectively. The other parameters are fixed at their true values.

Figure 3 Long description

The two panels share a common x-axis labeled r, ranging from 2 to 5. The y-axis represents the log approximation error.

* The left panel is titled (alpha sub 1, alpha sub 2) = (0.3, 0.4). The y-axis ranges from -25 to -10. A single solid black line shows a steep linear decrease, starting near -11 at r = 2 and ending near -24 at r = 5.

* The right panel is titled (alpha sub 1, alpha sub 2) = (0.95, 0.99). The y-axis ranges from -14 to -4. A single solid black line shows a linear decrease, starting near -5.5 at r = 2 and ending at -13 at r = 5.

Figure 4

Computation time and its logarithm, $n=20$ . Note: The figure displays the average computation time and its logarithm across 1,000 simulated datasets for $n=20$ and $(\alpha _1,\alpha _2)=(0.3,0.4)$ . The computation time for $(\alpha _1,\alpha _2)=(0.95,0.99)$ is almost identical and is omitted.

Figure 4 Long description

The image contains two panels arranged horizontally.

Left Panel: A line graph titled computation time, n equals 20. The y-axis is labeled seconds and scaled by 10 super minus 5, with tick marks from 0 to 1. The x-axis is labeled r and ranges from 2 to 7. The data shows a non-linear upward trend. The curve starts near 0 at r equals 2, rises slowly to approximately 0.24 at r equals 5, and then increases more steeply to approximately 0.8 at r equals 7.

Right Panel: A line graph titled log of computation time, n equals 20. The y-axis represents the logarithmic scale with tick marks from minus 15 to minus 11. The x-axis is labeled r and ranges from 2 to 7. The data shows a positive linear-like trend. The line starts at minus 15 for r equals 2, reaches approximately minus 14 at r equals 3, and continues a steady ascent to approximately minus 11.7 at r equals 7.

Next, we use the same simulated data to evaluate the method for bounding the approximation error. We consider two choices of $\overline {r}$ : $\overline {r}=r$ and $\overline {r}=5$ . Figure 2 plots average approximation errors and error bounds across the 1,000 datasets. The bounds with both $\overline {r}=r$ and $\overline {r}=5$ lie slightly above the errors, indicating that the bounding method performs well. Moreover, $\overline {r}=5$ provides a tighter upper bound.

7.3 Choice of r

We illustrate our procedure for selecting r. We simulate one dataset for each sample size $n=400, 1,200, 2,000$ under the baseline data-generating process, calculate the approximation error bounds, and record the computation time for $r=2,\ldots ,5$ and $\overline {r}=r$ and $\overline {r}=5$ . Figure 5 plots the computation time and the approximation error bounds evaluated at $(\alpha _1,\alpha _2)=(0.95,0.99)$ . Because the plots for $\overline {r}=r$ essentially overlap with those for $\overline {r}=5$ , we only plot the case of $\overline {r}=r$ . For reference, the log-likelihood values for $r=2$ when $n=400,1,200$ , and $2,000$ are $-1,165.5$ , $-3,592.5$ , and $-5,789.1$ , respectively. If the maximum tolerable error for a practitioner is 0.1, one may choose $r=2$ for $n=400$ and $r=3$ for $n=1,200$ and 2,000, provided that the computation times shown in the right panel of Figure 5 are acceptable.

Figure 5

Approximation error bounds and computation time. Note: The left and right panels display the error bounds and computation time, respectively, for $r=2,\ldots ,5$ , $\overline {r}=r$ , and $n=400, 1,200, 2,000$ . The error bounds are evaluated at $(\alpha _1,\alpha _2)=(0.95,0.99)$ and the remaining fixed at their true values.

Figure 5 Long description

The image consists of two side-by-side line graphs. Both graphs share a horizontal x-axis labeled r with integer values from 2 to 5. Each graph contains three lines representing different values of n: a solid line for n equals 400, a dashed line for n equals 1,200, and a dotted line for n equals 2,000.

* Left Panel: Titled Error bounds. The vertical y-axis ranges from 0 to 0.5. All three lines show a sharp non-linear decrease. At r equals 2, the error is highest, with the dotted line near 0.45 and the solid line near 0.1. By r equals 3, all lines drop significantly toward zero, remaining nearly flat and overlapping at the baseline for r values 4 and 5.

* Right Panel: Titled Computation time. The vertical y-axis ranges from 0 to 0.08. All three lines show a non-linear increase as r increases. The dotted line for n equals 2,000 shows the steepest upward trend, reaching approximately 0.07 at r equals 5. The dashed line for n equals 1,200 shows a moderate increase to approximately 0.045. The solid line for n equals 400 shows the slowest growth, ending near 0.015 at r equals 5.

7.4 Comparison of Kalman–Collapsing Approximation Algorithm and Auxiliary Particle Filter

This section uses simulation to compare the performance of the auxiliary particle filter and the Kalman–collapsing algorithm.

We generate 1,000 datasets, each of sample size $n=20$ using the baseline data-generating process. For this small sample size, exact likelihood functions remain computationally tractable. For each dataset, we compute: (i) approximation error $|\overline {\ell }_{n,\tilde {w}}^{(r)}(\theta )-\ell _{n,\tilde {w}}(\theta )|$ for $r=2,\ldots ,5$ and (ii) the auxiliary particle filter error $|\ell _{n,\tilde {w}}^{\text {APF}}(\theta )-\ell _{n,\tilde {w}}(\theta )|$ with particle size M and resample size R set to $M=R=10{,}000, 50{,}000, 100{,}000$ . All likelihood functions are evaluated at $(\alpha _1,\alpha _2)=(0.3,0.4)$ and $(0.95,0.99)$ , with the remaining parameters fixed at their true values. Computation times are also recorded for all methods.

For both algorithms, smaller errors are achieved at the cost of increased computation time. To facilitate comparison, Figure 6 plots the average log approximation errors against log computation times. The approximation algorithm reduces the error much faster as computation time increases than the auxiliary particle filter. Moreover, it achieves smaller errors with lower computation time than the auxiliary particle filter.

Figure 6

Log error and log computation time, $n=20$ . Note: This figure plots the log errors against the log computation time for the approximation algorithm with $r=2,\ldots ,5$ and the auxiliary particle filter (APF) with particle size M and resample size R set to $M=R=10,000, 50,000$ , and 100,000. Solid and dashed lines correspond to $(\alpha _1,\alpha _2)=(0.95,0.99)$ and (0.3,0.4), respectively. The other parameters are fixed at their true values.

Figure 6 Long description

The x-axis represents time in seconds on a log-scale ranging from minus 6 to 6. The y-axis represents absolute error on a log-scale ranging from minus 25 to 0.

Four distinct data series are plotted:

1. Approximation algorithm with alpha sub 1, alpha sub 2 equals 0.95, 0.99: Represented by a solid line with circular markers. It starts at the top-left at r equals 2 and descends diagonally to the right, ending at r equals 5 with an error of approximately minus 13.

2. Approximation algorithm with alpha sub 1, alpha sub 2 equals 0.3, 0.4: Represented by a dashed line with circular markers. It starts lower than the first series at r equals 2 and descends more steeply to the right, ending at r equals 5 with an error of approximately minus 24.

3. A P F error with alpha sub 1, alpha sub 2 equals 0.95, 0.99: Represented by a solid line with square markers located in the top-right quadrant. It shows a very slight downward slope as M increases from 10,000 to 100,000, maintaining a high error near minus 2.

4. A P F error with alpha sub 1, alpha sub 2 equals 0.3, 0.4: Represented by a dashed line with square markers located just below the previous A P F series. It also shows a slight downward slope as M increases, with error values between minus 3 and minus 4.

A legend in the bottom-right quadrant confirms these mappings, distinguishing between the approximation error and A P F error using solid and dashed lines for the different alpha parameter sets.

To assess how the two algorithms may perform when the sample size exceeds 20, Figure 7 displays the average errors at each time period, $|\log \overline {p}_{r,\tilde {w}}^{\theta }(y_t|\mathbf {y}_1^{t-1})-\log p_{\tilde {w}}^{\theta }(y_t|\mathbf {y}_1^{t-1})|$ and $|\log p_{\tilde {w}}^{\theta ,\text {APF}}(y_t|\mathbf {y}_1^{t-1})-\log p_{\tilde {w}}^{\theta }(y_t|\mathbf {y}_1^{t-1})|$ , $t=8,\ldots ,20$ , across 1,000 replications. The approximation errors remain largely stable over time, consistent with Proposition 1. In contrast, the auxiliary particle filter errors exhibit an upward trend as t increases. If this pattern persists for larger sample sizes, the Kalman–collapsing algorithm is expected to outperform the auxiliary particle filter as n grows.

Figure 7

APF and approximation errors at each period, $t=8,9,\ldots ,20$ . Note: The upper panel plots the average of the approximation errors of the Kalman-collapsing filter $|\log \overline {p}^{\theta }_{r,\tilde {w}}(y_t|\mathbf {y}_1^{t-1})-\log p^{\theta }_{\tilde {w}}(y_t|\mathbf {y}_1^{t-1})|$ across 1,000 datasets for $t = 8, 9, \ldots , 20$ . The lower panel plots the average of the computation errors of the auxiliary particle filter $|\log p^{\theta ,\text {APF}}_{\tilde {w}}(y_t|\mathbf {y}_1^{t-1})-\log p^{\theta }_{\tilde {w}}(y_t|\mathbf {y}_1^{t-1})|$ across 1,000 datasets for $t = 8, 9, \ldots , 20$ with the resampling size M and the particle size R being 10,000, 50,000, and 100,000.

Figure 7 Long description

The image consists of four line graphs arranged in a two-by-two grid. All graphs share an x-axis representing time periods from 8 to 20.

Top-left panel: Titled times 10 super minus 6 approximation error with alpha sub 1, alpha sub 2 equals 0.3, 0.4. The y-axis ranges from 0 to 4. A solid line for r equals 2 fluctuates between 2.5 and 3.2. Dashed and dotted lines for r equals 3 and r equals 4 remain nearly flat at zero.

Top-right panel: Titled times 10 super minus 3 approximation error with alpha sub 1, alpha sub 2 equals 0.95, 0.99. The y-axis ranges from 0 to 2.5. The solid line for r equals 2 fluctuates between 1.4 and 1.9. The dashed line for r equals 3 is flat near 0.1, and the dotted line for r equals 4 is flat at zero.

Bottom-left panel: Titled A P F error with alpha sub 1, alpha sub 2 equals 0.3, 0.4. The y-axis ranges from 0 to 0.2. Three lines show an upward trend. The solid line for M equals R equals 10,000 rises from 0.09 to 0.15. The dashed line for M equals R equals 50,000 rises from 0.05 to 0.08. The dotted line for M equals R equals 100,000 rises from 0.03 to 0.06.

Bottom-right panel: Titled A P F error with alpha sub 1, alpha sub 2 equals 0.95, 0.99. The y-axis ranges from 0.1 to 0.6. All lines show a steady linear increase. The solid line for M equals R equals 10,000 rises from 0.2 to 0.38. The dashed line for M equals R equals 50,000 rises from 0.14 to 0.29. The dotted line for M equals R equals 100,000 rises from 0.13 to 0.26.

7.5 Finite-Sample Distribution of Approximate MLE

This section uses simulation to examine how well the asymptotic distribution approximates the distribution of the MLE in finite samples. We generate 1,000 datasets from the baseline data-generating process for $n=400$ , $1,000$ , and $1,600$ . For each dataset, we estimate the approximate MLE with $r=2, 6$ . Using Theorem 10, we construct the confidence intervals from the Hessian matrix and the outer product of the score. Table 1 reports the frequency at which the 95% confidence interval contains the true parameter.Footnote ⁵

Table 1

Coverage frequency of the asymptotic 95% confidence intervals

Table 1 Long description

The table is divided into two sections. Panel A: Confidence interval constructed from the outer product of score. Panel B: Confidence interval constructed from the Hessian. Both panels share the same column headers: n, r, alpha sub 1, alpha sub 2, beta sub 1, beta sub 2, sigma sub 1, sigma sub 2, varsigma sub 2, q sub 1 1, and q sub 2 2.

In Panel A, for n = 400 and r = 2, values range from 0.930 to 0.981. For n = 1,000 and r = 2, values range from 0.933 to 0.975. For n = 1,600 and r = 2, values range from 0.936 to 0.967. Increasing r to 6 shows minimal variation in coverage frequency.

In Panel B, for n = 400 and r = 2, values show more variance, particularly for beta sub 1 at 0.774, while others range between 0.896 and 0.935. As n increases to 1,600, the coverage frequencies stabilize closer to the 0.95 target, with values ranging from 0.921 to 0.954.

The coverage frequencies in Panel A are closer to 95% than the ones in Panel B, especially for $n=400$ . The confidence intervals constructed from the outer product of score outperform those constructed from the Hessian matrix in small samples. The coverage frequencies with $r=2$ and $r=6$ are close. This might be because the approximate likelihood functions for $r=2$ and $r=6$ are close for this experiment. To further verify it, we compute $|\hat {\theta }_n^{(6)}-\hat {\theta }_n^{(2)}|$ for each replication and report the median in Table 2.

Table 2

Median of $|\hat {\theta }_n^{(6)}-\hat {\theta }_n^{(2)}|$

Table 2 Long description

The table consists of 10 columns and 4 rows including the header.

Header row labels from left to right: n, alpha sub 1, alpha sub 2, beta sub 1, beta sub 2, sigma sub 1, sigma sub 2, varsigma sub 2, q sub 11, and q sub 22.

Row 1: n is 400. Values are 0.0004, 0.0002, 0.0008, 0.0013, 0.0011, 0.0038, 0.0029, 0.00003, and 0.00003.

Row 2: n is 1,000. Values are 0.0002, 0.0001, 0.0006, 0.0010, 0.0007, 0.0028, 0.0013, 0.00002, and 0.00002.

Row 3: n is 1,600. Values are 0.0002, 0.0001, 0.0004, 0.0006, 0.0004, 0.0016, 0.0007, 0.00001, and 0.00001.

As n increases, the values for all parameters generally decrease.

Table 2 shows that $\hat {\theta }_n^{(6)}$ and $\hat {\theta }_n^{(2)}$ are close. Moreover, given that the coverage frequencies in Panel A for $r=2$ are mostly close to 95%, the estimator with $r=2$ and the confidence intervals constructed from the outer product of scores perform well in finite samples. Further increases in r do not lead to noticeable improvements for the sample sizes typically used in practice. These finite-sample findings are consistent with the theoretical requirement in Theorem 10, where asymptotic normality is established under the condition $n^3\alpha _+^{r_n}\to 0$ . This condition permits $r_n$ to diverge at a relatively slow rate and is satisfied, for example, by $r_n=C\log n$ with $C> 3/|\log \alpha _+|$ . Hence, a small value such as $r=2$ suffices in practice, which corroborates the good simulation results reported above.Footnote ⁶