Quantile hidden semi-Markov models for multivariate time series

Abstract

This paper develops a quantile hidden semi-Markov regression to jointly estimate multiple quantiles for the analysis of multivariate time series. The approach is based upon the Multivariate Asymmetric Laplace (MAL) distribution, which allows to model the quantiles of all univariate conditional distributions of a multivariate response simultaneously, incorporating the correlation structure among the outcomes. Unobserved serial heterogeneity across observations is modeled by introducing regime-dependent parameters that evolve according to a latent finite-state semi-Markov chain. Exploiting the hierarchical representation of the MAL, inference is carried out using an efficient Expectation-Maximization algorithm based on closed form updates for all model parameters, without parametric assumptions about the states’ sojourn distributions. The validity of the proposed methodology is analyzed both by a simulation study and through the empirical analysis of air pollutant concentrations in a small Italian city.

Appendix

Proof of Proposition 1

Firstly, we note that to ensure identifiability of the MAL density in (5) it suffices to apply Proposition 1 of Petrella and Raponi (2019). Secondly, for a general HSMM, identifiability has been proven up to label switching (Leroux 1992). Thus, to ensure identifiability, all one needs to prove is the identifiability of the marginal mixtures (Dannemann et al. 2014), which in our case are represented by the finite mixtures of MAL distributions. Based on the work of Holzmann et al. (2006), Browne and McNicholas (2015) prove identifiability of finite mixtures of multivariate generalized hyperbolic distributions. Since the MAL in (5) is a limiting case of the multivariate generalized hyperbolic distribution (see (3) and (4) in Browne and McNicholas (2015) with $\lambda = 1, \psi = 2$ and $\chi \rightarrow 0$), model identifiability follows by applying Corollary 2 of Browne and McNicholas (2015).

Proof of Proposition 2

The E-step of the EM algorithm considers the conditional expectation of the complete log-likelihood function in (8) given the observed data and the current parameter estimates $\hat{{\varvec{\Phi }}}^{(r-1)}_{\varvec{\tau }}$. At first, we recall that under the constraints imposed on ${\varvec{\tilde{\xi }}}$ and $\varvec{\Lambda }$, the representation in (7) implies that:

$$\begin{aligned} \mathbf{Y}\mid {\tilde{l}} = {\tilde{c}} \sim {{{\mathcal {N}}}}_p(\varvec{\mu } + \mathbf{D}\tilde{{\varvec{\xi }}} {\tilde{c}} \, , \, {\tilde{c}} \mathbf{D}{ {\varvec{\Sigma }}} \mathbf{D}), \qquad {\tilde{C}} \sim \text {Exp} (1). \end{aligned}$$

(35)

This means that the joint density function of $\mathbf{Y}$ and ${\tilde{C}}$ is:

