∎

11institutetext: Salvatore D. Tomarchio 22institutetext: Dipartimento di
Economia e Impresa, Università degli Studi di Catania, Catania, Italia
22email<EMAIL_ADDRESS>33institutetext: Antonio Punzo
44institutetext: Dipartimento di Economia e Impresa, Università degli Studi di
Catania, Catania, Italia
55institutetext: Antonello Maruotti 66institutetext: Dipartimento di
Giurisprudenza, Economia, Politica e Lingue Moderne, Libera Università Maria
Ss. Assunta, Roma, Italia
Department of Mathematics, University of Bergen, Bergen, Norway

#  Parsimonious Hidden Markov Models for Matrix-Variate Longitudinal Data

Salvatore D. Tomarchio
Antonio Punzo
Antonello Maruotti

###### Abstract

Hidden Markov models (HMMs) have been extensively used in the univariate and
multivariate literature. However, there has been an increased interest in the
analysis of matrix-variate data over the recent years. In this manuscript we
introduce HMMs for matrix-variate longitudinal data, by assuming a matrix
normal distribution in each hidden state. Such data are arranged in a four-way
array. To address for possible overparameterization issues, we consider the
spectral decomposition of the covariance matrices, leading to a total of 98
HMMs. An expectation-conditional maximization algorithm is discussed for
parameter estimation. The proposed models are firstly investigated on
simulated data, in terms of parameter recovery, computational times and model
selection. Then, they are fitted to a four-way real data set concerning the
unemployment rates of the Italian provinces, evaluated by gender and age
classes, over the last 16 years.

###### Keywords:

Hidden Markov models Matrix-variate Clustering Parsimonious models

##  1 Introduction

Multivariate longitudinal data have been widely analyzed in the literature
(Verbeke et al., 2014 and Verdam and Oort, 2019). Data are usually presented
in the standard three-way format, where units, times and variables are
arranged in software-ready manners. Because of their three-way structure,
multivariate longitudinal data have been recently arranged in a matrix-variate
fashion (Huang et al., 2019 and Viroli, 2011b): for each unit $i=1,\dots,I$,
we observe a $P\times T$ matrix, where $P$ and $T$ denote the number of
variables and times, respectively. Then, such data have been used for model-
based clustering via matrix-variate mixture models (see e.g. Wang and
Melnykov, 2020; Melnykov and Zhu, 2019; Tomarchio et al., 2020a, b). This
allows for both clustering units in homogeneous groups, defined according to
similarities between matrix-variate data, and separately modeling the
association between variables and times. Unfortunately, this procedure has two
side effects:

  1. (a)

using the time on either the rows or the columns of the matrices reduces the
types of longitudinal data structures that can be arranged in a matrix-variate
framework. For instance, spatio-temporal data are used to either to analyze
$P$ variables observed at $T$ times for $R$ different locations (Viroli,
2011b) or to evaluate one measurement on $R$ locations at $T$ times on a set
of $I$ units (Viroli, 2011a). However, it is not possible to jointly consider
$P$ variables at $R$ locations for $T$ times on $I$ units. A possible solution
could be to combine locations-times in one $RT$-dimension, as done by Viroli
(2011a), but this implies a loss in terms of interpretability as well as an
increase in the number of parameters of the models, given the higher
dimensionality of the matrices.

Another example consists of two-factor data, which have been commonly
considered in longitudinal settings (see e.g. Brunner and Puri, 2001,
Fitzmaurice and Ravichandran, 2008, Noguchi et al., 2012). Such data, have
been recently used in matrix-variate mixture models by (Sarkar et al., 2020)
in a not-longitudinal way, given that the factors fill the two dimensions of
the $P\times R$ matrices for the $I$ units, and an additional dimension for
the time is required. In other terms, in both examples it would be necessary
to move from three-way to four-way data structures.

  2. (b)

Currently available clustering approaches for matrix-variate data assume time-
constant clustering, i.e. it is not possible for the sample units to move
across clusters over time and the evolution over time of the clustering
structure is completely overlooked. Zhu and Melnykov (2021) recognize the
importance of allowing for time-dependence in a matrix-variate clustering
setting and propose a procedure capable of capturing the heterogeneity pattern
and estimating change points from all data groups simultaneously. Indeed,
time-varying heterogeneity is a specific important feature of longitudinal
data analysis, and as such appropriate modeling strategies should be
considered. Hidden Markov models (HMMs) have been extensively used to address
this longitudinal data peculiarity (Zucchini et al., 2017, Bartolucci et al.,
2012, Maruotti, 2011 and Altman, 2007). Being (dependent) mixtures, HMMs
simultaneously allow for clustering units and for modeling the evolution of
the clustering over time.

To jointly consider the aspects in (a) and (b), in this manuscript we
introduce and discuss HMMs for matrix-variate longitudinal data, with a
specific application on the two-factor longitudinal case. Such kind of data
can be arranged in a four-way array of dimension $P\times R\times I\times T$.
Under this four-way longitudinal data setting, HMMs simultaneously
investigates several further (null) hypotheses, beyond the recovering of the
time-varying clustering:

  * •

all means at different time points are the same;

  * •

there is an association between the factor levels.

The association is often the primary interest when researchers aim to study
whether levels differentially affect outcomes in different hidden states.

A side effect of working with four-way data is the large number of parameters
involved. This often occurs because of the (row- and column-specific)
covariance matrices, since $P(P+1)/2$ and $R(R+1)/2$ unique parameters must be
estimated. One of the most classical ways of addressing this
overparameterization issue involves the spectral decomposition of the
covariance matrices (Celeux and Govaert, 1995). This decomposition offers
remarkable flexibility and a geometric interpretation in terms of volume,
shape, and orientation of the hidden states (for other approaches available in
the HMMs literature see Maruotti et al., 2017 and Farcomeni et al., 2020). By
using the spectral decomposition of the covariance matrices, we obtain a
family of 98 parsimonious HMMs that will be described in Section 2.2, after
the presentation of the general model (Section 2.1). In this framework, model
parameters can be estimated by a full maximum likelihood method based on the
Expectation Conditional Maximization (ECM) algorithm (Meng and Rubin, 1993),
and recursions widely used in the HMM literature (Baum et al., 1970). An
iterative Minorization–Maximization (MM) algorithm (Browne and McNicholas,
2014) is also adopted to estimate a subset of the parsimonious HMMs.

We illustrate the proposal by a large-scale simulation study in order to
investigate the empirical behavior of the proposed approach with respect to
several aspects, such as the number of observed units and times, the number of
hidden states and the association structure between factor-levels. We focus on
goodness of clustering and parameters recovery, with a focus on computational
times and model selection procedures. Furthermore, we test the proposal by
analyzing a sample taken from the Italian National Institute of Statistics on
the unemployment rate in 98 Italian provinces recorded for 16 years, also
covering the 2008 crisis. We examine the unemployment rate arranged as a two-
factor design, i.e. taking into account gender and age classes, by allowing
some dynamics in the evolution of unemployment. We obtain a flexible model by
including different associations across levels, changing according to the
inferred dynamics, and by accounting for unobserved characteristics
influencing changes in the province’s unemployment patterns.

Finally, Section 5 summarizes the key aspects of our proposal along with
future possible extensions.

##  2 Methodology

###  2.1 The model

Let $\left\\{\mathcal{X}_{it};i=1,\ldots,I,t=1,\ldots,T\right\\}$ be a
sequence of matrix-variate longitudinal observations recorded on $I$ units
over $T$ times, with $\mathcal{X}_{it}\in\mathbb{R}^{R\times P}$, and let
$\left\\{S_{it};i=1,\ldots,I,t=1,\ldots,T\right\\}$ be a first-order Markov
chain defined on the state space $\left\\{1,\ldots,k,\ldots,K\right\\}$. As
mentioned in Section 1, a HMM is a particular type of dependent mixture model
consisting of two parts: an underlying unobserved process
$\left\\{S_{it}\right\\}$ that satisfies the Markov property, i.e.

$\text{Pr}\left(S_{it}=s_{it}|S_{i1}=s_{i1},\ldots,S_{it-1}=s_{it-1}\right)=\text{Pr}\left(S_{it}=s_{it}|S_{it-1}=s_{it-1}\right),$

and a state-dependent observation process $\left\\{\mathcal{X}_{it}\right\\}$
for which the conditional independence property holds, i.e.

$\displaystyle
f\left(\mathcal{X}_{it}=\mathbf{X}_{it}|\mathcal{X}_{i1}=\mathbf{X}_{i1},\ldots,\mathcal{X}_{it-1}=\mathbf{X}_{it-1},S_{i1}=s_{i1}\ldots,S_{it}=s_{it}\right)$
$\displaystyle=f\left(\mathcal{X}_{it}=\mathbf{X}_{it}|S_{it}=s_{it}\right),$

where $f(\cdot)$ is a generic probability density function (pdf). Therefore,
the unknown parameters in an HMM involve both the parameters of the Markov
chain and those of the state-dependent pdfs. In detail, the parameters of the
Markov chain are the initial probabilities
$\pi_{ik}=\text{Pr}\left(S_{i1}=k\right)$, $k=1,\ldots,K$, being $K$ the
number of states, and the transition probabilities

$\pi_{ik|j}=\text{Pr}\left(S_{it}=k|S_{it-1}=j\right),\hskip
5.0ptt=2,\ldots,T\hskip 5.0pt\text{and}\hskip 5.0ptj,k=1,\ldots,K,$

