Penalized estimation of flexible hidden Markov models for time series of counts

TL;DR

This paper introduces a nonparametric penalized estimation method for hidden Markov models tailored to count time series, enabling flexible, data-driven modeling of state-dependent distributions without pre-specifying distributional forms.

Contribution

It proposes a novel penalized likelihood approach that estimates state-dependent distributions nonparametrically, improving flexibility and avoiding distributional misspecification in HMMs for count data.

Findings

Method produces smooth, data-driven state-dependent distributions.

Simulation confirms effectiveness in avoiding overfitting.

Applications demonstrate practical utility in earthquake and animal movement data.

Abstract

Hidden Markov models are versatile tools for modeling sequential observations, where it is assumed that a hidden state process selects which of finitely many distributions generates any given observation. Specifically for time series of counts, the Poisson family often provides a natural choice for the state-dependent distributions, though more flexible distributions such as the negative binomial or distributions with a bounded range can also be used. However, in practice, choosing an adequate class of (parametric) distributions is often anything but straightforward, and an inadequate choice can have severe negative consequences on the model's predictive performance, on state classification, and generally on inference related to the system considered. To address this issue, we propose an effectively nonparametric approach to fitting hidden Markov models to time series of counts, where the state-dependent distributions are estimated in a completely data-driven way without the need to select a distributional family. To avoid overfitting, we add a roughness penalty based on higher-order differences between adjacent count probabilities to the likelihood, which is demonstrated to produce smooth probability mass functions of the state-dependent distributions. The feasibility of the suggested approach is assessed in a simulation experiment, and illustrated in two real-data applications, where we model the distribution of i) major earthquake counts and ii) acceleration counts of an oceanic whitetip shark (Carcharhinus longimanus) over time.

Tables

Table 1. Table 1: Average Kullback-Leibler divergences, mean absolute errors, and average state misclassification rates, as obtained under the different estimation approaches.

model specification	${\overline{\text{KLD}}}_1$	${\overline{\text{KLD}}}_2$	${\overline{\text{MAE}}}_{{\hat{\gamma }}_{12}}$	${\overline{\text{MAE}}}_{{\hat{\gamma }}_{21}}$	$\overline{\text{SMR}}$
(1a) true parametric	0.005	0.009	0.012	0.013	0.033
(1b) incorrect parametric	0.112	2.085	0.122	0.395	0.240
(2a) nonparam. w/o pen.	0.379	0.652	0.012	0.013	0.043
(2b) nonparam. w/ pen.	0.043	0.029	0.012	0.014	0.041

Equations

$${\gamma }_{ij}=\mathrm{Pr}\left(S_t=j|S_{t-1}=i\right),$$

$${\delta }_i=\mathrm{Pr}\left(S_1=i\right),$$

$${\pi }_{i,k}=\mathrm{Pr}(Y_t=k|S_t=i),$$

$$\mathcal{L}(\mathbf{\theta })=\mathrm{Pr}(y_1,\dots ,y_T|\mathbf{\theta })=\mathbf{\delta }\mathbf{P}(y_1)\prod _{t=2}^T\mathbf{\Gamma }\mathbf{P}(y_t)1,$$

$$\begin{array}{rl}{\mathbf{\alpha }}_1& =\mathbf{\delta }\mathbf{P}\left(y_1\right);\\ {\mathbf{\alpha }}_t& ={\mathbf{\alpha }}_{t-1}\mathbf{\Gamma }\mathbf{P}\left(y_t\right),\end{array}$$

$$\log \left({\mathcal{L}}_{\text{pen.}}(\mathbf{\theta })\right)=\log \left(\mathcal{L}(\mathbf{\theta })\right)-\sum _{i=1}^N{\lambda }_i\sum _{k=m}^K{\left({\Delta }^m{\pi }_{i,k}\right)}^2,$$

$${\gamma }_{ij}=\frac{\exp ({\gamma }_{ij}^{\ast })}{{\sum }_{l=1}^N\exp ({\gamma }_{il}^{\ast })},\text{leavevmode}{\delta }_i=\frac{\exp ({\delta }_i^{\ast })}{{\sum }_{l=1}^N\exp ({\delta }_l^{\ast })},\text{leavevmode}\text{and}\text{leavevmode}{\pi }_{i,k}=\frac{\exp ({\pi }_{i,k}^{\ast })}{{\sum }_{l=0}^K\exp ({\pi }_{i,l}^{\ast })},$$

$$\mathbf{\Gamma }=\left(\begin{array}{cc}0.95& 0.05\\ 0.05& 0.95\end{array}\right).$$

$${\overline{\text{KLD}}}_i$$

$${\overline{\text{MAE}}}_{{\hat{\gamma }}_{ij}}$$

$$\overline{\text{SMR}}=\frac{1}{200}\sum _{r=1}^{200}\frac{1}{500}\sum _{t=1}^{500}\text{mathds}{1}_{{\hat{s}}_t^{(r)}\ne s_t^{(r)}},$$

$$y_t=\#\text{ of earthquakes worldwide with magnitude}\ge \text{7 in year }t.$$

$$\hat{\mathbf{\Gamma }}=\left(\begin{array}{cc}0.934& 0.066\\ 0.128& 0.872\end{array}\right),$$