$$\begin{aligned} f_{\mathbf{Y},{\tilde{C}}}(\mathbf{y},{\tilde{c}})= & {} \frac{\exp {\left\{ (\mathbf{y}- \varvec{\mu }) ' \mathbf{D}^{-1}{\varvec{\Sigma }}^{-1} \tilde{{\varvec{\xi }}} \right\} }}{(2 \pi )^{p/2} \mid \mathbf{D}{ {\varvec{\Sigma }}} \mathbf{D}\mid ^{1/2}}\nonumber \\&\left( {\tilde{c}}^{-p/2} \exp {\left\{ -\frac{1}{2} \frac{\tilde{m}}{{\tilde{c}}} - \frac{1}{2} {\tilde{c}} ({\tilde{d}} +2) \right\} }\right) . \end{aligned}$$

(36)

By substituting (36) in (8) and taking the conditional expectation of the logarithm of (8), we obtain the expected complete log-likelihood function in (13).

To compute the conditional expectation of ${\tilde{c}}_{tk}$ and ${\tilde{z}}_{tk}$ in (13), ${\tilde{C}}$ is treated as an additional latent variable. Using the joint distribution of $\mathbf{Y}$ and ${\tilde{C}}$ derived in (36) and the MAL density of $\mathbf{Y}$ given in (5), we have that:

$$\begin{aligned} f_{{\tilde{C}}}({\tilde{C}} \mid \mathbf{Y}= & {} \mathbf{y}) = \frac{f_{{\tilde{C}},\mathbf{Y}}({\tilde{c}},\mathbf{y})}{f_{\mathbf{Y}}(\mathbf{y})} \nonumber \\= & {} \frac{{\tilde{c}}^{-p/2} \left( \frac{2+{\tilde{d}}}{{\tilde{m}}} \right) ^{\nu /2} \exp {\left\{ -\frac{{\tilde{m}}}{2 {\tilde{c}}}- \frac{{\tilde{c}} (2+{\tilde{d}})}{2} \right\} }}{2 K_{\nu }\left( \sqrt{(2+{\tilde{d}}){\tilde{m}}} \right) }, \end{aligned}$$

(37)

which corresponds to a Generalized Inverse Gaussian (GIG) distribution with parameters $\nu , {2+{\tilde{d}}}, \tilde{m_i}$, i.e.¹

$$\begin{aligned} f_{{\tilde{C}}}({\tilde{C}} \mid \mathbf{Y}= \mathbf{y}) \sim \text{ GIG }\left( \nu , {\tilde{d}}+2, \tilde{m}\right) . \end{aligned}$$

(38)

Then, it follows that

$$\begin{aligned} {\mathbb {E}}[{\tilde{C}} \mid \cdot ] = \left( \frac{\hat{\tilde{m}}}{2+\hat{{\tilde{d}}}} \right) ^{\frac{1}{2}} \frac{K_{\nu +1}\left( \sqrt{(2+\hat{{\tilde{d}}})\hat{\tilde{m}}} \right) }{K_{\nu }\left( \sqrt{(2+\hat{{\tilde{d}}}) \hat{\tilde{m}}}\right) } \end{aligned}$$

(39)

and

$$\begin{aligned} {\mathbb {E}}[{\tilde{C}}^{-1} \mid \cdot ] = \left( \frac{2+\hat{\tilde{d}}}{\hat{\tilde{m}}} \right) ^{\frac{1}{2}} \frac{K_{\nu +1} \left( \sqrt{(2+\hat{{\tilde{d}}})\hat{\tilde{m}}} \right) }{K_{\nu } \left( \sqrt{(2+\hat{{\tilde{d}}})\hat{\tilde{m}}} \right) } - \frac{2 \nu }{\hat{\tilde{m}}}. \end{aligned}$$

(40)

Denoting the two conditional expectations in (39) and (40) by $\hat{{\tilde{c}}}$ and $\hat{{\tilde{z}}}$ respectively, concludes the proof.

Proof of Proposition 3

Imposing the first order conditions on (13) with respect to each component of the set ${\varvec{\Phi }}_{\varvec{\tau }}$, gives the update estimates in (16), (17), (18) and (19). However, there is not closed formula solution to update the elements of the scale matrix $\mathbf{D}_j$; hence, the M-step update requires using numerical optimization techniques to maximize (13). A considerable disadvantage of this procedure is the necessary high computational effort which could be very time-consuming. For this reason, we utilize a simpler estimator for the scale parameters $\delta _{jk}, k = 1,\dots ,p$, which follows directly from the fact that all marginals of the MAL distribution are univariate AL distributions (see Yu and Zhang 2005):

$$\begin{aligned} \hat{\delta }_{jk} = \frac{1}{T} \sum _{t=1}^{T} {\hat{\gamma }}_{tj} \rho _\tau (Y_{t}^{(k)} - \hat{{\mu }}_{tjk}), \end{aligned}$$

(41)

where $\hat{{\mu }}_{tjk}$ is the k-th element of the vector $\hat{\varvec{\mu }}_{tj}$.