where $k$ refers to the current state and $j$ refers to the one previously
visited. To simplify the discussion, we will consider homogeneous HMMs, that
is $\pi_{ik|j}=\pi_{k|j}$ and $\pi_{ik}=\pi_{k},i=1,\ldots,I$. We collect the
initial probabilities in the $K$-dimensional vector $\bm{\pi}$, whereas the
time-homogenous transition probabilities are inserted in the $K\times K$
transition matrix $\mathbf{\Pi}$.

Regarding the conditional density for the observed process, it will be given
by a matrix-normal distribution, i.e.

$\phi\left(\mathbf{X}_{it}|S_{it}=k;\bm{\theta}_{k}\right)=\frac{\exp\left\\{-\frac{1}{2}\,\mbox{tr}\left[\mathbf{\Sigma}_{k}^{-1}(\mathbf{X}-\mathbf{M}_{k})\mathbf{\Psi}_{k}^{-1}(\mathbf{X}-\mathbf{M}_{k})^{\prime}\right]\right\\}}{(2\pi)^{\frac{PR}{2}}|\mathbf{\Sigma}_{k}|^{\frac{R}{2}}|\mathbf{\Psi}_{k}|^{\frac{P}{2}}},$
(1)

where $\mathbf{M}_{k}$ is the $P\times R$ matrix of means,
$\mathbf{\Sigma}_{k}$ is the $P\times P$ covariance matrix containing the
variances and covariances between the P rows, $\mathbf{\Psi}_{k}$ is the
$R\times R$ covariance matrix containing the variance and covariances of the
$R$ columns and
$\bm{\theta}_{k}=\left\\{\mathbf{M}_{k},\mathbf{\Sigma}_{k},\mathbf{\Psi}_{k}\right\\}$.
For an exhaustive description of the matrix-normal distribution and its
properties see Gupta and Nagar (2018).

###  2.2 Parsimonious models

As introduced in Section 1, a way to reduce the number of parameters of the
model is to introduce parsimony in the covariance matrices via the well-known
eigen decomposition. Specifically, a $Q\times Q$ covariance matrix can be
decomposed as

$\mathbf{\Phi}_{k}=\lambda_{k}\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\prime}$
(2)

where $\lambda_{k}=|\mathbf{\Phi}_{k}|^{1/Q}$, $\mathbf{\Gamma}_{k}$ is a
$Q\times Q$ orthogonal matrix of the eigenvectors of $\mathbf{\Phi}_{k}$ and
$\bm{\Delta}_{k}$ is a diagonal matrix with the eigenvalues of
$\mathbf{\Phi}_{k}$ located on the main diagonal. From a geometric point of
view, $\lambda_{k}$ determines the volume, $\mathbf{\Gamma}_{k}$ indicates the
orientation, and $\bm{\Delta}_{k}$ denotes the shape of the $k$th state. By
imposing constraints on the three components of (2), the fourteen parsimonious
models of Table 1 are obtained.

Table 1: Nomenclature, covariance matrix structure, and number of free
parameters in $\mathbf{\Phi}_{1},\ldots,\mathbf{\Phi}_{K}$ for the
parsimonious models obtained via the eigen decomposition of the state
covariance matrices.

Family | Model | Type | Volume | Shape | Orientation | # of free parameters
---|---|---|---|---|---|---
|  |  |  |  |  | in $\mathbf{\Phi}_{1},\ldots,\mathbf{\Phi}_{K}$
Spherical | EII | $\lambda\bm{I}$ | Equal | Spherical | - | 1
Spherical | VII | $\lambda_{k}\bm{I}$ | Variable | Spherical | - | $K$
Diagonal | EEI | $\lambda\bm{\Delta}$ | Equal | Equal | Axis-Aligned | $Q$
Diagonal | VEI | $\lambda_{k}\bm{\Delta}$ | Variable | Equal | Axis-Aligned | $K+Q-1$
Diagonal | EVI | $\lambda\bm{\Delta}_{k}$ | Equal | Variable | Axis-Aligned | $K(Q-1)+1$
Diagonal | VVI | $\lambda_{k}\bm{\Delta}_{k}$ | Variable | Variable | Axis-Aligned | $KQ$
General | EEE | $\lambda\mathbf{\Gamma}\bm{\Delta}\mathbf{\Gamma}^{\top}$ | Equal | Equal | Equal | $Q(Q+1)/2$
General | VEE | $\lambda_{k}\mathbf{\Gamma}\bm{\Delta}\mathbf{\Gamma}^{\top}$ | Variable | Equal | Equal | $Q(Q+1)/2+K-1$
General | EVE | $\lambda\mathbf{\Gamma}\bm{\Delta}_{k}\mathbf{\Gamma}^{\top}$ | Equal | Variable | Equal | $Q(Q-1)/2+K(Q-1)+1$
General | VVE | $\lambda_{k}\mathbf{\Gamma}\bm{\Delta}_{k}\mathbf{\Gamma}^{\top}$ | Variable | Variable | Equal | $Q(Q-1)/2+KQ$
General | EEV | $\lambda\mathbf{\Gamma}_{k}\bm{\Delta}\mathbf{\Gamma}_{k}^{\top}$ | Equal | Equal | Variable | $KQ(Q-1)/2+Q$
General | VEV | $\lambda_{k}\mathbf{\Gamma}_{k}\bm{\Delta}\mathbf{\Gamma}_{k}^{\top}$ | Variable | Equal | Variable | $KQ(Q-1)/2+K+Q-1$
General | EVV | $\lambda\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$ | Equal | Variable | Variable | $KQ(Q+1)/2-K+1$
General | VVV | $\lambda_{k}\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$ | Variable | Variable | Variable | $KQ(Q+1)/2$

Considering that we have two covariance matrices in (1), this would yield to
$14\times 14=196$ parsimonious HMMs. However, there is a non-identifiability
issue since
$\mathbf{\Psi}\otimes\mathbf{\Sigma}=\mathbf{\Psi}^{*}\otimes\mathbf{\Sigma}^{*}$
if $\mathbf{\Sigma}^{*}=a\mathbf{\Sigma}$ and
$\mathbf{\Psi}^{*}=a^{-1}\mathbf{\Psi}$. As a result, $\mathbf{\Sigma}$ and
$\mathbf{\Psi}$ are identifiable up to a multiplicative constant $a$ (Sarkar
et al., 2020). To avoid such issue, the column covariance matrix
$\mathbf{\Psi}$ is restricted to have $|\mathbf{\Psi}|=1$, implying that in
(2) the parameter $\lambda_{k}$ is unnecessary. This reduces the number of
models related to $\mathbf{\Psi}$ from 14 to 7, i.e.,
$\bm{I},\bm{\Delta},\bm{\Delta}_{k},\mathbf{\Gamma}\bm{\Delta}\mathbf{\Gamma}^{\top},\mathbf{\Gamma}\bm{\Delta}_{k}\mathbf{\Gamma}^{\top},\mathbf{\Gamma}_{k}\bm{\Delta}\mathbf{\Gamma}_{k}^{\top},\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$.
Therefore, we obtain $14\times 7=98$ parsimonious HMMs.

###  2.3 Maximum likelihood estimation

To fit our HMMs, we use the expectation-conditional maximization (ECM)
algorithm (Meng and Rubin, 1993). The ECM algorithm is a variant of the
classical expectation-maximization (EM) algorithm (Dempster et al., 1977),
from which it differs since the M-step is replaced by a sequence of simpler
and computationally convenient CM-steps.

Let $S=\left\\{\mathbf{X}_{it};i=1,\ldots,I,t=1,\ldots,T\right\\}$ be a sample
of matrix-variate longitudinal observations. Then, the incomplete-data
likelihood function is

$L\left(\bm{\Theta}|S\right)=\prod_{i=1}^{I}\bm{\pi}^{\prime}\bm{\phi}\left(\mathbf{X}_{i1}\right)\mathbf{\Pi}\bm{\phi}\left(\mathbf{X}_{i2}\right)\mathbf{\Pi}\ldots\bm{\phi}\left(\mathbf{X}_{iT-1}\right)\mathbf{\Pi}\bm{\phi}\left(\mathbf{X}_{iT}\right)\bm{1}_{K},$

where $\bm{\phi}\left(\mathbf{X}_{it}\right)$ is a $K\times K$ diagonal matrix
with conditional densities $\mathcal{X}_{it}=\mathbf{X}_{it}|S_{it}=k$ on the
main diagonal, $\bm{1}_{K}$ is a vector $K$ ones and $\bm{\Theta}$ contains
all the model parameters. In this setting, $\mathcal{S}$ is viewed as
incomplete because, for each observation, we do not know its state membership
and its evolution over time. For this reason, let us define the unobserved
state membership
$\bm{z}_{it}=\left(z_{it1},\ldots,z_{itk},\ldots,z_{itK}\right)^{\prime}$ and
the unobserved states transition

$\bm{z}\bm{z}_{it}=\begin{bmatrix}zz_{it11}&\ldots&zz_{it1k}&\ldots&zz_{it1K}\\\
\vdots&&\vdots&&\vdots\\\ zz_{itj1}&\ldots&zz_{itjk}&\ldots&zz_{itjK}\\\
\vdots&&\vdots&&\vdots\\\
zz_{itK1}&\ldots&zz_{itKk}&\ldots&zz_{itKK}\end{bmatrix}$

where

$z_{itk}=\begin{cases}1&\quad\text{if}\hskip 5.0ptS_{t}=k\\\
0&\quad\text{otherwise}\end{cases}\text{and}\hskip
5.0ptzz_{itjk}=\begin{cases}1&\quad\text{if}\hskip 5.0ptS_{it-1}=j\hskip
5.0pt\text{and}\hskip 5.0ptS_{it}=k\\\ 0&\quad\text{otherwise}\end{cases}.$