$$y_t=\sum _{s=1}^{60}\text{mathds}{1}_{{\text{ODBA}}_{s,t}>0.05},$$

$$\hat{\mathbf{\Gamma }}=\left(\begin{array}{cc}0.960& 0.040\\ 0.093& 0.907\end{array}\right),$$

Full text

Keywords: count data; nonparametric statistics; penalized likelihood; smoothing parameter selection; state-space model; time series modeling.

1 Introduction

Over the past half century, hidden Markov models (HMMs) have become increasingly popular tools for modeling time series data where, at each point in time, a hidden state process selects among a finite set of possible distributions for the observations (Zucchini et al., 2016). For example, in economic applications, the states of the Markov chain underlying the observations (which determines the state process) are often good proxies for market regimes such as periods of economic growth or recessions, while in ecology, they can regularly be linked to an animal’s behavioral modes such as resting, foraging, or traveling. Other areas of application include volcanology (Bebbington, 2007), psychology (Visser et al., 2002), and medicine (Jackson & Sharples, 2002), to name but a few examples. Depending on the application at hand, potential aims that can be addressed using HMMs are manifold, including the prediction of future values of a time series, decoding of the hidden states underlying the observations, and inference on the drivers for example on the state-switching dynamics.

HMMs constitute a versatile framework for modeling diverse types of time series data and can easily be tailored for example to binary data (Schliehe-Diecks et al., 2012), positive real-valued data (Langrock, 2012), circular data (Bulla et al., 2012), categorical data (Maruotti & Rocci, 2012), compositional data (Langrock et al., 2013), or count data (Lagona et al., 2015). Here we specifically focus on the latter type of data, i.e. sequences of non-negative integers. MacDonald & Zucchini (2000) provide an introduction to the use of HMMs for discrete-valued time series, including count data, while Weiß (2018) gives a more comprehensive overview of the various classes of models for time series of counts as well as further types of discrete-valued time series. A typical example for a count time series is the number of corporate defaults observed on a monthly, quarterly, or yearly basis: In periods of economic growth, these may be generated by some distribution with relatively small mean, whereas during recessions, another distribution with relatively higher mean may be active. Although the market regime is not directly observable, it clearly affects the observed counts (c.f. Li & Cheng, 2015). For an application that is similar in spirit, see Hambuckers et al. (2018), where both the number and the amount of operational losses of a bank are modeled using an HMM-type approach. Beyond economics, HMMs have been applied to time series of counts in a variety of different fields, including, inter alia, medicine (multiple sclerosis leisure counts; Altman & Petkau, 2005), geology (volcanic eruption counts; Bebbington, 2007), epidemiology (oliomyelitis counts; Le Strat & Carrat, 1999), ecology (pilot whale vocalization counts; Popov et al., 2017), and bioinformatics (T-lymphocyte counts; Marino et al., 2018).

Specifically for time series of counts, a Poisson distribution often provides a natural choice for the state-dependent distributions within an HMM, where one rate parameter is estimated for each state of the underlying Markov chain. While more flexible distributions such as the negative binomial can also be used, or models for bounded counts like the binomial, choosing an adequate family of parametric distributions generally remains a difficult task in practice, with potentially severe negative consequences in case an inadequate choice is made. For the case of continuous-valued time series, Langrock et al. (2015) propose a nonparametric approach to estimating the state-dependent distributions within an HMM, where the state-dependent densities are constructed using penalized B-splines (P-splines; Eilers & Marx, 1996), which yields flexible yet smooth functional shapes of the resulting densities without the need to make any distributional assumptions. For time series of counts that are either naturally bounded or considered to be effectively bounded (defining an upper threshold for the support on which the distribution is modeled), one can, in principle, easily avoid a distributional assumption by directly estimating the values of the state-dependent probability mass functions (p.m.f.s) on the support considered; this approach is in fact implemented in the R package hmm.discnp (Turner, 2018). However, without penalization such an approach will often lead to overfitting, which again limits the usefulness of corresponding models in particular for prediction and classification. Following Scott et al. (1980) and Simonoff (1983), and similar in spirit to Langrock et al. (2015), we here suggest to address this issue by considering a penalized likelihood function, where we add a roughness penalty based on higher-order differences between adjacent count probabilities. We further demonstrate how the penalty term can be adjusted in presence of (e.g. zero-) inflated observations, where small differences between corresponding count probabilities and their respective neighbors are not necessarily desired. This conceptually simple approach is demonstrated to be immensely effective in producing reliable estimates of both simple and complex state-dependent p.m.f.s within an HMM, with a smoothing parameter adjusting the required flexibility in a completely data-driven way.

The paper is structured as follows: In Section 2, we briefly recall the model formulation and dependence structure of HMMs, and introduce some notation specifically for the case of time series of counts. Furthermore, we provide an efficient algorithm for evaluating the likelihood and discuss how the model parameters can be estimated in a penalized maximum likelihood framework. In Section 3, we assess our approach in a simulation experiment, where we compare the performance of the penalized nonparametric approach to its parametric counterpart. In Section 4, we present two real-data applications, where we model the distribution of i) major earthquake counts and ii) acceleration counts of an oceanic whitetip shark (Carcharhinus longimanus) over time.

