<figure><img src="image_21.png" /><figcaption>**Fig. 12** Evolution of the estimated health indicator (HI) and corresponding operational states, for two training trajectories. Numbers 1, 2, 3, and 4 stand for *healthy*, *intermediate*, *faulty* and *failure* states, respectively. FD001 and FD003 are the two training datasets considered in our experiments (Color figure online)</figcaption></figure>

the doubt on the switching time-step location between states *i* and *j* was explicited. Thus, if $t \in [t_{i \to j}-5, t_{i \to j}+5]$ then $\sigma_t = \{i,j\}$; otherwise $\sigma_t$ equals the state provided by the segmentation.

Thus, PHMC[-MLAR] models were trained on CMAPSS training datasets with the partial annotations obtained previously. The number of states K was fixed at 4 and Gaussian white noises were considered. Moreover, amongst the 21 features (time series) that make up each trajectory, only those features {2, 3, 4, 7, 9, 11, 12, 14} were used, that show consistent monotonic degradation trends (Wang et al., 2008) and/or present “the highest content of domain-specific information relating to the influence of fault occurrences” (Aremu et al., 2020).

### 8.2.2 RUL prediction

Once PHMC[-MLAR] model parameters had been estimated, RUL prediction for each single testing trajectory was performed using the following three-step procedure.

(i) *Step 1: production of R completed partial degradation patterns, for each testing trajectory.*

Let *T* the length of the testing trajectory under consideration, which degradation pattern is known up to time-step *T*. We wished to produce several evolutions of the partial degradation pattern. We could not rely on the optimal point prediction (classical method), which provides only one value (mean). Instead, we sampled from the predictive density (defined in Sect. 6, last paragraph) which characterizes each time-step *T* + *h*, for *h* ranging from 1 to *H*. An iterative sampling of the predictive densities allowed us to generate one possible evolution of this degradation pattern (from time-step *T* + 1 to *T* + *H*). Besides, in order to obtain only reasonable completed patterns, the sampling was restricted to the values belonging to the interval [*Q*25, *Q*75], where *Q*25 and *Q*75 are the predictive-density first and third quartiles. These quantiles were estimated by Monte Carlo simulations. Note that the predictive density of PHMC is similar to that of PHMC-MLAR with the difference that in PHMC model, the dependency on past observations is removed.

Performing the previous operation *R* = 100 times, we generated *R* evolutions of the observed partial degradation pattern. In the sequel, the resulting completed trajectories, each of length *T* + *H*, are referred to as *completed degradation patterns*. The maximum
<figure><img src="image_17.png" /><figcaption>**Fig. 12** Evolution of the estimated health indicator (HI) and corresponding operational states, for two training trajectories. Numbers 1, 2, 3, and 4 stand for *healthy*, *intermediate*, *faulty* and *failure* states, respectively. FD001 and FD003 are the two training datasets considered in our experiments (Color figure online)</figcaption></figure>

the doubt on the switching time-step location between states *i* and *j* was explicited. Thus, if $t \in [t_{i \to j}-5, t_{i \to j}+5]$ then $\sigma_t = \{i,j\}$; otherwise $\sigma_t$ equals the state provided by the segmentation.

Thus, PHMC[-MLAR] models were trained on CMAPSS training datasets with the partial annotations obtained previously. The number of states K was fixed at 4 and Gaussian white noises were considered. Moreover, amongst the 21 features (time series) that make up each trajectory, only those features {2, 3, 4, 7, 9, 11, 12, 14} were used, that show consistent monotonic degradation trends (Wang et al., 2008) and/or present “the highest content of domain-specific information relating to the influence of fault occurrences” (Aremu et al., 2020).

### 8.2.2 RUL prediction

Once PHMC[-MLAR] model parameters had been estimated, RUL prediction for each single testing trajectory was performed using the following three-step procedure.

(i) *Step 1: production of R completed partial degradation patterns, for each testing trajectory.*

Let *T* the length of the testing trajectory under consideration, which degradation pattern is known up to time-step *T*. We wished to produce several evolutions of the partial degradation pattern. We could not rely on the optimal point prediction (classical method), which provides only one value (mean). Instead, we sampled from the predictive density (defined in Sect. 6, last paragraph) which characterizes each time-step *T* + *h*, for *h* ranging from 1 to *H*. An iterative sampling of the predictive densities allowed us to generate one possible evolution of this degradation pattern (from time-step *T* + 1 to *T* + *H*). Besides, in order to obtain only reasonable completed patterns, the sampling was restricted to the values belonging to the interval [*Q*₂₅, *Q*₇₅], where *Q*₂₅ and *Q*₇₅ are the predictive-density first and third quartiles. These quantiles were estimated by Monte Carlo simulations. Note that the predictive density of PHMC is similar to that of PHMC-MLAR with the difference that in PHMC model, the dependency on past observations is removed.

Performing the previous operation *R* = 100 times, we generated *R* evolutions of the observed partial degradation pattern. In the sequel, the resulting completed trajectories, each of length *T* + *H*, are referred to as *completed degradation patterns*. The maximum
<figure><img src="image_22.png" /><figcaption>(a) 11<sup>th</sup> testing trajectory, feature 3</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) 11<sup>th</sup> testing trajectory, feature 12</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) 47<sup>th</sup> training trajectory, feature 3</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) 47<sup>th</sup> training trajectory, feature 12</figcaption></figure>

**Fig. 13** Completion of partial degradation pattern for testing trajectories. **a** and **b** Illustration with 11th testing trajectory of dataset FD001. The solid line indicates the known degradation pattern - or partial degradation pattern -, that encompasses time-steps 1 to *T*, the length of the testing trajectory. The dotted line denotes one simulated evolution generated by PHMC-MLAR(K=4, p = 7) model. It encompasses time-steps *T* + 1 to *T* + *H*, where *H* is the maximum remaining useful life (RUL) value observed in the testing dataset. **c** and **d** Illustration with 47th training trajectory of dataset FD001. In our experiments, no completion of partial degradation pattern is achieved for training trajectories. Subfigures **c** and **d** are only provided here to show how simulated evolution generated by PHMC-MLAR(K = 4, p = 7) and true (ground truth) evolution compare. For this illustration, time-step *T*′ was set at *T* − *H* where *T* is the length of the training trajectory and *T* coincides with engine failure. The solid line indicates the true evolution encompassing time-steps 1 to *T*′. The dashed line denotes the true evolution as from *T*′ + 1 till failure at time-step *T*. The dotted line indicates one simulated evolution generated by PHMC-MLAR(K = 4, p = 7) model, in time interval [*T*′ + 1, *T*] (Color figure online)

forecast horizon *H* was set at 145, which is the maximum RUL value observed within the
testing datasets. An example of such completed degradation pattern is depicted in Fig. 13.

(ii) Step 2: segmentation of the R completed degradation patterns, production of R esti-
mates of the true RUL.

For each completed degradation pattern (previously generated), we know that up to
time-step *T*, the *failure* state has null probability since at time-step *T*, the true RUL is
not null. Namely, at each time-step *t* ≤ *T* of a testing trajectory, the annotation *σ*<sub>*t*</sub> is set
to {1, 2, 3} (We remind the reader that for each time-step of a training trajectory, *σ*<sub>*t*</sub> was
set to {*i*,*j*} around each switch *i* → *j*, and was set to the unique state obtained through
segmentation otherwise). We recall that such partial knowledge is taken into account by
our variant of the Viterbi algorithm. Afterwards, the R patterns produced for each test-
ing trajectory were each segmented into *healthy*, *intermediate*, *faulty* and *failure* states.
Let *ŷ*<sub>3→4</sub> the time-step at which the engine switches from *faulty* state to *failure* state. So,
*ŷ*<sub>3→4</sub> + *p* provides an estimation of engine end life time-step. To note, the term (+*p*) takes
into account the *p* initial observations whose states cannot be inferred because of the
autoregressive dynamics. Thus, the RUL estimate writes:
<figure><img src="image_18.png" /><figcaption>(a) 11<sup>th</sup> testing trajectory, feature 3</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) 11<sup>th</sup> testing trajectory, feature 12</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) 47<sup>th</sup> training trajectory, feature 3</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) 47<sup>th</sup> training trajectory, feature 12</figcaption></figure>

**Fig. 13** Completion of partial degradation pattern for testing trajectories. **a** and **b** Illustration with 11th testing trajectory of dataset FD001. The solid line indicates the known degradation pattern - or partial degradation pattern -, that encompasses time-steps 1 to *T*, the length of the testing trajectory. The dotted line denotes one simulated evolution generated by PHMC-MLAR(K=4, p = 7) model. It encompasses time-steps *T* + 1 to *T* + *H*, where *H* is the maximum remaining useful life (RUL) value observed in the testing dataset. **c** and **d** Illustration with 47th training trajectory of dataset FD001. In our experiments, no completion of partial degradation pattern is achieved for training trajectories. Subfigures **c** and **d** are only provided here to show how simulated evolution generated by PHMC-MLAR(K = 4, p = 7) and true (ground truth) evolution compare. For this illustration, time-step *T*′ was set at *T* − *H* where *T* is the length of the training trajectory and *T* coincides with engine failure. The solid line indicates the true evolution encompassing time-steps 1 to *T*′. The dashed line denotes the true evolution as from *T*′ + 1 till failure at time-step *T*. The dotted line indicates one simulated evolution generated by PHMC-MLAR(K = 4, p = 7) model, in time interval [*T*′ + 1, *T*] (Color figure online)

forecast horizon *H* was set at 145, which is the maximum RUL value observed within the
testing datasets. An example of such completed degradation pattern is depicted in Fig. 13.

(ii) Step 2: segmentation of the R completed degradation patterns, production of R esti-
mates of the true RUL.

For each completed degradation pattern (previously generated), we know that up to
time-step *T*, the *failure* state has null probability since at time-step *T*, the true RUL is
not null. Namely, at each time-step *t* ≤ *T* of a testing trajectory, the annotation *σ*<sub>*t*</sub> is set
to {1, 2, 3} (We remind the reader that for each time-step of a training trajectory, *σ*<sub>*t*</sub> was
set to {*i*,*j*} around each switch *i* → *j*, and was set to the unique state obtained through
segmentation otherwise). We recall that such partial knowledge is taken into account by
our variant of the Viterbi algorithm. Afterwards, the R patterns produced for each test-
ing trajectory were each segmented into *healthy*, *intermediate*, *faulty* and *failure* states.
Let *ŷ*<sub>3→4</sub> the time-step at which the engine switches from *faulty* state to *failure* state. So,
*ŷ*<sub>3→4</sub> + *p* provides an estimation of engine end life time-step. To note, the term (+*p*) takes
into account the *p* initial observations whose states cannot be inferred because of the
autoregressive dynamics. Thus, the RUL estimate writes:
$$ \hat{rul} = \hat{t}_{3 \to 4} + p - T. \qquad (31) $$

When the switching to the *failure state* was not achieved between the beginning and the end of the pattern, RUL estimate was set at *H*.

By doing so for the *R* patterns, we obtained *R* estimates of the true RUL denoted by $\{\hat{rul}_r\}_{r=1,...,R}$. Figure 14 presents the distribution of $\hat{rul}_r$'s obtained for two different testing trajectories.

(iii) Step 3: estimation of a final RUL from the R RUL estimates of a given testing trajectory.

To compute the final estimate of the RUL, the previously computed $\hat{rul}$'s were aggregated using a fusion rule $\mathcal{F}$:

$$ R\hat{UL} = \mathcal{F}(\hat{rul}_1, \hat{rul}_2, \dots, \hat{rul}_R). \qquad (32) $$

The mean and median functions were first considered: the corresponding fusion rules are referred to as $\mathcal{F}_{mean}$ and $\mathcal{F}_{median}$. Then, several combinations of the minimum and maximum of the $\hat{rul}$'s were also considered, as suggested by Ramasso (2014):

$$ \begin{gathered} \mathcal{F}_{a\_min\_max} = a \times \min(\hat{rul}\_s) + (1-a) \times \max(\hat{rul}\_s), \qquad (33) \\ \text{for } a \in \mathcal{A} = \{0, 0.1, 0.2, 0.3, 0.4, 0.5, 13/23, 0.6, 0.7, 0.8, 0.9, 1\}. \end{gathered} $$

The choice of $a = 13/23$ will be explained in Sect. 8.2.3.

The set of final RUL estimates, considering different fusion functions, is denoted

$$ \mathbf{R\hat{U}L} = \{\mathcal{F}_f, f \in \{\text{mean}, \text{median}, \text{max}, a\_min\_max_{a \in \mathcal{A}}\}\}. \qquad (34) $$

### 8.2.3 Performance metrics

To evaluate the performance of our model in RUL prediction, two performance metrics were used: the score function (SF) and the root mean square error (RMSE).

The score function defined by Saxena et al. (2008) has been largely used in the literature of CMAPSS datasets. This function writes:

<figure><img src="image_26.png" /><figcaption>**Fig. 14** Density of *R* = 100 estimates of the true RUL for testing trajectories. RUL: remaining useful life. **a** 17th trajectory from the FD001 dataset. **b** 26th trajectory from the FD001 dataset. Each estimate is obtained from a completed degradation pattern simulated from the PHMC-MLAR(K = 4, p = 7) model (Color figure online)</figcaption></figure>
$$ \hat{rul} = \hat{t}_{3 \to 4} + p - T. \qquad (31) $$

When the switching to the failure state was not achieved between the beginning and the end of the pattern, RUL estimate was set at $H$.

By doing so for the $R$ patterns, we obtained $R$ estimates of the true RUL denoted by $\{\hat{rul}_r\}_{r=1,...,R}$. Figure 14 presents the distribution of $\hat{rul}_r$'s obtained for two different testing trajectories.

(iii) Step 3: estimation of a final RUL from the R RUL estimates of a given testing trajectory.

To compute the final estimate of the RUL, the previously computed $\hat{rul}$'s were aggregated using a fusion rule $\mathcal{F}$:

$$ R\hat{UL} = \mathcal{F}(\hat{rul}_1, \hat{rul}_2, \dots, \hat{rul}_R). \qquad (32) $$

The mean and median functions were first considered: the corresponding fusion rules are referred to as $\mathcal{F}_{mean}$ and $\mathcal{F}_{median}$. Then, several combinations of the minimum and maximum of the $\hat{rul}$'s were also considered, as suggested by Ramasso (2014):

$$ \begin{gathered} \mathcal{F}_{a\_min\_max} = a \times \min(\hat{rul}\_s) + (1-a) \times \max(\hat{rul}\_s), \\ \text{for } a \in \mathcal{A} = \{0, 0.1, 0.2, 0.3, 0.4, 0.5, 13/23, 0.6, 0.7, 0.8, 0.9, 1\}. \end{gathered} \qquad (33) $$

The choice of $a = 13/23$ will be explained in Sect. 8.2.3.

The set of final RUL estimates, considering different fusion functions, is denoted

$$ \mathbf{R\hat{UL}} = \{\mathcal{F}_f, f \in \{\text{mean}, \text{median}, \text{max}, a\_min\_max_{a \in \mathcal{A}}\}\}. \qquad (34) $$

### 8.2.3 Performance metrics

To evaluate the performance of our model in RUL prediction, two performance metrics were used: the score function (SF) and the root mean square error (RMSE).

The score function defined by Saxena et al. (2008) has been largely used in the literature of CMAPSS datasets. This function writes:

<figure><img src="image_22.png" /><figcaption>Fig. 14 Density of $R = 100$ estimates of the true RUL for testing trajectories. RUL: remaining useful life. a 17th trajectory from the FD001 dataset. b 26th trajectory from the FD001 dataset. Each estimate is obtained from a completed degradation pattern simulated from the PHMC-MLAR($K = 4, p = 7$) model (Color figure online)</figcaption></figure>
$$SF = \sum_{i=1}^{N_{test}} SF_i, \quad SF_i = \begin{cases} \exp\left(-\frac{d_i}{13}\right) - 1, & \text{if } d_i < 0 \\ \exp\left(\frac{d_i}{10}\right) - 1, & \text{if } d_i \ge 0 \end{cases} \qquad (35)$$

where $N_{test}$ is the number of testing trajectories, $d_i = R\hat{U}L_i - RUL_i$ denotes the difference between the estimated RUL and true RUL for $i$th testing trajectory.

We point out that this score function assigns higher penalties on late predictions (over-estimations of RUL). Now function SF has been defined, we return to the choice of $a = 13/23$ in Eq. 33. This specific value was proposed by Wang et al. (2008). The motivation of these authors originates in the definition of the SF score function: this function penalizes over-estimations of RULs by 1/10 and penalizes under-estimations of RULs by 1/13. Over-estimations (late predictions) are therefore more penalized than under-estimations (early predictions). The authors transformed the 1/10 and 1/13 penalties into weights that sum to 1: 13/23 was assigned to the minimin in Eq. 33, since $min(\hat{r}ul's)$ favours RUL under-estimation, and 10/23 was assigned to the maximum in this same equation, since $max(\hat{r}ul's)$ favours over-estimation.

The RMSE of RUL prediction is computed as follows:

$$RMSE = \sqrt{\frac{1}{N_{test}} \sum_{i=1}^{N_{test}} d_i^2} \qquad (36)$$

Both performance metrics SF and RMSE have to be minimized.

## 8.3 Results and analysis

This subsection first describes the **best models selected** for PHMC-MLAR ($K = 4, p$) and MSAR($p$), on each of datasets FD001 and FD003. Then, we show and discuss the accuracies obtained for the **feature prediction** task performed with PHMC-MLAR, MSAR and PHMC, on both datasets. Third, we conduct a detailed analysis of the **RUL prediction performances** obtained for PHMC-MLAR and PHMC, on dataset FD001. The subsection ends with the comparison, on dataset FD001, of PHMC-MLAR against six RUL prediction methods recently reported in the literature.

### 8.3.1 Model selection

Remember that the number of states $K$ was fixed at 4. In this experiment, we identified the best PHMC-MLAR and MSAR models by varying the autoregressive order $p$ among values {1, 2, ..., 14}; the BIC score was used for the selection. Figure 15 displays the BIC scores of each model for both training datasets FD001 and FD003. For each model and each dataset, the value that provides the lowest BIC score was selected. Thus the best MSAR models are obtained for $p = 10$ in both datasets (see Fig. 15b and d). Regarding the PHMC-MLAR model, the BIC score allows to select $p = 13$ for FD003 dataset (see Fig. 15c). However, for FD001 dataset, models PHMC-MLAR($p = 6$), PHMC-MLAR($p = 7$) and PHMC-MLAR($p = 8$) obtained very close BIC scores (472.509; 472.433 and 472.501, respectively), making the selection difficult (see Fig. 15a). Therefore, using the testing trajectories, a second selection has been performed among these three models based on their 30-step ahead prediction
$$SF = \sum_{i=1}^{N_{test}} SF_i, \quad SF_i = \begin{cases} \exp\left(-\frac{d_i}{13}\right) - 1, & \text{if } d_i < 0 \\ \exp\left(\frac{d_i}{10}\right) - 1, & \text{if } d_i \ge 0 \end{cases} \quad (35)$$