Therefore, the complete data are
$\mathcal{S}_{c}=\left\\{\mathbf{X}_{it},\bm{z}_{it},\bm{z}\bm{z}_{it};i=1,\ldots,I,t=1,\ldots,T\right\\}$
and the corresponding complete-data log-likelihood is

$l_{c}\left(\bm{\Theta}|\mathcal{S}_{c}\right)=l_{c_{1}}\left(\bm{\pi}|\mathcal{S}_{c}\right)+l_{c_{2}}\left(\mathbf{\Pi}|\mathcal{S}_{c}\right)+l_{c_{3}}\left(\bm{\theta}|\mathcal{S}_{c}\right),$
(3)

with $\bm{\theta}=\left\\{\bm{\theta}_{k};k=1,\ldots,K\right\\}$ and

$\displaystyle l_{c_{1}}\left(\bm{\pi}|\mathcal{S}_{c}\right)$
$\displaystyle=\sum\limits_{i=1}^{I}\sum\limits_{k=1}^{K}z_{i1k}\log\left(\pi_{k}\right)$
$\displaystyle l_{c_{2}}\left(\mathbf{\Pi}|\mathcal{S}_{c}\right)$
$\displaystyle=\sum\limits_{i=1}^{I}\sum\limits_{t=2}^{T}\sum\limits_{k=1}^{K}\sum\limits_{j=1}^{K}zz_{itjk}\log\left(\pi_{k|j}\right)$
$\displaystyle l_{c_{3}}\left(\bm{\theta}|\mathcal{S}_{c}\right)$
$\displaystyle=\sum\limits_{i=1}^{I}\sum\limits_{t=1}^{T}\sum\limits_{k=1}^{K}z_{itk}\Bigg{[}\Bigg{.}-\frac{PR}{2}\log\left(2\pi\right)-\frac{R}{2}\log|\mathbf{\Sigma}_{k}|-\frac{P}{2}\log|\mathbf{\Psi}_{k}|$
$\displaystyle-\frac{1}{2}\,\mbox{tr}\left[\mathbf{\Sigma}_{k}^{-1}(\mathbf{X}-\mathbf{M}_{k})\mathbf{\Psi}_{k}^{-1}(\mathbf{X}-\mathbf{M}_{k})^{\prime}\right]\Bigg{.}\Bigg{]}.$

In the following, by adopting the notation used in Melnykov and Zhu (2019),
the parameters marked with one dot represent the updates at the previous
iteration and those marked with two dots are the updates at the current
iteration.

##### E-Step

The E-step requires calculation of the conditional expectation of (3), given
$\mathcal{S}_{c}$ and the current estimates of $\dot{\bm{\Theta}}$. Therefore,
we need to replace $z_{itk}$ and $z_{itjk}$ with their conditional
expectations, namely, $\ddot{z}_{itk}$ and $\ddot{zz}_{itjk}$. This can be
efficiently done by exploiting a forward recursion approach (Baum et al.,
1970; Zucchini et al., 2017).

Let us start by defining the forward probability

$\gamma_{itk}=\text{Pr}\left(\mathcal{X}_{i1}=\mathbf{X}_{i1},\ldots,\mathcal{X}_{it}=\mathbf{X}_{it},S_{it}=k\right),$

that is the probability of seeing the partial sequence finishing up in state
$k$ at time $t$, and the corresponding backward probability

$\beta_{itk}=\text{Pr}\left(\mathcal{X}_{it+1}=\mathbf{X}_{it+1},\ldots,\mathcal{X}_{iT}=\mathbf{X}_{iT}|S_{it}=k\right).$

It is known that the computation of the forward and backward probabilities is
susceptible to numerical overflow errors (Zucchini et al., 2017). To prevent,
or at least to decrease, the risk of such errors a scaling procedure is used.
Specifically, it will be convenient to work on the log-scale (Farcomeni,
2012). In order to do so, we use a simple computational device based on the
following equality:

$\log\left(a+b\right)=\log(a)+\log\left(1+\exp\left(\log\left(b\right)-\log\left(a\right)\right)\right).$

Therefore, if one has only the log of two quantities $\log\left(a\right)$ and
$\log\left(b\right)$, only their difference must be exponentiated to obtain
$\log\left(a+b\right)$, reducing the risks of underflow. By iterating this
reasoning, one can sum a vector of quantities on the log-scale. This operation
is called $\bigoplus$. Thus, when $t=1$, the forward recursion on the log-
scale is given by

$\log\left(\gamma_{i1k}\right)=\log\left[\phi\left(\mathbf{X}_{i1}|S_{i1}=k\right)\right]+\log\left(\pi_{k}\right)$

whereas, for $t=2,\ldots,T$, it is

$\log\left(\gamma_{itk}\right)=\log\left[\phi\left(\mathbf{X}_{it}|S_{it}=k\right)\right]+\log\left(\pi_{k}|j\right)+\bigoplus_{j=1}^{K}\log\left(\gamma_{it-1,j}\right).$

In a similar way, for the backward recursion on the log-scale is

$\log\left(\beta_{iTk}=0\right),$

and, for $t=T-1,\ldots,1$, we obtain

$\log\left(\beta_{itj}\right)=\bigoplus_{k=1}^{K}\log\left[\phi\left(\mathbf{X}_{it+1}|S_{it+1}=k\right)\right]+\log\left[\beta_{it+1,k}+\log\left(\pi_{k|j}\right)\right].$

Then, the updates required in the E-step can be computed as

$\ddot{z}_{itk}=\frac{\gamma_{itk}\beta_{itk}}{\sum\limits_{h=1}^{K}\gamma_{ith}\beta_{ith}},\quad\ddot{zz}_{itjk}=\frac{\gamma_{i\left(t-1\right)j}\pi_{k|j}\phi\left(\mathbf{X}_{it}|S_{it}=k\right)\beta_{itk}}{\sum\limits_{h=1}^{K}\gamma_{iTh}}.$

##### CM-Step 1

Consider $\bm{\Theta}=\left\\{\bm{\Theta}_{1},\bm{\Theta}_{2}\right\\}$, where
$\bm{\Theta}_{1}=\left\\{\pi_{k},\mathbf{\Pi},\mathbf{M}_{k},\mathbf{\Sigma}_{k}\right\\}$
and $\bm{\Theta}_{2}=\left\\{\mathbf{\Psi}_{k}\right\\}$. At the first CM-
step, we maximize the expectation of the complete-data log-likelihood with
respect to $\bm{\Theta}_{1}$, fixing $\bm{\Theta}_{2}$ at
$\dot{\bm{\Theta}_{2}}$. In particular, we obtain

$\ddot{\pi}_{k}=\frac{\sum_{i=1}^{I}\ddot{z}_{i1k}}{I},\quad\ddot{\pi}_{k|j}=\frac{\sum_{i=1}^{I}\sum_{t=2}^{T}\ddot{zz}_{itjk}}{\sum_{i=1}^{I}\sum_{t=2}^{T}\sum_{k=1}^{K}\ddot{zz}_{itjk}},$
$\ddot{\mathbf{M}}_{k}=\frac{\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}\mathbf{X}_{it}}{\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

The update for $\mathbf{\Sigma}_{k}$ depends on the parsimonious structure
considered. For notational simplicity, let
$\ddot{\mathbf{Y}}=\sum_{k=1}^{K}\ddot{\mathbf{Y}}_{k}$ be the update of the
within state row scatter matrix, where
$\ddot{\mathbf{Y}}_{k}=\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}\left(\mathbf{X}_{it}-\ddot{\mathbf{M}}_{k}\right)\dot{\mathbf{\Psi}}_{k}^{-1}\left(\mathbf{X}_{it}-\ddot{\mathbf{M}}_{k}\right)^{\prime}$
is the update of the row scatter matrix related to the $k$th state. The
updates for the 14 parsimonious structures of $\mathbf{\Sigma}_{k}$ are:

  * •

Model EII [$\mathbf{\Sigma}_{k}=\lambda\bm{I}$] In this setting, the row
covariance matrices of all states are spherical and have equal volume. We need
to estimate only $\lambda$ as

$\ddot{\lambda}=\frac{\,\mbox{tr}\left\\{\ddot{\mathbf{Y}}\right\\}}{PRTI}.$

  * •

Model VII [$\mathbf{\Sigma}_{k}=\lambda_{k}\bm{I}$] In this case, the row
covariance matrices are spherical but their volume is different. Thus, the
update for $\lambda_{k}$ is