2 Methodology

2.1 A nonparametric hidden Markov model for count data

Basic HMMs comprise two stochastic processes, only one of which is observed, namely the time series of interest, $\{Y_{t}\}_{t=1,\dots,T}$. The observed process is driven by another, hidden process, which we denote by $\{S_{t}\}_{t=1,\dots,T}$. The latter process is the so-called state process, which is usually modeled as an $N$-state Markov chain. Throughout this paper we consider a first-order Markov chain, i.e. we assume the state process to satisfy the Markov property, $\Pr(S_{t}|S_{1},\dots,S_{t-1})=\Pr(S_{t}|S_{t-1})$, $t=2,\dots,T$. This simplifying dependence assumption is exploited in the likelihood calculations provided in Section 2.2 and can be relaxed to higher-order Markov chains or semi-Markov chains if deemed necessary (Zucchini et al., 2016). Assuming the first-order Markov chain to be time-homogeneous, the state transition probabilities are summarized in the $N\times N$ transition probability matrix (t.p.m.) $\boldsymbol{\Gamma}$, with elements

[TABLE]

$i,j=1,\dots,N$. The initial state probabilities, i.e. the probabilities of the state process being in the different states at time $t=1$, are summarized in the $N$-dimensional row vector $\boldsymbol{\delta}$, with elements

[TABLE]

$i=1,\dots,N$. If the Markov chain is assumed to be stationary, which is reasonable in many applications, then the initial distribution is the stationary distribution, i.e. the solution to the equation system $\boldsymbol{\delta}\boldsymbol{\Gamma}=\boldsymbol{\delta}$ subject to $\sum_{i=1}^{N}\delta_{i}=1$ (Zucchini et al., 2016). Otherwise, the initial state probabilities are parameters which need to be estimated. The state process is completely specified by the initial state probabilities and the state transition probabilities.

The basic dependence structure is such that the observations are assumed to be conditionally independent of each other, given the states. The states then directly select which of $N$ possible distributions generates the observation at any time point. This dependence structure is illustrated in Figure 1. For time series of counts, the observed (state-dependent) process $\{Y_{t}\}_{t=1,\dots,T}$ is determined by its state-dependent p.m.f.s, $\Pr(Y_{t}=k|S_{t}=i)$, $k=0,\dots,K$ and $i=1,\dots,N$. It is common to use some distributional family such as the class of Poisson distributions to model these state-dependent p.m.f.s (see, e.g., Altman & Petkau, 2005; Bebbington, 2007). Here we do not make any such assumption, and instead assign a state-specific probability to each possible count on the bounded support $\{0,\dots,K\}$, i.e. we consider the model parameters

[TABLE]

$i=1,\dots,N$ and $k=0,\dots,K$. While count data can be unbounded, we here consider an upper threshold as to obtain a fixed number of parameters. The support should be bounded in a reasonable way; specifically, it should at least cover all observed counts. The state-dependent process is completely specified by the count probabilities ($NK$ free model parameters). Thus, with this model formulation, we consider a possibly large number of parameters rather than say only one (as e.g. in the case of a Poisson or binomial distribution) or two (as e.g. in the case of a negative binomial distribution).

Although the parameter space in this model formulation is still finite-dimensional, it will usually have a fairly high dimension, with the individual parameters $\pi_{i,k}$ not being of direct interest themselves. As a consequence, we follow Turner (2018) and call our approach “nonparametric”. In particular, this label emphasizes that the state-dependent p.m.f.s are not determined by a small number of parameters, as would be the case when a distributional family such as the class of Poisson distributions would be used. In addition, with the given model formulation we are not restricted to any particular functional shape of the state-dependent p.m.f.s, and instead have full flexibility to let the data speak for themselves, just like with other methods for which the label “nonparametric” is commonly used in the literature.

2.2 Likelihood evaluation

For a given parameter vector $\boldsymbol{\theta}$, in our case comprising the state transition probabilities $\gamma_{ij}$, $i,j=1,\ldots,N$ with $i\neq j$, and the count probabilities $\pi_{i,k}$, $i=1,\dots,N$ and $k=0,\dots,K-1$, and a time series of counts $y_{1},\dots,y_{T}$, the likelihood function of the HMM as formulated above is given by

[TABLE]

where $\mathbf{P}(k)=\text{diag}(\pi_{1,k},\dots,\pi_{N,k})$, and where $\mathbf{1}\in\mathbb{R}^{N}$ is a column vector of ones. The evaluation of (1) corresponds to the application of the so-called forward algorithm. Specifically, defining the forward probabilities $\alpha_{t}(i)=\Pr(y_{1},\dots,y_{t},S_{t}=i)$, which are summarized in the row vectors \boldsymbol{\alpha}_{t}=\bigl{(}\alpha_{t}(1),\dots,\alpha_{t}(N)\bigr{)}, the recursive scheme

[TABLE]