where $N_{test}$ is the number of testing trajectories, $d_i = R\hat{U}L_i - RUL_i$ denotes the difference between the estimated RUL and true RUL for $i$th testing trajectory.

We point out that this score function assigns higher penalties on late predictions (over-estimations of RUL). Now function SF has been defined, we return to the choice of $a = 13/23$ in Eq. 33. This specific value was proposed by Wang et al. (2008). The motivation of these authors originates in the definition of the SF score function: this function penalizes over-estimations of RULs by 1/10 and penalizes under-estimations of RULs by 1/13. Over-estimations (late predictions) are therefore more penalized than under-estimations (early predictions). The authors transformed the 1/10 and 1/13 penalties into weights that sum to 1: 13/23 was assigned to the minimin in Eq. 33, since $min(\hat{r}ul's)$ favours RUL under-estimation, and 10/23 was assigned to the maximum in this same equation, since $max(\hat{r}ul's)$ favours over-estimation.

The RMSE of RUL prediction is computed as follows:

$$RMSE = \sqrt{\frac{1}{N_{test}} \sum_{i=1}^{N_{test}} d_i^2} \quad (36)$$

Both performance metrics SF and RMSE have to be minimized.

## 8.3 Results and analysis

This subsection first describes the **best models selected** for PHMC-MLAR ($K = 4, p$) and MSAR($p$), on each of datasets FD001 and FD003. Then, we show and discuss the accuracies obtained for the **feature prediction** task performed with PHMC-MLAR, MSAR and PHMC, on both datasets. Third, we conduct a detailed analysis of the **RUL prediction performances** obtained for PHMC-MLAR and PHMC, on dataset FD001. The subsection ends with the comparison, on dataset FD001, of PHMC-MLAR against six RUL prediction methods recently reported in the literature.

### 8.3.1 Model selection

Remember that the number of states $K$ was fixed at 4. In this experiment, we identified the best PHMC-MLAR and MSAR models by varying the autoregressive order $p$ among values {1, 2, ..., 14}; the BIC score was used for the selection. Figure 15 displays the BIC scores of each model for both training datasets FD001 and FD003. For each model and each dataset, the value that provides the lowest BIC score was selected. Thus the best MSAR models are obtained for $p = 10$ in both datasets (see Fig. 15b and d). Regarding the PHMC-MLAR model, the BIC score allows to select $p = 13$ for FD003 dataset (see Fig. 15c). However, for FD001 dataset, models PHMC-MLAR($p = 6$), PHMC-MLAR($p = 7$) and PHMC-MLAR($p = 8$) obtained very close BIC scores (472.509; 472.433 and 472.501, respectively), making the selection difficult (see Fig. 15a). Therefore, using the testing trajectories, a second selection has been performed among these three models based on their 30-step ahead prediction
<figure><img src="image_27.png" /><figcaption>**Fig. 15** Model selection for the PHMC-MLAR(K = 4, p) and MSAR(p) models, for datasets FD001 and FD003, when the autoregressive order *p* varies from 1 to 14. The BIC score is used. Mind the different scales in the four subfigures (Color figure online)</figcaption></figure>

performances. This time, the model PHMC-MLAR(p = 7) was selected. We observe that this model also obtained the lowest BIC score.

### 8.3.2 Feature prediction performance

For each dataset, we compared the feature prediction accuracies of the previously selected PHMC-MLAR and MSAR models against that of the PHMC model. To this end, we considered short, medium and long-term forecasts: *h* = 5, 10, 20, 30. For each testing trajectory, we performed rolling *h*-step ahead forecasts, starting from *t* = 15 and considering increments of 5 time-steps, up to *t* + 5*i*\* with *i*\* the highest integer *i* such that *t* + 5*i* ≤ *T* − *h*. Thus, for each trajectory, we obtained one estimate of the *h*-step ahead forecast RMSE. We recall that for feature prediction, the RMSE is computed as shown in Eq. 29, with *N*<sub>rep</sub> = *i*\*.
To produce reliable estimates of RMSE, very short trajectories were removed (that is the ones for which less than 10 *h*-step ahead projections can be performed).

For the three models PHMC-MLAR, PHMC and MSAR, Table 6 (resp. Table 7) in Appendix B presents the mean and 95% confidence interval of RMSE for dataset 1 (respectively dataset 3), and for each pair (feature, prediction horizon).

The results show that PHMC-MLAR (respectively MSAR) obtains the lowest RMSE for features F9 and F14 (respectively features F2, F3, F7 and F12) on both datasets FD001 and FD003. For these two models, Table 4 displays the sum of the RMSE means computed over all features, at each forecast horizon. It can be highlighted that the two models have the same performance at short forecast horizon (that is *h* = 5). However, our model outperforms MSAR when medium and long-term forecasts are considered (that is for *h* = 10, 20, 30).
<figure><img src="image_23.png" /><figcaption>**Fig. 15** Model selection for the PHMC-MLAR(K = 4, p) and MSAR(p) models, for datasets FD001 and FD003, when the autoregressive order *p* varies from 1 to 14. The BIC score is used. Mind the different scales in the four subfigures (Color figure online)</figcaption></figure>

performances. This time, the model PHMC-MLAR(p = 7) was selected. We observe that this model also obtained the lowest BIC score.

### 8.3.2 Feature prediction performance

For each dataset, we compared the feature prediction accuracies of the previously selected PHMC-MLAR and MSAR models against that of the PHMC model. To this end, we considered short, medium and long-term forecasts: *h* = 5, 10, 20, 30. For each testing trajectory, we performed rolling *h*-step ahead forecasts, starting from *t* = 15 and considering increments of 5 time-steps, up to *t* + 5*i*\* with *i*\* the highest integer *i* such that *t* + 5*i* ≤ *T* − *h*. Thus, for each trajectory, we obtained one estimate of the *h*-step ahead forecast RMSE. We recall that for feature prediction, the RMSE is computed as shown in Eq. 29, with *N*<sub>rep</sub> = *i*\*.
To produce reliable estimates of RMSE, very short trajectories were removed (that is the ones for which less than 10 *h*-step ahead projections can be performed).

For the three models PHMC-MLAR, PHMC and MSAR, Table 6 (resp. Table 7) in Appendix B presents the mean and 95% confidence interval of RMSE for dataset 1 (respectively dataset 3), and for each pair (feature, prediction horizon).

The results show that PHMC-MLAR (respectively MSAR) obtains the lowest RMSE for features F9 and F14 (respectively features F2, F3, F7 and F12) on both datasets FD001 and FD003. For these two models, Table 4 displays the sum of the RMSE means computed over all features, at each forecast horizon. It can be highlighted that the two models have the same performance at short forecast horizon (that is *h* = 5). However, our model outperforms MSAR when medium and long-term forecasts are considered (that is for *h* = 10, 20, 30).
<table><caption>**Table 4** Global comparison of feature prediction performances for PHMC-MLAR and MSAR models</caption><thead><tr><th>Dataset</th><th>Model</th><th>h = 5</th><th>h = 10</th><th>h = 20</th><th>h = 30</th></tr></thead><tbody><tr><td rowspan="2">FD001</td><td>PHMC-MLAR(K = 4, p = 7)</td><td><strong>16.904</strong></td><td><strong>17.239</strong></td><td><strong>18.697</strong></td><td><strong>21.284</strong></td></tr><tr><td>MSAR(K = 4, p = 10)</td><td>16.949</td><td>17.486</td><td>19.356</td><td>22.126</td></tr><tr><td rowspan="2">FD003</td><td>PHMC-MLAR(K = 4, p = 13)</td><td><strong>17.316</strong></td><td><strong>17.651</strong></td><td><strong>18.544</strong></td><td><strong>20.008</strong></td></tr><tr><td>MSAR(K = 4, p = 10)</td><td>17.397</td><td>17.863</td><td>19.134</td><td>21.267</td></tr></tbody></table>

Sum of RMSE means over all features, at each forecast horizon *h*. The figures in bold highlight the minimum

Compared to the PHMC model, our model obtains much better prediction performance on both datasets. Figure 16 displays the percentage of RMSE improvement (that is RMSE reduction) obtained using our model instead of the PHMC model. This figure shows that our model improves the prediction performance over PHMC for all features whatever the forecast horizon. For dataset FD001 (respectively FD003), the improvement is greater than 40% for features F9 and F14 (respectively features F7, F9, F12 and F14), whatever the forecast horizon. From these results, we can conclude that, on CMAPSS datasets, the autoregressive dynamics included in our model results in better prediction performance. Thus, for these datasets, more information is obtained from the history of past observations compared with the multivariate data point sampled at a single time-step.

### 8.3.3 RUL prediction performance

We compared the RUL prediction performance of our proposal with those of the PHMC and several existing methods, on testing dataset FD001. We point out that, on this task, a comparison with MSAR would have been biased since the MSAR framework is fully unsupervised. Figure 17 depicts the performance metrics (SF, RMSE) of PHMC(K = 4) model and PHMC-MLAR(K = 4, p = 7) model, for the different fusion functions considered in this work. The results show that the fusion function $F_{0.7\_min\_max}$ yields the best performance for PHMC-MLAR(K = 4, p = 7) model, which obtains the scores $(SF, RMSE) = (472, 17.49)$, whereas PHMC(K = 4) model shows the smallest scores $(SF, RMSE) = (17066, 34.32)$, when considering fusion function $F_{0.9\_min\_max}$.

<figure><img src="image_28.png" /><figcaption>**Fig. 16** Comparison of feature prediction performances for PHMC-MLAR and PHMC models. Percentages of RMSE improvement of PHMC-MLAR model over PHMC model compared along various pairs (feature, prediction horizon). **a** PHMC-MLAR(K = 4, p = 7) versus PHMC(K = 4), experiment with dataset FD001. **b** PHMC-MLAR(K = 4, p = 13) versus PHMC(K = 4), experiment with dataset FD003. An RMSE improvement corresponds to a decrease (Color figure online)</figcaption></figure>
<table><caption>**Table 4** Global comparison of feature prediction performances for PHMC-MLAR and MSAR models</caption><thead><tr><th>Dataset</th><th>Model</th><th>h = 5</th><th>h = 10</th><th>h = 20</th><th>h = 30</th></tr></thead><tbody><tr><td rowspan="2">FD001</td><td>PHMC-MLAR(K = 4, p = 7)</td><td><strong>16.904</strong></td><td><strong>17.239</strong></td><td><strong>18.697</strong></td><td><strong>21.284</strong></td></tr><tr><td>MSAR(K = 4, p = 10)</td><td>16.949</td><td>17.486</td><td>19.356</td><td>22.126</td></tr><tr><td rowspan="2">FD003</td><td>PHMC-MLAR(K = 4, p = 13)</td><td><strong>17.316</strong></td><td><strong>17.651</strong></td><td><strong>18.544</strong></td><td><strong>20.008</strong></td></tr><tr><td>MSAR(K = 4, p = 10)</td><td>17.397</td><td>17.863</td><td>19.134</td><td>21.267</td></tr></tbody></table>

Sum of RMSE means over all features, at each forecast horizon *h*. The figures in bold highlight the minimum

Compared to the PHMC model, our model obtains much better prediction performance on both datasets. Figure 16 displays the percentage of RMSE improvement (that is RMSE reduction) obtained using our model instead of the PHMC model. This figure shows that our model improves the prediction performance over PHMC for all features whatever the forecast horizon. For dataset FD001 (respectively FD003), the improvement is greater than 40% for features F9 and F14 (respectively features F7, F9, F12 and F14), whatever the forecast horizon. From these results, we can conclude that, on CMAPSS datasets, the autoregressive dynamics included in our model results in better prediction performance. Thus, for these datasets, more information is obtained from the history of past observations compared with the multivariate data point sampled at a single time-step.

### 8.3.3 RUL prediction performance

We compared the RUL prediction performance of our proposal with those of the PHMC and several existing methods, on testing dataset FD001. We point out that, on this task, a comparison with MSAR would have been biased since the MSAR framework is fully unsupervised. Figure 17 depicts the performance metrics (SF, RMSE) of PHMC(K = 4) model and PHMC-MLAR(K = 4, p = 7) model, for the different fusion functions considered in this work. The results show that the fusion function $F_{0.7\_min\_max}$ yields the best performance for PHMC-MLAR(K = 4, p = 7) model, which obtains the scores (SF, RMSE) = (472, 17.49), whereas PHMC(K = 4) model shows the smallest scores (SF, RMSE) = (17066, 34.32), when considering fusion function $F_{0.9\_min\_max}$.

<figure><img src="image_24.png" /><figcaption>**Fig. 16** Comparison of feature prediction performances for PHMC-MLAR and PHMC models. Percentages of RMSE improvement of PHMC-MLAR model over PHMC model compared along various pairs (feature, prediction horizon). **a** PHMC-MLAR(K = 4, p = 7) versus PHMC(K = 4), experiment with dataset FD001. **b** PHMC-MLAR(K = 4, p = 13) versus PHMC(K = 4), experiment with dataset FD003. An RMSE improvement corresponds to a decrease (Color figure online)</figcaption></figure>
<table><caption>**Table 5** RUL prediction</caption><thead><tr><th>Method</th><th>RMSE</th><th>SF</th></tr></thead><tbody><tr><td>MLP Switching Kalman Filter Ensemble (Lim et al., 2015)</td><td>[18.4, 18.9]</td><td>[560, 580]</td></tr><tr><td>Extreme Learning Machine Fuzzy Clustering (Javed et al., 2015)</td><td>---</td><td>1046</td></tr><tr><td>Vanilla DLSTM (Wu et al., 2018)</td><td>19.74</td><td>---</td></tr><tr><td>Deep CNN (Sateesh Babu et al., 2016)</td><td>18.45</td><td>1287</td></tr><tr><td>DLSTM (Zhao et al., 2019)</td><td>[19.3, 23.8]</td><td>[750, 1200]</td></tr><tr><td>CEEMD DLSTM (Zhao et al., 2019)</td><td><strong>14.72</strong></td><td><strong>262</strong></td></tr><tr><td>PHMC(K = 4)</td><td>34.33 ± 0.15</td><td>17066 ± 590</td></tr><tr><td>The proposed method - PHMC-MLAR(K = 4, p = 7)</td><td><strong>17.49 ± 5.2</strong></td><td><strong>472 ± 42</strong></td></tr></tbody></table>

Comparison of PHMC-MLAR against six recent methods of the state-of-the-art. The figures in bold highlight the minima for RMSE values and SF values. The second minimum values are underlined

<figure><img src="image_29.png" /><figcaption>**Fig. 17** Comparison of RUL prediction performances for PHMC-MLAR and PHMC models. Mean and standard deviation of performance metrics (SF, RMSE) computed from ten executions of the methods based on PHMC-MLAR and PHMC. Performance metrics are presented for fusion functions $F_{mean}$, $F_{median}$ and $F_{a_min_max}$ for $a = \{0.5, 13/23, 0.6, 0.7, 0.8, 0.9\}$. The scores obtained by the other fusion functions are too high; therefore they are not displayed in this figure. Mind the different scales in the two subfigures (Color figure online)</figcaption></figure>

Overall, our model largely outperforms PHMC since its SF and RMSE values are much smaller: the RMSE was reduced by half and the SF was divided by 34, with respect to PHMC. When we go deeper into details and analyze RUL estimates, we found out that PHMC is more subject to late predictions of RUL (over-estimations), with 70% of testing trajectories for which the difference between the estimated and actual RULs is greater than 10, against only 24% for our model. This explains the very large difference between the SF of the two models (in comparison with RMSE) because the SF metric assigns higher penalties to late predictions of the RUL (see Eq. 35).

Finally, Table 5 compares the method proposed in this paper with six recent state-of-the-art methods used to predict the remaining useful life of a machine: Multi-Layer Perceptron (MLP) Switching Kalman filter ensemble (Lim et al., 2015), extreme learning machine fuzzy clustering (Javed et al., 2015), deep convolutional neural network (CNN) (Sateesh Babu et al., 2016), Vanilla long short-term memory (LSTM) (Wu et al., 2018), double LSTM (DLSTM) (Zhao et al., 2019), and complete ensemble empirical mode decomposition DLSTM (CEEMD DLSTM) (Zhao et al., 2019).
<table><caption>**Table 5** RUL prediction</caption><thead><tr><th>Method</th><th>RMSE</th><th>SF</th></tr></thead><tbody><tr><td>MLP Switching Kalman Filter Ensemble (Lim et al., 2015)</td><td>[18.4, 18.9]</td><td>[560, 580]</td></tr><tr><td>Extreme Learning Machine Fuzzy Clustering (Javed et al., 2015)</td><td>---</td><td>1046</td></tr><tr><td>Vanilla DLSTM (Wu et al., 2018)</td><td>19.74</td><td>---</td></tr><tr><td>Deep CNN (Sateesh Babu et al., 2016)</td><td>18.45</td><td>1287</td></tr><tr><td>DLSTM (Zhao et al., 2019)</td><td>[19.3, 23.8]</td><td>[750, 1200]</td></tr><tr><td>CEEMD DLSTM (Zhao et al., 2019)</td><td><b>14.72</b></td><td><b>262</b></td></tr><tr><td>PHMC(K = 4)</td><td>34.33 ± 0.15</td><td>17066 ± 590</td></tr><tr><td>The proposed method - PHMC-MLAR(K = 4, p = 7)</td><td><b>17.49 ± 5.2</b></td><td><b>472 ± 42</b></td></tr></tbody></table>

Comparison of PHMC-MLAR against six recent methods of the state-of-the-art. The figures in bold highlight the minima for RMSE values and SF values. The second minimum values are underlined

<figure><img src="image_25.png" /><figcaption>**Fig. 17** Comparison of RUL prediction performances for PHMC-MLAR and PHMC models. Mean and standard deviation of performance metrics (SF, RMSE) computed from ten executions of the methods based on PHMC-MLAR and PHMC. Performance metrics are presented for fusion functions $F_{mean}$, $F_{median}$ and $F_{a_min_max}$ for $a = \{0.5, 13/23, 0.6, 0.7, 0.8, 0.9\}$. The scores obtained by the other fusion functions are too high; therefore they are not displayed in this figure. Mind the different scales in the two subfigures (Color figure online)</figcaption></figure>

Overall, our model largely outperforms PHMC since its SF and RMSE values are much smaller: the RMSE was reduced by half and the SF was divided by 34, with respect to PHMC. When we go deeper into details and analyze RUL estimates, we found out that PHMC is more subject to late predictions of RUL (over-estimations), with 70% of testing trajectories for which the difference between the estimated and actual RULs is greater than 10, against only 24% for our model. This explains the very large difference between the SF of the two models (in comparison with RMSE) because the SF metric assigns higher penalties to late predictions of the RUL (see Eq. 35).