$\ddot{\lambda}_{k}=\frac{\,\mbox{tr}\left\\{\ddot{\mathbf{Y}}_{k}\right\\}}{PR\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

  * •

Model EEI [$\mathbf{\Sigma}_{k}=\lambda\bm{\Delta}$] Here, the row covariance
matrices of all states have equal volume, shape and are axis-aligned. The
updates for $\bm{\Delta}$ and $\lambda$ are

$\ddot{\bm{\Delta}}=\frac{\text{diag}\left(\ddot{\mathbf{Y}}\right)}{\left|\text{diag}\left(\ddot{\mathbf{Y}}\right)\right|^{\frac{1}{P}}},\quad\ddot{\lambda}=\frac{\left|\text{diag}\left(\ddot{\mathbf{Y}}\right)\right|^{\frac{1}{P}}}{RTI}.$

  * •

Model VEI [$\mathbf{\Sigma}_{k}=\lambda_{k}\bm{\Delta}$] In this setting, the
row covariance matrices of all states have equal shape and are axis-aligned,
but they are allowed to have different volumes. We need to update
$\bm{\Delta}$ and $\lambda$ as

$\ddot{\bm{\Delta}}=\frac{\text{diag}\left(\sum\limits_{k=1}^{K}\dot{\lambda}_{k}^{-1}\ddot{\mathbf{Y}}_{k}\right)}{\left|\text{diag}\left(\sum\limits_{k=1}^{K}\dot{\lambda}_{k}^{-1}\ddot{\mathbf{Y}}_{k}\right)\right|^{\frac{1}{P}}},\quad\ddot{\lambda}_{k}=\frac{\,\mbox{tr}\left\\{\ddot{\mathbf{Y}}_{k}\ddot{\bm{\Delta}}^{-1}\right\\}}{PR\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

  * •

Model EVI [$\mathbf{\Sigma}_{k}=\lambda\bm{\Delta}_{k}$] In this case, the row
covariance matrices have equal volume and are axis-aligned, but they have
different shapes. The updates for $\bm{\Delta}_{k}$ and $\lambda$ are

$\ddot{\bm{\Delta}}_{k}=\frac{\text{diag}\left(\ddot{\mathbf{Y}}_{k}\right)}{\left|\text{diag}\left(\ddot{\mathbf{Y}}_{k}\right)\right|^{\frac{1}{P}}},\quad\ddot{\lambda}=\frac{\sum\limits_{k=1}^{K}\left|\text{diag}\left(\ddot{\mathbf{Y}}_{k}\right)\right|^{\frac{1}{P}}}{RTI}.$

  * •

Model VVI [$\mathbf{\Sigma}_{k}=\lambda_{k}\bm{\Delta}_{k}$] The most generic
case in the diagonal family has row covariance matrices with different volume
and shape, in addition to being axis-aligned. Then, the updates for
$\bm{\Delta}_{k}$ and $\lambda_{k}$ are

$\ddot{\bm{\Delta}}_{k}=\frac{\text{diag}\left(\ddot{\mathbf{Y}}_{k}\right)}{\left|\text{diag}\left(\ddot{\mathbf{Y}}_{k}\right)\right|^{\frac{1}{P}}},\quad\ddot{\lambda}_{k}=\frac{\left|\text{diag}\left(\ddot{\mathbf{Y}}_{k}\right)\right|^{\frac{1}{P}}}{D\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

  * •

Model EEE
[$\mathbf{\Sigma}_{k}=\lambda\mathbf{\Gamma}\bm{\Delta}\mathbf{\Gamma}^{\top}$]
The most constrained member of the general family has row covariance matrices
with same volume, shape and orientation. Thus, the update for
$\mathbf{\Sigma}$ is given by

$\ddot{\mathbf{\Sigma}}=\frac{\ddot{\mathbf{Y}}}{RTI}.$

  * •

Model VEE
[$\mathbf{\Sigma}_{k}=\lambda_{k}\mathbf{\Gamma}\bm{\Delta}\mathbf{\Gamma}^{\top}$]
This model assumes row covariances matrices with same shape and orientation
but different volumes. Let
$\mathbf{C}=\ddot{\mathbf{\Gamma}}\ddot{\bm{\Delta}}\ddot{\mathbf{\Gamma}}^{\top}$.
The updates for $\mathbf{C}$ and $\lambda_{k}$ are

$\ddot{\mathbf{C}}=\frac{\sum\limits_{k=1}^{K}\dot{\lambda}_{k}^{-1}\ddot{\mathbf{Y}}_{k}}{\left|\sum\limits_{k=1}^{K}\dot{\lambda}_{k}^{-1}\ddot{\mathbf{Y}}_{k}\right|^{\frac{1}{P}}},\quad\ddot{\lambda}_{k}=\frac{\,\mbox{tr}\left\\{\ddot{\mathbf{C}}^{-1}\ddot{\mathbf{Y}}_{k}\right\\}}{PR\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

  * •

Model EVE
[$\mathbf{\Sigma}_{k}=\lambda\mathbf{\Gamma}\bm{\Delta}_{k}\mathbf{\Gamma}^{\top}$]
Here, the row covariance matrices have equal volume and orientation but
different shapes. Given that there is no analytical solution for
$\mathbf{\Gamma}$, while keeping fixed the other parameters, an iterative
Minorization–Maximization (MM) algorithm (Browne and McNicholas, 2014) is
employed. In detail, a surrogate function can be constructed as

$f\left(\mathbf{\Gamma}\right)=\sum\limits_{k=1}^{K}\,\mbox{tr}\left\\{\mathbf{Y}_{k}\mathbf{\Gamma}\bm{\Delta}_{k}^{-1}\mathbf{\Gamma}^{\top}\right\\}\leq
S+\,\mbox{tr}\left\\{\bm{F}\mathbf{\Gamma}\right\\},$

where $S$ is a constant and
$\bm{F}=\sum_{k=1}^{K}\left(\bm{\Delta}_{k}^{-1}\dot{\mathbf{\Gamma}}^{\top}\mathbf{Y}_{k}-e_{k}\bm{\Delta}_{k}^{-1}\dot{\mathbf{\Gamma}}^{\top}\right)$,
with $e_{k}$ being the largest eigenvalue of $\mathbf{Y}_{k}$. The update of
$\mathbf{\Gamma}$ is given by
$\ddot{\mathbf{\Gamma}}=\dot{\bm{G}}\dot{\bm{H}}^{\top}$, where $\dot{\bm{G}}$
and $\dot{\bm{H}}$ are obtained from the singular value decomposition of
$\bm{F}$. This process is repeated until a specified convergence criterion is
met and the estimate $\ddot{\mathbf{\Gamma}}$ is obtained from the last
iteration. Then, we obtain the update for $\bm{\Delta}_{k}$ and $\lambda$ as

$\ddot{\bm{\Delta}}_{k}=\frac{\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{Y}}_{k}\ddot{\mathbf{\Gamma}}\right)}{\left|\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{Y}}_{k}\ddot{\mathbf{\Gamma}}\right)\right|^{\frac{1}{P}}},\quad\ddot{\lambda}=\frac{\sum\limits_{k=1}^{K}\,\mbox{tr}\left(\ddot{\mathbf{\Gamma}}\ddot{\bm{\Delta}}_{k}^{-1}\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{Y}}_{k}\right)}{PRTI}.$

  * •

Model VVE
[$\mathbf{\Sigma}_{k}=\lambda_{k}\mathbf{\Gamma}\bm{\Delta}_{k}\mathbf{\Gamma}^{\top}$]
In this case, the row covariance matrices have the same orientation, but
varying volumes and shapes. Again, there is no analytical solution for
$\mathbf{\Gamma}$, and its update is obtained by employing the MM algorithm as
described for the EVE model. Then, the updates for $\bm{\Delta}_{k}$ and
$\lambda_{k}$ are