$t=2,\dots,T$, can be applied to arrive at $\boldsymbol{\alpha}_{T}$, from which we obtain $\mathcal{L}(\boldsymbol{\theta})=\Pr(y_{1},\dots,y_{T}|\boldsymbol{\theta})$ $=\sum_{i=1}^{N}\alpha_{T}(i)=\boldsymbol{\alpha}_{T}\mathbf{1}$ by the law of total probability (cf. Zucchini et al., 2016).

Using (2), evaluating the likelihood requires $\mathcal{O}(TN^{2})$ operations, which makes an estimation of the model parameters by numerically maximizing the likelihood (or log-likelihood in case of numerical underflow) practically feasible even for relatively long time series and a moderately high number of states. Alternatively, the expectation-maximization algorithm, which also arrives at a (local) maximum of the likelihood, can be used. For the model formulation considered here, the latter approach is implemented in the R package hmm.discnp (Turner, 2018).

2.3 Roughness penalization

The downside of the above nonparametric and hence very flexible approach to estimating the state-dependent p.m.f.s is its propensity to overfit any given data. Especially in cases where $T$ is small relative to the number of model parameters, $N(N-1)+NK$, the fitted state-dependent p.m.f.s will often be anything but smooth, and may even involve isolated spikes with implausible gaps in between (corresponding to a lack of data in regions where observations would in fact be expected to occur in the long run). For short time series, it can in fact easily happen that for specific values well within the plausible range of observations to occur in future, each state-dependent probability is estimated to be zero, namely if no such observations are present in the training data (see Section 4.1 for an example of this problem). The consequence of this would be that the model deems the corresponding values to be impossible to occur in future, which could be problematic for forecasting. To avoid such kind of overfitting, we add a roughness penalty to the logarithm of the likelihood given in (1), which leads to the penalized log-likelihood

[TABLE]

where $\lambda_{i},i=1,\dots,N$, denotes a smoothing parameter associated with the $i$-th state-dependent distribution, and where $\Delta^{m}\pi_{i,k}=\Delta^{m-1}(\Delta\pi_{i,k})$, $\Delta\pi_{i,k}=\pi_{i,k}-\pi_{i,k-1}$, denotes the $m$-th order differences between adjacent count probabilities.

The inclusion of the penalty term, together with the associated smoothing parameters, allows to control the variance of the otherwise unrestricted and hence highly variable estimation of the state-dependent p.m.f.s. Optimizing the penalized log-likelihood (3) then amounts to finding a good compromise between goodness of fit, as measured by the logarithm of the likelihood given in (1), and smoothness of the state-dependent p.m.f.s, as measured by the $m$-th order differences. The order of the differences, $m$, should be chosen pragmatically based on the data at hand. With $m=1$, a uniform distribution is obtained as $\lambda_{i}\rightarrow\infty$ (the penalty term vanishes if all count probabilities are equal). For $m=2$, we obtain a triangular distribution in the limit (the penalty term vanishes if all count probabilities lie on a straight line with arbitrary slope). Based on our experience, $m=3$ or an even higher order produces the most reliable estimates across a range of scenarios of varying complexity. In case of zero-inflation, which in practice often occurs when dealing with count data, it may make sense not to penalize differences between the probability mass on 0 and the adjacent count probabilities (i.e. those on $1,\dots,m$), as otherwise the penalization will shrink the estimate of $\pi_{i,0}$ and increase its neighboring count probabilities as to ensure smoothness of the resulting state-dependent p.m.f.s, which in case of a genuine excess of zeros can be undesirable (see Section 4.2 for an example of this problem). The penalty term in (3) can then be replaced by the inflation-adjusted penalty term $\sum_{i=1}^{N}\lambda_{i}\sum_{k=m+1}^{K}(\Delta^{m}\pi_{i,k})^{2}$, such that the $\pi_{i,0}$’s are estimated without any constraints related to the smoothness of the resulting state-dependent p.m.f.s.

Regarding the choice of $K$, i.e. the size of the support on which the conditional p.m.f.s are estimated, it is again necessary to use a $K$ greater than or equal to the highest count observed. However, with unbounded counts, it does in fact make sense to choose a somewhat larger $K$, as the absence of values greater than $K$ in the training data does not guarantee that such values will not occur in future. While without penalization, a nonparametric approach would estimate the probability of such future events to be 0, the penalized approach will place positive probability on counts slightly larger than the maximal value observed due to the enforced smoothing.

2.4 Model fitting

For a given set of smoothing parameters, maximum penalized likelihood estimates for the model parameters can be obtained by numerically maximizing the penalized log-likelihood (3) using some Newton-Raphson-type optimization routine, e.g. the R function nlm (R Core Team, 2017). To increase the chance of having found the global rather than a local maximum, a multiple-start-point strategy can be applied, where the penalized log-likelihood (3) is maximized from different, possibly randomly selected starting values, then choosing the estimate corresponding to the highest penalized log-likelihood (Zucchini et al., 2016).

In our model, all parameters to be estimated are probabilities. The corresponding constraints to be satisfied are as follows:

1. (1)

All model parameters to be estimated lie in the interval $[0,1]$; 2. (2)

the rows of the t.p.m. sum to one, i.e. $\sum_{j=1}^{N}\gamma_{ij}=1$, $i=1,\dots,N$; 3. (3)