Finally, Table 5 compares the method proposed in this paper with six recent state-of-the-art methods used to predict the remaining useful life of a machine: Multi-Layer Perceptron (MLP) Switching Kalman filter ensemble (Lim et al., 2015), extreme learning machine fuzzy clustering (Javed et al., 2015), deep convolutional neural network (CNN) (Sateesh Babu et al., 2016), Vanilla long short-term memory (LSTM) (Wu et al., 2018), double LSTM (DLSTM) (Zhao et al., 2019), and complete ensemble empirical mode decomposition DLSTM (CEEMD DLSTM) (Zhao et al., 2019).
The results show that the method proposed in this work presents better performance than the first five comparison methods with a smaller RMSE and SF metrics. However, thanks to the extensive feature extraction procedure (CEEMD) used to build the input features of DLSTM, CEEMD DLSTM method presents better performance than our proposal. Note that our model is directly trained on noisy features recorded from sensors. Moreover, in contrast to CEEMD DLSTM method, our proposal provides an accurate (short, medium and long-term) prediction of features which can be used in practice to extract useful information about system operational states.

# 9 Conclusion

In this work, we have introduced the PHMC-MLAR model to analyze time series subject to switches in regimes. Our model is a generalization of the well-known Hidden Regime-Switching AutoRegressive (HRSAR) and Observed Regime-Switching AutoRegressive (ORSAR) models when regime-switching is modelled by a Markov Chain. Our model allows to handle the intermediate case where the state process is partially observed.

In the evaluation, we conducted experiments on simulated data and considered both inference performance and prediction accuracy. The results show that the partially observed states (when they represent a reasonable proportion) allow a better characterization of training data (reflected by greater log-likelihood), in comparison with the unsupervised case. An interesting characteristic of the PHMC-MLAR model is that the partially observed states allow faster convergence for the learning algorithm. This performance is obtained with no or practically no impact on the quality of hidden state inference, as from labelling percentages around 20–30%; the prediction accuracy is also preserved above such percentage thresholds. Furthermore, faster EM convergence is also verified in a fully supervised scheme where part of the observations is ill-labelled. Model selection strategies can therefore rely on an approximate labelling function (provided by an expert or by a supervised algorithm learnt on a small subset of data for which the true labels are known), to explore larger grids of hyper-parameter values. In addition, complementary experimental studies have revealed the robustness of our model to labelling errors, particularly when large training datasets and moderate labelling error rates are considered. Finally, we showed the ability of our variant of the Viterbi algorithm to infer partially-labelled sequences.

We also conducted experiments on two realistic machine condition data (CMAPSS datasets FD001 and FD003) and considered both short/medium/long-term feature forecast accuracy and prediction accuracy for machine remaining useful life (RUL). Regarding feature prediction, we compared PHMC-MLAR with the two models it extends, MSAR and PHMC. The results show that at medium and long-term forecast horizons ($h = 10, 20, 30$) our model presents higher forecast accuracies than the MSAR model, whereas both models obtain comparable accuracies at short-term forecast horizon ($h = 5$). In comparison with the PHMC model, our model achieves much better performance (whatever the horizon $h = 5, 10, 20, 30$). These results show the relevance of including an autoregressive model within each regime (as we suggested in this work). Regarding RUL prediction, our proposal outperforms PHMC and five out of six recent state-of-the-art RUL prediction methods, including four artificial intelligence-based methods.

A natural extension of the PHMC-MLAR model consists in putting uncertainty on partial knowledge: for instance instead of states observed with no doubt, a subset of possible states with various occurrence probabilities can be considered at each time-step. On the
The results show that the method proposed in this work presents better performance than the first five comparison methods with a smaller RMSE and SF metrics. However, thanks to the extensive feature extraction procedure (CEEMD) used to build the input features of DLSTM, CEEMD DLSTM method presents better performance than our proposal. Note that our model is directly trained on noisy features recorded from sensors. Moreover, in contrast to CEEMD DLSTM method, our proposal provides an accurate (short, medium and long-term) prediction of features which can be used in practice to extract useful information about system operational states.

# 9 Conclusion

In this work, we have introduced the PHMC-MLAR model to analyze time series subject to switches in regimes. Our model is a generalization of the well-known Hidden Regime-Switching AutoRegressive (HRSAR) and Observed Regime-Switching AutoRegressive (ORSAR) models when regime-switching is modelled by a Markov Chain. Our model allows to handle the intermediate case where the state process is partially observed.

In the evaluation, we conducted experiments on simulated data and considered both inference performance and prediction accuracy. The results show that the partially observed states (when they represent a reasonable proportion) allow a better characterization of training data (reflected by greater log-likelihood), in comparison with the unsupervised case. An interesting characteristic of the PHMC-MLAR model is that the partially observed states allow faster convergence for the learning algorithm. This performance is obtained with no or practically no impact on the quality of hidden state inference, as from labelling percentages around 20–30%; the prediction accuracy is also preserved above such percentage thresholds. Furthermore, faster EM convergence is also verified in a fully supervised scheme where part of the observations is ill-labelled. Model selection strategies can therefore rely on an approximate labelling function (provided by an expert or by a supervised algorithm learnt on a small subset of data for which the true labels are known), to explore larger grids of hyper-parameter values. In addition, complementary experimental studies have revealed the robustness of our model to labelling errors, particularly when large training datasets and moderate labelling error rates are considered. Finally, we showed the ability of our variant of the Viterbi algorithm to infer partially-labelled sequences.

We also conducted experiments on two realistic machine condition data (CMAPSS datasets FD001 and FD003) and considered both short/medium/long-term feature forecast accuracy and prediction accuracy for machine remaining useful life (RUL). Regarding feature prediction, we compared PHMC-MLAR with the two models it extends, MSAR and PHMC. The results show that at medium and long-term forecast horizons ($h = 10, 20, 30$) our model presents higher forecast accuracies than the MSAR model, whereas both models obtain comparable accuracies at short-term forecast horizon ($h = 5$). In comparison with the PHMC model, our model achieves much better performance (whatever the horizon $h = 5, 10, 20, 30$). These results show the relevance of including an autoregressive model within each regime (as we suggested in this work). Regarding RUL prediction, our proposal outperforms PHMC and five out of six recent state-of-the-art RUL prediction methods, including four artificial intelligence-based methods.

A natural extension of the PHMC-MLAR model consists in putting uncertainty on partial knowledge: for instance instead of states observed with no doubt, a subset of possible states with various occurrence probabilities can be considered at each time-step. On the
other hand, it is more realistic to consider time-dependent state processes, especially when large time series are analyzed. These directions will be investigated in future work.

## A Appendix: backward-forward-backward algorithm

The Backward-forward-backward algorithm introduced by Scheffer and Wrobel (2001) for PHMC model learning has been adapted to the PHMC-MLAR framework. This algorithm makes it possible to compute the probabilities

$$ \xi_t(k, \ell) = P(S_{t-1} = k, S_t = \ell | X_{1-p}^T = x_{1-p}^T, \Sigma; \hat{\theta}), \quad \text{for } t = 2, \dots, T, \quad 1 \le k, \ell \le K \tag{37} $$

in $\mathcal{O}(TK^2)$ operations. The analytical development for the above quantity involves three additional probabilities:

$$ \xi_t(k, \ell) = \frac{\beta_t(\ell) P(S_t = \ell | S_{t-1} = k; \hat{\theta}) P(X_t = x_t | X_{t-p}^{T-1}, S_t = \ell; \hat{\theta}) \alpha_{t-1}(k) \tau_t(\ell)}{P(X_1^T = x_1^T | X_{1-p}^0, \Sigma; \hat{\theta}) \tau_{t-1}(k)} \tag{38} $$
$$ \times \mathbf{1}_{\{\ell \in \sigma_t, k \in \sigma_{t-1}\}}, $$

with

$$ \tau_t(s) = P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, \hat{\theta}), $$

$$ \alpha_t(s) = P(S_t = s, X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta}), $$

$$ \beta_t(s) = P(X_{t+1}^T = x_{t+1}^T | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta}). $$

The algorithm operates recursively in three steps: two backward steps chained through a forward step. The first backward step computes the set of probabilities $\tau_t(s)$ (Sect. A.1); the forward step computes probabilities $\alpha_t(s)$ (Sect. A.2); the second backward step computes probabilities $\beta_t(s)$. In Sect. A.4, we describe a scaling method that is necessary to prevent floating point underflow when running the algorithm, especially when large sequences are considered.

Note that the $\xi_t(k, \ell)$ quantities can be used to compute the $\gamma_t(\ell)$ probabilities

$$ \gamma_1(\ell) = \sum_{j=1}^{K} \xi_2(\ell, j) \quad \gamma_t(\ell) = \sum_{j=1}^{K} \xi_t(j, \ell), \text{ for } t = 2, \dots, T_t. $$

In practice $\gamma_t(\ell)$ have to be normalized by dividing them by the sum $\sum_{\ell=1}^K \gamma_t(\ell)$.

**Proof** First, in Eq. 39, the conditional probability is transformed into a joint probability. Then, in Eq. 40, we successively maginalize $X_{t+1}^T$, $X_t$, $S_t$ and $(S_{t-1}, X_1^{T-1})$. According to the conditional independence graph of the PHMC-MLAR model, the marginalization of $X_{t+1}^T$ gives $\beta_t(\ell)$, that of $X_t$ yields $P(X_t = x_t | S_t = \ell, X_{t-p}^{T-1}, \Sigma; \hat{\theta})$, that of $S_t$ gives $P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta})$ and that of $(S_{t-1}, X_1^{T-1})$ provides $\alpha_{t-1}(k)$. Finally, in Eqs. 42–44, the probability $P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta})$ is developed using Bayes’ rule. Note that in Eq. 44, the probability $P(S_t = \ell, \sigma_t | S_{t-1} = k, \sigma_{t-1}; \hat{\theta})$ is null for $\ell \notin \sigma_t$ and is not defined for $k \notin \sigma_{t-1}$.
other hand, it is more realistic to consider time-dependent state processes, especially when large time series are analyzed. These directions will be investigated in future work.

## A Appendix: backward-forward-backward algorithm

The Backward-forward-backward algorithm introduced by Scheffer and Wrobel (2001) for PHMC model learning has been adapted to the PHMC-MLAR framework. This algorithm makes it possible to compute the probabilities

$$ \xi_t(k, \ell) = P(S_{t-1} = k, S_t = \ell | X_{1-p}^T = x_{1-p}^T, \Sigma; \hat{\theta}), \quad \text{for } t = 2, \dots, T, \quad 1 \le k, \ell \le K \tag{37} $$

in $\mathcal{O}(TK^2)$ operations. The analytical development for the above quantity involves three additional probabilities:

$$ \xi_t(k, \ell) = \frac{\beta_t(\ell) P(S_t = \ell | S_{t-1} = k; \hat{\theta}) P(X_t = x_t | X_{t-p}^{T-1}, S_t = \ell; \hat{\theta}) \alpha_{t-1}(k) \tau_t(\ell)}{P(X_1^T = x_1^T | X_{1-p}^0, \Sigma; \hat{\theta}) \tau_{t-1}(k)} \tag{38} $$
$$ \times \mathbf{1}_{\{\ell \in \sigma_t, k \in \sigma_{t-1}\}}, $$

with

$$ \tau_t(s) = P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, \hat{\theta}), $$

$$ \alpha_t(s) = P(S_t = s, X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta}), $$

$$ \beta_t(s) = P(X_{t+1}^T = x_{t+1}^T | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta}). $$

The algorithm operates recursively in three steps: two backward steps chained through a forward step. The first backward step computes the set of probabilities $\tau_t(s)$ (Sect. A.1); the forward step computes probabilities $\alpha_t(s)$ (Sect. A.2); the second backward step computes probabilities $\beta_t(s)$. In Sect. A.4, we describe a scaling method that is necessary to prevent floating point underflow when running the algorithm, especially when large sequences are considered.

Note that the $\xi_t(k, \ell)$ quantities can be used to compute the $\gamma_t(\ell)$ probabilities

$$ \gamma_1(\ell) = \sum_{j=1}^{K} \xi_2(\ell, j) \quad \gamma_t(\ell) = \sum_{j=1}^{K} \xi_t(j, \ell), \text{ for } t = 2, \dots, T_i. $$

In practice $\gamma_t(\ell)$ have to be normalized by dividing them by the sum $\sum_{\ell=1}^K \gamma_t(\ell)$.

**Proof** First, in Eq. 39, the conditional probability is transformed into a joint probability. Then, in Eq. 40, we successively maginalize $X_{t+1}^T$, $X_t$, $S_t$ and $(S_{t-1}, X_1^{T-1})$. According to the conditional independence graph of the PHMC-MLAR model, the marginalization of $X_{t+1}^T$ gives $\beta_t(\ell)$, that of $X_t$ yields $P(X_t = x_t | S_t = \ell, X_{t-p}^{T-1}, \Sigma; \hat{\theta})$, that of $S_t$ gives $P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta})$ and that of $(S_{t-1}, X_1^{T-1})$ provides $\alpha_{t-1}(k)$. Finally, in Eqs. 42–44, the probability $P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta})$ is developed using Bayes’ rule. Note that in Eq. 44, the probability $P(S_t = \ell, \sigma_t | S_{t-1} = k, \sigma_{t-1}; \hat{\theta})$ is null for $\ell \notin \sigma_t$ and is not defined for $k \notin \sigma_{t-1}$.
$$ \begin{aligned} \xi_t(k, \ell) &= P(S_{t-1} = k, S_t = \ell | X_{1-p}^T = x_{1-p}^T, \Sigma; \hat{\theta}) \\ &= \frac{P(S_{t-1} = k, S_t = \ell, X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})}{P(X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})} \end{aligned} \quad (39) $$

$$ \begin{aligned} &= P(X_{t+1}^T = x_{t+1}^T | S_{t-1} = k, S_t = \ell, X_1^t, X_{1-p}^0, \Sigma; \hat{\theta}) \\ &\quad \times P(X_t = x_t | S_{t-1} = k, S_t = \ell, X_1^{t-1}, X_{1-p}^0, \Sigma; \hat{\theta}) \\ &\quad \times P(S_t = \ell | S_{t-1} = k, X_1^{t-1}, X_{1-p}^0, \Sigma; \hat{\theta}) \frac{P(S_{t-1} = k, X_1^{t-1} = x_1^{t-1} | X_{1-p}^0, \Sigma; \hat{\theta})}{P(X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})} \end{aligned} \quad (40) $$

$$ = \frac{\beta_t(\ell) P(X_t = x_t | S_t = \ell, X_{t-p}^{t-1}, \Sigma; \hat{\theta}) \alpha_{t-1}(k)}{P(X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})} \times P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta}) \quad (41) $$

with

$$ P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta}) = \frac{P(S_t = \ell, S_{t-1} = k, \sigma_1^T; \hat{\theta})}{P(S_{t-1} = k, \sigma_1^T; \hat{\theta})} \quad (42) $$

$$ = \frac{P(\sigma_{t+1}^T | S_t = \ell, S_{t-1} = k, \sigma_1^t; \hat{\theta}) P(S_t = \ell, S_{t-1} = k, \sigma_1^t; \hat{\theta})}{P(\sigma_t^T | S_{t-1} = k, \sigma_1^{t-1}; \hat{\theta}) P(S_{t-1} = k, \sigma_1^{t-1}; \hat{\theta})} \quad (43) $$

$$ = \frac{\tau_t(\ell)}{\tau_{t-1}(k)} \times P(S_t = \ell, \sigma_t | S_{t-1} = k, \sigma_{t-1}; \hat{\theta}) \quad (44) $$

$$ = \begin{cases} \frac{\tau_t(\ell)}{\tau_{t-1}(k)} \times P(S_t = \ell | S_{t-1} = k; \hat{\theta}) & \text{if } k \in \sigma_{t-1}, \ell \in \sigma_t \\ 0 & \text{otherwise.} \end{cases} $$

□

## A.1 First backward step

The first backward step computes probabilities $\tau_t(s)$, the probabilities of the remaining possible states given that state $s \in \{1, ..., K\}$ is observed at time-step $t \in \{1, ..., T\}$: $\tau_t(s) = P(\sigma_{t+1}, ..., \sigma_T | S_t = s, \hat{\theta}) = P(S_{t+1} \in \sigma_{t+1}, ..., S_T \in \sigma_T | S_t = s, \hat{\theta})$. This set of probabilities is computed recursively as follows:
$$ 
\begin{aligned}
\xi_t(k, \ell) &= P(S_{t-1} = k, S_t = \ell | X_{1-p}^T = x_{1-p}^T, \Sigma; \hat{\theta}) \\
&= \frac{P(S_{t-1} = k, S_t = \ell, X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})}{P(X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})} \quad (39)
\end{aligned}
$$

$$ 
\begin{aligned}
&= P(X_{t+1}^T = x_{t+1}^T | S_{t-1} = k, S_t = \ell, X_1^t, X_{1-p}^0, \Sigma; \hat{\theta}) \\
&\quad \times P(X_t = x_t | S_{t-1} = k, S_t = \ell, X_1^{t-1}, X_{1-p}^0, \Sigma; \hat{\theta}) \\
&\quad \times P(S_t = \ell | S_{t-1} = k, X_1^{t-1}, X_{1-p}^0, \Sigma; \hat{\theta}) \frac{P(S_{t-1} = k, X_1^{t-1} = x_1^{t-1} | X_{1-p}^0, \Sigma; \hat{\theta})}{P(X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})} \quad (40)
\end{aligned}
$$

$$ = \frac{\beta_t(\ell) P(X_t = x_t | S_t = \ell, X_{t-p}^{t-1}, \Sigma; \hat{\theta}) \alpha_{t-1}(k)}{P(X_1^T | X_{1-p}^0, \Sigma; \hat{\theta})} \times P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta}) \quad (41) $$

with

$$ P(S_t = \ell | S_{t-1} = k, \Sigma; \hat{\theta}) = \frac{P(S_t = \ell, S_{t-1} = k, \sigma_1^T; \hat{\theta})}{P(S_{t-1} = k, \sigma_1^T; \hat{\theta})} \quad (42) $$

$$ = \frac{P(\sigma_{t+1}^T | S_t = \ell, S_{t-1} = k, \sigma_1^t; \hat{\theta}) P(S_t = \ell, S_{t-1} = k, \sigma_1^t; \hat{\theta})}{P(\sigma_t^T | S_{t-1} = k, \sigma_1^{t-1}; \hat{\theta}) P(S_{t-1} = k, \sigma_1^{t-1}; \hat{\theta})} \quad (43) $$

$$ = \frac{\tau_t(\ell)}{\tau_{t-1}(k)} \times P(S_t = \ell, \sigma_t | S_{t-1} = k, \sigma_{t-1}; \hat{\theta}) \quad (44) $$

$$ = \begin{cases} \frac{\tau_t(\ell)}{\tau_{t-1}(k)} \times P(S_t = \ell | S_{t-1} = k; \hat{\theta}) & \text{if } k \in \sigma_{t-1}, \ell \in \sigma_t \\ 0 & \text{otherwise.} \end{cases} $$

□

## A.1 First backward step

