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.
76,826 characters · 22 sections · 52 citation commands
Sequential Bayesian Learning for Hidden Semi-Markov Models
\newacronym{ssm}{SSM}{State Space Model}
\newacronym{mc}{MC}{Monte Carlo} \newacronym{ibis}{IBIS}{Iterated Batch Importance Sampling} \newacronym{smc}{SMC}{Sequential Monte Carlo} \newacronym{smc2}{SMC2}{Sequential Monte Carlo Squared}
\newacronym{mcmc}{MCMC}{Markov Chain Monte Carlo} \newacronym{hmc}{HMC}{Hamiltonian Monte Carlo} \newacronym{nuts}{NUTS}{No U-Turn Sampling} \newacronym{pmh}{PMH}{Particle Metropolis Hastings} \newacronym{pgibbs}{PGIBBS}{Particle Gibbs} \newacronym{pgas}{PGAS}{Particle Gibbs with ancestors sampling} \newacronym{pmcmc}{PMCMC}{Particle MCMC}
\newacronym{pf}{PF}{Particle Filter} \newacronym{cpf}{CPF}{Conditional Particle Filter} \newacronym{bpf}{BPF}{Bootstrap Particle Filter} \newacronym{apf}{APF}{Auxiliary Particle Filter}
\newacronym{hmm}{HMM}{Hidden Markov Model} \newacronym{hsmm}{HSMM}{Hidden Semi-Markov Model} \newacronym{edhmm}{EDHMM}{Explicit-duration Hidden Markov Model} \newacronym{arhsmm}{ARHSMM}{Autoregressive Hidden Semi-Markov Model} \newacronym{sv}{SV}{Stochastic Volatility}
\newacronym{hsmm-egarch}{HSMM-eGARCH}{Hidden Semi-Markov Model with eGarch volatility dynamics}
\newacronym{svm}{SVM}{Stochastic Volatility Model} \newacronym{mjm}{MJM}{Markov Jump Model} \newacronym{garch}{GARCH}{Generalized Autoregressive Conditional Heteroskedasticity} \newacronym{egarch}{eGARCH}{Exponential Generalized Autoregressive Conditional Heteroskedasticity}
\newacronym{crps}{CRPS}{Continuous Ranked Probability Score} \newacronym{bmi}{BMI}{Basic Marginal Likelihood Identity} \newacronym{is}{IS}{Importance Sampling} \newacronym{sis}{SIS}{Sequential Importance Sampling} \newacronym{sir}{SIR}{Sequential Importance Resampling} \newacronym{ci}{CI}{Credible Interval} \newacronym{pl}{PL}{Predictive Likelihood} \newacronym{bf}{BF}{Bayes Factor} \newacronym{clpbf}{CLPBF}{Cumulative Log Predictive Bayes Factor} \newacronym{ar}{AR}{Autoregressive} \newacronym{ess}{ESS}{Effective Sample Size} \newacronym{pacf}{PACF}{Partial Autocorrelation Function} \glsaddall
Discrete State Space Models (SSM) provide a flexible class of models with applications in ecology, economics, finance, robotics and signal processing Bulla06, Lindsten13, Papaspiliopoulos20, corenflos21, among others. They can handle structural breaks, shifts, or time-varying parameters and still have an interpretable structure. Moreover, such models are generative and allow for multi-step forecasting. However, analytical forms of the likelihood function are only available in special cases, and standard parameter optimization routines are often challenging to implement. Given the observed data $e_{1:T}=(e_1,\dots,e_T)$ and parameter $\theta$, the major challenge in the estimation of \acrshort{ssm}s is thus the generally intractable likelihood function $p_{\theta}(e_{1:T})$, which integrates over the latent state trajectory $s_{1:T}$ such that $p_{\theta}(e_{1:T}) = \int p_{\theta}(e_{1:T}, s_{1:T}) d s_{1:T}$.
A flexible discrete \acrshort{ssm} on which we focus in this paper is known as \acrlong{hsmm}. \acrshort{hsmm}s have a flexible state duration distribution, well suited for processes that remain in any particular state for an extended period of time, and can be considered as generalizations of the well-known basic \acrlong{hmm} introduced in Baum66. \acrshort{hsmm}s have been employed in ecology, epidemiology, finance Bulla06, pohle21, visani21 and many other fields Yu16, but are typically used more sporadically than their standard \acrshort{hmm} counterparts because the likelihood function is significantly more costly to evaluate. In the \acrshort{hmm} case, the likelihood has computational complexity of $\mathcal{O}(K^2T)$, where K = number of latent states, T = number of data points, see Baum66. For the \acrshort{hsmm}, this is a much more expensive operation of order $\mathcal{O}(K^2(d_{max} - d_{min})^2T)$ Murphy02, Wood12, where $d_{min}$ and $d_{max}$ denote the minimal and maximal state duration in a latent regime. In practice, $(d_{max} - d_{min}) >> K$, as described in more detail in Section (ref), which often leads to computationally expensive inference algorithms. Such considerations have led to the use of approximate methods in applications where \acrshort{hsmm}s provide valuable models, see for example HadjAmar22 and Xiao18.
In this paper, we follow an alternative route aiming to construct efficient computational schemes that operate on the joint space of latent states and parameters using \acrfull{smc} methods that are exact, in the sense that they are only amenable to Monte Carlo error. The fundamental building block of the proposed schemes is the \acrfull{pf}, see for example Doucet11 and the references therein. Traditionally, \acrshort{smc} samplers such as \acrshort{pf}s have been used to estimate the underlying state sequence of \acrshort{ssm}s, while standard \acrfull{mcmc} samplers facilitate Bayesian inference for the model parameters. More recently, combining these methods is becoming increasingly popular, see daviet18 and buchholz2020. A natural computational framework that jointly infers the latent state sequence and model parameter is known as \acrfull{pmcmc} Andrieu10, andrieu2009. To our knowledge, \acrshort{pmcmc} has not been used for \acrshort{hsmm}s, so we work within this framework aiming to construct an efficient implementation. In particular, we focus on the \acrfull{pgibbs} version to implement parameter updates, conditional on the latent state trajectory, via Hamiltonian \acrshort{mcmc} neal12 variants. Given that \acrshort{ssm}s are typically used in applications with data of sequential nature, it is essential to explore techniques where previous parameter estimates can be reused once the data is updated. An example for that is provided by the \acrfull{smc2} algorithm, introduced in chopin2012, which can be viewed as an extension of the main \acrshort{smc} framework of Chopin02 and DelMoral06; see also Dai20 for some recent work that includes a survey of applications in different contexts. Other similar approaches include fearnhead10 and crisan17. More information on these methods is provided in the Section (ref).
The major motivation of this paper is thus to develop methods to estimate \acrshort{hsmm}s (1) in a computationally feasible time and (2) in a sequential manner. The contribution of this paper is two-fold: First, we offer \acrlong{smc} schemes on \acrlong{hsmm}s by tailoring ideas from Andrieu10 and chopin2013smc2 for batch and sequential estimation. This offers several benefits over standard deterministic filtering techniques, including computational efficiency. The developed \acrshort{smc} schemes can also facilitate Bayesian model choice and assessment of predictive performance in an efficient manner. Second, we propose a novel class of models by linking \acrlong{ar}-type models with \acrshort{hsmm}s to better describe data consisting of financial and econometric time series. Sequential estimation of such models is particularly important as \acrshort{ar} \acrshort{hsmm}s have the potential to detect substantial changes in the data, which we illustrate in a case study on data that evolves rapidly during the Covid-19 pandemic.
The paper is organized as follows: Section (ref) formally introduces \acrshort{hsmm}s via a suitable formulation to apply sequential Monte Carlo methods such as \acrlong{pf}ing. It also provides justification for the use of \acrshort{pf}s instead of deterministic filtering techniques. In Section (ref), the developed methodology of this paper is presented, which includes the model choice criteria available from by-products of the estimation process. Section (ref) explores the performance of the developed methods via simulation based experiments. In Section (ref), we focus on the performance of \acrshort{hsmm}s, estimated with the developed methodology of this paper, on real-world applications such as financial time series of the VIX index. Comparisons of different \acrshort{hsmm}s as well as benchmark \acrshort{hmm}s are conducted. Model selection, and in particular choice of the number of states using \acrshort{smc2}, is also put into test. Finally, Section (ref) concludes with some relevant discussion.
A standard \acrlong{hmm} may be specified via a bivariate stochastic process $ \{ e_t, s_t \}_{t = 1,2,\ldots} $, where $s_t$ is an unobserved Markov chain and $e_t$ is an observed sequence of independent random variables, conditional on $s_t$. The model is fully specified by the transition distribution $f_\theta$, $s_t \sim f_\theta(s_t \mid s_{t-1})$, $t \geq 2$, the corresponding initial distribution $\pi_\theta$, $s_1 \sim \pi_\theta( s_1 )$, and the observation distribution $g_\theta$, $e_t \sim g_\theta(e_t \mid s_{t})$, $t \geq 1$. Directly computing the likelihood function of this model involves summing up over all possible state sequences,
Hence, various filtering techniques have been proposed that take into account the memory of the latent state variable to reduce the computational costs to $\mathcal{O}(K^2T)$. One shortcoming of \acrshort{hmm}s is their explicit distributional assumption regarding the duration in any particular state. To give insight into this issue, we denote $p( s_{t+k} = j,~ s_{t+1:t+k-1} = i \mid s_{t} = i )$, the probability a state remains in any current state until it switches, as state duration distribution. In the \acrshort{hmm} case, this probability is implicitly geometric. Set $p(s_t = i \mid s_{t-1} = i ) = \mathcal{T}_{ii}$, and assume there are only 2 states, then for a homogeneous Markov chain, using the chain rule and the Markov assumption, it holds:
In general, for $t+k$ steps, it holds that
where the geometric distribution has to be interpreted as the length of state duration up to and including the transition to the other state. For processes that tend to stay in any particular state for a long-time horizon, this may be a poor modelling choice. Alternatively, the state duration could be explicitly modelled. A \acrlong{hsmm}, see Murphy02, Yu10, Yu16, is a generalization of an \acrshort{hmm}, which may be viewed as a \acrshort{hsmm} with Geometric state duration distribution. A graph structure and a comparison to the standard \acrshort{hmm} can be seen in figures (ref) and (ref). A specific formulation of the \acrshort{hsmm} that explicitly defines the duration distribution is known as \acrfull{edhmm}. Transitions are allowed only at the end of each state, resulting in the following definition:
Popular choices for the duration distribution $h_\theta$ are the Poisson or the Negative Binomial distribution, for greater flexibility at the cost of an additional model parameter per state. The observation distribution $g_{\theta}$ can be set according to the specifics of the application at hand, which includes higher order data dependency such as $g_{\theta}(e_t \mid s_t, e_{t-k:t-1})$, for $t \ge 1$. An example for an excellent data dependency use case is provided in Section (ref), where \acrshort{ar}(1) models are used in each latent regime. The joint distribution of an \acrshort{edhmm} given the parameter can be stated as
The likelihood can be obtained by integrating out both $s_{1:T}$ and $d_{1:T}$. Due to the additional latent variables $d_{1:T}$, this is a much more computationally expensive operation than in the standard HMM case. In order to gain more insight on this, we can shrink the graphical model structure of a \acrshort{hmm} and \acrshort{hsmm} to a single time step. In order to compute the likelihood of observation $e_{t+1}$ given the current state $s_t$, a sum over all possible state transitions has to be taken, as can be seen in equation (ref). This resembles a standard mixture model computation, with the transition matrix of the \acrshort{hmm} replacing the mixture component weights.
For the \acrshort{hsmm}, however, the duration variable means an additional sum over a random variable that has at worst an infinite number of terms in the case of duration distributions with countably infinite support, shown in equation (ref).
Note that $\sum_{z_{t+1}} = \sum_{d_{t+1}}\sum_{s_{t+1}}$, which sums up all possible durations over all states, has at worst an infinite number of terms if the duration distributions have countably infinite support, and at best a large number of terms for long sequences, see Wood12. The standard approach to tackle this problem is to set up a minimum and maximum duration $d_{min}$ and $d_{max}$, where the computational complexity of the forward-backward algorithm reduces to $O(T(K(d_{max} - d_{min})^2) $, compared to the original $O(TK^2)$ in the \acrshort{hmm}, see Murphy02. Choosing an appropriate maximum duration varies depending on the underlying data. If the truncation is too small, then inference will typically fail, if is too large then calculations might become infeasible. Hence, $(d_{max} - d_{min})$ may increase the computational complexity to burdensome levels, which requires the modeler to set $d_{max}$ too small. Other approaches include Johnson13, Johnson14, who decrease computational complexity by censoring the initial or end time. To give a numerical example, in order to appropriately estimate a \acrshort{hsmm}, we assume that at least a single state transition has to occur, hence $d_{max} < T$, but $(d_{max} - d_{min}) >> K$. The computational costs for a $5$-state \acrshort{hmm} and $1000$ data points would be $\mathcal{O}(5^2 \times 1000)$, while for the \acrshort{hsmm}, assuming $d_{max} = 500$ and $d_{min} = 0$, $\mathcal{O}(5^2 \times 500^2 \times 1000)$ for a single likelihood call.
Alternatively, a particle filter can be used for the likelihood computation in $\mathcal{O}(NT)$ operations, even if the model has \acrshort{hsmm} dynamics. Exact inference is retained, subject to Monte Carlo error, using the \acrlong{pmcmc} algorithm Andrieu10. $N$ denotes the number of particles used in the filter, and we observed that it is sufficient to set $N = \frac{T}{2}$ for sequential Monte Carlo schemes on \acrshort{hsmm}s, see Section (ref) for more detail.
In a Bayesian framework, the typical goal is to infer the posterior distribution of the model parameter $\theta$ given the observed data $e_{1:T}$, $p(\theta \mid e_{1:T}) = \frac{ p_\theta(e_{1:T}) ~ p(\theta)}{ p( e_{1:T} )}$. This is a challenging task for \acrshort{ssm}s as it involves integrating $s_{1:T}$ over the likelihood $p_\theta(e_{1:t} ) = \int p_\theta(e_{1:T}, s_{1:T})~d s_{1:T}$, which is typically intractable or costly to evaluate. Hence, usually the full posterior distribution
is inferred. If $s_{1:T}$ is continuous, the target distribution $p( s_{1:t}, \theta \mid e_{1:t}) \propto p(e_{1:t} \mid s_{1:t}, \theta ) \times p( s_{1:t}\mid \theta ) \times p( \theta ) $ can theoretically be estimated via \acrshort{mcmc}. However, in this case, a state trajectory of ${p( s_{1:T} \mid {\theta})}$ has to be sampled while evaluating the target function, usually resulting in very poor outcomes of this strategy. Alternatively, classic Gibbs sampling strategies could be applied, which iterate the estimation process between sampling the latent states given the continuous model parameter and vice versa. However, for the latent state trajectory proposal step, a forward-backward algorithm would have to be employed again, as other choices such as one-at-a-time updates or overlapping blocks are known to cause slow mixing Kalogeropoulos10, Golightly09 of the Markov chain.
A more general attempt to jointly target the full joint posterior $p( s_{1:T}, \theta \mid e_{1:T}) $ can be shown as follows:
The last term in equation (ref) has been simplified by using the \acrfull{bmi} of Chib95. This framework allows to jointly sample $\theta$ and $s_{1:t}$, but the in general intractable likelihood function is still contained in step 3. While this term can be computed analytically for the discrete \acrshort{hsmm} via the so-called forward-backward algorithms, they are prohibitively expensive to run, as described in Section (ref). Going forward, we introduce the algorithmic machinery known as \acrlong{pmcmc} Andrieu10, which replaces the likelihood evaluation with an estimate $\hat{p}_\theta(e_{1:T})$ from a \acrlong{pf}.
\acrlong{pf}s are often used to solve filtering equations in the form of $\pi_t( x_{1:t} ) = \frac{ \tau_t(x_{1:t}) }{ z_t }$. The goal is to sequentially sample a sequence of random variables, $x_t, t \in (1,..., T)$ that come from a sequence of target probabilities $\pi_t( x_{1:t} )$ with the same computational complexity at each time step. If it impossible to directly sample from $\pi_t$, a similar proposal distribution $q_t$ can be used, s.t. $\pi_t(x_{1:t}) > 0 \Rightarrow q_t(x_{1:t}) > 0 $. The fraction of $\tau_t(x_{1:t})$ and $q_t(x_{1:t})$ is known as un-normalized weight function $ w_t( x_{1:t} ) = \frac{ \tau_t(x_{1:t}) }{ q_t(x_{1:t}) }$, s.t. the target distribution can be rewritten as $ \pi_t(x_{1:t}) = \frac{w_t(x_{1:t}) q_t(x_{1:t})}{z_t}$. This method is recognized as \acrfull{is}. Often, the variable of interest is the normalizing constant $z_t= \int \tau_t(x_{1:t}) d x_{1:t} = \int w_t(x_{1:t}) q_t(x_{1:t}) d x_{1:t} $, which can be approximated via the un-normalized weight functions
This technique does, unfortunately, rapidly degenerate as $t$ becomes larger. A technique to sequentially sample from such distributions is called \acrfull{sis}, which keeps the computational costs fixed given additional time steps by decomposing the joint distribution as $\tau_t(x_{1:t}) = \tau_{t-1}(x_{1:t-1}) \tau_t( x_t \mid x_{1:t-1})$. Similarly, the importance distribution can be decomposed as $q_t( x_{1:t} ) = q_1(x_1) \prod_{n=2}^{t} q_n(x_n \mid x_{1:n-1})$. The associated un-normalized weights can then be computed recursively via $w_t( x_{1:t} ) = w_1(x_1) \prod_{k=1}^{t} \alpha_k(x_{1:k})$, where the incremental importance weights $\alpha_t(x_{1:t})$ are given as $ \alpha_t(x_{1:t}) = \frac{ \tau_t(x_{t} \mid x_{1:t-1}) }{ q_t(x_t \mid x_{1:t-1}) }$.
In most state space models, the memory for $q_t(x_t \mid x_{1:t-1})$ is limited, so the target distribution can be evaluated and sampled from via the methods above at fixed computational costs. The only freedom in this framework is choosing an appropriate $q_t$. Unfortunately, it can be shown that the variance of the corresponding weights grows with $t$, and we refer to this problem as weight degeneracy. Due to weights degeneracy, the variance for the estimator of the normalizing constant $\hat{z}$ is also increasing with $t$. To alleviate this obstacle, a resampling step for the particle trajectories $x_{1:t}$ that normalizes the corresponding weights can be applied. Algorithms that do so at each iteration are known as \acrfull{sir}. However, this creates a new challenge known as sample path degeneracy, which refers to the problem that continuously resampling particle paths ultimately ends with very few unique trajectories. Balancing weight and sample path degeneracy is an ongoing research topic, and the most common method is to resample trajectories only at specific iterations, for example, if the \acrfull{ess}
of the particles is less than an a priori set threshold. For a discussion on several different resampling techniques, see Douc05. Adaptive resampling mitigates the exploding variance of the particle weights and keeps sample path degeneracy in check. Algorithms that apply this machinery are commonly referred as \acrlong{pf}s. They have a fixed computational complexity that is both linear in time $T$ and in number of particles $N$, $\mathcal{O}(NT)$, and return an estimate of the normalizing constant $\hat{z}_t$ and a particle path of $\hat{\pi}(x_{1:t})$. In the \acrshort{ssm} case, the joint distribution and the normalizing constant are of the form $\tau_t(x_{1:t}) = p_\theta(s_{1:t}, e_{1:t})$ and $z_t = p_\theta(e_{1:t})$. An approximation for the likelihood can be computed via the weights $\hat{z_t} = \frac{1}{N} \sum_{n=1}^{N} w_t( s_{1:t}^n, e_{1:t} )$ , which can be decomposed in the following recursive form:
The incremental weight $\alpha$ is defined as $\alpha_t(s_{1:t}, e_{1:t}) = \frac{ p_\theta( e_{t} \mid s_{1:t}, e_{1:t-1} ) ~p_\theta( s_{t} \mid s_{1:t-1}, e_{1:t-1} ) }{ q(s_t \mid s_{1:t-1}, e_{1:t}) }$, and the full likelihood estimate can be expressed as $\hat{z_T} = \frac{1}{N} \sum_{n=1}^{N} \prod_{k=1}^{T} \alpha_k(s^n_{1:k}, e_{1:k})$, which is usually the preferred method in particle filter software implementations as this avoids memory allocations. This permits the estimation of the incremental likelihood $p_{\theta}(e_t \mid e_{1:t-1}) \approx \frac{1}{N} \sum_{n=1}^{N} \alpha_t(s^n_{1:t}, e_{1:t})$ as well, which becomes relevant in the model selection Section (ref). $p_\theta( e_{t} \mid s_{1:t}, e_{1:t-1} )$ and $p_\theta( s_{t} \mid s_{1:t-1}, e_{1:t-1} )$ are model distributions, and have usually limited memory. The only free distribution to choose is $q(s_t \mid s_{1:t-1}, e_{1:t})$, which should ideally look like $p_\theta(s_t \mid s_{1:t-1}, e_{1:t})$. A common and simple choice is known as \acrfull{bpf}, which takes $q(s_t \mid s_{1:t-1}, e_{1:t}) = p_\theta( s_{t} \mid s_{1:t-1}, e_{1:t-1} )$, reducing the incremental weight to $\alpha_t(s_{1:t}, e_{1:t}) = p_\theta( e_{t} \mid s_{1:t}, e_{1:t-1} )$. Another popular approach is the so called \acrfull{apf} Pitt99, which assumes $\alpha_t(s_{1:t}, e_{1:t})$ to be independent of $s_t$ (in the Bootstrap \acrlong{pf}, this does not hold!). See Kantas15, Doucet11 for a more in-depth review. A pseudo algorithm implementation for a standard \acrshort{pf} can be found in Algorithm (ref), where the auxiliary variable $a^i_{t}$ refers to the ancestor path of a particular particle $s^i$ at time t. Hence a particle trajectory can be recursively defined as $s^i_{1:t} = (s^{a^i_{t}}_{1:t-1}, s^i_t)$. Resampling the whole particle trajectory is equivalent to sampling a new ancestor path. Note that it is usually much faster to sample ancestors one at a time and then recursively recover the resampled particle path than to resample the whole particle trajectory at each iteration. A variant of this algorithm is known as \acrfull{cpf}, where a single particle path, $s'_{1:T}$, is chosen a priori as reference trajectory. This implementation tracks a slightly different target distribution, $\hat{p}(s_{1:T} \mid s'_{1:T}, e_{1:T})$, and is used in Section (ref) and (ref). A pseudo algorithm for this \acrfull{cpf} with ancestor sampling can be found in Algorithm (ref), and a more thorough review can be read in lindsten14, lindsten15.
Once a particle filter has been run, observed and latent data can be forecasted by first sampling a new state $s_{T+1} \sim p_{\theta}(\cdot \mid s_{1:T}, e_{1:T}) $ and then a new data point given this state $e_{T+1} \sim p_{\theta}(\cdot \mid s_{1:T+1}, e_{1:T})$. $s_{T+1}$ can be sampled by forward propagating algorithm (ref) or (ref) from $T$ to $T+1$, thereby reusing particles from $1$ to $T$. This can be repeated for multiple time steps as well, resulting in a very fast procedure to sample from predictive distributions.
The most common Bayesian inference technique for model parameter $\theta$ is known as \acrlong{mcmc}. Basic familiarly with this concept is assumed, and we refer to Rosenthal14 for a more detailed review about standard \acrshort{mcmc} techniques. A pseudo algorithm for a basic Metropolis step can be found in Algorithm (ref). The major difficulty in algorithm (ref) is finding a good proposal distribution $f$, which often results in slow mixing. A \acrshort{mcmc} kernel that automatically tunes its proposal distribution at the cost of additional tuning hyper-parameter is known as \acrfull{hmc}, see neal12 for an introduction and betancourt18 for a review on this topic. A pseudo algorithm can be seen in Algorithm (ref). The additional tuning hyper-parameter can be configured on the fly in the famous extension \acrfull{nuts}, which was proposed by Gelman14. Notable improvements for this algorithm are suggested in betancourt16.
Note that both the \acrshort{hmc} and \acrshort{nuts} kernel require the target density to be fully differentiable with respect to the model parameter $\theta$. In the \acrshort{ssm} case, targetting the marginal posterior distribution $p(\theta \mid e_{1:T}) \propto p_\theta(e_{1:T}) ~p(\theta)$ is difficult or impossible via \acrshort{mcmc}, as the latent variables in $p_\theta(e_{1:T}) = \int_{s_{1:T} } p_\theta(e_{1:T}, s_{1:T})$ have to be integrated out in the proposal ratio. However, $p_\theta(e_{1:T}, s_{1:T})$ is usually computable pointwise, so $p(s_{1:T}, \theta \mid e_{1:T})$ can be targeted. In the \acrfull{pmh} case, formally introduced in Andrieu10 and shown in pseudo Algorithm (ref), a particle filter is used to obtain approximations for $p_{\theta}(e_{1:T} )$ and $p_{\theta}(s_{1:T} \mid e_{1:T})$ as substitutes for the analytical solutions to target $p(s_{1:T}, \theta \mid e_{1:T})$ jointly. In this setting, andrieu2009 have shown the puzzling result that one can do so and still target the exact posterior distribution of interest. A major difficulty for this method is finding a good \acrshort{mcmc} kernel $ K_{mcmc}(e_{1:T}, \theta )$, because the $\theta$ proposal will be accepted based on the particle filter likelihood estimate, so tuning might be very noisy. Moreover, gradient based \acrshort{mcmc} sampler do not work in this case as $p_\theta(e_{1:T})$ typically cannot be evaluated pointwise. A common critique on \acrshort{pmh} is thus that this algorithm is ill-suited for a higher dimensional model parameter $\theta$. A \acrshort{pmcmc} variant that can mitigate this is known as \acrfull{pgas}, see lindsten14, lindsten15. A pseudo algorithm is shown in Algorithm (ref). To account for sampling from an approximation via a particle filter and to preserve the invariance principle, a slightly adjusted $\hat{p}_{\theta^{\star}}(s^{\star}_{1:T} \mid s_{1:T}, e_{1:T})$ distribution is used to sample from the state trajectory. This method does not jointly estimate the state sequence and model parameter, but the state sequence is fixed when the new model parameter $\theta^{\star} \sim p_\theta(\theta^{\star} \mid s_{1:T}, e_{1:T})$ are sampled. In this step, more advanced \acrshort{hmc} style \acrshort{mcmc} kernels can be used as $p_{\theta^{\star}}(e_{1:T}\mid s_{1:T})$, which is easy to evaluate, replaces $p_{\theta^{\star}}(e_{1:T} )$ in the acceptance ratio. Hence, more advanced \acrshort{mcmc} kernels, such as the \acrshort{nuts} sampler, can be used to estimate model parameter $\theta$ in the \acrshort{pgas} setting.
In a times series setting, forecasting is of major relevance. The standard way for prediction in a Bayesian setting is simple: obtain the posterior predictive distribution by integrating out the model parameter $\theta$ and, in the \acrshort{ssm} case, the state trajectories ${s_{1:T}}$,
Once a sample for $p(s_{1:T}, \theta \mid e_{1:T})$ is obtained, the predictive distributions for $ s_{T+1} \mid s_{1:T}, e_{1:T}, \theta $ and $ e_{T+1} \mid s_{T+1}, s_{1:T}, e_{1:T}, \theta $ are trivial to sample from. A major drawback of the \acrshort{pmcmc} machinery is that, even though this algorithm primarily works for models suited to times series settings, it only works as batch estimation. Once additional data is observed, the algorithm needs to be run again to target $p(s_{1:T+1}, \theta \mid e_{1:T+1})$. A method that uses \acrshort{pmcmc} in a sequential setting is known as \acrlong{smc2}, see chopin2012. \acrshort{smc} based algorithm often expand in the density region, known as density tempering or annealing, or the data dimension, known as data tempering or annealing, see gunawan21. \acrshort{smc2} is a data tempering algorithm and moves a collection of particles that consist of the model parameters and latent states by incrementally adding data to the estimation process. At the beginning, particles are drawn from the prior. Subsequent particles are explored iteratively by using multiple \acrlong{pf} and particle \acrshort{mcmc} sampler. At each iteration, N \acrlong{pf} are used to obtain the incremental likelihood estimates $\hat{p}_{\theta^n}(e_{1:t} \mid e_{1:t-1}) = \hat{Z}_t = \frac{1}{N} \sum_{n=1}^{N} \alpha_{\theta^n}(s^n_{1:t}, e_{1:t})$ and state trajectories $s^{n}_{1:t} \sim \hat{p}_{\theta^{n}}(s^{n}_{1:t} \mid e_{1:t})$ for $n = 1,\dots, N$. If the estimates $\hat{p}_{\theta^n}(e_{1:t} \mid e_{1:t-1})$ are getting too noisy, the particles are jittered via \acrshort{pmcmc}. Note that either the Particle Gibbs or the Particle Metropolis Hastings variant can be chosen for the \acrshort{pmcmc} kernel. The density at the final iteration is the posterior distribution of interest. A pseudo algorithm can be seen in Algorithm (ref).
A powerful feature of the \acrshort{smc2} algorithm is that at each iteration, an unbiased estimate of the incremental marginal likelihood $\hat{p}(e_t \mid e_{1:t-1})$ is obtained at practically no extra costs, see equation (ref). From here on, it is straight forward to obtain an estimate for the marginal likelihood, $\hat{p}(e_{1:T} ) = \prod_{t=1}^T \hat{p}(e_t \mid e_{1:t-1})$, for model comparison or \acrfull{bf} calculations. Moreover, this machinery is highly parallelizable, a feature that is typically difficult to include in standard \acrshort{mcmc} techniques, and most \acrshort{smc2} iterations can be performed online as the particle filter can be propagated forward if no resampling step has been taken at the previous iteration. Additionally, during the propagation step, samples for $s_{t+1}$ and $e_{t+1}$ for posterior predictive distribution analysis can be obtained at each time index.
In our setup, we used a bootstrap particle filter with transition distribution equal to the model dynamics as defined in Section (ref). The particle resampling method was chosen to be systematic, and the resampling threshold for the ESS calculation was set to $75\%$. The only free tuning parameter in this case is the number of particles $N$. The higher $N$, the lower the variance for the log likelihood estimate, but the higher the computational costs per \acrshort{pmcmc} iteration. We note that ultimately, the \acrlong{pmcmc} samples are drawn from the correct target distribution, only subject to \acrlong{mc} error, independent of $N$. However, the mixing of the \acrshort{mcmc} chains might be slower if less particles are used. In this case, a larger $N$ might lead to better results at a fixed computational time horizon. Often, people set $N$ equal to the number of data points received, but in practice, significantly less particles may be used. As a sanity check, we performed two experiments to provide insight. First, we computed a likelihood estimate for a range of parameter values for $ \theta = \{ \mu, \sigma, r, \phi \}$ of a 2-state \acrshort{hsmm} for 1000 data points. $\mu$ and $\sigma$ represent the parameter for a Normal distribution given a latent state, $r$ and $\phi$ are the parameter of a Negative Binomial duration distribution. The transition distribution has no unknown parameter in the 2-state \acrshort{hsmm} case, as the diagonal elements of the transition matrix are separately modeled by the duration. The true parameter can be seen as vertical grey lines on each subplot in Figure (ref), which shows the log likelihood estimate for a range of each parameter that was chosen based on the 95% \acrfull{ci} of a \acrshort{pmcmc} run in Section (ref), keeping all other parameter fixed. It can be seen that the variance is reasonably similar for $N = 500$, $1000$ or $2000$. Second, we directly examined the particle filter performance during a \acrshort{pmcmc} run for the model defined above. The \acrfull{pacf} plot of the log likelihood \acrshort{pf} estimates for a varying number of particles can be seen in Figure (ref). As there is little difference in the likelihood estimate variance between \acrshort{pf}s with $500$ and $2000$ particles and the \acrshort{pacf} look reasonably similar for the same particle range, we will use $N = \frac{T}{2}$ for our analysis going forward.
Once a particle filter is designed, only a suitable MCMC kernel has to be chosen. As discussed in Section (ref), more advanced gradient based \acrshort{mcmc} sampler do not work in the \acrlong{pmh} case, as $p_\theta(e_{1:T})$ typically cannot be evaluated pointwise. Thus, we will use \acrlong{pgas} and target $\theta^{\star} \sim p_\theta(\theta^{\star} \mid s_{1:T}, e_{1:T})$ in the MCMC step. As gradients for $p_{\theta^{\star}}(e_{1:T}\mid s_{1:T})$ can easily be calculated in this case, we choose the NUTS MCMC variant Gelman14, betancourt16 as MCMC kernel, as other kernels often take significantly more proposal steps to move toward the typical set. This is especially relevant for \acrshort{pmcmc} on our model, as we will run a comparatively expensive particle filter after each \acrshort{mcmc} proposal.
As described in pseudo-algorithm (ref), \acrshort{smc2} has a particle filter and a \acrshort{pmcmc} algorithm assigned for each particle. These particles may be propagated in parallel, so the number of \acrshort{smc2} are typically chosen to be a multiple of the available computer cores. Initial model parameter drawn from the prior distributions, and the associated \acrshort{pf} and \acrshort{pmcmc} algorithm will be initiated based on the starting data of length $t_0<<T$. Another tuning parameter is the number of jittering steps in the resampling step. As a guideline, we will continue jittering until the maximum parameter correlation is below $75\%$.
Once parameter are estimated, how should the performance of a model be evaluated? A powerful method for comparison is to validate models based on their marginal likelihood $p(e_{1:T} ) = \int p(e_{1:T}, \theta)~d \theta = \prod_{t=1}^T p(e_t \mid e_{1:t-1})$. This distribution is typically intractable or very costly to evaluate, but during an \acrshort{smc2} run, an estimate for $p(e_t \mid e_{1:t-1})$ can be obtained at practically no extra cost at each time step. In the particle filter propagation step in algorithm (ref), the incremental weight given current model parameter $\theta_m$ is computed via
where $\alpha_{t, \theta_n}$ is defined as in Section (ref) and $M$ is the number of particles that are used for the particle filter associated to continuous parameter vector $\theta_n$. After all particles have been propagated forward, a Monte Carlo estimator for $p(e_t \mid e_{1:t-1})$ is obtained by weighting these likelihood increments from equation (ref) with the corresponding normalized particle weight $w_n \propto w_{n-1} \hat{p}_{\theta_n}(e_t \mid e_{1:t-1})$ associated to parameter $\theta_n$,
Moving forward, we refer to $p(e_{t+1} \mid e_{1:t})$ as (one step ahead) \acrfull{pl} at $t+1$, $PL_{t+1}$, see Kastner16. After the final iteration, the marginal likelihood estimate can then be computed as $\hat{p}(e_{1:t}) = \prod_{i=1}^t \hat{PL}_t$. Based on the cumulative sums of log \acrshort{pl}s, model choice can be performed via the so called \acrfull{clpbf}. To compare model $A$ and $B$, for $u > 0$, the \acrshort{clpbf} is defined as
A positive \acrshort{clpbf} indicates evidence in favor for model A and, if $t=0$ and $u=T$, this factor is known as log \acrlong{bf}. Note that in Section (ref), we mentioned that we usually use $t_0 > 1$ data points to initialize the jitter kernels and then record all future \acrshort{pl} increments going forward, hence the resulting estimate will be the slightly different $\hat{p}(e_{t_{0}+1:T} \mid e_{1:t_0}) = \prod_{t=t_0+1}^T \hat{PL}_{t}$. The $t_0$ data points can be seen as training data for the jitter kernels to be initialized in a reasonable parameter region.
This section consists of simulation experiments conducted to study the performance of the \acrshort{pf}, the \acrshort{pmcmc} and the \acrshort{smc2} algorithm on \acrlong{hsmm}s. We first generate 1000 data points from a \acrshort{hsmm} with Negative Binomial state duration distributions. A plot with sampled observed and latent data is shown in figure (ref), which also depicts a \acrshort{pf} estimate of the latent parameter for known continuous model parameter. The advantages of explicitly modelling duration are visible on the third subplot, in which durations in each state vary drastically. The last sub-section shows how the number of hidden states in the data can be estimated using \acrshort{smc2}.
The model consists of parameter: $\theta = \{\mu, \sigma, p, r, \phi \}$, where the data $e_t \sim N(\mu_{s_t}, \sigma_{s_t})$ has a normal distribution given the latent state. The latent states have the same dynamics as explained in Section (ref), latent state and duration \\ $s_t \sim
$, $d_t \sim
$.
The $\mu$ parameter have truncated Normal priors with equal variance and different means, $\mu_1 \sim Normal_{(-100, 0)}(\mu = -2, \sigma = 10^5)$, $\mu_2 \sim Normal_{(0, 100)}(\mu = 2, \sigma = 10^5)$. We assigned a truncated Normal prior for the variances $\sigma \sim Normal_{(0, 10)}(\mu = 2, \sigma = 10^5)$ and for the Negative Binomial parameter, $r \sim Normal_{(0, 100]}(\mu = 10, \sigma = 10^5)$. The second duration distribution parameter, $\phi$, has equal mass from 0 to 1, $\phi \sim Beta(\alpha = 1, \beta = 1)$. Similarly, we assigned a Dirichlet prior for the transition probabilities, $p$, that favors equal weights, $p \sim Dirichlet( \alpha_1 = \alpha_2 = ... = \alpha_k = k),$ where $k =$ number of latent states.
If parameter $\theta$ are unknown, both the \acrshort{pmcmc} and \acrshort{smc2} machinery can be used. As the latter builds on the former, we first show estimation results for the \acrlong{pmcmc} run. The traceplots of four chains for the continuous model parameter can be seen in figure (ref), and for the latent state sequence in figure (ref). After parameter are initialized from the prior distributions, they rapidly converge to the values used to generate the sample data. Common MCMC output statistics are summarized in Table (ref). \acrshort{smc2} results can be seen in figure (ref) for the model parameter and in figure (ref) for the latent variables. Just as in the \acrshort{pmcmc} case, samples converge fast toward the typical set. Common MCMC output statistics of the final \acrshort{smc2} iteration for the continuous model parameter are displayed in Table (ref).
Before estimating model parameter for a given data set, the researcher has to choose the number of hidden states in a discrete \acrshort{ssm}, which is often challenging a priori. To tackle this problem, we choose methods from the overfitting mixtures literature, see rousseau11, frue06, frue18. In particular, rousseau11 show that the posterior distribution of a mixture model is much more stable if the prior weights of the mixture components are concentrated on the boundary regions of the parameter space. We borrow this concept for \acrshort{ssm}s and show that by assigning appropriate prior weights on the transition matrices of a \acrshort{hsmm}, we can choose more states than necessarily describe the data and still have an interpretable structure, as superfluous latent states will never be visited during the estimation process. As an experiment, we fit a 2-, 3-, and 5-state HSMM to data that was generated by a 3-state HSMM. The parameter used to generate the data can be seen in Table (ref), and the corresponding data in plot (ref). In Section (ref), we mentioned that the \acrlong{pl} can be estimated as a by-product of the inference procedure in \acrshort{smc2}. Hence, we tracked the cumulative log \acrshort{pl} for each model at each iteration, and plotted the results in figure (ref). It can be seen that the 3- and 5-state \acrshort{hsmm}s \acrlong{pl} perfectly align, and we conclude that it is possible to determine the number of states via \acrshort{smc2}, and even check for differences during the data propagation. All estimates for the model parameter can be seen in Table (ref). In this table, state $2$ and $4$ are almost never visited and act as superfluous regimes. The remaining state estimates contain the correct parameter in the 95% \acrlong{ci}s.
This section consists of two interesting applications for \acrlong{ssm}s. The first contains a case study on model selection and prediction for a financial data set, the other displays how \acrshort{ssm}s can be used for clustering purposes.
In this Section, we are modeling VIX index data directly via \acrshort{hsmm}s and other popular benchmark models. The VIX index is derived from options with near-term expiration dates on the US major equity S&P 500 index, and is a popular indicator for future short term volatility expectations. Volatility is often modeled indirectly based on (log) stock prices via \acrfull{sv} models and its autoregressive nature is well recognized, see e.g. hull:1987, kim:1998, Kastner16. We incorporate this behaviour in our model by assigning autoregressive weights to the location parameters in each latent regime. The corresponding model can be viewed as a separate \acrshort{ar}(1) model in each state, where a latent variable following a semi-Markov chain governs the regime changes and durations. Multiple duration distribution choices, such as the negative Binomial and Poisson distribution, are tested and benchmarked against \acrshort{hmm} and AR(1) models.
For real data collection, we first obtain 1000 daily end-of-day data points of this instrument, which can be seen in figure (ref), as of January 1st 2022 from the Thomson Reuters database. Notably, the Covid-19 epidemic is included in the data, which causes the evidence to be much more volatile from the start of 2020 and to transition to a different regime from there onward. Common ideas to address this issue are to assign change points across the times series before estimating the model parameter. Such methods have an easily interpretable structure, but lack the information from a stochastic process governing the model dynamics, which enhances inference capabilities for the data. For example, due to the discrete state space formulation of our proposed model, parameter interpretation is straightforward, while the \acrshort{hsmm} dynamics allow to incorporate significantly more decision-making tools, such as regime change and duration forecasting.
The AR(1) HSMM with Negative Binomial duration consists of parameter: $\theta = \{\mu, \sigma, w, p, r, \phi \}$, where the data $e_t \sim N(w_{s_t} \times e_{t-1} + \mu_{s_t}, \sigma_{s_t})$ has a normal distribution given the latent state and the previous data point.
The latent states have the same dynamics as explained in Section (ref), latent state and duration \\ $s_t \sim
$, $d_t \sim
$.
The $\mu$ parameter have truncated Normal priors with equal variance and different means, $\mu_1 \sim Normal_{(0, 10)}(\mu = 3, \sigma = 10^5)$, $\mu_2 \sim Normal_{(0, 10)}(\mu = 1, \sigma = 10^5)$. We assigned a truncated Normal prior for the variances $\sigma_{1,2} \sim Normal_{(0, 10)}(\mu = 0.2, \sigma = 10^5)$ and for the Negative Binomial parameter, $r_1 \sim Normal_{(0, 100]}(\mu = 5, \sigma = 10^5)$, $r_2 \sim Normal_{(0, 100]}(\mu = 2.5, \sigma = 10^5)$. The second duration distribution parameter $\phi$ has equal mass from 0 to 1, $\phi_{1,2} \sim Beta(\alpha = 1, \beta = 1)$. Similarly, we assigned a Dirichlet prior for the transition probabilities $p$ that favors equal weights, $p \sim Dirichlet( \alpha_1 = \alpha_2 = ... = \alpha_k = k),$ where $k =$ number of latent states. The autoregressive parameter, $w$, is bounded between $-1$ and $1$ by assigning a truncated prior, $w_{1,2} \sim Normal_{(-1, 1)}(\mu = 0, \sigma = 10^5)$. The intuition behind the prior assignment is that there is a single state with a large constant that mimics sudden jumps in the index, while the other state has a much smaller constant, but possibly a higher autoregressive weight $w$. This assumption is confirmed in the actual estimates, visible in figure (ref). The AR(1) HSMM with Poisson duration only differentiates with respect to the duration parameter, $\theta = \{\mu, \sigma, w, p, \lambda \}$. In this case, the latent duration has a Poisson distribution, \\ $d_t \sim
$, where $\lambda \sim Normal_{(0, 100]}(\mu = 20, \sigma = 10^5)$. The AR(1) HMM with parameter $\theta = \{\mu, \sigma, w, p\}$, has the same data dynamics and latent state dynamics $s_t \sim Categorical(p_{s_{t-1}})$. All prior configurations are assigned from the previous models. The AR(1) Model with parameter $\theta = \{\mu, \sigma, w \}$ does not have latent variables, and data dynamics $e_t \sim Normal(w \times e_{t-1} + \mu, \sigma)$. The priors for $\mu$ are set to $\mu \sim Normal_{(0, 10)}(\mu = 2.0, \sigma = 10^5)$, where the location parameter has been set as the mean of the location parameter for the 2 states in the state space models defined above. $\sigma$ and $w$ have the same prior as the other models.
A detailed methodology for model comparison can be found in section (ref). The \acrshort{smc2} machinery with the tuning configurations discussed in Section (ref) is used to estimate model parameter. Results for our proposed model can be seen in Figure (ref) and Table (ref) for model parameter. The traceplots show all model parameter estimates at each point in time with the corresponding 95% credible interval. There is a clear distinction between a lower and higher volatility state, in which the latter has volatility levels of roughly four times the normal state. Interestingly, there seems to be a regime shift at the time COVID-19 makes significant headlines in the global markets, which is captured very fast during the estimation process, and would be impossible for standard batch estimation methods. At the start of March 2020, the volatility parameter in the high-volatility state changes from $0.15$ to $0.20$ within days, and continues to be around this level until the end of the times series, making the regime change easily recognizable. Figure (ref) shows the filtered state trajectory of the latent variables at each time step. A re-scaled posterior mean of the latent variable at each time index is shown against the real data in the bottom sub-plot, which displays the clear distinction between a volatile and a stable state in the times series. There is more variation at the initial stages before more data is added to the algorithm. Figure (ref) depicts model predictions against the realized data at each time step. The \acrlong{clpbf} as discussed in Section (ref) is shown in figure (ref). The top plot in this figure shows the \acrshort{clpbf} of the AR(1) \acrshort{hsmm} with negative Binomial duration distribution against the other models, which is positive and in favor of the proposed model consistently over time against all other models. Diagnostics at the final time period are summarized in Table (ref).
Once parameter have been estimated, the latent state trajectory estimates can be used to cluster the log VIX data into different regimes. A re-scaled state posterior mean of the latent trajectory at the last \acrshort{smc2} against the actual data can be seen at the top plot in Figure (ref). State 1 corresponds to a short duration state with jumps and drastic changes in levels, while state 2 corresponds to a regime with more normal volatility levels that is observed for the majority of times. Based on the posterior mean, we clustered the changes in the index for each state, which can be seen in the bottom plot of Figure (ref). Here, state 1 is clearly associated with drastic movements in either direction, while state 2 represents a more normal market environment.
In this paper, we discussed sequential parameter estimation techniques for \acrlong{ssm}s with a focus on \acrlong{hsmm}s. We compared the forecasting accuracy of various models on financial data and concluded that the \acrshort{hsmm} has a superior predictive performance against other popular discrete \acrshort{ssm}s. Moreover, we demonstrated how by-products emerging from a \acrshort{smc2} estimation run can be used to determine the number of latent regimes governing such models. While the additional duration variable in the \acrshort{hsmm} typically causes parameter inference to be more challenging due to the increased computational costs of the likelihood function, it adds significantly more flexibility in modelling the latent process. Our proposed inference technique has the same computational costs for both the basic \acrshort{hmm} and the \acrshort{hsmm} and is particularly suitable for sequential data. As for future research topics, more optimized techniques to adaptively select the number of particles in a \acrshort{pf} may lead to faster runs and improved mixing for both the \acrshort{pmcmc} and \acrshort{smc2} algorithm. Furthermore, as tuning the individual \acrshort{mcmc} and \acrshort{pf} kernels was handled independently during the \acrshort{smc2} runs, adding information from all chains may drastically increase the tuning process for the individual jitter kernels.
The data and code used to run the algorithms in this paper can be be seen in \url{https://github.com/paschermayr/Publish_SequentialHSMM}. For more detailed information about the implementations for running all algorithms and computing all tables can be found in \url{https://github.com/paschermayr/Baytes.jl} and its sub-libraries. The corresponding plots are defined in \url{https://github.com/paschermayr/BaytesInference.jl}.
\printglossary[type=\acronymtype]