the initial state probabilities sum to one, i.e. $\sum_{i=1}^{N}\delta_{i}=1$; 4. (4)

the count probabilities of each state-dependent distribution sum to one, i.e. $\sum_{k=0}^{K}\pi_{i,k}$ $=1$, $i=1,\dots,N$.

To ensure that these restrictions are satisfied, the constrained model parameters can be transformed into unconstrained ones using multinomial logit links, i.e.

[TABLE]

then maximizing the penalized log-likelihood (3) with respect to the unconstrained parameters $\gamma_{ij}^{*},\delta_{i}^{*}$, and $\pi_{i,k}^{*}$, fixing $\gamma_{ii}^{*}=\delta_{1}^{*}=\pi_{i,0}^{*}=0$, $i=1\ldots,N$, for identifiability. The constrained parameters are then obtained by applying the transformation (4). Regarding identifiability in general, and specifically in case of the very flexible model formulation considered here, Alexandrovich et al. (2016) show that for the HMM to be identifiable it is sufficient if the t.p.m. has full rank and the state-dependent distributions are distinct, conditions that can be expected to be met in most practical settings where HMMs seem natural candidate models.

An adequate choice of the smoothing parameters is crucial for finding a good balance between model fit and estimator variance. We here adopt the cross-validation approach as presented in Langrock et al. (2015), where we try to find the optimal vector $(\lambda_{1}^{*},\ldots,\lambda_{N}^{*})$ from a pre-specified, $N$-dimensional grid using a greedy search algorithm:

1. (1)

Choose an initial vector $(\lambda^{(0)*}_{1},\dots,\lambda^{(0)*}_{N})$ from the grid and set $r=0$; 2. (2)

calculate the average out-of-sample log-likelihood for $(\lambda^{(r)*}_{1},\dots,\lambda^{(r)*}_{N})$ and each direct neighbor on the grid; 3. (3)

from these values choose the vector $(\lambda^{(r+1)*}_{1},\dots,\lambda^{(r+1)*}_{N})$ as the one that yielded the highest out-of-sample log-likelihood averaged across folds; 4. (4)

repeat steps 2. and 3. until $(\lambda^{(r+1)*}_{1},\dots,\lambda^{(r+1)*}_{N})=(\lambda^{(r)*}_{1},\dots,\lambda^{(r)*}_{N})$.

In step (2), the out-of-sample log-likelihood is evaluated as follows: First, for any given smoothing parameter vector, the out-of-sample observations — e.g. 5% of all data points in each of 20 folds — are treated as missing data for model training using maximum penalized likelihood estimation, hence replacing the corresponding $\mathbf{P}(y_{t})$-matrices in the likelihood function given in (1) by identity matrices. Then the out-of-sample log-likelihood can be calculated in the same way, now treating the in-sample observations as missing data and using the estimated model parameters to calculate the out-of-sample (unpenalized) log-likelihood.

3 Simulation experiment

To demonstrate the feasibility of the suggested approach, and to illustrate its potential usefulness, we conducted the following simulation experiment: In each of $200$ simulation runs, we first simulated $T=500$ realizations from a $2$-state Markov chain with initial (stationary) distribution $\boldsymbol{\delta}=(0.5,0.5)$ and t.p.m.

[TABLE]

Conditional on the simulated states, the observations were then drawn from either of the following distributions: A Conway-Maxwell-Poisson (CMP) distribution (when the state process was in state 1), or a 2-component mixture of a Poisson and a CMP distribution (when in state 2). Compared to the Poisson distribution, the CMP distribution has an additional parameter, which allows to account for under- or overdispersion relative to the Poisson. The shape of the two state-dependent p.m.f.s, and additionally an example time series as simulated in one of the simulation runs, are displayed in panels (1) and (3) in Figure 2. The marginal distribution could be fairly well captured by a mixture of two Poisson or two negative binomial distributions (c.f. panels (2a)–(2b) in Figure 2), such that a 2-state Poisson or negative binomial HMM would seem to be a natural choice for modeling. However, the underlying state-dependent distributions do in fact substantially deviate from a Poisson, exhibiting underdispersion in state 1, and overdispersion as well as bimodality in state 2. This fairly complex model formulation was selected as to show the full potential of the new estimation framework, but also to highlight potential pitfalls that may occur when choosing simplistic parametric model formulations based in particular on an inspection of the marginal distribution of the data.

We consider four different model specifications: (1a) The true parametric model, as a benchmark only (noting that, in practice, a model specification as complex as the given one effectively could not be guessed based on an inspection of the marginal distribution of the data; cf. Figure 2); (1b) an incorrectly specified parametric model with Poisson state-dependent distributions (which as discussed above would seem to be a reasonable model specification based on the marginal distribution); (2a) the unpenalized nonparametric estimation approach; and (2b) the nonparametric estimation approach with third-order difference penalty. For the latter two, the support of the state-dependent p.m.f.s was taken to be $\max(40,\max(y_{1},\dots,y_{T}))$ in each simulation run. The smoothing parameters required within model (2b) were selected by $20$-fold cross validation over the grid $(10^{1},10^{2},\dots,10^{6})\times(10^{1},10^{2},\dots,10^{6})$. The means of the smoothing parameters selected over the 200 runs were $8173.5$ (state 1) and $570610.0$ (state 2), respectively.