The first backward step computes probabilities $\tau_t(s)$, the probabilities of the remaining possible states given that state $s \in \{1, ..., K\}$ is observed at time-step $t \in \{1, ..., T\}$: $\tau_t(s) = P(\sigma_{t+1}, ..., \sigma_T | S_t = s, \hat{\theta}) = P(S_{t+1} \in \sigma_{t+1}, ..., S_T \in \sigma_T | S_t = s, \hat{\theta})$. This set of probabilities is computed recursively as follows:
$$ \left\{ \begin{array}{l} \tau_T(s) := 1 \\ \tau_t(s) = \sum_{i \in \sigma_{t+1}} \tau_{t+1}(i) P(S_{t+1} = i | S_t = s; \hat{\theta}). \end{array} \right. \qquad (45) $$

**Proof** Base case: $t = T - 1$

By applying the definition of $\tau_{T-1}$, we obtain:

$$ \tau_{T-1}(s) = P(\sigma_T | S_{T-1} = s, \hat{\theta}) = P(S_T \in \sigma_T | S_{T-1} = s; \hat{\theta}) \qquad (46) $$

$$ = \sum_{i \in \sigma_T} \tau_T(i) P(S_T = i | S_{T-1} = s; \hat{\theta}). \qquad (47) $$

**Recursive case:** $t = T - 2, \dots, 1$

We first use the law of total probabilities (Eq. 48), followed by Bayes' rule (Eq. 49). Note that in Eq. 49, the probability $P(\sigma_{t+1}, \dots, \sigma_T | S_{t+1} = i, S_t = s, \hat{\theta})$ is null for $i \notin \sigma_{t+1}$ (since $\sigma_{t+1}$ is the set of possible states at time-step $t+1$); otherwise it equals $P(\sigma_{t+2}, \dots, \sigma_T | S_{t+1} = i, \hat{\theta}) = \tau_{t+1}(i)$ (Eq. 50). Thus, we obtain the recursive formula presented in Eq. 45.

$$ \begin{aligned} \tau_t(s) &= P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, \hat{\theta}) \\ &= \sum_{i=1}^{K} P(\sigma_{t+1}, \dots, \sigma_T, S_{t+1} = i | S_t = s, \hat{\theta}) \end{aligned} \qquad (48) $$

$$ = \sum_{i=1}^{K} P(\sigma_{t+1}, \dots, \sigma_T | S_{t+1} = i, S_t = s, \hat{\theta}) P(S_{t+1} = i | S_t = s, \hat{\theta}) \qquad (49) $$

$$ = \sum_{i \in \sigma_{t+1}} P(\sigma_{t+2}, \dots, \sigma_T | S_{t+1} = i, \hat{\theta}) P(S_{t+1} = i | S_t = s, \hat{\theta}) \qquad (50) $$

$$ = \sum_{i \in \sigma_{t+1}} \tau_{t+1}(i) P(S_{t+1} = i | S_t = s, \hat{\theta}). \qquad (51) $$

□

## A.2 Forward step

This step allows to compute the probabilities of being in regime $s$ at time-step $t$ while observing sequence $x_1, \dots, x_t$. These probabilities, denoted by $\alpha_t(s)$, are defined as $\alpha_t(s) = P(S_t = s, X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta})$ for $1 \le t \le T$, $1 \le s \le K$. They are computed as follows:
$$ \left\{ \begin{array}{l} \tau_T(s) := 1 \\ \tau_t(s) = \sum_{i \in \sigma_{t+1}} \tau_{t+1}(i) P(S_{t+1} = i | S_t = s; \hat{\theta}). \end{array} \right. \qquad (45) $$

**Proof** Base case: $t = T - 1$

By applying the definition of $\tau_{T-1}$, we obtain:

$$ \tau_{T-1}(s) = P(\sigma_T | S_{T-1} = s, \hat{\theta}) = P(S_T \in \sigma_T | S_{T-1} = s; \hat{\theta}) \quad (46) $$

$$ = \sum_{i \in \sigma_T} \tau_T(i) P(S_T = i | S_{T-1} = s; \hat{\theta}). \qquad (47) $$

**Recursive case:** $t = T - 2, \dots, 1$

We first use the law of total probabilities (Eq. 48), followed by Bayes' rule (Eq. 49). Note that in Eq. 49, the probability $P(\sigma_{t+1}, \dots, \sigma_T | S_{t+1} = i, S_t = s, \hat{\theta})$ is null for $i \notin \sigma_{t+1}$ (since $\sigma_{t+1}$ is the set of possible states at time-step $t+1$); otherwise it equals $P(\sigma_{t+2}, \dots, \sigma_T | S_{t+1} = i, \hat{\theta}) = \tau_{t+1}(i)$ (Eq. 50). Thus, we obtain the recursive formula presented in Eq. 45.

$$ \begin{aligned} \tau_t(s) &= P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, \hat{\theta}) \\ &= \sum_{i=1}^{K} P(\sigma_{t+1}, \dots, \sigma_T, S_{t+1} = i | S_t = s, \hat{\theta}) \end{aligned} \qquad (48) $$

$$ = \sum_{i=1}^{K} P(\sigma_{t+1}, \dots, \sigma_T | S_{t+1} = i, S_t = s, \hat{\theta}) P(S_{t+1} = i | S_t = s, \hat{\theta}) \quad (49) $$

$$ = \sum_{i \in \sigma_{t+1}} P(\sigma_{t+2}, \dots, \sigma_T | S_{t+1} = i, \hat{\theta}) P(S_{t+1} = i | S_t = s, \hat{\theta}) \quad (50) $$

$$ = \sum_{i \in \sigma_{t+1}} \tau_{t+1}(i) P(S_{t+1} = i | S_t = s, \hat{\theta}). \qquad (51) $$

□

## A.2 Forward step

This step allows to compute the probabilities of being in regime $s$ at time-step $t$ while observing sequence $x_1, \dots, x_t$. These probabilities, denoted by $\alpha_t(s)$, are defined as $\alpha_t(s) = P(S_t = s, X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta})$ for $1 \le t \le T$, $1 \le s \le K$. They are computed as follows:
$$
\left\{
\begin{aligned}
\alpha_1(s) &= P(X_1 = x_1 | X_{1-p}^0, S_1 = s; \hat{\theta}) P(S_1 = s; \hat{\theta}) \frac{\tau_1(s)}{\sum_{i \in \sigma_1} \tau_1(i) P(S_1 = i; \hat{\theta})} \\
\alpha_t(s) &= P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta}) \sum_{i \in \sigma_{t-1}} \alpha_{t-1}(i) P(S_t = s | S_{t-1} = i; \hat{\theta}) \frac{\tau_t(s)}{\tau_{t-1}(i)} \times \mathbf{1}_{\{s \in \sigma_t\}}
\end{aligned}
\right.
\qquad (52)
$$

To note, the likelihood of sequence $x_1^T$ can be easily computed by integrating out $S_t$ in $\alpha_T$:

$$ P(X_1^T = x_1^T | X_{1-p}^0, \Sigma; \hat{\theta}) = \sum_{s=1}^{K} \alpha_T(s). \qquad (53) $$

The likelihood of $N$ independent sequences is therefore calculated by multiplying the individual likelihoods across the sequences.

**Proof** Base case: $t=1$

In Eq. 54, using the conditional independence graph of the PHMC-MLAR model, we transform the joint probability into two conditional probabilities, $P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta})$ and $P(S_1 = s | \Sigma; \hat{\theta})$. Then, in Eq. 55, Bayes' rule is applied to the latter conditional probability. It can be easily shown that $P(\Sigma; \hat{\theta}) = \sum_{i \in \sigma_1} \tau_1(i) P(S_1 = i \hat{\theta})$. Thus we obtain Eq. 56.

$$
\begin{aligned}
\alpha_1(s) &= P(S_1 = s, X_1 = x_1 | X_{1-p}^0, \Sigma; \hat{\theta}) \\
&= P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta}) \times P(S_1 = s | \Sigma; \hat{\theta})
\end{aligned}
\qquad (54)
$$

$$ = P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta}) \times \frac{P(\sigma_1, \dots, \sigma_T | S_1 = s; \hat{\theta}) P(S_1 = s; \hat{\theta})}{P(\Sigma; \hat{\theta})} \quad (55) $$

$$ = P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta}) P(S_1 = s; \hat{\theta}) \frac{\tau_1(s)}{\sum_{i \in \sigma_1} \tau_1(i) P(S_1 = i \hat{\theta})}. \quad (56) $$

**Recursive case: $t=2, \dots, T$**

As previously, the joint probability is split into two conditional probabilities (Eq. 57). We use the law of total probabilities to introduce $S_{t-1}$ in Eq. 58. From Eqs. 58 to 59, Bayes' rule is applied on the terms within the sum. Then, in Eq. 60, recursive terms $\alpha_{t-1}$ weighted by probabilities $P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta})$ appear within the sum. Finally, probabilities $P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta})$ are computed through the calculations presented in Eqs. 61–64. Thus, by substituting Eq. 64 in Eq. 60, we obtain the recursive case (Eq. 52).

$$
\begin{aligned}
\alpha_t(s) &= P(S_t = s, X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta}) \\
&= P(X_t = x_t | X_1^{t-1}, S_t = s, X_{1-p}^0, \Sigma; \hat{\theta}) \times P(X_1^{t-1} = x_1^{t-1}, S_t = s | X_{1-p}^0, \Sigma; \hat{\theta}) \quad (57) \\
&= P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta})
\end{aligned}
$$

$$
\begin{aligned}
& \sum_{i=1}^{K} P(X_1^{i-1} = x_1^{i-1}, S_t = s, S_{t-1} = i | X_{1-p}^0, \Sigma; \hat{\theta}) \quad (58) \\
& = P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta})
\end{aligned}
$$
$$
\left\{
\begin{aligned}
\alpha_1(s) &= P(X_1 = x_1 | X_{1-p}^0, S_1 = s; \hat{\theta}) P(S_1 = s; \hat{\theta}) \frac{\tau_1(s)}{\sum_{i \in \sigma_1} \tau_1(i) P(S_1 = i; \hat{\theta})} \\
\alpha_t(s) &= P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta}) \sum_{i \in \sigma_{t-1}} \alpha_{t-1}(i) P(S_t = s | S_{t-1} = i; \hat{\theta}) \frac{\tau_t(s)}{\tau_{t-1}(i)} \times \mathbf{1}_{\{s \in \sigma_t\}}
\end{aligned}
\right.
\tag{52}
$$

To note, the likelihood of sequence $x_1^T$ can be easily computed by integrating out $S_t$ in $\alpha_T$:

$$ P(X_1^T = x_1^T | X_{1-p}^0, \Sigma; \hat{\theta}) = \sum_{s=1}^{K} \alpha_T(s). \tag{53} $$

The likelihood of $N$ independent sequences is therefore calculated by multiplying the individual likelihoods across the sequences.

**Proof** Base case: $t=1$

In Eq. 54, using the conditional independence graph of the PHMC-MLAR model, we transform the joint probability into two conditional probabilities, $P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta})$ and $P(S_1 = s | \Sigma; \hat{\theta})$. Then, in Eq. 55, Bayes' rule is applied to the latter conditional probability. It can be easily shown that $P(\Sigma; \hat{\theta}) = \sum_{i \in \sigma_1} \tau_1(i) P(S_1 = i; \hat{\theta})$. Thus we obtain Eq. 56.

$$
\begin{aligned}
\alpha_1(s) &= P(S_1 = s, X_1 = x_1 | X_{1-p}^0, \Sigma; \hat{\theta}) \\
&= P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta}) \times P(S_1 = s | \Sigma; \hat{\theta})
\end{aligned}
\tag{54}
$$

$$ = P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta}) \times \frac{P(\sigma_1, \dots, \sigma_T | S_1 = s; \hat{\theta}) P(S_1 = s; \hat{\theta})}{P(\Sigma; \hat{\theta})} \tag{55} $$

$$ = P(X_1 = x_1 | S_1 = s, X_{1-p}^0; \hat{\theta}) P(S_1 = s; \hat{\theta}) \frac{\tau_1(s)}{\sum_{i \in \sigma_1} \tau_1(i) P(S_1 = i; \hat{\theta})} \tag{56} $$

**Recursive case:** $t=2, \dots, T$

As previously, the joint probability is split into two conditional probabilities (Eq. 57). We use the law of total probabilities to introduce $S_{t-1}$ in Eq. 58. From Eqs. 58 to 59, Bayes' rule is applied on the terms within the sum. Then, in Eq. 60, recursive terms $\alpha_{t-1}$ weighted by probabilities $P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta})$ appear within the sum. Finally, probabilities $P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta})$ are computed through the calculations presented in Eqs. 61–64. Thus, by substituting Eq. 64 in Eq. 60, we obtain the recursive case (Eq. 52).

$$
\begin{aligned}
\alpha_t(s) &= P(S_t = s, X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta}) \\
&= P(X_t = x_t | X_1^{t-1}, S_t = s, X_{1-p}^0, \Sigma; \hat{\theta}) \times P(X_1^{t-1} = x_1^{t-1}, S_t = s | X_{1-p}^0, \Sigma; \hat{\theta}) \tag{57} \\
&= P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta})
\end{aligned}
$$

$$
\begin{aligned}
\sum_{i=1}^{K} P(X_1^{i-1} = x_1^{i-1}, S_t = s, S_{t-1} = i | X_{1-p}^0, \Sigma; \hat{\theta}) \tag{58} \\
= P(X_t = x_t | X_{t-p}^{i-1}, S_t = s; \hat{\theta})
\end{aligned}
$$
$$ \sum_{i=1}^{K} P(X_{1}^{t-1} = x_{1}^{t-1}, S_{t-1} = i | X_{1-p}^{0}, \Sigma; \hat{\theta}) P(S_{t} = s | S_{t-1} = i, \Sigma; \hat{\theta}) \quad (59) $$

$$ = P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta}) \sum_{i=1}^{K} \alpha_{t-1}(i) P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta}) \quad (60) $$

where

$$ P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta}) = \frac{P(S_t = s, S_{t-1} = i, \Sigma; \hat{\theta})}{P(S_{t-1} = i, \Sigma; \hat{\theta})} \quad (61) $$

$$ = \frac{P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, S_{t-1} = i, \sigma_1^t; \hat{\theta})}{P(\sigma_t, \dots, \sigma_T | S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta}) P(S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta})} \times P(S_t = s, \sigma_t | S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta}) P(S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta}) \quad (62) $$

$$ = \frac{\tau_t(s)}{\tau_{t-1}(i)} \times P(S_t = s, \sigma_t | S_{t-1} = i, \sigma_{t-1}; \hat{\theta}) \quad (63) $$

$$ = \frac{\tau_t(s)}{\tau_{t-1}(i)} \times \begin{cases} P(S_t = s | S_{t-1} = i; \hat{\theta}) & \text{if } i \in \sigma_{t-1}, s \in \sigma_t \\ 0 & \text{otherwise.} \end{cases} \quad (64) $$

□

### A.3 Second backward step

In this second backward step, quantities $\beta_t(s) = P(X_{t+1}^T = x_{t+1}^T | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta})$ are computed. $\beta_t(s)$ denotes the probability to observe sequence $x_{t+1}, \dots, x_T$ given that state $s$ has been observed at time-step $t$. These probabilities are recursively computed as follows:

$$ \left\{ \begin{aligned} \beta_T(s) &:= 1 \\ \beta_t(s) &= \sum_{i \in \sigma_{t+1}} \beta_{t+1}(i) P(S_{t+1} = i | S_t = s; \hat{\theta}) \frac{\tau_{t+1}(i)}{\tau_t(s)} \quad (65) \\ P(X_{t+1} = x_{t+1} | X_{t+1-p}^T, S_{t+1} = i; \hat{\theta}) \times \mathbf{1}_{\{s \in \sigma_t\}} & \end{aligned} \right. $$

**Proof** Base case: $t = T - 1$

Equation 66 is obtained by applying the law of total probabilities. In Eq. 67, Bayes’ rule is applied to $P(S_T = i | S_{T-1} = s, \Sigma; \hat{\theta})$ and a quotient of probabilities appears. Then, the numerator and denominator of this quotient are transformed into products of conditional probabilities (Eq. 68). In Eq. 69, we introduce backward propagation terms
$$ \sum_{i=1}^{K} P(X_{1}^{t-1} = x_{1}^{t-1}, S_{t-1} = i | X_{1-p}^{0}, \Sigma; \hat{\theta}) P(S_{t} = s | S_{t-1} = i, \Sigma; \hat{\theta}) \quad (59) $$

$$ = P(X_t = x_t | X_{t-p}^{t-1}, S_t = s; \hat{\theta}) \sum_{i=1}^{K} \alpha_{t-1}(i) P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta}) \quad (60) $$

where

$$ P(S_t = s | S_{t-1} = i, \Sigma; \hat{\theta}) = \frac{P(S_t = s, S_{t-1} = i, \Sigma; \hat{\theta})}{P(S_{t-1} = i, \Sigma; \hat{\theta})} \quad (61) $$

$$ = \frac{P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, S_{t-1} = i, \sigma_1^t; \hat{\theta})}{P(\sigma_t, \dots, \sigma_T | S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta}) P(S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta})} \times P(S_t = s, \sigma_t | S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta}) P(S_{t-1} = i, \sigma_1^{t-1}; \hat{\theta}) \quad (62) $$

$$ = \frac{\tau_t(s)}{\tau_{t-1}(i)} \times P(S_t = s, \sigma_t | S_{t-1} = i, \sigma_{t-1}; \hat{\theta}) \quad (63) $$

$$ = \frac{\tau_t(s)}{\tau_{t-1}(i)} \times \begin{cases} P(S_t = s | S_{t-1} = i; \hat{\theta}) & \text{if } i \in \sigma_{t-1}, s \in \sigma_t \\ 0 & \text{otherwise.} \end{cases} \quad (64) $$

□

### A.3 Second backward step

In this second backward step, quantities $\beta_t(s) = P(X_{t+1}^T = x_{t+1}^T | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta})$ are computed. $\beta_t(s)$ denotes the probability to observe sequence $x_{t+1}, \dots, x_T$ given that state $s$ has been observed at time-step $t$. These probabilities are recursively computed as follows:

$$ \left\{ \begin{aligned} \beta_T(s) &:= 1 \\ \beta_t(s) &= \sum_{i \in \sigma_{t+1}} \beta_{t+1}(i) P(S_{t+1} = i | S_t = s; \hat{\theta}) \frac{\tau_{t+1}(i)}{\tau_t(s)} \\ P(X_{t+1} = x_{t+1} | X_{t+1-p}^T, S_{t+1} = i; \hat{\theta}) \times \mathbf{1}_{\{s \in \sigma_t\}} & \end{aligned} \right. \quad (65) $$

**Proof** Base case: $t = T - 1$

Equation 66 is obtained by applying the law of total probabilities. In Eq. 67, Bayes’ rule is applied to $P(S_T = i | S_{T-1} = s, \Sigma; \hat{\theta})$ and a quotient of probabilities appears. Then, the numerator and denominator of this quotient are transformed into products of conditional probabilities (Eq. 68). In Eq. 69, we introduce backward propagation terms
$\beta_T(i)$ and $\tau_T(i)$, which each equal one (by definition); thanks to Markov property, probability $P(S_T = i, \sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})$ is equal to $P(S_T = i | S_{T-1} = s; \hat{\theta})$ if $i \in \sigma_T$ and $s \in \sigma_{T-1}$, and this probability is null if $i \notin \sigma_T$ and is undefined if $s \notin \sigma_{T-1}$ (hence the indicator function $1_{\{s \in \sigma_{T-1}, i \in \sigma_T\}}$). Besides, in Eq. 68, a common term appears at numerator and denominator, which entails a simplification. Finally, probability $P(\sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})$ appearing at denominator equals $\tau_{T-1}(s)$ thanks to Markov property.