$\ddot{\bm{\Delta}}_{k}=\frac{\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{Y}}_{k}\ddot{\mathbf{\Gamma}}\right)}{\left|\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{Y}}_{k}\ddot{\mathbf{\Gamma}}\right)\right|^{\frac{1}{P}}},\quad\ddot{\lambda}_{k}=\frac{\left|\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{Y}}_{k}\ddot{\mathbf{\Gamma}}\right)\right|^{\frac{1}{P}}}{R\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

  * •

Model EEV
[$\mathbf{\Sigma}_{k}=\lambda\mathbf{\Gamma}_{k}\bm{\Delta}\mathbf{\Gamma}_{k}^{\top}$]
Here, row covariance matrices have the same volume and shape, but different
orientations. An algorithm similar to the one proposed by Celeux and Govaert
(1995) can here employed. In detail, the eigen-decomposition
$\mathbf{Y}_{k}=\bm{L}_{k}\bm{\Omega}_{k}\bm{L}_{k}^{\top}$ is firstly
considered, with eigenvalues in the diagonal matrix $\bm{\Omega}_{k}$
following descending order and orthogonal matrix $\bm{L}_{k}$ composed of the
corresponding eigenvectors. Then, we obtain the update for
$\mathbf{\Gamma}_{k}$, $\bm{\Delta}$ and $\lambda$ as

$\ddot{\mathbf{\Gamma}}_{k}=\ddot{\bm{L}}_{k},\quad\ddot{\bm{\Delta}}=\frac{\sum\limits_{k=1}^{K}\ddot{\bm{\Omega}}_{k}}{\left|\sum\limits_{k=1}^{K}\ddot{\bm{\Omega}}_{k}\right|^{\frac{1}{P}}},\quad\ddot{\lambda}=\frac{\left|\sum\limits_{k=1}^{K}\ddot{\bm{\Omega}}_{k}\right|^{\frac{1}{P}}}{RTI}.$

  * •

Model VEV
[$\mathbf{\Sigma}_{k}=\lambda_{k}\mathbf{\Gamma}_{k}\bm{\Delta}\mathbf{\Gamma}_{k}^{\top}$]
In this setting, row covariance matrices have the same shape, but different
volumes and orientations. By using the same algorithm applied for the EEV
model, the update for $\mathbf{\Gamma}_{k}$, $\bm{\Delta}_{k}$ and
$\lambda_{k}$ are

$\ddot{\mathbf{\Gamma}}_{k}=\ddot{\bm{L}}_{k},\quad\ddot{\bm{\Delta}}=\frac{\sum\limits_{k=1}^{K}\lambda_{k}^{-1}\ddot{\bm{\Omega}}_{k}}{\left|\sum\limits_{k=1}^{K}\lambda_{k}^{-1}\ddot{\bm{\Omega}}_{k}\right|^{\frac{1}{P}}},\quad\ddot{\lambda}_{k}=\frac{\,\mbox{tr}\left\\{\ddot{\bm{\Omega}}_{k}\ddot{\bm{\Delta}}^{-1}\right\\}}{PR\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

  * •

Model EVV
[$\mathbf{\Sigma}_{k}=\lambda\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$]
In this model, row covariance matrices have varying shapes and orientations
but equal volume. The updates of this model can be obtained in a similar
fashion of the EVI model. Thus, by considering
$\mathbf{C}_{k}=\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$,
we estimate $\mathbf{\Gamma}_{k}$, $\bm{\Delta}_{k}$ and $\lambda$ as

$\ddot{\mathbf{C}}_{k}=\frac{\ddot{\mathbf{Y}}_{k}}{\left|\ddot{\mathbf{Y}}_{k}\right|^{\frac{1}{P}}},\quad\ddot{\lambda}=\frac{\sum\limits_{k=1}^{K}\left|\ddot{\mathbf{Y}}_{k}\right|^{\frac{1}{P}}}{RTI}.$

  * •

Model VVV
[$\mathbf{\Sigma}_{k}=\lambda_{k}\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$]
In the full unconstrained case, we obtain

$\ddot{\mathbf{\Sigma}}_{k}=\frac{\ddot{\mathbf{Y}}_{k}}{R\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}}.$

##### CM-Step 2

At the second CM-step, we maximize the expectation of the complete-data log-
likelihood with respect to $\bm{\Theta}_{2}$, keeping $\bm{\Theta}_{1}$ fixed
at $\ddot{\bm{\Theta}}_{1}$. The update for $\mathbf{\Psi}_{k}$ depends on
which of the 7 parsimonious structure is considered. For notational
simplicity, let $\ddot{\mathbf{W}}=\sum_{k=1}^{K}\ddot{\mathbf{W}}_{k}$ be the
update of the within state column scatter matrix, where
$\ddot{\mathbf{W}}_{k}=\sum_{i=1}^{I}\sum_{t=1}^{T}\ddot{z}_{itk}\left(\mathbf{X}_{it}-\ddot{\mathbf{M}}_{k}\right)^{\prime}\ddot{\mathbf{\Sigma}}_{k}^{-1}\left(\mathbf{X}_{it}-\ddot{\mathbf{M}}_{k}\right)$
is the update of the column scatter matrix related to the $k$th state. In
detail, we have:

  * •

Model II [$\mathbf{\Psi}_{k}=\bm{I}$] This is the simpler model, since the
column covariance matrices are spherical and assumed to be $R\times R$
identity matrices. Therefore, there are no parameters to be estimated.

  * •

Model EI [$\mathbf{\Psi}_{k}=\bm{\Delta}$] In this setting, the column
covariance matrices have the same shape and are axis-aligned. The update for
$\bm{\Delta}$ is

$\ddot{\bm{\Delta}}=\frac{\text{diag}\left(\ddot{\mathbf{W}}\right)}{\left|\text{diag}\left(\ddot{\mathbf{W}}\right)\right|^{\frac{1}{R}}}.$

  * •

Model VI [$\mathbf{\Psi}_{k}=\bm{\Delta}_{k}$] Here, the column covariance
matrices have different shapes and are axis-aligned. Thus, we update
$\bm{\Delta}_{k}$ as

$\ddot{\bm{\Delta}}_{k}=\frac{\text{diag}\left(\ddot{\mathbf{W}}_{k}\right)}{\left|\text{diag}\left(\ddot{\mathbf{W}}_{k}\right)\right|^{\frac{1}{R}}}.$

  * •

Model EE
[$\mathbf{\Psi}_{k}=\mathbf{\Gamma}\bm{\Delta}\mathbf{\Gamma}^{\top}$] In this
case, the column covariance matrices have equal shapes and orientations.
Therefore, we can directly obtain

$\ddot{\mathbf{\Psi}}=\frac{\ddot{\mathbf{W}}}{\left|\ddot{\mathbf{W}}\right|^{\frac{1}{R}}}.$

  * •

Model VE
[$\mathbf{\Psi}_{k}=\mathbf{\Gamma}\bm{\Delta}_{k}\mathbf{\Gamma}^{\top}$] In
this setting, the column covariance matrices have equal orientation but
different shapes. Similarly to the EVE and VVE models in CM-Step 1, there is
no analytical solution for $\mathbf{\Gamma}$, while keeping fixed the other
parameters. Therefore, the MM algorithm is implemented by following the same
procedure explained for the EVE model and by replacing $\mathbf{Y}$ with
$\mathbf{W}$. Then, the update of $\bm{\Delta}_{k}$ is

$\ddot{\bm{\Delta}}_{k}=\frac{\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{W}}_{k}\ddot{\mathbf{\Gamma}}\right)}{\left|\text{diag}\left(\ddot{\mathbf{\Gamma}}^{\top}\ddot{\mathbf{W}}_{k}\ddot{\mathbf{\Gamma}}\right)\right|^{\frac{1}{R}}}$

  * •

Model EV
[$\mathbf{\Psi}_{k}=\mathbf{\Gamma}_{k}\bm{\Delta}\mathbf{\Gamma}_{k}^{\top}$]
Here, the column covariance matrices have common shape but varying
orientations. By using the same approach of the EEV and VEV models, and by
changing $\ddot{\mathbf{Y}}$ with $\ddot{\mathbf{W}}$, we obtain the updates
of $\mathbf{\Gamma}_{k}$ and $\bm{\Delta}$ as

$\ddot{\mathbf{\Gamma}}_{k}=\ddot{\bm{L}}_{k},\quad\ddot{\bm{\Delta}}=\frac{\sum\limits_{k=1}^{K}\ddot{\bm{\Omega}}_{k}}{\left|\sum\limits_{k=1}^{K}\ddot{\bm{\Omega}}_{k}\right|^{\frac{1}{R}}}.$

  * •

Model VV
[$\mathbf{\Psi}_{k}=\mathbf{\Gamma}_{k}\bm{\Delta}_{k}\mathbf{\Gamma}_{k}^{\top}$]
In the full unconstrained case, we obtain

$\ddot{\mathbf{\Psi}}_{k}=\frac{\ddot{\mathbf{W}}_{k}}{\left|\ddot{\mathbf{W}}_{k}\right|^{\frac{1}{R}}}.$

####  2.3.1 A note on the initialization strategy

To start our ECM algorithm, we followed the approach of Tomarchio et al.
(2020b), where a generalization of the short-EM initialization strategy
proposed by Biernacki et al. (2003) has been implemented. It consists in $H$
short runs of the algorithm from several random positions. The term “short”
means that the algorithm is run for a few iterations $s$, without waiting for
convergence. In this manuscript, we set $H=100$ and $s=1$. Then, the parameter
set producing the largest log-likelihood is used to initialize the ECM
algorithm. In both simulated and real data analyses this procedure has shown
stable results after multiple runs.

##  3 Simulated analyses

In this section, we examine different aspects of our HMMs through large-scale
simulation studies. Considering the high number of models, we will only focus
on two of them for the sake of simplicity. In detail, the EII-II HMM (which is
the most parsimonious model) and the VVE-VE HMM (which is one of the two
models for which an MM algorithm is used both for $\mathbf{\Sigma}_{k}$ and
$\mathbf{\Psi}_{k}$) are considered. For each model, several experimental
conditions are evaluated. Specifically, we set $P=R=2$, $I=100$,
$T\in\left\\{5,10,15\right\\}$, $K\in\left\\{2,4\right\\}$ and two levels of
overlap, that will be labeled as “Overlap 1” and “Overlap 2”. Therefore,
$3\times 2\times 2=12$ scenarios are analyzed and, for each of them, 50 data
sets are generated by the considered HMMs.

About the parameters used to generate the data, when $K=2$ we set

  * •

EII-II Model

$\mathbf{\Sigma}_{1}=\mathbf{\Sigma}_{2}=\begin{bmatrix}1.50&0.00\\\
0.00&1.50\end{bmatrix},\quad\mathbf{\Psi}_{1}=\mathbf{\Psi}_{2}=\begin{bmatrix}1.00&0.00\\\
0.00&1.00\end{bmatrix},$

  * •

VVE-VE Model

$\mathbf{\Sigma}_{1}=\begin{bmatrix}0.85&0.29\\\ 0.29&0.85\end{bmatrix},\
\mathbf{\Sigma}_{2}=\begin{bmatrix}0.50&0.30\\\ 0.30&0.50\end{bmatrix},\
\mathbf{\Psi}_{1}=\begin{bmatrix}1.06&0.36\\\ 0.36&1.06\end{bmatrix},\
\mathbf{\Psi}_{2}=\begin{bmatrix}1.25&0.75\\\ 0.75&1.25\end{bmatrix},$

while for both HMMs we set $\bm{\pi}=\left(0.5,0.5\right)$,

$\mathbf{\Pi}=\begin{bmatrix}0.60&0.40\\\
0.20&0.80\end{bmatrix},\quad\mathbf{M}_{1}=\begin{bmatrix}1.00&1.50\\\
0.50&1.00\end{bmatrix}.$

The mean matrix of the second state ($\mathbf{M}_{2}$) is obtained by adding a
constant $c$ to each element of $\mathbf{M}_{1}$, which depends on the level
of overlap. Specifically, we consider $c=2$ under the “Overlap 1” scenarios,
whereas $c=5$ under the “Overlap 2” scenarios.

When $K=4$, the first two hidden states have the same
$\left\\{\mathbf{\Sigma}_{k},\mathbf{\Psi}_{k},\mathbf{M}_{k};k=1,2\right\\}$
as before. Clearly, the covariance matrices of the third and fourth hidden
states for the EII-II Model are still equal to those of the first two states.
On the contrary, for the VVE-VE model we have