The fitted state-dependent p.m.f.s under the four different models are illustrated in Figure 3. It can be seen that the proposed penalized nonparametric approach, model (2b), produced estimates very similar to those obtained when using the true parametric model, i.e. model (1a) — estimator variance is slightly higher when using model (2b), and also there is some underestimation of peaks and overestimation of troughs, as would be expected when using this effectively nonparametric estimation approach. In any case, given that the true parametric model specification is unknown in practice, this first impression regarding the performance of the penalized nonparametric approach is encouraging. Regarding the other two competitors, the incorrectly specified parametric model, model (1b), clearly lacks the flexibility to capture the functional shape of the true state-dependent p.m.f.s (thus exhibiting a substantial bias), while the unpenalized nonparametric approach, model (2a), leads to a very much higher estimator variance due to overfitting.

To more formally assess the performance of the suggested approach, we calculated three different performance measures as obtained for each of the four different estimation methods considered:

• (1)

The average Kullback-Leibler divergences (KLDs) between the true and the estimated state-dependent distributions, averaged across all $200$ simulation runs, i.e.

[TABLE]

where $\hat{\pi}_{i,k}^{(r)}$ denotes the estimated state-dependent count probabilities in state $i$ as obtained in the $r$-th simulation run;

• (2)

the mean absolute errors (MAEs) of the estimated state transition probabilities,

[TABLE]

where $\hat{\gamma}_{ij}^{(r)}$, $i\neq j$, denotes the estimate obtained in the $r$-th simulation run;

• (3)

the average state misclassification rates (SMRs),

[TABLE]

where $\hat{s}_{t}^{(r)}$ denotes the Viterbi-decoded state at time $t$, and $s_{t}^{(r)}$ denotes the true state at time $t$, in the $r$-th simulation run.

The different values of these performance measures, as obtained under the four different estimation methods considered, are stated in Table 1. The average time required to fit the penalized nonparametric model (2b) was $9.1$ seconds, which is comparable to the average time required to fit the nonparametric model without penalization (2a) ($6.7$ seconds). The penalized log-likelihood (3) was evaluated in C++, and maximized using the R function nlm (R Core Team, 2017), on a 3.6 GHz Intel® Core™ i7 CPU.

Unsurprisingly, the incorrectly specified parametric model (1b) shows the worst performance, with the largest average KLDs (due to the lack of flexibility to capture the shape of the conditional p.m.f.s, in particular within state 2), which notably results in a very high average state misclassification rate ($\overline{\text{SMR}}=0.240$; most observations lower than $15$ are assigned to state $1$, although a considerable number of these observations were actually generated from state $2$) as well as relatively large average MAEs of the estimated transition probabilities ($0.122$ and $0.395$; an obvious consequence of a considerable proportion of observations in $[5,15]$ incorrectly allocated to state 1). The unpenalized nonparametric approach, (2a), as implemented in the R package hmm.discnp (Turner, 2018), shows a much better performance thanks to its flexibility in particular to capture the bimodality within state 2. Due to the substantial reduction in the variance, the inclusion of the penalty term within the suggested penalized nonparametric approach, (2b), further considerably improves in particular the average deviation of the fitted model from the true model, as measured using the average KLDs. As expected, the true parametric model, (1a), shows the best performance — in practice, it is however by no means guaranteed that a parametric model formulation close to the true data-generating process can easily be identified.

4 Real-data case studies

4.1 Major earthquake counts

To illustrate the suggested approach in a real-data application, we model the distribution of major earthquake counts over time. In particular, we consider the observed variable

[TABLE]

The data, which are used as a running example in Zucchini et al. (2016), cover the period from 1900 to 2006, thus comprising $T=107$ observations in total. As a benchmark for the approach developed in the present paper, we consider a 2-state Poisson HMM as presented in Zucchini et al. (2016). We additionally fit 2-state HMMs using the unpenalized as well as the novel penalized nonparametric estimation approach. For the latter, we used third-order difference penalties and selected the smoothing parameters via 20-fold cross validation over the grid $(10^{1},10^{2},\dots,10^{10})\times(10^{1},10^{2},\dots,10^{10})$, which led to the optimal values $\lambda_{1}^{*}=10^{8}$ and $\lambda_{2}^{*}=10^{9}$.

With the penalized nonparametric approach, (2b), the t.p.m. of the underlying Markov chain was estimated as

[TABLE]

such that the associated stationary distribution is $\boldsymbol{\delta}=(0.660,0.340)$, indicating that about $66$% and $34$% of the observations were generated in state $1$ and $2$, respectively. These estimates are very close to those obtained under the 2-state Poisson HMM (cf. Zucchini et al., 2016). However, as can be seen in Figure 4, the associated state-dependent p.m.f.s as obtained under these two approaches do clearly differ. In particular, the parametric model clearly lacks the flexibility to account for the overdispersion present in the data, particularly within state 2, as captured by the nonparametric approaches. In contrast, due to the short length of the time series, the unpenalized nonparametric approach heavily overfits the data, which demonstrates the need for roughness penalization. In particular, the values 9 and 33 are both assigned a conditional probability of exactly zero in either of the two states — simply because these counts did not occur between 1900 and 2006 — such that according to the model these values also cannot occur in future years, which is implausible. Despite the differences in the estimated conditional p.m.f.s, each of the approaches identifies essentially the same pattern, where state 1 may be interpreted as a calm regime with relatively low seismic activity, whereas state 2 corresponds to periods exhibiting relatively high seismic activity. This is further illustrated by means of the globally decoded state sequence under the penalized nonparametric approach, displayed in panel (3) of Figure 4.