$$ \beta_{T-1}(s) = P(X_T = x_T | X_{T-p}^{T-1}, S_{T-1} = s, \Sigma; \hat{\theta}) \\ = \sum_{i=1}^{K} P(X_T = x_T | X_{T-p}^{T-1}, S_T = i, \Sigma; \hat{\theta}) P(S_T = i | S_{T-1} = s, \Sigma; \hat{\theta}) \quad (66) $$

$$ = \sum_{i=1}^{K} P(X_T = x_T | X_{T-p}^{T-1}, S_T = i; \hat{\theta}) \frac{P(S_T = i, S_{T-1} = s, \Sigma; \hat{\theta})}{P(S_{T-1} = s, \Sigma; \hat{\theta})} \quad (67) $$

$$ = \sum_{i=1}^{K} P(X_T = x_T | X_{T-p}^{T-1}, S_T = i; \hat{\theta}) \times \frac{P(S_T = i, \sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta}) P(S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})}{P(\sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta}) P(S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})} \quad (68) $$

$$ = \sum_{i=1}^{K} \beta_T(i) P(X_T = x_T | X_{T-p}^{T-1}, S_T = i; \hat{\theta}) \frac{\tau_T(i)}{\tau_{T-1}(s)} \times P(S_T = i | S_{T-1} = s; \hat{\theta}) \times 1_{\{s \in \sigma_{T-1}, i \in \sigma_T\}} \quad (69) $$

Recursive case: $t = T - 2, \dots, 1$

The application of the law of total probabilities yields Eq. 70. In Eq. 71, $X_{t+2}^T$ then $X_{t+1}$ are marginalized, which allows to make appear the recursive term $\beta_{t+1}$ together with the conditional probability of $X_{t+1}$ given $S_{t+1}$ and past values in Eq. 72. As in the base case, probability $P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta})$ is computed using Bayes’ rule (Eqs. 73–76).

$$ \beta_t(s) = P(X_{t+1}^T = x_{t+1}^T | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta}) \\ = \sum_{i=1}^{K} P(X_{t+1} = x_{t+1}, X_{t+2}^T = x_{t+2}^T, S_{t+1} = i | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta}) \quad (70) $$

$$ = \sum_{i=1}^{K} P(X_{t+2}^T = x_{t+2}^T | X_{t+1-p}^{t+1}, S_{t+1} = i, \Sigma; \hat{\theta}) \times P(X_{t+1} = x_{t+1} | X_{t+1-p}^t, S_{t+1} = i, \Sigma; \hat{\theta}) P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta}) \quad (71) $$

$$ = \sum_{i=1}^{K} \beta_{t+1}(i) P(X_{t+1} = x_{t+1} | X_{t+1-p}^t, S_{t+1} = i, \Sigma; \hat{\theta}) P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta}) \quad (72) $$
$\beta_T(i)$ and $\tau_T(i)$, which each equal one (by definition); thanks to Markov property, probability $P(S_T = i, \sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})$ is equal to $P(S_T = i | S_{T-1} = s; \hat{\theta})$ if $i \in \sigma_T$ and $s \in \sigma_{T-1}$, and this probability is null if $i \notin \sigma_T$ and is undefined if $s \notin \sigma_{T-1}$ (hence the indicator function $1_{\{s \in \sigma_{T-1}, i \in \sigma_T\}}$). Besides, in Eq. 68, a common term appears at numerator and denominator, which entails a simplification. Finally, probability $P(\sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})$ appearing at denominator equals $\tau_{T-1}(s)$ thanks to Markov property.

$$ \beta_{T-1}(s) = P(X_T = x_T | X_{T-p}^{T-1}, S_{T-1} = s, \Sigma; \hat{\theta}) = \sum_{i=1}^{K} P(X_T = x_T | X_{T-p}^{T-1}, S_T = i, \Sigma; \hat{\theta}) P(S_T = i | S_{T-1} = s, \Sigma; \hat{\theta}) \quad (66) $$

$$ = \sum_{i=1}^{K} P(X_T = x_T | X_{T-p}^{T-1}, S_T = i; \hat{\theta}) \frac{P(S_T = i, S_{T-1} = s, \Sigma; \hat{\theta})}{P(S_{T-1} = s, \Sigma; \hat{\theta})} \quad (67) $$

$$ = \sum_{i=1}^{K} P(X_T = x_T | X_{T-p}^{T-1}, S_T = i; \hat{\theta}) \times \frac{P(S_T = i, \sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta}) P(S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})}{P(\sigma_T | S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta}) P(S_{T-1} = s, \sigma_1, \dots, \sigma_{T-1}; \hat{\theta})} \quad (68) $$

$$ = \sum_{i=1}^{K} \beta_T(i) P(X_T = x_T | X_{T-p}^{T-1}, S_T = i; \hat{\theta}) \frac{\tau_T(i)}{\tau_{T-1}(s)} \times P(S_T = i | S_{T-1} = s; \hat{\theta}) \times 1_{\{s \in \sigma_{T-1}, i \in \sigma_T\}} \quad (69) $$

Recursive case: $t = T - 2, \dots, 1$

The application of the law of total probabilities yields Eq. 70. In Eq. 71, $X_{t+2}^T$ then $X_{t+1}$ are marginalized, which allows to make appear the recursive term $\beta_{t+1}$ together with the conditional probability of $X_{t+1}$ given $S_{t+1}$ and past values in Eq. 72. As in the base case, probability $P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta})$ is computed using Bayes’ rule (Eqs. 73–76).

$$ \beta_t(s) = P(X_{t+1}^T = x_{t+1}^T | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta}) = \sum_{i=1}^{K} P(X_{t+1} = x_{t+1}, X_{t+2}^T = x_{t+2}^T, S_{t+1} = i | X_{t+1-p}^T, S_t = s, \Sigma; \hat{\theta}) \quad (70) $$

$$ = \sum_{i=1}^{K} P(X_{t+2}^T = x_{t+2}^T | X_{t+1-p}^{t+1}, S_{t+1} = i, \Sigma; \hat{\theta}) \times P(X_{t+1} = x_{t+1} | X_{t+1-p}^t, S_{t+1} = i, \Sigma; \hat{\theta}) P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta}) \quad (71) $$

$$ = \sum_{i=1}^{K} \beta_{t+1}(i) P(X_{t+1} = x_{t+1} | X_{t+1-p}^t, S_{t+1} = i, \Sigma; \hat{\theta}) P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta}) \quad (72) $$
where

$$P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta}) = \frac{P(S_{t+1} = i, S_t = s, \Sigma; \hat{\theta})}{P(S_t = s, \Sigma; \hat{\theta})} \quad (73)$$

$$= (\sigma_{t+2}, \dots, \sigma_T | S_{t+1} = i, S_t = s, \sigma_1, \dots, \sigma_{t+1}; \hat{\theta}) \\ \times \frac{P(S_{t+1} = i, \sigma_{t+1} | S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta}) P(S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta})}{P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta}) P(S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta})} \quad (74)$$

$$= \frac{\tau_{t+1}(i)}{\tau_t(s)} \times P(S_{t+1} = i, \sigma_{t+1} | S_t = s, \sigma_t; \hat{\theta}) \quad (75)$$

$$= \frac{\tau_{t+1}(i)}{\tau_t(s)} \times P(S_{t+1} = i | S_t = s; \hat{\theta}) \times \mathbf{1}_{\{s \in \sigma_t, i \in \sigma_{t+1}\}}. \quad (76)$$

□

## A.4 Scaling of backward-forward-backward algorithm

For large sequences, i.e. large value of T, the quantities $\tau_t(s)$, $\alpha_t(s)$ and $\beta_t(s)$ tend to zero as products of probabilities. Thus, the computations will require beyond the precision range of machine and PHMC-MLAR parameter estimate will be inaccurate. Generally, this problem is solved by normalizing $\tau_t(s)$, $\alpha_t(s)$ and $\beta_t(s)$ by a term of same order of magnitude (Florez-Larrahondo, 2020; Koenig & Simmons, 1996). Thus, we propose the following normalization:

$$\tilde{\tau}_t(s) = \frac{\tau_t(s)}{P(\sigma_t, \dots, \sigma_T | \sigma_{t-1}; \hat{\theta})}, \quad (77)$$

$$\tilde{\alpha}_t(s) = \frac{\alpha_t(s)}{P(X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta})}, \quad (78)$$

$$\tilde{\beta}_t(s) = \frac{\beta_t(s)}{P(X_t^T = x_t^T | X_{1-p}^{T-1}, \Sigma; \hat{\theta})}. \quad (79)$$

As previously, $\tilde{\tau}_t(s)$, $\tilde{\alpha}_t(s)$ and $\tilde{\beta}_t(s)$ can be computed recursively. The recursive formula for these quantities can be deduced from those of $\tau_t(s)$ (Eq. 45), $\alpha_t(s)$ (Eq. 52) and $\beta_t(s)$
where

$$P(S_{t+1} = i | S_t = s, \Sigma; \hat{\theta}) = \frac{P(S_{t+1} = i, S_t = s, \Sigma; \hat{\theta})}{P(S_t = s, \Sigma; \hat{\theta})} \quad (73)$$

$$= (\sigma_{t+2}, \dots, \sigma_T | S_{t+1} = i, S_t = s, \sigma_1, \dots, \sigma_{t+1}; \hat{\theta}) \\ \times \frac{P(S_{t+1} = i, \sigma_{t+1} | S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta}) P(S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta})}{P(\sigma_{t+1}, \dots, \sigma_T | S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta}) P(S_t = s, \sigma_1, \dots, \sigma_T; \hat{\theta})} \quad (74)$$

$$= \frac{\tau_{t+1}(i)}{\tau_t(s)} \times P(S_{t+1} = i, \sigma_{t+1} | S_t = s, \sigma_t; \hat{\theta}) \quad (75)$$

$$= \frac{\tau_{t+1}(i)}{\tau_t(s)} \times P(S_{t+1} = i | S_t = s; \hat{\theta}) \times \mathbf{1}_{\{s \in \sigma_t, i \in \sigma_{t+1}\}}. \quad (76)$$

<figure><img src="image_26.png" /></figure>

## A.4 Scaling of backward-forward-backward algorithm

For large sequences, i.e. large value of T, the quantities $\tau_t(s)$, $\alpha_t(s)$ and $\beta_t(s)$ tend to zero as products of probabilities. Thus, the computations will require beyond the precision range of machine and PHMC-MLAR parameter estimate will be inaccurate. Generally, this problem is solved by normalizing $\tau_t(s)$, $\alpha_t(s)$ and $\beta_t(s)$ by a term of same order of magnitude (Florez-Larrahondo, 2020; Koenig & Simmons, 1996). Thus, we propose the following normalization:

$$\tilde{\tau}_t(s) = \frac{\tau_t(s)}{P(\sigma_t, \dots, \sigma_T | \sigma_{t-1}; \hat{\theta})}, \quad (77)$$

$$\tilde{\alpha}_t(s) = \frac{\alpha_t(s)}{P(X_1^t = x_1^t | X_{1-p}^0, \Sigma; \hat{\theta})}, \quad (78)$$

$$\tilde{\beta}_t(s) = \frac{\beta_t(s)}{P(X_t^T = x_t^T | X_{1-p}^{T-1}, \Sigma; \hat{\theta})}. \quad (79)$$

As previously, $\tilde{\tau}_t(s)$, $\tilde{\alpha}_t(s)$ and $\tilde{\beta}_t(s)$ can be computed recursively. The recursive formula for these quantities can be deduced from those of $\tau_t(s)$ (Eq. 45), $\alpha_t(s)$ (Eq. 52) and $\beta_t(s)$
(Eq. 65). To do so, Eqs. 45, 52 and 65 are respectively divided by the normalization terms $P(\sigma_t, \dots, \sigma_T | \sigma_{t-1}; \hat{\theta})$, $P(X_t^1 = x_t^1 | X_{1-p}^0, \Sigma; \hat{\theta})$ and $P(X_t^T = x_t^T | X_{1-p}^{T-1}, \Sigma; \hat{\theta})$. After decomposing the formula obtained and after some calculations, we obtain the subsequent recursive formulas for $\tilde{\tau}_t$, $\tilde{\alpha}_t$ and $\tilde{\beta}_t$.

First backward propagation