$\mathbf{\Sigma}_{3}=\begin{bmatrix}1.45&1.05\\\ 1.05&1.45\end{bmatrix},\
\mathbf{\Sigma}_{4}=\begin{bmatrix}1.33&0.29\\\ 0.29&1.33\end{bmatrix},\
\mathbf{\Psi}_{3}=\begin{bmatrix}1.45&1.00\\\ 1.00&1.45\end{bmatrix},\
\mathbf{\Psi}_{4}=\begin{bmatrix}1.03&0.23\\\ 0.23&1.03\end{bmatrix}.$

Then, for both HMMs we set $\bm{\pi}=\left(0.25,0.25,0.25,0.25\right)$ and

$\mathbf{\Pi}=\begin{bmatrix}0.55&0.00&0.21&0.24\\\ 0.03&0.52&0.18&0.27\\\
0.06&0.15&0.49&0.30\\\ 0.09&0.12&0.33&0.46\end{bmatrix}.$

To obtain $\mathbf{M}_{3}$ and $\mathbf{M}_{4}$ we add $c=4$ and $c=-2$ to
each element of $\mathbf{M}_{1}$, respectively.

###  3.1 Discussion

First of all, we evaluate the recovery and the consistency of the estimated
parameters by computing the mean square errors (MSEs). Considering the high
number of parameters that should be reported, we follow an approach similar to
the one used by Farcomeni and Punzo (2020), i.e. we calculate the average of
the MSEs of each parameter of the model over the $K$ states, allowing us to
summarize in a single number the MSE of each parameter. Furthermore, before
showing the obtained results, it is important to underline the well-known
label switching issue, caused by the invariance of the likelihood function
under relabeling the model states (Frühwirth-Schnatter, 2006). There are no
generally accepted labeling methods, and we simply attribute the labels by
looking at the estimated $\mathbf{M}_{k}$.

Table 2 and Table 3 report the average MSEs, computed after fitting the EII-II
and VVE-VE HMMs, with the corresponding $K$, to the respective data sets. Note
that the $\mathbf{\Psi}$ covariance matrix is not reported in Table 2 since it
is not estimated in the EII-II HMM.

Table 2: Average of the MSEs of the parameter estimates, over the $K$ states and 50 data sets, for the EII-II HMM under each scenario. $K$ | Parameter | Overlap 1 | Overlap 2
---|---|---|---
|  | $T=5$ | $T=10$ | $T=15$ | $T=5$ | $T=10$ | $T=15$
2 | $\mathbf{M}$ | 0.0083 | 0.0040 | 0.0024 | 0.0063 | 0.0033 | 0.0023
| $\mathbf{\Sigma}$ | 0.0020 | 0.0016 | 0.0007 | 0.0016 | 0.0012 | 0.0007
| $\bm{\pi}$ | 0.0026 | 0.0029 | 0.0023 | 0.0024 | 0.0029 | 0.0023
| $\mathbf{\Pi}$ | 0.0013 | 0.0004 | 0.0004 | 0.0010 | 0.0005 | 0.0004
4 | $\mathbf{M}$ | 0.0164 | 0.0084 | 0.0055 | 0.0135 | 0.0069 | 0.0044
| $\mathbf{\Sigma}$ | 0.0029 | 0.0010 | 0.0007 | 0.0024 | 0.0013 | 0.0010
| $\bm{\pi}$ | 0.0022 | 0.0017 | 0.0022 | 0.0017 | 0.0021 | 0.0020
| $\mathbf{\Pi}$ | 0.0009 | 0.0006 | 0.0004 | 0.0009 | 0.0004 | 0.0003
Table 3: Average of the MSEs of the parameter estimates, over the $K$ states and 50 data sets, for the VVE-VE HMM under each scenario. $K$ | Parameter | Overlap 1 | Overlap 2
---|---|---|---
|  | $T=5$ | $T=10$ | $T=15$ | $T=5$ | $T=10$ | $T=15$
2 | $\mathbf{M}$ | 0.0052 | 0.0024 | 0.0015 | 0.0033 | 0.0017 | 0.0012
| $\mathbf{\Sigma}$ | 0.0021 | 0.0010 | 0.0006 | 0.0014 | 0.0008 | 0.0005
| $\mathbf{\Psi}$ | 0.0028 | 0.0014 | 0.0009 | 0.0020 | 0.0008 | 0.0007
| $\bm{\pi}$ | 0.0042 | 0.0032 | 0.0030 | 0.0029 | 0.0030 | 0.0020
| $\mathbf{\Pi}$ | 0.0010 | 0.0007 | 0.0003 | 0.0009 | 0.0005 | 0.0004
4 | $\mathbf{M}$ | 0.0176 | 0.0092 | 0.0050 | 0.0094 | 0.0050 | 0.0035
| $\mathbf{\Sigma}$ | 0.0110 | 0.0049 | 0.0032 | 0.0058 | 0.0030 | 0.0022
| $\mathbf{\Psi}$ | 0.0076 | 0.0038 | 0.0023 | 0.0041 | 0.0026 | 0.0016
| $\bm{\pi}$ | 0.0030 | 0.0021 | 0.0024 | 0.0020 | 0.0024 | 0.0021
| $\mathbf{\Pi}$ | 0.0011 | 0.0006 | 0.0005 | 0.0007 | 0.0005 | 0.0003

As we can see, the MSEs can be considered negligible in all the considered
scenarios. It is interesting to note that, for a fixed overlap, their values
become better with the increase of $T$ and that, fixed $T$, their values
roughly improve as we move from “Overlap 1” to “Overlap 2”, thus indicating a
decrease in the level of overlap. Additionally, when the VVE-VE HMM is
considered, it seems that the MM algorithms used for estimating the covariance
matrices produce reliable values.

Another aspect that is interesting to evaluate, is the computational time
required for fitting the HMMs. In detail, on each of the above data sets, all
the 98 HMMs are now fitted for the corresponding $K$, and their computational
times (in seconds) are illustrated by using the heat maps of Figures 1, 2, 3
and 4.

(a)

(b)

Figure 1: Heat maps of the average computational time for the 98 HMMs,
computed over 50 data sets, when the data are generated by a EII-II HMM with
$K=2$ and “Overlap 1” (a) or “Overlap 2” (b).

(a)

(b)

Figure 2: Heat maps of the average computational time for the 98 HMMs,
computed over 50 data sets, when the data are generated by a EII-II HMM with
$K=4$ and “Overlap 1” (a) or “Overlap 2” (b).

(a)

(b)

Figure 3: Heat maps of the average computational time for the 98 HMMs,
computed over 50 data sets, when the data are generated by a VVE-VE HMM with
$K=2$ and “Overlap 1” (a) or “Overlap 2” (b).

(a)

(b)

Figure 4: Heat maps of the average computational time for the 98 HMMs,
computed over 50 data sets, when the data are generated by a VVE-VE HMM with
$K=4$ and “Overlap 1” (a) or “Overlap 2” (b).

Computation is performed on a Windows 10 PC, with AMD Ryzen 7 3700x CPU, 16.0
GB RAM, using the R 64-bit statistical software (R Core Team, 2019), and the
proc.time() function of the base package is used to measure the time. As it is
reasonable to expect, the computational time grows as $T$ increases on each
scenario, and it approximately halves when we pass from “Overlap 1” to
“Overlap 2”, highlighting the easier of estimation in this case. Furthermore,
with the exclusion of Figure LABEL:fig:vve_c_4, the computational time
approximately triplicates when we move from fitting HMMs with $K=2$ to HMMs
with $K=4$ hidden states. It is interesting to note that the EVE-VE and VVE-VE
HMMs, which are the two models for which we use a MM algorithm for estimating
both covariance matrices, are the most time consuming, with a computational
burden that seems to double with respect to the other models. This is
particularly evident in the “Overlap 2” scenarios.

The total computational time can be strongly reduced by exploiting parallel
computing. In detail, Table 4 shows the overall time taken by fitting the 98
HMMs sequentially (default in R) and by parallelizing them on 14 cores. As we
can see, the computational burden is decreased by about 10 times, and all the
models can be fitted in a reasonable fast way (with some exceptions in the
“Overlap 1” scenarios).

Table 4: Computational times (in seconds) for running the algorithm sequentially or via parallel computing. They refer to the fit of all the 98 HMMs with $K$ states, averaged over the 50 data sets, and generated by the two HMMs on each scenario. Type | HMM | $K$ | Overlap 1 | Overlap 2
---|---|---|---|---
|  |  | $T=5$ | $T=10$ | $T=15$ | $T=5$ | $T=10$ | $T=15$
Sequential | EII-II | 2 | 82.38 | 159.96 | 252.27 | 41.37 | 78.03 | 103.66
4 | 346.51 | 753.71 | 1236.52 | 89.32 | 176.54 | 259.24
VVE-VE | 2 | 103.31 | 199.63 | 301.53 | 43.19 | 84.66 | 123.98
| 4 | 473.74 | 1018.19 | 1551.68 | 116.69 | 270.18 | 380.48
Parallel | EII-II | 2 | 9.65 | 15.53 | 22.68 | 6.63 | 9.42 | 11.33
4 | 29.87 | 60.01 | 95.54 | 10.39 | 17.34 | 23.82
VVE-VE | 2 | 11.12 | 18.43 | 26.06 | 6.70 | 9.75 | 12.79
| 4 | 43.11 | 86.26 | 130.60 | 12.19 | 23.77 | 32.46