For the penalized nonparametric approach, (2b), the quantile-quantile plot as well as the autocorrelation function of the normal ordinary pseudo-residuals are displayed in Figure 5. Ordinary pseudo-residuals indicate whether an observation is extreme relative to the conditional distribution under the fitted model, given all other observations (Zucchini et al., 2016). If the model fits the data well, then the normal ordinary pseudo-residuals approximately follow a standard normal distribution. In that respect, the quantile-quantile plot in Figure 5 does not reveal any problematic lack of fit, and the sample autocorrelation function of the pseudo-residuals reveals only little residual autocorrelation. Overall, the model fitted using the penalized nonparametric approach, (2b), shows a satisfactory goodness of fit. In particular, there is no indication that a third state needs to be included in the model, as is the case when using Poisson HMMs (see Chapter 6 in Zucchini et al., 2016).

4.2 Oceanic whitetip shark acceleration counts

Oceanic whitetip sharks (Carcharhinus longimanus) are circumtropical pelagic apex predators, which are currently listed as “critically endangered” in the Northwest Atlantic and “vulnerable” globally by the International Union for the Conservation of Nature (IUCN; Baum et al., 2006). Despite an increasing interest in conservation, knowledge about the movement behaviors of this species is very limited (Howey-Jordan et al., 2013). Here we are particularly interested in activity patterns. From a tri-axial accelerometer, deployed from 8–12 May 2013, we obtained roughly $10.1$ million observations of sway, surge, and heave accelerations at a $30$ Hz resolution. These measures were combined into measures of overall dynamic body acceleration (ODBA) by evaluating the sum over the absolute values of the three different variables, and subsequently averaged over $1$ second intervals (see Leos-Barajas et al. (2017) and Langrock et al. (2013) for further background on acceleration data).

In an experimental setting, Lear et al. (2017) show for another species from the Carchiarhinus-family, namely lemon sharks (Negaprion brevirostris), that cruising behavior (low activity) corresponds to ODBA values $\leq 0.05$, while bursts in activity correspond to ODBA values $>0.05$. To reduce the very large size of the data set while maintaining effectively all biologically relevant information on the shark’s behavioral dynamics, we thus defined the count variable

[TABLE]

where $\text{ODBA}_{s,t}$ denotes the ODBA value observed during the $s$-th second of the $t$-th minute. The observations $y_{t}$ — which are bounded counts with range $\{0,\ldots,60\}$ — thus simply indicate for how many seconds from any given minute the shark exhibited bursts in activity. The downsampling yielded a total of $T=5560$ observations to be modeled, corresponding to an observation period of about $3.9$ days.

Due to the way the count variable $y_{t}$ was built, there is an inflation of both 0s and 60s in the data, which within the penalized nonparametric estimation approach, (1c), we account for by not penalizing the differences between the corresponding probabilities and their respective neighbors (i.e. 1 & 2, and 58 & 59). For comparison purposes, we further analyzed the data using the unpenalized nonparametric approach, (1a), as well as the penalized nonparametric approach without taking the inflation of 0s and 60s into account, (1b). The smoothing parameters for (1c) were selected via 20-fold cross validation over the grid $(10^{1},10^{2},\dots,10^{10})\times(10^{1},10^{2},\dots,10^{10})$, which led to the optimal values $\lambda_{1}^{*}=10^{4}$ and $\lambda_{2}^{*}=10^{7}$ (for comparison, these values were also used for (1b)).

With the penalized nonparametric approach, (1c), i.e. accounting for the inflation of 0s and 60s, the t.p.m. of the underlying Markov chain was estimated as

[TABLE]

such that the associated stationary distribution is $\boldsymbol{\delta}=(0.697,0.303)$, indicating that about $70$% and $30$% of the observations were generated in state $1$ and $2$, respectively. The state-dependent p.m.f.s fitted with the different estimation approaches, and the globally decoded time series of acceleration counts as obtained under the penalized nonparametric approach, (1c), are displayed in Figure 6. Since it does not take the excess 0s and 60s into account when penalizing, using the standard penalty term, (1b), results in an underestimation of the probabilities of these counts, and an overestimation of the probabilities of their respective neighbors. Without any penalization, i.e. with (1a), as well as with the inflation-adjusted penalty term, i.e. with (1c), this bias is avoided. The latter approach yields slightly smoother state-dependent p.m.f.s., but overall the models estimated with (1a) and (1c) are very similar — this is due to the relatively large sample size, which reduces the importance of penalization. Focusing on the penalized nonparametric approach, (1c), we find that state $1$ captures low activity counts, with a high probability mass on 0 ($\hat{\pi}_{1,0}=0.546$), which may therefore be interpreted as a resting/cruising state. State $2$ relates to moderate and high activity counts, with a (relatively) high probability mass on 60 ($\hat{\pi}_{2,60}=0.052$), which may therefore be linked to a traveling/hunting behavior.