$$ 
\left\{ 
\begin{aligned}
\tilde{\tau}_T(s) &= \frac{1}{P(\sigma_T | \sigma_{T-1}; \hat{\theta})} \\
\tilde{\tau}_t(s) &= \sum_{i \in \sigma_{t+1}} \tilde{\tau}_{t+1}(i) \frac{P(S_{t+1} = i | S_t = s; \hat{\theta})}{P(\sigma_t | \sigma_{t-1}; \hat{\theta})}, \quad \text{for } t = T-1, \dots, 1,
\end{aligned}
\right. 
\qquad (80) 
$$

with

$$ 
\begin{aligned}
P(\sigma_t | \sigma_{t-1}; \hat{\theta}) &= P(S_t \in \sigma_t | S_{t-1} \in \sigma_{t-1} \hat{\theta}) = \sum_{i \in \sigma_{t-1}} \sum_{j \in \sigma_t} P(S_t = j | S_{t-1} = i; \hat{\theta}). \\
P(\sigma_1; \hat{\theta}) &= P(S_1 \in \sigma_1 \hat{\theta}) = \sum_{i \in \sigma_1} P(S_1 = i; \hat{\theta}).
\end{aligned} 
\qquad (81) 
$$

Forward propagation

$$ 
\left\{ 
\begin{aligned}
\tilde{\alpha}_1(s) &= \frac{P(X_1 = x_1 | X_{1-p}^0, S_1 = s; \hat{\theta}) P(S_1 = s; \hat{\theta})}{C_1} \times \frac{\tilde{\tau}_1(s)}{\sum_{i \in \sigma_1} \tilde{\tau}_1(i) P(S_1 = i; \hat{\theta})} \\
\tilde{\alpha}_t(s) &= P(X_t = x_t | X_{t-p}^{T-1}, S_t = s; \hat{\theta}) \left[ \sum_{i \in \sigma_{t-1}} \tilde{\alpha}_{t-1}(i) P(S_t = s | S_{t-1} = i; \hat{\theta}) \frac{\tilde{\tau}_t(s)}{\tilde{\tau}_{t-1}(i)} \right] \\
&\quad \times \frac{1}{P(\sigma_{t-1} | \sigma_{t-2}; \hat{\theta}) C_t} \times \mathbf{1}_{\{s \in \sigma_t\}}
\end{aligned}
\right. 
\qquad (82) 
$$

with $P(\sigma_{t-1} | \sigma_{t-2}; \hat{\theta})$ defined in Eq. 81 and $C_t$ the scaling term defined and computed as follows:

$$ C_1 = P(X_1 = x_1 | X_{1-p}^0, \Sigma; \hat{\theta}) = \sum_{i \in \sigma_1} P(X_1 = x_1 | X_{1-p}^0, S_1 = i; \hat{\theta}) P(S_1 = i; \hat{\theta}) \qquad (83) $$

$$ C_t = P(X_t = x_t | X_{1-p}^{T-1}, \Sigma; \hat{\theta}) \quad \text{for } t = 2, \dots, T \qquad (84) $$

$$ = \sum_{s \in \sigma_t} P(X_t = x_t | X_{t-p}^{T-1}, S_t = s; \hat{\theta}) \times \left[ \sum_{i \in \sigma_{t-1}} \tilde{\alpha}_{t-1}(i) P(S_t = s | S_{t-1} = i; \hat{\theta}) \right]. \qquad (85) $$
(Eq. 65). To do so, Eqs. 45, 52 and 65 are respectively divided by the normalization terms $P(\sigma_t, \dots, \sigma_T | \sigma_{t-1}; \hat{\theta})$, $P(X_t^1 = x_t^1 | X_{1-p}^0, \Sigma; \hat{\theta})$ and $P(X_t^T = x_t^T | X_{1-p}^{T-1}, \Sigma; \hat{\theta})$. After decomposing the formula obtained and after some calculations, we obtain the subsequent recursive formulas for $\tilde{\tau}_t$, $\tilde{\alpha}_t$ and $\tilde{\beta}_t$.

First backward propagation

$$ 
\left\{ 
\begin{aligned}
\tilde{\tau}_T(s) &= \frac{1}{P(\sigma_T | \sigma_{T-1}; \hat{\theta})} \\
\tilde{\tau}_t(s) &= \sum_{i \in \sigma_{t+1}} \tilde{\tau}_{t+1}(i) \frac{P(S_{t+1} = i | S_t = s; \hat{\theta})}{P(\sigma_t | \sigma_{t-1}; \hat{\theta})}, \quad \text{for } t = T-1, \dots, 1,
\end{aligned}
\right. 
\quad (80) 
$$

with

$$ 
\begin{aligned}
P(\sigma_t | \sigma_{t-1}; \hat{\theta}) &= P(S_t \in \sigma_t | S_{t-1} \in \sigma_{t-1} \hat{\theta}) = \sum_{i \in \sigma_{t-1}} \sum_{j \in \sigma_t} P(S_t = j | S_{t-1} = i; \hat{\theta}). \\
P(\sigma_1; \hat{\theta}) &= P(S_1 \in \sigma_1 \hat{\theta}) = \sum_{i \in \sigma_1} P(S_1 = i; \hat{\theta}).
\end{aligned} 
\quad (81) 
$$

Forward propagation

$$ 
\left\{ 
\begin{aligned}
\tilde{\alpha}_1(s) &= \frac{P(X_1 = x_1 | X_{1-p}^0, S_1 = s; \hat{\theta}) P(S_1 = s; \hat{\theta})}{C_1} \times \frac{\tilde{\tau}_1(s)}{\sum_{i \in \sigma_1} \tilde{\tau}_1(i) P(S_1 = i; \hat{\theta})} \\
\tilde{\alpha}_t(s) &= P(X_t = x_t | X_{t-p}^{T-1}, S_t = s; \hat{\theta}) \left[ \sum_{i \in \sigma_{t-1}} \tilde{\alpha}_{t-1}(i) P(S_t = s | S_{t-1} = i; \hat{\theta}) \frac{\tilde{\tau}_t(s)}{\tilde{\tau}_{t-1}(i)} \right] \\
&\quad \times \frac{1}{P(\sigma_{t-1} | \sigma_{t-2}; \hat{\theta}) C_t} \times \mathbf{1}_{\{s \in \sigma_t\}}
\end{aligned}
\right. 
\quad (82) 
$$

with $P(\sigma_{t-1} | \sigma_{t-2}; \hat{\theta})$ defined in Eq. 81 and $C_t$ the scaling term defined and computed as follows:

$$ C_1 = P(X_1 = x_1 | X_{1-p}^0, \Sigma; \hat{\theta}) = \sum_{i \in \sigma_1} P(X_1 = x_1 | X_{1-p}^0, S_1 = i; \hat{\theta}) P(S_1 = i; \hat{\theta}) \quad (83) $$

$$ C_t = P(X_t = x_t | X_{1-p}^{T-1}, \Sigma; \hat{\theta}) \quad \text{for } t = 2, \dots, T \quad (84) $$

$$ = \sum_{s \in \sigma_t} P(X_t = x_t | X_{t-p}^{T-1}, S_t = s; \hat{\theta}) \times \left[ \sum_{i \in \sigma_{t-1}} \tilde{\alpha}_{t-1}(i) P(S_t = s | S_{t-1} = i; \hat{\theta}) \right]. \quad (85) $$
The proof is straightforward and is left to the reader. Note that
$P(X_1^T = x_1^T | X_{1-p}^0, \Sigma: \hat{\theta}) = \prod_{t=1}^T C_t$.

Second backward propagation

$$
\left\{
\begin{aligned}
\tilde{\beta}_T(s) &= \frac{1}{C_T} \\
\tilde{\beta}_t(s) &= \sum_{i \in \sigma_{t+1}} \left[ \tilde{\beta}_{t+1}(i) P(S_{t+1} = i | S_t = s; \hat{\theta}) \frac{\tilde{\tau}_{t+1}(i)}{\tilde{\tau}_t(s)} P(X_{t+1} = x_{t+1} | X_{t+1-p}^\prime, S_{t+1} = i; \hat{\theta}) \right] \\
&\quad \times \frac{1}{P(\sigma_t | \sigma_{t-1}; \hat{\theta}) C_t} \times \mathbf{1}_{\{s \in \sigma_t\}}
\end{aligned}
\right.
\qquad (86)
$$

where $C_t$ and $P(\sigma_t | \sigma_{t-1}; \hat{\theta})$ are defined in Eqs. 83–84 and Eq. 81 respectively.

$\xi_t(k, \ell)$ computation

In Eq. 37 probabilities $\xi_t(k, \ell)$ are defined in function of quantities $\tau_t$, $\tau_{t-1}$, $\alpha_{t-1}$ and $\beta_t$. These quantities can be easily expressed in function of their normalized versions $\tilde{\tau}_t$, $\tilde{\tau}_{t-1}$, $\tilde{\alpha}_{t-1}$ and $\tilde{\beta}_t$ using Eqs. 77, 78 and 79. After substituting $\tau_t$, $\tau_{t-1}$, $\alpha_{t-1}$ and $\beta_t$ by the resulting expressions and after some simplifications, we obtain the following formula:

$$
\xi_t(k, \ell) = \frac{\tilde{\beta}_t(\ell) P(S_t = \ell | S_{t-1} = k; \hat{\theta}) P(X_t = x_t | X_{t-p}^{t-1}, S_t = \ell; \hat{\theta}) \tilde{\alpha}_{t-1}(k) \tilde{\tau}_t(\ell)}{P(\sigma_{t-1} | \sigma_{t-2}; \hat{\theta}) \tilde{\tau}_{t-1}(k)} \times \mathbf{1}_{\{\ell \in \sigma_t, k \in \sigma_{t-1}\}} \quad (87)
$$

B Appendix: evaluation of feature prediction performance
on the CMAPSS real-world data

See Tables 6 and 7.
The proof is straightforward and is left to the reader. Note that
$P(X_1^T = x_1^T | X_{1-p}^0, \Sigma: \hat{\theta}) = \prod_{t=1}^T C_t$.

Second backward propagation

$$
\left\{
\begin{aligned}
\tilde{\beta}_T(s) &= \frac{1}{C_T} \\
\tilde{\beta}_t(s) &= \sum_{i \in \sigma_{t+1}} \left[ \tilde{\beta}_{t+1}(i) P(S_{t+1} = i | S_t = s; \hat{\theta}) \frac{\tilde{\tau}_{t+1}(i)}{\tilde{\tau}_t(s)} P(X_{t+1} = x_{t+1} | X_{t+1-p}^\prime, S_{t+1} = i; \hat{\theta}) \right] \\
&\quad \times \frac{1}{P(\sigma_t | \sigma_{t-1}; \hat{\theta}) C_t} \times \mathbf{1}_{\{s \in \sigma_t\}}
\end{aligned}
\right.
\qquad (86)
$$

where $C_t$ and $P(\sigma_t | \sigma_{t-1}; \hat{\theta})$ are defined in Eqs. 83–84 and Eq. 81 respectively.

$\xi_t(k, \ell)$ computation

In Eq. 37 probabilities $\xi_t(k, \ell)$ are defined in function of quantities $\tau_t$, $\tau_{t-1}$, $\alpha_{t-1}$ and $\beta_t$. These quantities can be easily expressed in function of their normalized versions $\tilde{\tau}_t$, $\tilde{\tau}_{t-1}$, $\tilde{\alpha}_{t-1}$ and $\tilde{\beta}_t$ using Eqs. 77, 78 and 79. After substituting $\tau_t$, $\tau_{t-1}$, $\alpha_{t-1}$ and $\beta_t$ by the resulting expressions and after some simplifications, we obtain the following formula:

$$
\xi_t(k, \ell) = \frac{\tilde{\beta}_t(\ell) P(S_t = \ell | S_{t-1} = k; \hat{\theta}) P(X_t = x_t | X_{t-p}^{t-1}, S_t = \ell; \hat{\theta}) \tilde{\alpha}_{t-1}(k) \tilde{\tau}_t(\ell)}{P(\sigma_{t-1} | \sigma_{t-2}; \hat{\theta}) \tilde{\tau}_{t-1}(k)} \times \mathbf{1}_{\{\ell \in \sigma_t, k \in \sigma_{t-1}\}} \quad (87)
$$

B Appendix: evaluation of feature prediction performance
on the CMAPSS real-world data

See Tables 6 and 7.
<table><caption>**Table 6** Comparison of feature prediction performances for PHMC-MLAR, PHMC and MSAR models on dataset FD001</caption><thead><tr><th>Features</th><th>Models</th><th>h = 5</th><th>h = 10</th><th>h = 20</th><th>h = 30</th></tr></thead><tbody><tr><td rowspan="3">F2</td><td>PHMC-MLAR</td><td>0.296 [0.214, 0.386]</td><td>0.304 [0.199, 0.413]</td><td>0.315 [0.217, 0.405]</td><td>0.337 [0.218, 0.459]</td></tr><tr><td>PHMC</td><td>0.340 [0.233, 0.452]</td><td>0.345 [0.247, 0.467]</td><td>0.349 [0.243, 0.487]</td><td>0.353 [0.244, 0.482]</td></tr><tr><td>MSAR</td><td><strong>0.295</strong> [0.208, 0.394]</td><td><strong>0.300</strong> [0.206, 0.397]</td><td><strong>0.303</strong> [0.216, 0.387]</td><td><strong>0.311</strong> [0.230, 0.400]</td></tr><tr><td rowspan="3">F3</td><td>PHMC-MLAR</td><td>4.023 [2.698, 5.674]</td><td>4.043 [2.570, 5.674]</td><td>4.195 [2.663, 6.153]</td><td>4.379 [2.811, 5.848]</td></tr><tr><td>PHMC</td><td>4.385 [3.052, 6.074]</td><td>4.390 [2.965, 6.345]</td><td>4.390 [2.753, 6.157]</td><td>4.385 [2.854, 6.030]</td></tr><tr><td>MSAR</td><td><strong>4.003</strong> [2.677, 5.597]</td><td><strong>4.003</strong> [2.626, 5.667]</td><td><strong>4.024</strong> [2.439, 5.879]</td><td><strong>4.006</strong> [2.645, 5.475]</td></tr><tr><td rowspan="3">F4</td><td>PHMC-MLAR</td><td>4.108 [2.831, 5.255]</td><td>4.148 [2.789, 5.441]</td><td>4.489 [3.112, 5.908]</td><td>5.168 [3.164, 7.013]</td></tr><tr><td>PHMC</td><td>5.336 [3.391, 7.929]</td><td>5.327 [3.175, 8.089]</td><td>5.409 [3.696, 7.966]</td><td>5.548 [3.908, 8.416]</td></tr><tr><td>MSAR</td><td><strong>4.084</strong> [2.865, 5.232]</td><td><strong>4.089</strong> [2.780, 5.311]</td><td><strong>4.225</strong> [3.007, 5.543]</td><td><strong>4.491</strong> [3.049, 6.280]</td></tr><tr><td rowspan="3">F7</td><td>PHMC-MLAR</td><td>0.418 [0.272, 0.555]</td><td>0.429 [0.303, 0.609]</td><td>0.478 [0.351, 0.712]</td><td>0.553 [0.370, 0.775]</td></tr><tr><td>PHMC</td><td>0.561 [0.346, 0.872]</td><td>0.561 [0.362, 0.868]</td><td>0.576 [0.367, 0.898]</td><td>0.592 [0.382, 0.936]</td></tr><tr><td>MSAR</td><td><strong>0.415</strong> [0.262, 0.547]</td><td><strong>0.420</strong> [0.296, 0.584]</td><td><strong>0.440</strong> [0.306, 0.646]</td><td><strong>0.461</strong> [0.295, 0.636]</td></tr><tr><td rowspan="3">F9</td><td>PHMC-MLAR</td><td><strong>4.303</strong> [2.791, 5.616]</td><td><strong>4.487</strong> [2.840, 6.063]</td><td><strong>4.971</strong> [3.237, 6.987]</td><td><strong>5.886</strong> [3.720, 9.157]</td></tr><tr><td>PHMC</td><td>9.581 [4.373, 24.60]</td><td>9.629 [4.244, 24.82]</td><td>10.18 [4.390, 25.98]</td><td>10.45 [4.381, 25.36]</td></tr><tr><td>MSAR</td><td>4.363 [2.599, 5.826]</td><td>4.674 [2.772, 6.310]</td><td>5.533 [3.526, 7.300]</td><td>6.868 [4.023, 9.986]</td></tr><tr><td rowspan="3">F11</td><td>PHMC-MLAR</td><td><strong>0.104</strong> [0.083, 0.128]</td><td>0.108 [0.080, 0.137]</td><td>0.119 [0.080, 0.167]</td><td>0.142 [0.093, 0.209]</td></tr><tr><td>PHMC</td><td>0.146 [0.091, 0.216]</td><td>0.148 [0.087, 0.216]</td><td>0.151 [0.089, 0.225]</td><td>0.156 [0.093, 0.245]</td></tr><tr><td>MSAR</td><td><strong>0.104</strong> [0.078, 0.129]</td><td><strong>0.106</strong> [0.080, 0.133]</td><td><strong>0.111</strong> [0.081, 0.142]</td><td><strong>0.118</strong> [0.080, 0.164]</td></tr><tr><td rowspan="3">F12</td><td>PHMC-MLAR</td><td>0.309 [0.199, 0.412]</td><td>0.318 [0.193, 0.428]</td><td>0.356 [0.232, 0.482]</td><td>0.420 [0.255, 0.593]</td></tr><tr><td>PHMC</td><td>0.437 [0.278, 0.684]</td><td>0.444 [0.255, 0.725]</td><td>0.453 [0.279, 0.699]</td><td>0.465 [0.273, 0.742]</td></tr><tr><td>MSAR</td><td><strong>0.307</strong> [0.196, 0.399]</td><td><strong>0.312</strong> [0.192, 0.412]</td><td><strong>0.325</strong> [0.199, 0.413]</td><td><strong>0.349</strong> [0.228, 0.461]</td></tr><tr><td rowspan="3">F14</td><td>PHMC-MLAR</td><td><strong>3.343</strong> [2.447, 4.375]</td><td><strong>3.402</strong> [2.295, 4.337]</td><td><strong>3.774</strong> [2.296, 5.549]</td><td><strong>4.399</strong> [2.395, 7.693]</td></tr><tr><td>PHMC</td><td>8.679 [3.435, 23.03]</td><td>8.677 [3.400, 23.68]</td><td>9.047 [3.403, 24.39]</td><td>9.251 [3.515, 23.56]</td></tr><tr><td>MSAR</td><td>3.378 [2.521, 4.343]</td><td>3.582 [2.371, 4.687]</td><td>4.395 [2.511, 6.928]</td><td>5.522 [3.097, 9.282]</td></tr></tbody></table>

Mean and 95% confidence interval of RMSE at different forecast horizons *h* for the eight features (2, 3, 4, 7, 9, 11, 12, 14) used in the training stage. The figures in bold highlight the minimum mean RMSE across the three models PHMC-MLAR(K = 4, p = 7), PHMC(K = 4) and MSAR(p = 10) for each pair (feature, forecast horizon)
<table><caption>**Table 6** Comparison of feature prediction performances for PHMC-MLAR, PHMC and MSAR models on dataset FD001</caption><thead><tr><th>Features</th><th>Models</th><th>h = 5</th><th>h = 10</th><th>h = 20</th><th>h = 30</th></tr></thead><tbody><tr><td rowspan="3">F2</td><td>PHMC-MLAR</td><td>0.296 [0.214, 0.386]</td><td>0.304 [0.199, 0.413]</td><td>0.315 [0.217, 0.405]</td><td>0.337 [0.218, 0.459]</td></tr><tr><td>PHMC</td><td>0.340 [0.233, 0.452]</td><td>0.345 [0.247, 0.467]</td><td>0.349 [0.243, 0.487]</td><td>0.353 [0.244, 0.482]</td></tr><tr><td>MSAR</td><td><b>0.295</b> [0.208, 0.394]</td><td><b>0.300</b> [0.206, 0.397]</td><td><b>0.303</b> [0.216, 0.387]</td><td><b>0.311</b> [0.230, 0.400]</td></tr><tr><td rowspan="3">F3</td><td>PHMC-MLAR</td><td>4.023 [2.698, 5.674]</td><td>4.043 [2.570, 5.674]</td><td>4.195 [2.663, 6.153]</td><td>4.379 [2.811, 5.848]</td></tr><tr><td>PHMC</td><td>4.385 [3.052, 6.074]</td><td>4.390 [2.965, 6.345]</td><td>4.390 [2.753, 6.157]</td><td>4.385 [2.854, 6.030]</td></tr><tr><td>MSAR</td><td><b>4.003</b> [2.677, 5.597]</td><td><b>4.003</b> [2.626, 5.667]</td><td><b>4.024</b> [2.439, 5.879]</td><td><b>4.006</b> [2.645, 5.475]</td></tr><tr><td rowspan="3">F4</td><td>PHMC-MLAR</td><td>4.108 [2.831, 5.255]</td><td>4.148 [2.789, 5.441]</td><td>4.489 [3.112, 5.908]</td><td>5.168 [3.164, 7.013]</td></tr><tr><td>PHMC</td><td>5.336 [3.391, 7.929]</td><td>5.327 [3.175, 8.089]</td><td>5.409 [3.696, 7.966]</td><td>5.548 [3.908, 8.416]</td></tr><tr><td>MSAR</td><td><b>4.084</b> [2.865, 5.232]</td><td><b>4.089</b> [2.780, 5.311]</td><td><b>4.225</b> [3.007, 5.543]</td><td><b>4.491</b> [3.049, 6.280]</td></tr><tr><td rowspan="3">F7</td><td>PHMC-MLAR</td><td>0.418 [0.272, 0.555]</td><td>0.429 [0.303, 0.609]</td><td>0.478 [0.351, 0.712]</td><td>0.553 [0.370, 0.775]</td></tr><tr><td>PHMC</td><td>0.561 [0.346, 0.872]</td><td>0.561 [0.362, 0.868]</td><td>0.576 [0.367, 0.898]</td><td>0.592 [0.382, 0.936]</td></tr><tr><td>MSAR</td><td><b>0.415</b> [0.262, 0.547]</td><td><b>0.420</b> [0.296, 0.584]</td><td><b>0.440</b> [0.306, 0.646]</td><td><b>0.461</b> [0.295, 0.636]</td></tr><tr><td rowspan="3">F9</td><td>PHMC-MLAR</td><td><b>4.303</b> [2.791, 5.616]</td><td><b>4.487</b> [2.840, 6.063]</td><td><b>4.971</b> [3.237, 6.987]</td><td><b>5.886</b> [3.720, 9.157]</td></tr><tr><td>PHMC</td><td>9.581 [4.373, 24.60]</td><td>9.629 [4.244, 24.82]</td><td>10.18 [4.390, 25.98]</td><td>10.45 [4.381, 25.36]</td></tr><tr><td>MSAR</td><td>4.363 [2.599, 5.826]</td><td>4.674 [2.772, 6.310]</td><td>5.533 [3.526, 7.300]</td><td>6.868 [4.023, 9.986]</td></tr><tr><td rowspan="3">F11</td><td>PHMC-MLAR</td><td><b>0.104</b> [0.083, 0.128]</td><td>0.108 [0.080, 0.137]</td><td>0.119 [0.080, 0.167]</td><td>0.142 [0.093, 0.209]</td></tr><tr><td>PHMC</td><td>0.146 [0.091, 0.216]</td><td>0.148 [0.087, 0.216]</td><td>0.151 [0.089, 0.225]</td><td>0.156 [0.093, 0.245]</td></tr><tr><td>MSAR</td><td><b>0.104</b> [0.078, 0.129]</td><td><b>0.106</b> [0.080, 0.133]</td><td><b>0.111</b> [0.081, 0.142]</td><td><b>0.118</b> [0.080, 0.164]</td></tr><tr><td rowspan="3">F12</td><td>PHMC-MLAR</td><td>0.309 [0.199, 0.412]</td><td>0.318 [0.193, 0.428]</td><td>0.356 [0.232, 0.482]</td><td>0.420 [0.255, 0.593]</td></tr><tr><td>PHMC</td><td>0.437 [0.278, 0.684]</td><td>0.444 [0.255, 0.725]</td><td>0.453 [0.279, 0.699]</td><td>0.465 [0.273, 0.742]</td></tr><tr><td>MSAR</td><td><b>0.307</b> [0.196, 0.399]</td><td><b>0.312</b> [0.192, 0.412]</td><td><b>0.325</b> [0.199, 0.413]</td><td><b>0.349</b> [0.228, 0.461]</td></tr><tr><td rowspan="3">F14</td><td>PHMC-MLAR</td><td><b>3.343</b> [2.447, 4.375]</td><td><b>3.402</b> [2.295, 4.337]</td><td><b>3.774</b> [2.296, 5.549]</td><td><b>4.399</b> [2.395, 7.693]</td></tr><tr><td>PHMC</td><td>8.679 [3.435, 23.03]</td><td>8.677 [3.400, 23.68]</td><td>9.047 [3.403, 24.39]</td><td>9.251 [3.515, 23.56]</td></tr><tr><td>MSAR</td><td>3.378 [2.521, 4.343]</td><td>3.582 [2.371, 4.687]</td><td>4.395 [2.511, 6.928]</td><td>5.522 [3.097, 9.282]</td></tr></tbody></table>

Mean and 95% confidence interval of RMSE at different forecast horizons *h* for the eight features (2, 3, 4, 7, 9, 11, 12, 14) used in the training stage. The figures in bold highlight the minimum mean RMSE across the three models PHMC-MLAR(K = 4, p = 7), PHMC(K = 4) and MSAR(p = 10) for each pair (feature, forecast horizon)
<table><caption>**Table 7** Comparison of feature prediction performances for PHMC-MLAR, PHMC and MSAR models on dataset FD003</caption><thead><tr><th>Features</th><th>Models</th><th>h=5</th><th>h=10</th><th>h=20</th><th>h=30</th></tr></thead><tbody><tr><td rowspan="3">F2</td><td>PHMC-MLAR</td><td>0.306 [0.230, 0.376]</td><td>0.307 [0.240, 0.377]</td><td><strong>0.306</strong> [0.240, 0.408]</td><td><strong>0.308</strong> [0.233, 0.386]</td></tr><tr><td>PHMC</td><td>0.361 [0.480, 0.274]</td><td>0.363 [0.268, 0.478]</td><td>0.363 [0.265, 0.482]</td><td>0.369 [0.240, 0.481]</td></tr><tr><td>MSAR</td><td><strong>0.305</strong> [0.230, 0.379]</td><td><strong>0.306</strong> [0.242, 0.377]</td><td><strong>0.306</strong> [0.234, 0.402]</td><td><strong>0.308</strong> [0.234, 0.382]</td></tr><tr><td rowspan="3">F3</td><td>PHMC-MLAR</td><td>4.175 [2.989, 5.384]</td><td>4.161 [3.117, 5.298]</td><td>4.145 [3.039, 5.292]</td><td>4.142 [2.700, 5.322]</td></tr><tr><td>PHMC</td><td>4.790 [3.276, 6.781]</td><td>4.793 [3.302, 6.882]</td><td>4.868 [3.369, 7.342]</td><td>4.940 [3.544, 7.632]</td></tr><tr><td>MSAR</td><td><strong>4.116</strong> [2.994, 5.235]</td><td><strong>4.106</strong> [3.106, 5.283]</td><td><strong>4.080</strong> [2.976, 5.255]</td><td><strong>4.097</strong> [2.526, 5.302]</td></tr><tr><td rowspan="3">F4</td><td>PHMC-MLAR</td><td><strong>4.089</strong> [2.790, 5.174]</td><td>4.134 [2.746, 5.329]</td><td><strong>4.229</strong> [2.545, 5.686]</td><td><strong>4.343</strong> [3.090, 5.976]</td></tr><tr><td>PHMC</td><td>5.796 [4.088, 8.876]</td><td>5.818 [3.844, 8.958]</td><td>5.951 [3.559, 8.897]</td><td>6.097 [3.662, 9.492]</td></tr><tr><td>MSAR</td><td>4.101 [2.771, 5.227]</td><td><strong>4.118</strong> [2.777, 5.217]</td><td>4.264 [2.527, 5.833]</td><td>4.390 [2.509, 5.973]</td></tr><tr><td rowspan="3">F7</td><td>PHMC-MLAR</td><td>0.485 [0.315, 0.802]</td><td>0.517 [0.330, 0.878]</td><td>0.613 [0.307, 1.409]</td><td>0.768 [0.356, 2.063]</td></tr><tr><td>PHMC</td><td>1.601 [0.454, 3.101]</td><td>1.653 [0.447, 3.118]</td><td>1.754 [0.403, 3.334]</td><td>1.850 [0.416, 3.461]</td></tr><tr><td>MSAR</td><td><strong>0.480</strong> [0.316, 0.837]</td><td><strong>0.506</strong> [0.313, 0.871]</td><td><strong>0.567</strong> [0.304, 0.970]</td><td><strong>0.689</strong> [0.348, 1.474]</td></tr><tr><td rowspan="3">F9</td><td>PHMC-MLAR</td><td><strong>4.347</strong> [3.105, 5.674]</td><td><strong>4.472</strong> [3.220, 5.932]</td><td><strong>4.748</strong> [3.213, 6.906]</td><td><strong>5.349</strong> [3.415, 9.431]</td></tr><tr><td>PHMC</td><td>8.609 [3.637, 19.27]</td><td>8.809 [3.747, 20.30]</td><td>9.234 [3.354, 23.72]</td><td>9.635 [3.690, 19.35]</td></tr><tr><td>MSAR</td><td>4.415 [3.222, 5.854]</td><td>4.601 [3.465, 6.077]</td><td>5.099 [3.210, 7.041]</td><td>6.069 [3.906, 9.781]</td></tr><tr><td rowspan="3">F11</td><td>PHMC-MLAR</td><td>0.106 [0.080, 0.145]</td><td><strong>0.107</strong> [0.082, 0.142]</td><td><strong>0.112</strong> [0.085, 0.150]</td><td><strong>0.117</strong> [0.077, 0.158]</td></tr><tr><td>PHMC</td><td>0.168 [0.098, 0.240]</td><td>0.167 [0.101, 0.233]</td><td>0.173 [0.105, 0.250]</td><td>0.176 [0.098, 0.275]</td></tr><tr><td>MSAR</td><td><strong>0.105</strong> [0.081, 0.144]</td><td>0.108 [0.080, 0.142]</td><td>0.113 [0.084, 0.153]</td><td>0.118 [0.070, 0.160]</td></tr><tr><td rowspan="3">F12</td><td>PHMC-MLAR</td><td>0.407 [0.251, 0.704]</td><td>0.442 [0.256, 0.802]</td><td>0.554 [0.261, 1.312]</td><td>0.709 [0.275, 1.959]</td></tr><tr><td>PHMC</td><td>1.478 [0.357, 2.928]</td><td>1.531 [0.338, 2.882]</td><td>1.634 [0.308, 3.048]</td><td>1.733 [0.316, 3.254]</td></tr><tr><td>MSAR</td><td><strong>0.395</strong> [0.241, 0.712]</td><td><strong>0.421</strong> [0.250, 0.722]</td><td><strong>0.493</strong> [0.252, 0.973]</td><td><strong>0.616</strong> [0.278, 1.359]</td></tr><tr><td rowspan="3">F14</td><td>PHMC-MLAR</td><td><strong>3.401</strong> [2.065, 4.880]</td><td><strong>3.511</strong> [2.160, 5.268]</td><td><strong>3.837</strong> [2.352, 6.051]</td><td><strong>4.272</strong> [2.587, 7.920]</td></tr><tr><td>PHMC</td><td>7.496 [2.758, 19.05]</td><td>7.693 [3.053, 20.07]</td><td>8.039 [2.882, 22.58]</td><td>8.244 [3.114, 19.40]</td></tr><tr><td>MSAR</td><td>3.480 [2.094, 4.761]</td><td>3.697 [2.356, 5.249]</td><td>4.212 [2.520, 6.507]</td><td>4.980 [2.892, 8.583]</td></tr></tbody></table>

Mean and 95% confidence interval of RMSE at different forecast horizons *h* for the eight features (2, 3, 4, 7, 9, 11, 12, 14) used in training stage. The figures in bold highlight the minimum mean RMSE across the three models PHMC-MLAR(K = 4, p = 13), PHMC(K = 4) and MSAR(p = 10) for each pair (feature, forecast horizon)
<table><caption>**Table 7** Comparison of feature prediction performances for PHMC-MLAR, PHMC and MSAR models on dataset FD003</caption><thead><tr><th>Features</th><th>Models</th><th>h=5</th><th>h=10</th><th>h=20</th><th>h=30</th></tr></thead><tbody><tr><td rowspan="3">F2</td><td>PHMC-MLAR</td><td>0.306 [0.230, 0.376]</td><td>0.307 [0.240, 0.377]</td><td><strong>0.306</strong> [0.240, 0.408]</td><td><strong>0.308</strong> [0.233, 0.386]</td></tr><tr><td>PHMC</td><td>0.361 [0.480, 0.274]</td><td>0.363 [0.268, 0.478]</td><td>0.363 [0.265, 0.482]</td><td>0.369 [0.240, 0.481]</td></tr><tr><td>MSAR</td><td><strong>0.305</strong> [0.230, 0.379]</td><td><strong>0.306</strong> [0.242, 0.377]</td><td><strong>0.306</strong> [0.234, 0.402]</td><td><strong>0.308</strong> [0.234, 0.382]</td></tr><tr><td rowspan="3">F3</td><td>PHMC-MLAR</td><td>4.175 [2.989, 5.384]</td><td>4.161 [3.117, 5.298]</td><td>4.145 [3.039, 5.292]</td><td>4.142 [2.700, 5.322]</td></tr><tr><td>PHMC</td><td>4.790 [3.276, 6.781]</td><td>4.793 [3.302, 6.882]</td><td>4.868 [3.369, 7.342]</td><td>4.940 [3.544, 7.632]</td></tr><tr><td>MSAR</td><td><strong>4.116</strong> [2.994, 5.235]</td><td><strong>4.106</strong> [3.106, 5.283]</td><td><strong>4.080</strong> [2.976, 5.255]</td><td><strong>4.097</strong> [2.526, 5.302]</td></tr><tr><td rowspan="3">F4</td><td>PHMC-MLAR</td><td><strong>4.089</strong> [2.790, 5.174]</td><td>4.134 [2.746, 5.329]</td><td><strong>4.229</strong> [2.545, 5.686]</td><td><strong>4.343</strong> [3.090, 5.976]</td></tr><tr><td>PHMC</td><td>5.796 [4.088, 8.876]</td><td>5.818 [3.844, 8.958]</td><td>5.951 [3.559, 8.897]</td><td>6.097 [3.662, 9.492]</td></tr><tr><td>MSAR</td><td>4.101 [2.771, 5.227]</td><td><strong>4.118</strong> [2.777, 5.217]</td><td>4.264 [2.527, 5.833]</td><td>4.390 [2.509, 5.973]</td></tr><tr><td rowspan="3">F7</td><td>PHMC-MLAR</td><td>0.485 [0.315, 0.802]</td><td>0.517 [0.330, 0.878]</td><td>0.613 [0.307, 1.409]</td><td>0.768 [0.356, 2.063]</td></tr><tr><td>PHMC</td><td>1.601 [0.454, 3.101]</td><td>1.653 [0.447, 3.118]</td><td>1.754 [0.403, 3.334]</td><td>1.850 [0.416, 3.461]</td></tr><tr><td>MSAR</td><td><strong>0.480</strong> [0.316, 0.837]</td><td><strong>0.506</strong> [0.313, 0.871]</td><td><strong>0.567</strong> [0.304, 0.970]</td><td><strong>0.689</strong> [0.348, 1.474]</td></tr><tr><td rowspan="3">F9</td><td>PHMC-MLAR</td><td><strong>4.347</strong> [3.105, 5.674]</td><td><strong>4.472</strong> [3.220, 5.932]</td><td><strong>4.748</strong> [3.213, 6.906]</td><td><strong>5.349</strong> [3.415, 9.431]</td></tr><tr><td>PHMC</td><td>8.609 [3.637, 19.27]</td><td>8.809 [3.747, 20.30]</td><td>9.234 [3.354, 23.72]</td><td>9.635 [3.690, 19.35]</td></tr><tr><td>MSAR</td><td>4.415 [3.222, 5.854]</td><td>4.601 [3.465, 6.077]</td><td>5.099 [3.210, 7.041]</td><td>6.069 [3.906, 9.781]</td></tr><tr><td rowspan="3">F11</td><td>PHMC-MLAR</td><td>0.106 [0.080, 0.145]</td><td><strong>0.107</strong> [0.082, 0.142]</td><td><strong>0.112</strong> [0.085, 0.150]</td><td><strong>0.117</strong> [0.077, 0.158]</td></tr><tr><td>PHMC</td><td>0.168 [0.098, 0.240]</td><td>0.167 [0.101, 0.233]</td><td>0.173 [0.105, 0.250]</td><td>0.176 [0.098, 0.275]</td></tr><tr><td>MSAR</td><td><strong>0.105</strong> [0.081, 0.144]</td><td>0.108 [0.080, 0.142]</td><td>0.113 [0.084, 0.153]</td><td>0.118 [0.070, 0.160]</td></tr><tr><td rowspan="3">F12</td><td>PHMC-MLAR</td><td>0.407 [0.251, 0.704]</td><td>0.442 [0.256, 0.802]</td><td>0.554 [0.261, 1.312]</td><td>0.709 [0.275, 1.959]</td></tr><tr><td>PHMC</td><td>1.478 [0.357, 2.928]</td><td>1.531 [0.338, 2.882]</td><td>1.634 [0.308, 3.048]</td><td>1.733 [0.316, 3.254]</td></tr><tr><td>MSAR</td><td><strong>0.395</strong> [0.241, 0.712]</td><td><strong>0.421</strong> [0.250, 0.722]</td><td><strong>0.493</strong> [0.252, 0.973]</td><td><strong>0.616</strong> [0.278, 1.359]</td></tr><tr><td rowspan="3">F14</td><td>PHMC-MLAR</td><td><strong>3.401</strong> [2.065, 4.880]</td><td><strong>3.511</strong> [2.160, 5.268]</td><td><strong>3.837</strong> [2.352, 6.051]</td><td><strong>4.272</strong> [2.587, 7.920]</td></tr><tr><td>PHMC</td><td>7.496 [2.758, 19.05]</td><td>7.693 [3.053, 20.07]</td><td>8.039 [2.882, 22.58]</td><td>8.244 [3.114, 19.40]</td></tr><tr><td>MSAR</td><td>3.480 [2.094, 4.761]</td><td>3.697 [2.356, 5.249]</td><td>4.212 [2.520, 6.507]</td><td>4.980 [2.892, 8.583]</td></tr></tbody></table>

Mean and 95% confidence interval of RMSE at different forecast horizons *h* for the eight features (2, 3, 4, 7, 9, 11, 12, 14) used in training stage. The figures in bold highlight the minimum mean RMSE across the three models PHMC-MLAR(K = 4, p = 13), PHMC(K = 4) and MSAR(p = 10) for each pair (feature, forecast horizon)
**Acknowledgements** The software development and the realization of the experiments were performed at the CCIPL (Centre de Calcul Intensif des Pays de la Loire, Nantes, France).

**Author Contributions** Not Applicable.

**Funding** Fatoumata Dama is supported by a PhD scholarship granted by the French Ministry for Higher Education, Research and Innovation.

**Data availability** The CMAPSS datasets are available from https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repo-sitory/#turbofan.

**Code availability** Not Applicable.

## Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

**Consent to participate** Not Applicable.

**Consent for publication** Not Applicable.

**Ethical approval** Not Applicable.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Ailliot, P., & Monbet, V. (2012). Markov-switching autoregressive models for wind time series. *Environmental Modelling & Software, 30*, 92–101.

Ailliot, P., Bessac, J., Monbet, V., & Pene, F. (2015). Non-homogeneous hidden Markov-switching models for wind time series. *Journal of Statistical Planning and Inference, 160*, 75–88.

Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control, 19*, 716–723.

Aremu, O. O., Cody, R. A., Hyland-Wood, D., & McAree, P. R. (2020). A relative entropy based feature selection framework for asset data in predictive maintenance. *Computers & Industrial Engineering, 145*, 106536.

Baum, L. E., Petrie, T., Soules, G., & Weiss, N. (1970). A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. *The Annals of Mathematical Statistics, 41*(1), 164–171.

Bauwens, L., Carpantier, J. F., & Dufays, A. (2017). Autoregressive moving average infinite hidden Markov-switching models. *Journal of Business and Economic Statistics, 35*(2), 162–182.

Berg, J., Reckordt, T., Richter, C., & Reinhart, G. (2018). Action recognition in assembly for human-robot-cooperation using Hidden Markov models. *Procedia CIRP, 76*, 205–210.

Bergmeir, C., Hyndman, R. J., & Benítez, J. M. (2016). Bagging exponential smoothing methods using STL decomposition and Box-Cox transformation. *International Journal of Forecasting, 32*(2), 303–312.

Bessac, J., Ailliot, P., Cattiaux, J., & Monbet, V. (2016). Comparison of hidden and observed regime-switching autoregressive models for (u, v)-components of wind fields in the Northeast Atlantic. *Advances in Statistical Climatology, Meteorology and Oceanography, 2*(1), 1–16.
**Acknowledgements** The software development and the realization of the experiments were performed at the CCIPL (Centre de Calcul Intensif des Pays de la Loire, Nantes, France).

**Author Contributions** Not Applicable.

**Funding** Fatoumata Dama is supported by a PhD scholarship granted by the French Ministry for Higher Education, Research and Innovation.

**Data availability** The CMAPSS datasets are available from https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repo-sitory/#turbofan.

**Code availability** Not Applicable.

## Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

**Consent to participate** Not Applicable.

**Consent for publication** Not Applicable.

**Ethical approval** Not Applicable.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Ailliot, P., & Monbet, V. (2012). Markov-switching autoregressive models for wind time series. *Environmental Modelling & Software*, 30, 92–101.

Ailliot, P., Bessac, J., Monbet, V., & Pene, F. (2015). Non-homogeneous hidden Markov-switching models for wind time series. *Journal of Statistical Planning and Inference*, 160, 75–88.

Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, 19, 716–723.

Aremu, O. O., Cody, R. A., Hyland-Wood, D., & McAree, P. R. (2020). A relative entropy based feature selection framework for asset data in predictive maintenance. *Computers & Industrial Engineering*, 145, 106536.

Baum, L. E., Petrie, T., Soules, G., & Weiss, N. (1970). A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. *The Annals of Mathematical Statistics*, 41(1), 164–171.

Bauwens, L., Carpantier, J. F., & Dufays, A. (2017). Autoregressive moving average infinite hidden Markov-switching models. *Journal of Business and Economic Statistics*, 35(2), 162–182.

Berg, J., Reckordt, T., Richter, C., & Reinhart, G. (2018). Action recognition in assembly for human-robot-cooperation using Hidden Markov models. *Procedia CIRP*, 76, 205–210.

Bergmeir, C., Hyndman, R. J., & Benítez, J. M. (2016). Bagging exponential smoothing methods using STL decomposition and Box-Cox transformation. *International Journal of Forecasting*, 32(2), 303–312.

Bessac, J., Ailliot, P., Cattiaux, J., & Monbet, V. (2016). Comparison of hidden and observed regime-switching autoregressive models for (u, v)-components of wind fields in the Northeast Atlantic. *Advances in Statistical Climatology, Meteorology and Oceanography*, 2(1), 1–16.
Bharathi, R., & Selvarani, R. (2020). Hidden Markov model approach for software reliability estimation with logic error. *International Journal of Automation and Computing*, 17, 305.

Box, G. E., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). *Time series analysis: Forecasting and control* (5th ed.). Wiley.

Cardenas-Gallo, I., Sanchez-Silva, M., Akhavan-Tabatabaei, R., & Bastidas-Arteaga, E. (2016). A Markov regime-switching framework application for describing El Niño Southern Oscillation (ENSO) patterns. *Natural Hazards*, 81(2), 829–843.

Clements, M. P., & Krolzig, H. M. (1998). A comparison of the forecast performance of Markov-switching and threshold autoregressive models of US GNP. *The Econometrics Journal*, 1(1), 47–75.

Degtyarev, A. B., & Gankevich, I. (2019). Evaluation of hydrodynamic pressures for autoregressive model of irregular waves. In *Contemporary ideas on ship stability* (pp. 37–47). Springer.

Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. *Journal of the Royal Statistical Society: Series B (Methodological)*, 39(1), 1–22.

Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. *Journal of the American Statistical Association*, 74(366), 427–431.

Durand, J. B. (2003). Modèles à structure cachée : Inférence, sélection de modèles et applications. PhD thesis, Grenoble 1. (in French)

Filardo, A. J. (1994). Business-cycle phases and their transitional dynamics. *Journal of Business & Economic Statistics*, 12(3), 299–308.

Flecher, C., Naveau, P., Allard, D., & Brisson, N. (2010). A stochastic daily weather generator for skewed data. *Water Resources Research*, 46, W07519.

Florez-Larrahondo, G. (2020). Incremental learning of discrete hidden Markov models. PhD thesis, Mississippi State University.

Forney, G. D. (1973). The Viterbi algorithm. *Proceedings of the IEEE*, 61(3), 268–278.

Fox, E., Sudderth, E., Jordan, M., & Willsky, A. (2010). Bayesian nonparametric learning of Markov switching processes. *IEEE Signal Processing Magazine*, 27(6), 43–54.

Fox, E., Sudderth, E., Jordan, M., & Willsky, A. (2011). A sticky HDP-HMM with application to speaker diarization. *Annals of Applied Statistics*, 5(2A), 1020–1056.

Gardner, E., Jr., & Everette, S. (2006). Exponential smoothing: The state of the art—Part II. *International Journal of Forecasting*, 22(4), 637–666.

Ghahramani, Z. (2001). An introduction to hidden Markov models and Bayesian networks. *International Journal of Pattern Recognition and Artificial Intelligence*, 15(1), 9–42.

Ghasvarian Jahromi, K., Gharavian, D., & Mahdiani, H. (2020). A novel method for day-ahead solar power prediction based on hidden Markov model and cosine similarity. *Soft Computing*, 24(7), 4991–5004.

Hamilton, J. D. (1989). A new approach to the economic analysis of nonstationary time series and the business cycle. *Econometrica*, 57, 357–384.

Hamilton, J. D. (1990). Analysis of time series subject to changes in regime. *Journal of Econometrics*, 45(1–2), 39–70.

Javed, K., Gouriveau, R., & Zerhouni, N. (2015). A new multivariate approach for prognostics based on extreme learning machine and fuzzy clustering. *IEEE Transactions on Cybernetics*, 45(12), 2626–2639.

Juesas, P., & Ramasso, E. (2016). Ascertainment-adjusted parameter estimation approach to improve robustness against misspecification of health monitoring methods. *Mechanical Systems and Signal Processing*, 81, 387–401.

Kim, C. J. (1994). Dynamic linear models with Markov-switching. *Journal of Econometrics*, 60, 1–22.

Koenig, S., & Simmons, R. G. (1996). Unsupervised learning of probabilistic models for robot navigation. *Proceedings of IEEE International Conference on Robotics and Automation*, IEEE, 3, 2301–2308.

Kuck, K., & Schweikert, K. (2017). A Markov regime-switching model of crude oil market integration. *Journal of Commodity Markets*, 6, 16–31.

Kwiatkowski, D., Phillips, P. C., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root? *Journal of Econometrics*, 54(1–3), 159–178.

Lhuissier, S. (2019). Bayesian inference for Markov-switching skewed autoregressive models. Banque de France working paper #726.

Li, K., & Fu, Y. (2012). ARMA-HMM: A new approach for early recognition of human activity. In *21st International Conference on Pattern Recognition (ICPR)* (pp 1779–1782).

Lim, P., Goh, C. K., Tan, K. C., & Dutta, P. (2015). Multimodal degradation prognostics based on switching Kalman filter ensemble. *IEEE Transactions on Neural Networks and Learning Systems*, 28(1), 136–148.
Bharathi, R., & Selvarani, R. (2020). Hidden Markov model approach for software reliability estimation with logic error. *International Journal of Automation and Computing*, 17, 305.

Box, G. E., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). *Time series analysis: Forecasting and control* (5th ed.). Wiley.

Cardenas-Gallo, I., Sanchez-Silva, M., Akhavan-Tabatabaei, R., & Bastidas-Arteaga, E. (2016). A Markov regime-switching framework application for describing El Niño Southern Oscillation (ENSO) patterns. *Natural Hazards*, 81(2), 829–843.

Clements, M. P., & Krolzig, H. M. (1998). A comparison of the forecast performance of Markov-switching and threshold autoregressive models of US GNP. *The Econometrics Journal*, 1(1), 47–75.

Degtyarev, A. B., & Gankevich, I. (2019). Evaluation of hydrodynamic pressures for autoregressive model of irregular waves. In *Contemporary ideas on ship stability* (pp. 37–47). Springer.

Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. *Journal of the Royal Statistical Society: Series B (Methodological)*, 39(1), 1–22.

Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. *Journal of the American Statistical Association*, 74(366), 427–431.

Durand, J. B. (2003). Modèles à structure cachée : Inférence, sélection de modèles et applications. PhD thesis, Grenoble 1. (in French)

Filardo, A. J. (1994). Business-cycle phases and their transitional dynamics. *Journal of Business & Economic Statistics*, 12(3), 299–308.

Flecher, C., Naveau, P., Allard, D., & Brisson, N. (2010). A stochastic daily weather generator for skewed data. *Water Resources Research*, 46, W07519.

Florez-Larrahondo, G. (2020). Incremental learning of discrete hidden Markov models. PhD thesis, Mississippi State University.

Forney, G. D. (1973). The Viterbi algorithm. *Proceedings of the IEEE*, 61(3), 268–278.

Fox, E., Sudderth, E., Jordan, M., & Willsky, A. (2010). Bayesian nonparametric learning of Markov switching processes. *IEEE Signal Processing Magazine*, 27(6), 43–54.

Fox, E., Sudderth, E., Jordan, M., & Willsky, A. (2011). A sticky HDP-HMM with application to speaker diarization. *Annals of Applied Statistics*, 5(2A), 1020–1056.

Gardner, E., Jr., & Everette, S. (2006). Exponential smoothing: The state of the art—Part II. *International Journal of Forecasting*, 22(4), 637–666.

Ghahramani, Z. (2001). An introduction to hidden Markov models and Bayesian networks. *International Journal of Pattern Recognition and Artificial Intelligence*, 15(1), 9–42.

Ghasvarian Jahromi, K., Gharavian, D., & Mahdiani, H. (2020). A novel method for day-ahead solar power prediction based on hidden Markov model and cosine similarity. *Soft Computing*, 24(7), 4991–5004.

Hamilton, J. D. (1989). A new approach to the economic analysis of nonstationary time series and the business cycle. *Econometrica*, 57, 357–384.

Hamilton, J. D. (1990). Analysis of time series subject to changes in regime. *Journal of Econometrics*, 45(1–2), 39–70.

Javed, K., Gouriveau, R., & Zerhouni, N. (2015). A new multivariate approach for prognostics based on extreme learning machine and fuzzy clustering. *IEEE Transactions on Cybernetics*, 45(12), 2626–2639.

Juesas, P., & Ramasso, E. (2016). Ascertainment-adjusted parameter estimation approach to improve robustness against misspecification of health monitoring methods. *Mechanical Systems and Signal Processing*, 81, 387–401.

Kim, C. J. (1994). Dynamic linear models with Markov-switching. *Journal of Econometrics*, 60, 1–22.

Koenig, S., & Simmons, R. G. (1996). Unsupervised learning of probabilistic models for robot navigation. *Proceedings of IEEE International Conference on Robotics and Automation*, IEEE, 3, 2301–2308.

Kuck, K., & Schweikert, K. (2017). A Markov regime-switching model of crude oil market integration. *Journal of Commodity Markets*, 6, 16–31.

Kwiatkowski, D., Phillips, P. C., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root? *Journal of Econometrics*, 54(1–3), 159–178.

Lhuissier, S. (2019). Bayesian inference for Markov-switching skewed autoregressive models. Banque de France working paper #726.

Li, K., & Fu, Y. (2012). ARMA-HMM: A new approach for early recognition of human activity. In *21st International Conference on Pattern Recognition (ICPR)* (pp 1779–1782).

Lim, P., Goh, C. K., Tan, K. C., & Dutta, P. (2015). Multimodal degradation prognostics based on switching Kalman filter ensemble. *IEEE Transactions on Neural Networks and Learning Systems*, 28(1), 136–148.
Michalek, S., Wagner, M., & Timmer, J. (2000). A new approximate likelihood estimator for ARMA-filtered Hidden Markov Models. *IEEE Transactions on Signal Processing*, **48**(6), 1537–1547.

Morwal, S., Jahan, N., & Chopra, D. (2012). Named entity recognition using hidden Markov model (HMM). *International Journal on Natural Language Computing* (IJNLC), **1**(4), 15–23.

Mouhcine, R., Mustapha, A., & Zouhir, M. (2018). Recognition of cursive Arabic handwritten text using embedded training based on HMMs. *Journal of Electrical Systems and Information Technology*, **5**(2), 245–251.

Noman, F., Alkawsi, G., Alkahtani, A. A., Al-Shetwi, A. Q., Tiong, S. K., Alalwan, N., et al. (2020). Multistep short-term wind speed prediction using nonlinear auto-regressive neural network with exogenous variable selection. *Alexandria Engineering Journal*, **60**, 1221–1229.

Phillips, P. C., & Perron, P. (1988). Testing for a unit root in time series regression. *Biometrika*, **75**(2), 335–346.

Pinto, C., & Spezia, L. (2015). Markov switching autoregressive models for interpreting vertical movement data with application to an endangered marine apex predator. *Methods in Ecology and Evolution*. https://doi.org/10.1111/2041-210X12494.

Pohle, J., Langrock, R., van Beest, F., & Schmidt, N. (2017). Selecting the number of states in hidden Markov models: Pragmatic solutions illustrated using animal movement. *Journal of Agricultural, Biological and Environmental Statistics*, **22**(3), 270–293.

Psaradakis, Z., & Spagnolo, N. (2003). On the determination of the number of regimes in Markov-switching autoregressive models. *Journal of Time Series Analysis*, **24**(2), 237–252.

Psaradakis, Z., & Spagnolo, N. (2006). Joint determination of the state dimension and autoregressive order for models with Markov regime switching. *Journal of Time Series Analysis*, **27**(5), 753–766.

Ramasso, E. (2014). Investigating computational geometry for failure prognostics. *International Journal of Prognostics and Health Management*, **5**(1), 005.

Ramasso, E. (2016). Segmentation of CMAPSS health indicators into discrete states for sequence-based classification and prediction purposes. Tech. rep., 6839, FEMTO-ST Institute.

Ramasso, E., & Denoeux, T. (2013). Making use of partial knowledge about hidden states in HMMs: An approach based on belief functions. *IEEE Transactions on Fuzzy Systems*, **22**(2), 395–405.

Sateesh Babu, G., Zhao, P., & Li, X. (2016). Deep convolutional neural network based regression approach for estimation of remaining useful life. In *International conference on database systems for advanced applications* (pp. 214–228).

Saxena, A., Goebel, K., Simon, D., & Eklund, N. (2008). Damage propagation modeling for aircraft engine run-to-failure simulation. In *International conference on prognostics and health management* (pp. 1–9).

Scheffer, T., & Wrobel, S. (2001). Active learning of partially hidden Markov models. In *Proceedings of the ECML/PKDD workshop on instance selection*.

Schuller, B., Rigoll, G., & Lang, M. (2003). Hidden Markov model-based speech emotion recognition. In *IEEE International Conference on Multimedia and Expo (ICME)* (pp. 401–404).

Schwarz, G. (1978). Estimating the dimension of a model. *Annals of Statistics*, **6**, 461–464.

Smith, A., Naik, P., & Tsai, C. L. (2006). Markov-switching model selection using Kullback-Leibler divergence. *Journal of Econometrics*, **134**(2), 553–577.

Ubilava, D., & Helmers, C. G. (2013). Forecasting ENSO with a smooth transition autoregressive model. *Environmental Modelling & Software*, **40**, 181–190.

Wang, P., Wang, H., & Yan, R. (2019). Bearing degradation evaluation using improved cross recurrence quantification analysis and nonlinear auto-regressive neural network. *IEEE Access*, **7**, 38937–38946.

Wang, T., Yu, J., Siegel, D., & Lee, J. (2008). A similarity-based prognostics approach for remaining useful life estimation of engineered systems. In *International Conference on Prognostics and Health Management* (pp. 1–6).

Wold, H. (1954). *A study in the analysis of stationary time series* (2nd ed.). Almqvist and Wiksell Book Co.

Wu, Y., Yuan, M., Dong, S., Lin, L., & Liu, Y. (2018). Remaining useful life estimation of engineered systems using vanilla LSTM neural networks. *Neurocomputing*, **275**, 167–179.

Yu, L., Zhou, L., Tan, L., Jiang, H., Wang, Y., Wei, S., & Nie, S. (2014). Application of a new hybrid model with seasonal auto-regressive integrated moving average (ARIMA) and nonlinear auto-regressive neural network (NARNN) in forecasting incidence cases of HFMD in Shenzhen, China. *PLoS ONE*, **9**(6), e98241.

Zhao, S., Zhang, Y., Wang, S., Zhou, B., & Cheng, C. (2019). A recurrent neural network approach for remaining useful life prediction utilizing a novel trend features construction method. *Measurement*, **146**, 279–288.
Michalek, S., Wagner, M., & Timmer, J. (2000). A new approximate likelihood estimator for ARMA-filtered Hidden Markov Models. *IEEE Transactions on Signal Processing*, **48**(6), 1537–1547.

Morwal, S., Jahan, N., & Chopra, D. (2012). Named entity recognition using hidden Markov model (HMM). *International Journal on Natural Language Computing* (IJNLC), **1**(4), 15–23.

Mouhcine, R., Mustapha, A., & Zouhir, M. (2018). Recognition of cursive Arabic handwritten text using embedded training based on HMMs. *Journal of Electrical Systems and Information Technology*, **5**(2), 245–251.

Noman, F., Alkawsi, G., Alkahtani, A. A., Al-Shetwi, A. Q., Tiong, S. K., Alalwan, N., et al. (2020). Multistep short-term wind speed prediction using nonlinear auto-regressive neural network with exogenous variable selection. *Alexandria Engineering Journal*, **60**, 1221–1229.

Phillips, P. C., & Perron, P. (1988). Testing for a unit root in time series regression. *Biometrika*, **75**(2), 335–346.

Pinto, C., & Spezia, L. (2015). Markov switching autoregressive models for interpreting vertical movement data with application to an endangered marine apex predator. *Methods in Ecology and Evolution*. https://doi.org/10.1111/2041-210X12494.

Pohle, J., Langrock, R., van Beest, F., & Schmidt, N. (2017). Selecting the number of states in hidden Markov models: Pragmatic solutions illustrated using animal movement. *Journal of Agricultural, Biological and Environmental Statistics*, **22**(3), 270–293.

Psaradakis, Z., & Spagnolo, N. (2003). On the determination of the number of regimes in Markov-switching autoregressive models. *Journal of Time Series Analysis*, **24**(2), 237–252.

Psaradakis, Z., & Spagnolo, N. (2006). Joint determination of the state dimension and autoregressive order for models with Markov regime switching. *Journal of Time Series Analysis*, **27**(5), 753–766.

Ramasso, E. (2014). Investigating computational geometry for failure prognostics. *International Journal of Prognostics and Health Management*, **5**(1), 005.

Ramasso, E. (2016). Segmentation of CMAPSS health indicators into discrete states for sequence-based classification and prediction purposes. Tech. rep., 6839, FEMTO-ST Institute.

Ramasso, E., & Denoeux, T. (2013). Making use of partial knowledge about hidden states in HMMs: An approach based on belief functions. *IEEE Transactions on Fuzzy Systems*, **22**(2), 395–405.

Sateesh Babu, G., Zhao, P., & Li, X. (2016). Deep convolutional neural network based regression approach for estimation of remaining useful life. In *International conference on database systems for advanced applications* (pp. 214–228).

Saxena, A., Goebel, K., Simon, D., & Eklund, N. (2008). Damage propagation modeling for aircraft engine run-to-failure simulation. In *International conference on prognostics and health management* (pp. 1–9).

Scheffer, T., & Wrobel, S. (2001). Active learning of partially hidden Markov models. In *Proceedings of the ECML/PKDD workshop on instance selection*.

Schuller, B., Rigoll, G., & Lang, M. (2003). Hidden Markov model-based speech emotion recognition. In *IEEE International Conference on Multimedia and Expo (ICME)* (pp. 401–404).

Schwarz, G. (1978). Estimating the dimension of a model. *Annals of Statistics*, **6**, 461–464.

Smith, A., Naik, P., & Tsai, C. L. (2006). Markov-switching model selection using Kullback-Leibler divergence. *Journal of Econometrics*, **134**(2), 553–577.

Ubilava, D., & Helmers, C. G. (2013). Forecasting ENSO with a smooth transition autoregressive model. *Environmental Modelling & Software*, **40**, 181–190.

Wang, P., Wang, H., & Yan, R. (2019). Bearing degradation evaluation using improved cross recurrence quantification analysis and nonlinear auto-regressive neural network. *IEEE Access*, **7**, 38937–38946.

Wang, T., Yu, J., Siegel, D., & Lee, J. (2008). A similarity-based prognostics approach for remaining useful life estimation of engineered systems. In *International Conference on Prognostics and Health Management* (pp. 1–6).

Wold, H. (1954). *A study in the analysis of stationary time series* (2nd ed.). Almqvist and Wiksell Book Co.

Wu, Y., Yuan, M., Dong, S., Lin, L., & Liu, Y. (2018). Remaining useful life estimation of engineered systems using vanilla LSTM neural networks. *Neurocomputing*, **275**, 167–179.

Yu, L., Zhou, L., Tan, L., Jiang, H., Wang, Y., Wei, S., & Nie, S. (2014). Application of a new hybrid model with seasonal auto-regressive integrated moving average (ARIMA) and nonlinear auto-regressive neural network (NARNN) in forecasting incidence cases of HFMD in Shenzhen, China. *PLoS ONE*, **9**(6), e98241.

Zhao, S., Zhang, Y., Wang, S., Zhou, B., & Cheng, C. (2019). A recurrent neural network approach for remaining useful life prediction utilizing a novel trend features construction method. *Measurement*, **146**, 279–288.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.