Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
34,280 characters · 7 sections · 16 citation commands
Recursive Score and Hessian Computation in Regime-Switching Models
Keywords: Regime-switching model, Hessian, Score, Inference
JEL codes: C22, C12
Regime-switching models have been widely applied since hamilton1989new to analyze how time-series patterns shift across different latent economic states. These models are typically estimated via maximum likelihood, with the computation of the score vector and the Hessian matrix playing a central role in producing reliable standard errors and enabling hypothesis tests. Although hamilton1996specification provides an algorithm for computing the score using smoothed probabilities, a direct and generally accepted method for calculating the Hessian matrix remains lacking in the literature. At the same time, inference based on the Hessian matrix is common practice. Consequently, applied research often resorts to numerical approximations for inference, as seen in studies by psaradakis1998finite and pouzo2022maximum.
The algorithm for smoothed functionals in cappe2005inference can, in principle, be used to compute both the score and the Hessian. However, it remains unfamiliar to most researchers for several reasons. First, its illustration has been limited to a simple one-parameter model, and although generalizable, the complete algorithm has not been laid out in a general setting. Implementing it requires a preliminary step of applying Fisher’s and Louis’ identities (cappe2005inference) to express the score as the conditional expectation of the complete score and to relate the Hessian to that of the complete-data log-likelihood. These expressions must then be reformulated into the specific form required by the smoothing functional algorithm (cappe2005inference). For a score vector of dimension $k\times 1$ and a Hessian of size $k\times k$, researchers must derive smoothing functionals for each element, a process that can become cumbersome and tedious.
The contributions of this study are twofold. First, we propose a unified recursive framework to compute the score and Hessian matrix simultaneously. This framework is designed for a general regime-switching specification that encompasses many useful models in the empirical literature as special cases. Unlike the approach in hamilton1996specification, our algorithm does not require computing or storing smoothed probabilities and can be integrated directly into the prediction and update steps of likelihood evaluation. It is straightforward to implement, as it avoids the need for preliminary transformations via Fisher’s or Louis’ identities and works directly with derivatives of the period likelihood function. Researchers need only supply the first and second derivatives of certain quantities, and our algorithm returns the score and Hessian, significantly reducing pre-computation effort compared to the smoothed functional approach. Second, to our knowledge, no thorough investigation has been conducted to compare inference based on the exact score versus the exact Hessian in this context. We supplement our methodological contribution with simulation evidence to examine their finite-sample properties.
{\bf Notations: }let $\mathbf{Y}_m^n\triangleq(Y_{m},\ldots,Y_n)'$ for $m\leq n$, and let $\mathbf{Y}_m^n$ be an empty set for $m>n$. Similarly, define $\mathbf{S}_m^n$ and $\mathbf{X}_m^n$. We use $\mathbb{N}(\mu,\sigma^2)$ to denote the density of the normal distribution with mean $\mu$ and variance $\sigma^2$.
The conditional distribution of the observable process $Y_t$ is governed by a latent regime process $S_t$, and it admits a density
where $p,q\geq0$, $X_t$ is a predetermined variable (vector), $g_\theta$ is a family of density functions indexed by $\theta$. The regime $S_t$ is a first-order Markov chain on $\mathbb{S}=\{1,2,\ldots,J\}$ with transition probability $q_{ij}=q_\theta(S_t=j|S_{t-1}=i)$ for $i,j=1,\ldots,J$. Given $(\mathbf{S}_{t-\max\{p,1\}}^{t-1},\mathbf{Y}_{t-q}^{t-1},X_t)$, $(Y_{t},S_{t})$ is independent of $(\{Y_k\}_{k\leq t-q-1},\{S_k\}_{k\leq t-\max\{p,q\}},\{X_k\}_{k\leq t-1})$, and the conditional density can be written as
For short notation, let $\overline{\mathbf{S}}_t\triangleq \mathbf{S}_{t-\max\{p,1\}+1}^{t}$. We abbreviate the dependence on $Y_{t}$, $\overline{\mathbf{Y}}_{t-1}$, and $X_t$ into the subscript $t$. Then ((ref)) can be expressed as $\mathrm{f}_{\theta,t}(S_t,\overline{\mathbf{S}}_{t-1})$. Its first and second derivatives are expressed as $\nabla_\theta \mathrm{f}_{\theta,t}(S_t,\overline{\mathbf{S}}_{t-1})$ and $\nabla_\theta^2 \mathrm{f}_{\theta,t}(S_t,\overline{\mathbf{S}}_{t-1})$.
{\bf Example.} (Autoregressive model with regime switching) Consider the following model:
where $u_t$ is an independent and identically distributed (i.i.d.) sequence of random variables. The parameters $\mu_{s_t}$ and $\sigma_{s_t}$ depend on a latent regime variable $s_t$. $s_t=1$ or 2 and follows a first-order Markov chain. In this example, $p=q=2$, and
This model was employed in psaradakis1998finite to study the finite-sample properties of the maximum likelihood estimator.
This study provides the algorithms to compute the score and hessian matrix of the log-likelihood function given the initial observations $\mathbf{Y}^0_{-q+1}$, $\ell_{n,\nu}(\theta)=\log p_{\theta,\nu}(Y_1,\ldots,Y_n|\mathbf{Y}^0_{-q+1},\mathbf{X}_1^n)$, where $\nu$ is the short notation for the initial distribution of regimes $\nu_\theta(\overline{\mathbf{S}}_0)\triangleq p_\theta(\overline{\mathbf{S}}_0|\mathbf{Y}_{-q+1}^0)$. Examples of the initial distributions include a fixed starting point $\overline{\mathbf{s}}_0$, $\nu(\overline{\mathbf{S}}_0)=\mathbbm{1}\{\overline{\mathbf{S}}_0=\overline{\mathbf{s}}_0\}$. If the Markov chain is ergodic, the initial distribution can also be chosen to be the unconditional distribution. The log-likelihood function can be expressed as
For notational simplicity, we slightly abuse notation and let $\mathrm{f}_{\theta,0}(S_0,\overline{\mathbf{S}}_{-1})$ denote $\nu(\overline{\mathbf{S}}_{0})$, where $\overline{\mathbf{S}}_{-1}\triangleq \mathbf{S}^{-1}_{-\max\{p,1\}+1}$. Then ((ref)) can be expressed as
We use $\mathrm{p}_{\theta,n}$ as a short notation for the likelihood $p_{\theta,\nu}(Y_1,\ldots,Y_n|\overline{\mathbf{Y}}_0,\mathbf{X}_1^n)$. Then $\ell_{n,\nu}(\theta)=\log \mathrm{p}_{\theta,n}$. The score function can be expressed as
where
The Hessian matrix can be expressed as
where $\nabla_\theta^2\mathrm{p}_{\theta,n}=\mathrm{H}_{\theta,n}+\mathrm{h}_{\theta,n}+\mathrm{h}_{\theta,n}^T$ and
The next section presents the algorithm used to compute ((ref)) and ((ref)). Note that the algorithm requires researchers to supply only ((ref)), the initial distribution $\nu_\theta(\overline{\mathbf{S}}_t)$, and their first and second derivatives. Unlike the smoothed functional approach of hamilton1996specification and cappe2005inference, our method does not require any preliminary transformations using Fisher's or Louis' identities. This significantly reduces the pre-computation effort.
We can compute the score and Hessian matrix if we know how to compute ((ref))--((ref)) and ((ref))--((ref)). ((ref)) can be computed as ((ref)) by replacing $\nabla_\theta \mathrm{f}_{\theta,t}(S_t,\overline{\mathbf{S}}_{t-1})$ with $\nabla_\theta^2 \mathrm{f}_{\theta,t}(S_t,\overline{\mathbf{S}}_{t-1})$. Thus, the main challenge is to compute ((ref)), ((ref)), and ((ref)) efficiently. To explain the idea of the algorithm, we first consider the simple model where there is no regime switching in Section (ref).
In the simple model without regime switching, ((ref)), ((ref)), and ((ref)) are simplified to
They can be computed according to the following algorithm.
To see why we update as ((ref)), notice that each summand in ((ref)) is the product of $t+1$ terms, with one being the the derivative and the rest being the period likelihood. The summand in $\mathrm{s}_{\theta,t-1}$ already contains one derivative, so we multiply it with the period likelihood function $\mathrm{f}_{\theta,t}$, which brings the first term in ((ref)). We add the second term in ((ref)) to include the term with derivative at time $t$. To write it more specifically,
((ref)) is updated with the same reasoning. Each summand in ((ref)) is the product of $t+1$ terms, with two being the derivatives and the rest being the period likelihood. The summand in $\mathrm{h}_{\theta,t-1}$ already contains two derivatives, so we multiply it with $\mathrm{f}_{\theta,t}$, which is the first term in ((ref)). For the second term in ((ref)), the summand in $\mathrm{s}_{\theta,t-1}$ contains one derivative, so we multiply it with $\nabla_\theta^T\mathrm{f}_{\theta,t}$. Next subsection explains the algorithm to compute the score and Hessian matrix for regime-switching models.
Section (ref) explains the algorithm to compute ((ref)), ((ref)), and ((ref)) in the model without regime switching. To obtain the complete algorithm for the model with regime switching, we need to make two modifications. First, for models with regime switching, we need to take summation over regimes at appropriate steps. For each regime $S_t$, it is involved in $\mathrm{f}_{\theta,k}(S_{k},\overline{\mathbf{S}}_{k-1})$, $t\leq k\leq t+p$, but not in $\mathrm{f}_{\theta,\ell}(S_{\ell},\overline{\mathbf{S}}_{\ell-1})$, $\ell\geq t+p+1$. Thus, we can take summation over $S_t$ at the recursion step of $t+p$. Second, in ((ref)) and ((ref)), the computation of the denominator $\mathrm{p}_{\theta,n}=p_{\theta,\nu}(Y_1,\ldots,Y_n|\overline{\mathbf{Y}}_0,\mathbf{X}_1^n)$ might suffer from numerical overflow or underflow.\footnote{To see this more clearly, note that under suitable regularity conditions, the normalized log-likelihood satisfies $n^{-1}\log p_{\theta,\nu}(Y_1,\ldots,Y_n\mid\overline{\mathbf{Y}}_0,\mathbf{X}_1^n)\overset{p}{\longrightarrow}\ell(\theta)$ in the consistency proof; see douc2004asymptotic and kasahara2019asymptotic. Consequently, for large $n$, the likelihood itself behaves as $p_{\theta,\nu}(Y_1,\ldots,Y_n\mid\overline{\mathbf{Y}}_0,\mathbf{X}_1^n)\approx\exp\big(n\ell(\theta)\big)$. If $\ell(\theta)<0$, then $\exp(n\ell(\theta))$ decays exponentially to zero. Once its magnitude falls below the machine epsilon, floating‑point software will round it to zero, resulting in numerical underflow. Conversely, if $\ell(\theta)>0$, the term grows exponentially and may exceed the largest representable floating‑point number, causing numerical overflow (stored as INF).} Thus, instead of making a division in the last step, we re-scale each value with the period likelihood at each recursion step. From $\mathrm{p}_{\theta,t}=\prod_{k=1}^tp_{\theta,\nu}(Y_k|\overline{\mathbf{Y}}_0^{k-1},\mathbf{X}_1^k)$, we can re-scale at each recursion with $p_{\theta,\nu}(Y_k|\overline{\mathbf{Y}}_0^{k-1},\mathbf{X}_1^k)$. The complete algorithm can be described as follows.
Unlike the smoothed functional approach in hamilton1996specification, our algorithm is recursive and does not require the computation and storage of the smoothed probabilities. The connection between our algorithm and the standard prediction-and-update steps for evaluating the likelihood function can be seen as follows. For $t\geq 2$, $\mathpzc{p}_{\theta,t}(\overline{\mathbf{S}}_{t}) =p_{\theta,\nu}(Y_t,\overline{\mathbf{S}}_{t}|\overline{\mathbf{Y}}_0^{t-1},\mathbf{X}_1^t)$ $\sum_{\overline{\mathbf{S}}_{t}}\mathpzc{p}_{\theta,t}(\overline{\mathbf{S}}_{t}) =p_{\theta,\nu}(Y_t|\overline{\mathbf{Y}}_0^{t-1},\mathbf{X}_1^t)$ are the terms computed in the prediction step. $\widetilde{\mathpzc{p}}_{\theta,t}(\overline{\mathbf{S}}_{t})=p_{\theta,\nu}(\overline{\mathbf{S}}_{t}|\overline{\mathbf{Y}}_0^{t},\mathbf{X}_1^t)$ is the term computed in the updating step. Thus, the recursion in this algorithm extends the standard prediction and update algorithm of likelihood evaluation to compute the score and the Hessian matrix simultaneously.
kasahara2019asymptotic established the asymptotic normality of the MLE in regime-switching models and showed that the asymptotic covariance matrix of $\sqrt{n}(\hat{\theta}_{n,\nu}-\theta_0)$ can be consistently estimated using either the Hessian-based estimator
or the outer product of the score
This section compares these two methods for estimating the asymptotic covariance matrix. We consider a comprehensive list of autoregressive models with regime switching. The data-generating process is the fully regime-switching autoregression
and its seven nested variants obtained by restricting any subset of $(\mu,\phi,\sigma)$ remain regime-invariant. The resulting models are denoted as $M_{\mu\phi\sigma}$, $M_{\phi\sigma}$, $M_{\mu\sigma}$, $M_{\mu\phi}$, $M_{\mu}$, $M_{\phi}$ and $M_{\sigma}$, where the subscripts indicate which parameters switch. For instance, the model $M_{\mu\phi\sigma}$ is ((ref)). The model $M_{\mu\sigma}$ is ((ref)). psaradakis1998finite examined the performance of the inference based on the numerical Hessian matrix for model $M_{\mu\sigma}$. According to their simulation results, the asymptotic approximation is worse if $\phi,q_{11},$ and $q_{22}$ are close to unity. Thus, we consider the following parameters when the asymptotic distribution is not favored. $q_{11}=q_{22}=0.95$. For the parameters that are allowed to switch,
when the parameters do not switch,
For each model, we generate 1000 data sets of sample sizes $n = $100, 200, 400, and 800.\footnote{We simulate $(800 + n)$ periods and use the last $n$ observations as our sample, so that the initial value of the simulated data set is approximately drawn from the stationary distribution.} For each data set, we compute the MLE using the Broyden-Fletcher-Goldfarb-Shanno optimization method and compute the standard errors based on ((ref)) and ((ref)).
We report the coverage frequencies of the 95% confidence intervals in Table (ref) and the ratio of the standard deviation of the MLE to the median of the standard errors across 1000 datasets in Table (ref) for $M_\mu$ and $M_\sigma$. The asymptotic normality provides a better approximation in $M_\mu$ than $M_{\sigma}$. More specifically, the coverage frequencies are closer to the nominal level of 95% and the ratios of sampling standard deviations to the estimated standard errors are closer to unity in the model $M_\mu$ than $M_\sigma$. The performance in the remaining models falls between these two cases. To save space, we present only results for $M_\mu$ and $M_\sigma$. Full results are available upon request. When $n=100$, the performance of the asymptotic normal distribution is somewhat discouraging. The coverage frequencies are lower than the nominal level of 95%, and the estimated standard errors are smaller than the sampling standard deviations of the MLE. When $n$ gets larger to 800, the performance improves. Moreover, inference based on the outer product of the score ((ref)) outperforms that based on ((ref)), in the sense that coverage frequencies are closer to 95% and the ratio of standard deviations to standard errors is closer to unity. This pattern holds across all models examined.
This study proposes a recursive algorithm for computing the score vector and Hessian matrix in regime-switching models, without relying on pre-computation via Fisher’s or Louis’ identities. Simulation results show that asymptotic approximations based on the outer product of the score are favored over those derived from the Hessian matrix. While the presented framework is general, it does not encompass more complex specifications such as regime-switching GARCH models or regime-switching state-space models, where likelihood evaluation requires approximation. Extending the algorithm to compute score and Hessian matrix for such models will be addressed in future work.