In terms of model checking, pseudo-residuals are not very helpful in this example due to the excess of zeros in the count data (the quantile-quantile plot, which is not shown here, is dominated by these). Instead, we simply simulate data from the model fitted using the penalized nonparametric approach, (1c), and compare the patterns found in the simulated data, in particular regarding the marginal distribution and the correlation structure, to those found in the empirical data (see Figure 7). Though slightly less formal than a residual analysis, such a simple check, which is similar in spirit to a posterior predictive check as commonly done in Bayesian analyses, is often very helpful in revealing any potential lack of fit; see also Section 2.4 in Weiß (2018). Indeed, we find that our model captures the marginal distribution of the data very well (cf. panels (1a)–(1b) in Figure 7) — which is no surprise given its nonparametric nature — but fails to replicate the strong serial correlation present in the real data (cf. panels (2a)–(2b) in Figure 7). In particular, the sequence of the first 500 acceleration counts, as shown in panel (2) of Figure 6, clearly indicates that our assumption of the observations being conditionally independent, given the states, is violated: activity levels seem to sometimes gradually decrease or increase over time, a pattern that is not accounted for within our model formulation. Furthermore panel (2) in Figure 6 indicates periodicity in the shark’s behavioral patterns, which could potentially be further investigated by incorporating corresponding covariates (e.g. the time of day) into the state process, or alternatively using hidden semi-Markov models, which would allow us to more explicitly model the duration of a stay in any given state (Langrock & Zucchini, 2012). However, given that this example is supposed to only illustrate the new methodology, we here refrain from considering more complex model formulations that could help to address these nonstandard patterns.

5 Discussion

We introduced a penalized nonparametric estimation framework for flexibly tailoring an HMM to a time series of counts. The proposed methodology was shown to be a powerful alternative to parametric HMMs, and may be superior when simple parametric state-dependent distributions are not able to capture some of the features in the data. Specifically, the increased flexibility to capture complex distributional shapes can improve the accuracy of time series forecasts, can reduce state misclassification rates, and can help to avoid making biased inference related, in particular, to the dynamics of the state process. In any case, the approach can also be regarded as an exploratory tool, which can be applied when it is unclear which distributional family should be used for the state-dependent distributions. However, a simple parametric family is to be preferred whenever appropriate, as it will usually be easier to implement and to interpret, and also the computational effort will be much lower than when using the penalized nonparametric approach.

A notorious difficulty with HMMs, which can partly be addressed using the suggested penalized nonparametric approach, is the choice of the number of states. When using model selection criteria to choose an adequate number of states, then these tend to point to models with more states than can plausibly be interpreted, especially (though not only) in ecological applications (Pohle et al., 2017). Inflexibility of (parametric) state-dependent distributions to capture particular features in the data is a common cause of this problem — such inflexibility can always be compensated for by including additional states that will then usually have no meaningful interpretation anymore (Langrock et al., 2015). In this respect, the effectively unlimited flexibility of the nonparametric approach can help to reduce the required number of states, which will often substantially improve interpretability. In fact, with the new approach only a single state is required to capture the marginal distribution, and potential additional states only need to be included if they help to capture the dependence structure. On the other hand, model selection using information criteria is in fact more difficult within the nonparametric estimation framework, as it is necessary to somehow derive the effective number of parameters used to fit the model (Langrock et al., 2013), a challenge we did not address in the present paper.

In our simulation experiment and real-data case studies, the size of the support on which the state-dependent p.m.f.s were to be modeled was moderate — the largest support was $\{0,1,\ldots,60\}$ (in the acceleration data example). If p.m.f.s on much larger supports, say on $\{1,\ldots,5000\}$, are to be modeled, then the high dimensionality of the parameter space can become problematic, especially with regard to the computing time. In those instances, one option would be to simply treat the data as stemming from a continuous distribution, building the state-dependent distributions as linear combinations of B-splines (Langrock et al., 2015). Alternatively, the data could be binned, say in intervals $I_{1}=\{1,\ldots,10\}$, $I_{2}=\{11,\ldots,20\}$, $\ldots$, $I_{500}=\{4991,\ldots,5000\}$, then estimating conditional p.m.f.s defined on the intervals instead of directly on the counts (in the given example thus reducing the size of the support by a factor 10). In fact, these considerations highlight that the approach considered is by no means restricted to time series of counts. Instead, any time series data where observations are at least of ordinal scale can, in principle, be modeled using the approach developed in this work (in the same spirit as presented in Simonoff, 1983). This we believe could be relevant in particular for data on Likert-type scales. Finally, our penalization approach could potentially also be adapted to different types of models for count time series; as an example, one may think of combining it with the semiparametric estimation approach for integer-valued autoregressive (INAR) models by Drost et al. (2009).

Acknowledgements

The authors wish to thank Yannis Papastamatiou and Yuuki Watanabe for providing the oceanic whitetip shark data.