Lastly, the capability of the Bayesian information criterion (BIC; Schwarz et
al., 1978) in identifying the true parsimonious structure and the correct
number of groups is investigated. This is because, so far, we have fitted
models with $K$ equal to the true number of states present in the data, and we
need to assess if the BIC, which is one of the most famous and used tools in
model-based clustering, accurately works. Therefore, on each of the above data
sets, the 98 HMMs are fitted for $K\in\left\\{1,\ldots,K+1\right\\}$, and the
results are reported in Table 5. First of all, in each scenario, the true $K$
has been always selected by the best fitting model according to the BIC (for
this reason this information is not reported in Table 5). Additionally, we
notice that in almost all the cases the true data generating model has been
identified by the BIC. In those few cases where the BIC selects a wrong model,
this is because of an incorrect choice of the parsimonious structure for one
of the two covariance matrices $\mathbf{\Sigma}$ or $\mathbf{\Psi}$.

Table 5: Number of times, over the 50 data sets generated by the two HMMs on each scenario, for which the true parsimonious structure is selected by the BIC when all the 98 HMMs are fitted for $k\in\left\\{1,\ldots,K+1\right\\}$. HMM | $K$ | Overlap 1 | Overlap 2
---|---|---|---
|  | $T=5$ | $T=10$ | $T=15$ | $T=5$ | $T=10$ | $T=15$
EII-II | 2 | 47 | 48 | 48 | 49 | 49 | 47
| 4 | 46 | 50 | 50 | 50 | 50 | 49
VVE-VE | 2 | 45 | 49 | 48 | 48 | 49 | 48
| 4 | 47 | 49 | 50 | 50 | 50 | 50

##  4 Real data example

In this section, we analyze data concerning the unemployment rate in the
Italian provinces (NUTS3, according to the European Nomenclature of
Territorial Units for Statistics). The data comes from the Italian National
Institute of Statistics (ISTAT), a public research organization and the main
producer of official statistics in the service of citizens and policy-makers
in Italy, and are freely accessible at http://dati.istat.it/#. In detail, we
investigate the $I=98$ Italian provinces for which the unemployment rate is
available from the beginning of the data collection at the provincial level
(2004) to the most recent year (2019). This implies that we are considering
$T=16$ years of data. Note that some provinces are not included in the
analysis since the data were available for only few years.

For each province, the unemployment rate is recorded in a two-factor format.
The first factor, gender, has two levels (i.e. $P=2$): males and females. The
second factor, age, has three levels (i.e. $R=3$) driven by the age category:
15–24, 25–34 and 35–older. Therefore, the whole data set is presented in a
four-way array having dimensions $2\times 3\times 98\times 16$. The
unemployment rates are then mapped to the real line by using the logit
transformation, as commonly done in this branch of literature for overcoming
boundary bias problems (see, e.g. Wallis, 1987; Koop and Potter, 1999;
Hudomiet, 2015).

In analyzing this data set, several questions arise concerning the existence
of areas with similar unemployment levels among the Italian provinces.
Historically Southern Italy have always shown worse economic performance with
respect to the rest of the Country (Daniele et al., 2007). However, within
each region (NUTS2, according to the European Nomenclature of Territorial
Units for Statistics) there can be considerable differences among the
provinces that are part of it. Also of interest is the strength of time
dependence as measured by the transition probability matrix, as well as how
the provinces move between the hidden states. This latter aspect can be
particularly of interest in light of the two main recessions that the Italian
economy faced in 2008 and 2011. Relatedly, an overview of which provinces have
best withstood the crises or which have been able to recover from the crises
can be easily obtained.

###  4.1 Discussion

Our 98 HMMs are fitted to the data for $K\in\left\\{1,\ldots,10\right\\}$ and,
according to the BIC, the best model is the VEV-EE with $K=7$ hidden states.
Despite the logit scale, we can easily interpret the $K=7$ states by looking
at the estimated mean matrices

$\mathbf{M}_{1}=\begin{bmatrix}-2.34&-3.50&-4.14\\\
-1.84&-2.82&-3.33\end{bmatrix},\quad\mathbf{M}_{2}=\begin{bmatrix}-1.61&-2.88&-3.57\\\
-1.27&-2.39&-3.05\end{bmatrix},$
$\mathbf{M}_{3}=\begin{bmatrix}-1.67&-2.87&-3.38\\\
-1.27&-2.08&-2.69\end{bmatrix},\quad\mathbf{M}_{4}=\begin{bmatrix}-1.05&-2.33&-3.02\\\
-0.80&-1.93&-2.68\end{bmatrix},$
$\mathbf{M}_{5}=\begin{bmatrix}-0.80&-1.86&-2.83\\\
-0.44&-1.36&-2.35\end{bmatrix},\quad\mathbf{M}_{6}=\begin{bmatrix}-0.37&-1.38&-2.31\\\
-0.09&-1.03&-2.10\end{bmatrix},$
$\mathbf{M}_{7}=\begin{bmatrix}0.14&-0.85&-1.81\\\
0.39&-0.53&-1.70\end{bmatrix}.$

As we can note, it is possible to sort the states according to growing
unemployment levels, both in the gender and ages factors. Specifically, as we
move from the first to the seventh state the unemployment rises, and each
state becomes worse than the previous ones under each point of view. The only
occasion where this happens partially concerns the states two and three.
Indeed, regardless of gender, state two shows better rates for people over 25
and conversely state three is preferred for people under 25. However, since 4
times out of 6 state two is preferred, and considering that this involves the
majority of people, we could globally consider it better than state three. We
can also observe that, regardless of the considered state, the unemployment
levels are higher for females and get lower as the age increases. It might be
also interesting to report that, regardless of the age class, the lowest
relative differences between the two genders are in the seventh state, the
worst.

Useful insights can also be obtained from the analysis of the estimated
covariance matrices. The best fitting model suggests various volumes and
orientations but equal shapes for the gender-related covariance matrices,
whereas all shapes and orientations are found to be the same for the age-based
covariance matrices. In detail, the estimated gender-related covariance
matrices are

$\mathbf{\Sigma}_{1}=\begin{bmatrix}0.27&0.03\\\ 0.03&0.20\end{bmatrix},\hskip
1.70717pt\mathbf{\Sigma}_{2}=\begin{bmatrix}0.13&0.01\\\
0.01&0.10\end{bmatrix},\hskip
1.70717pt\mathbf{\Sigma}_{3}=\begin{bmatrix}0.36&0.02\\\
0.02&0.25\end{bmatrix},$ $\mathbf{\Sigma}_{4}=\begin{bmatrix}0.11&0.02\\\
0.02&0.10\end{bmatrix},\hskip
1.70717pt\mathbf{\Sigma}_{5}=\begin{bmatrix}0.09&0.02\\\
0.02&0.12\end{bmatrix},\hskip
1.70717pt\mathbf{\Sigma}_{6}=\begin{bmatrix}0.12&0.02\\\
0.02&0.16\end{bmatrix},\hskip
1.70717pt\mathbf{\Sigma}_{7}=\begin{bmatrix}0.09&0.01\\\
0.01&0.11\end{bmatrix}.$

We notice that, excluding the first and third, all the states have similar
covariance matrices. Furthermore, while the first four states have higher
variances for men, the last three show greater variances for women.
Additionally, it is interesting to see that the third state has the largest
variances both for both men and women. With respect to the age-based
covariance matrices, showed in following,

$\mathbf{\Psi}_{1},\ldots,\mathbf{\Psi}_{7}=\begin{bmatrix}1.66&0.18&0.19\\\
0.18&1.02&0.11\\\ 0.19&0.11&0.63\end{bmatrix},$

we can note that the variances become lower as the age categories grows.

Lastly, before showing how these states cluster the Italian provinces, it is
worth analyzing the estimated transition probability matrix

$\mathbf{\Pi}=\begin{bmatrix}0.80&0.20&0.00&0.00&0.00&0.00&0.00\\\
0.03&0.72&0.01&0.24&0.00&0.00&0.00\\\ 0.00&0.14&0.72&0.05&0.06&0.03&0.00\\\
0.00&0.03&0.08&0.83&0.06&0.00&0.00\\\ 0.00&0.00&0.00&0.02&0.84&0.14&0.00\\\
0.00&0.00&0.01&0.00&0.04&0.86&0.09\\\
0.00&0.00&0.00&0.00&0.00&0.04&0.96\end{bmatrix},$

and how many provinces have changed their state in each year (from the second
year onwards). The latter aspect is illustrated in the bar plot of Figure 5.

Figure 5: Number of provinces that changed state in each year (from the second
year onwards).

As we can note by the estimated transition probability matrix, transitions
between states are not uncommon between adjacent states, whereas are null
among distant states. Furthermore, it seems that the persistence of staying in
a state, increases as we move from the fourth to the seventh, i.e. it appears
increasingly difficult for the provinces clustered in the troubled states to
improve their position. In particular, the last state shows a high
persistence.

From the analysis of the bar plot, the highest number of switches between the
states occurs over the years 2009–2013, as a consequence of the two
aforementioned economic recessions. Indeed, all the provinces showed an
increase in the unemployment rates in those years, causing a change towards
worse states. This can be better understood by looking at the Italian
provinces maps of Figures 6 and 7, that are colored according to state
memberships. Note that the provinces not included in the analysis are colored
in gray. For simplicity, we avoid to plot a map for each of the 16-years of
data, and we limit to report some key years. It is interesting to note that,
in almost all the cases, not all the provinces belonging to the same region
are clustered in the same state. Therefore, working at province-level may
offer additional insight than working with regional-level.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Italian provinces map colored according to the estimated state
memberships

(a)

(b)

Figure 7: Italian provinces map colored according to the estimated state
memberships (cont.)

Starting from the first year of analysis, i.e. 2004, in Figure LABEL:fig:2004
we can recognize several clusters of provinces that, as we move towards the
south, belong to states with higher unemployment rates. After some years
(2005–2008) characterized by relatively few changes among the states, the
first economic recession began to produce its effects in 2009, where a lot of
provinces started to perform badly (see Figure LABEL:fig:2009). Such event,
further strengthened by the second economic recession of 2011, lead to a
continuous rise of the unemployment levels, that (see LABEL:fig:2010,
LABEL:fig:2011, LABEL:fig:2012 and LABEL:fig:2013) lead the majority of
provinces to the worst states, reducing the pre-existing differences between
them. After two years of small changes (2014–2015), the provinces located in
the upper part of the Country started to recover, whereas the rest of the
provinces failed to start again, remaining anchored in the post-recessions
difficulties (see LABEL:fig:2018 and LABEL:fig:2019 that show the most two
recent years). In any case, these signs of recovery are going to be
dramatically arrested by the COVID-19 pandemic, and its effects will have
serious repercussions in the next years.

It is also interesting to report the behavior of some specific provinces. For
example, there is only one province that, over the 16 years, has never changed
its state. Specifically, the province of Bolzano, the northernmost in Italy,
has been always in state one, the best. This means that, regardless of the
economic crises, the unemployment issue has been virtuously managed in this
province. We can also mention the two provinces that have worsened their
position the most over the 16 years. In detail, the province of Ancona and
Ferrara were in state one in the first years of our analysis, but they are
actually located in the fifth state (since 2013 and 2016, respectively).
Conversely, none of the provinces is actually in a state that can be
considered better than the one had in 2004. This means that all the provinces
have kept (Bolzano), lost (most of the provinces) and at best regained (only
some provinces) the state they had at the beginning of our analysis.

##  5 Conclusions

In this manuscript we introduced parsimonious hidden Markov models for matrix-
variate longitudinal data. Being (dependent) mixture models, they allow the
recovery of homogenous latent subgroups and, simultaneously, provide
meaningful interpretation on how the sample units move between the hidden
states over time. The parsimony has been introduced via the eigen
decomposition of the state covariance matrices, producing a family of 98 HMMs.
An ECM algorithm has been illustrated for parameter estimation. At first, the
parameter recovery of our algorithm has been evaluated under different
scenarios, providing good results. This can be particularly interesting for
those HMMs that use a MM algorithm at each step of the ECM algorithm.
Relatedly, we have analyzed the computational times for fitting all the 98
HMMs. The computational burden of the HMMs using MM algorithm is definitely
higher, even if we are able to fit all the HMMs in a quite fast way when
parallel computing is considered. The BIC has proven to be effective in
detecting the true number of states in the data as well as the parsimonious
structure. The real data example has shown the usefulness of our HMMs. Indeed,
other than identifying different states, they have provided a tool for easily
analyzing the evolution over time of the unemployment at province level, and
for obtaining useful insight on the behavior of some specific provinces.

There are different possibilities for further work, some of which are worth
mentioning. First of all, we can extend our HMMs by using skewed or heavy
tailed state dependent probability density functions (Gallaugher and
McNicholas, 2017, 2019; Tomarchio et al., 2020b, a), in order to model
possible features commonly present in the data. Another extension would deal
with the regression setting (Viroli, 2012), where covariates shared by all
units in the same hidden state are used. This can be done both in a fixed and
in random covariates framework.

## References

  * Altman (2007) Altman RM (2007) Mixed hidden Markov models. Journal of the American Statistical Association 102(477):201–210
  * Bartolucci et al. (2012) Bartolucci F, Farcomeni A, Pennoni F (2012) Latent Markov models for longitudinal data. CRC Press
  * Baum et al. (1970) Baum LE, 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
  * Biernacki et al. (2003) Biernacki C, Celeux G, Govaert G (2003) Choosing starting values for the em algorithm for getting the highest likelihood in multivariate Gaussian mixture models. Computational Statistics & Data Analysis 41(3-4):561–575
  * Browne and McNicholas (2014) Browne RP, McNicholas PD (2014) Estimating common principal components in high dimensions. Advances in Data Analysis and Classification 8(2):217–226
  * Brunner and Puri (2001) Brunner E, Puri ML (2001) Nonparametric methods in factorial designs. Statistical papers 42(1):1–52
  * Celeux and Govaert (1995) Celeux G, Govaert G (1995) Gaussian parsimonious clustering models. Pattern recognition 28(5):781–793
  * Daniele et al. (2007) Daniele V, Malanima P, et al. (2007) Il prodotto delle regioni e il divario nord-sud in italia (1861-2004). Rivista di Politica Economica 97(2):267–316
  * Dempster et al. (1977) Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 39(1):1–38
  * Farcomeni (2012) Farcomeni A (2012) Quantile regression for longitudinal data based on latent markov subject-specific parameters. Statistics and Computing 22(1):141–152
  * Farcomeni and Punzo (2020) Farcomeni A, Punzo A (2020) Robust model-based clustering with mild and gross outliers. TEST 29:989–1007
  * Farcomeni et al. (2020) Farcomeni A, Ranalli M, Viviani S (2020) Dimension reduction for longitudinal multivariate data by optimizing class separation of projected latent markov models. TEST pp 1–19
  * Fitzmaurice and Ravichandran (2008) Fitzmaurice GM, Ravichandran C (2008) A primer in longitudinal data analysis. Circulation 118(19):2005–2010
  * Frühwirth-Schnatter (2006) Frühwirth-Schnatter S (2006) Finite mixture and Markov switching models. Springer Science & Business Media
  * Gallaugher and McNicholas (2017) Gallaugher MP, McNicholas PD (2017) A matrix variate skew-t distribution. Stat 6(1):160–170
  * Gallaugher and McNicholas (2019) Gallaugher MP, McNicholas PD (2019) Three skewed matrix variate distributions. Statistics & Probability Letters 145:103–109
  * Gupta and Nagar (2018) Gupta AK, Nagar DK (2018) Matrix variate distributions, vol 104. CRC Press
  * Huang et al. (2019) Huang L, Bai J, Ivanescu A, Harris T, Maurer M, Green P, Zipunnikov V (2019) Multilevel matrix-variate analysis and its application to accelerometry-measured physical activity in clinical populations. Journal of the American Statistical Association 114:553–564
  * Hudomiet (2015) Hudomiet P (2015) The role of occupation specific adaptation costs in explaining the educational gap in unemployment. Tech. rep., Working paper
  * Koop and Potter (1999) Koop G, Potter SM (1999) Dynamic asymmetries in us unemployment. Journal of Business & Economic Statistics 17(3):298–312
  * Maruotti (2011) Maruotti A (2011) Mixed hidden markov models for longitudinal data: An overview. International Statistical Review 79(3):427–454
  * Maruotti et al. (2017) Maruotti A, Bulla J, Lagona F, Picone M, Martella F (2017) Dynamic mixtures of factor analyzers to characterize multivariate air pollutant exposures. Annals of Applied Statistics 11(3):1617–1648, DOI 10.1214/17-AOAS1049, URL https://doi.org/10.1214/17-AOAS1049
  * Melnykov and Zhu (2019) Melnykov V, Zhu X (2019) Studying crime trends in the usa over the years 2000–2012. Advances in Data Analysis and Classification 13(1):325–341
  * Meng and Rubin (1993) Meng XL, Rubin DB (1993) Maximum likelihood estimation via the ecm algorithm: A general framework. Biometrika 80(2):267–278
  * Noguchi et al. (2012) Noguchi K, Gel YR, Brunner E, Konietschke F (2012) nparLD: an R software package for the nonparametric analysis of longitudinal data in factorial experiments. Journal of Statistical Software 50(12):1–23
  * R Core Team (2019) R Core Team (2019) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, URL https://www.R-project.org/
  * Sarkar et al. (2020) Sarkar S, Zhu X, Melnykov V, Ingrassia S (2020) On parsimonious models for modeling matrix data. Computational Statistics & Data Analysis 142:106822
  * Schwarz et al. (1978) Schwarz G, et al. (1978) Estimating the dimension of a model. Annals of Statistics 6(2):461–464
  * Tomarchio et al. (2020a) Tomarchio SD, Gallaugher MP, Punzo A, McNicholas PD (2020a) Mixtures of contaminated matrix variate normal distributions. arXiv preprint arXiv:200503861
  * Tomarchio et al. (2020b) Tomarchio SD, Punzo A, Bagnato L (2020b) Two new matrix-variate distributions with application in model-based clustering. Computational Statistics & Data Analysis 152:107050
  * Verbeke et al. (2014) Verbeke G, Fieuws S, Molenberghs G, Davidian M (2014) The analysis of multivariate longitudinal data: A review. Statistical Methods in Medical Research 23(1):42–59
  * Verdam and Oort (2019) Verdam M, Oort F (2019) The analysis of multivariate longitudinal data: An instructive application of the longitudinal three-mode model. Multivariate behavioral research 54(4):457–474
  * Viroli (2011a) Viroli C (2011a) Finite mixtures of matrix normal distributions for classifying three-way data. Statistics and Computing 21(4):511–522
  * Viroli (2011b) Viroli C (2011b) Model based clustering for three-way data structures. Bayesian Analysis 6(4):573–602
  * Viroli (2012) Viroli C (2012) On matrix-variate regression analysis. Journal of Multivariate Analysis 111:296–309
  * Wallis (1987) Wallis KF (1987) Time series analysis of bounded economic variables. Journal of Time Series Analysis 8(1):115–123
  * Wang and Melnykov (2020) Wang Y, Melnykov V (2020) On variable selection in matrix mixture modelling. Stat 9(1):e278
  * Zhu and Melnykov (2021) Zhu X, Melnykov Y (2021) On finite mixture modeling of change-point processes. Journal of Classification pp 1–20
  * Zucchini et al. (2017) Zucchini W, MacDonald IL, Langrock R (2017) Hidden Markov models for time series: an introduction using R. CRC press