# Copula Modelling of Serially Correlated Multivariate Data with Hidden
Structures

Robert Zimmerman, Radu V. Craiu and Vianey Leos-Barajas

###### Abstract

We propose a copula-based extension of the hidden Markov model (HMM) which
applies when the observations recorded at each time in the sample are
multivariate. The joint model produced by the copula extension allows decoding
of the hidden states based on information from multiple observations. However,
unlike the case of independent marginals, the copula dependence structure
embedded into the likelihood poses additional computational challenges. We
tackle the latter using a theoretically-justified variation of the EM
algorithm developed within the framework of inference functions for margins.
We illustrate the method using numerical experiments and an analysis of house
occupancy.

Keywords: copulas, EM algorithm, hidden Markov model, inference functions for
margins.

##  1 Introduction

Statistical models aim to capture the generating mechanism of the observed
data. In many instances, one must consider various forms of dark data [Hand,
2020], such as missing data that could not be observed at all or that has been
modified during the sampling stage, or data that is unobservable but can be
injected into the model as latent variables to add meaning and to enhance the
interpretability of the statistical model. The latter approach is used to set
up a hidden Markov model (HMM) in which the observed data consists of serially
correlated observations on each item. The latent aspect of the model is a
Markov chain with a discrete state space that evolves on the same time scale
as the observed processes. The HMM postulates that the distribution of the
observed data depends on the state of the Markov chain at that same time.
Estimating the hidden structure can often illuminate the underlying workings
of the system [Rabiner and Juang, 1986, Zucchini et al., 2017]; in some cases,
the Markov latent structure can be associated to a real hidden mechanism with
some certainty. Examples of HMMs abound in numerous domains including
economics [Du et al., 2020], medicine [Williams et al., 2020], ecology
[Ailliot et al., 2009], and sports [Ötting et al., 2021].

Several distinct problems can be tackled using the HMM formulation. If one
aims to identify the number of states that could have generated the data, and
subsequently to _decode_ the state sequence, an unsupervised approach relies
on clustering the observed features and assigning each to the state most
likely to have generated it. On the other hand, the supervised version of the
analysis implies that enough items have been monitored and their hidden states
labeled to constitute a complete set of training data. In this case, one aims
to predict the hidden states for observations outside the training sample. In
either the supervised or unsupervised setting, one might also aim to impute or
predict missing values from the observed features in a dataset that exhibits
missing patterns.

Our work is motivated by situations in which the observations recorded at each
time for each item are multivariate, and one must decide how to _jointly_
integrate the information they contain about the hidden states. When the
observed data in an HMM setup is vector-valued, the state-dependent
distributions can be constructed by assuming either contemporaneous
conditional independence or longitudinal conditional independence [Zucchini et
al., 2017]. In the former case, the multivariate state-dependent distribution
for the vector-valued observations is constructed as a product of marginal
distributions, while in the latter case it is typically assumed to be
multivariate Gaussian, which is highly tractable and easily implemented within
the standard array of HMM algorithms [Cappé et al., 2005, Zucchini et al.,
2017]. However, the data which the multivariate Gaussian distribution can
model is limited; in the bivariate case, for example, it cannot capture
dependence in the extremes of the upper-right or lower-left tails (a property
known as _tail dependence_ [Embrechts et al., 2001]). One of our main aims
here is to propose a general approach for integration of the continuous-valued
information provided by multivariate observations in HMM models.

Copulas have become a ubiquitous tool in modelling complex dependence
structures. Sklar [1959] provided a theoretical foundation and demonstrated
that any multivariate distribution can be represented by its marginal
distributions and a copula that fully describes their interdependence.
Moreover, this decomposition is unique whenever the marginal distributions are
continuous. Copulas have been widely used in analysis of survival data
[Goetghebeur and Ryan, 1995, Hougaard, 2000] and in various statistical
applications, including actuarial science [Frees and Valdez, 1998], hydrology
[Genest et al., 2007, Dupuis, 2007] and finance [Cherubini et al., 2004, da
Silva Filho et al., 2012], to name only a few. The developments of conditional
copulas [Acar et al., 2011, Veraverbeke et al., 2011, Chavez-Demoulin and
Vatter, 2015, Craiu and Sabeti, 2012, Levi and Craiu, 2018] and vine copulas
[Aas et al., 2009, Min and Czado, 2010] have further expanded the range of
dependence structures available for statistical modelling.

The paper makes several contributions. First, we extend the modelling toolbox
for HMMs with multivariate observations by considering a copula model for the
distribution of observed data. Our numerical experiments show that when the
copula is allowed to vary with the hidden state variable, the accuracy of
hidden state identification increases. Because using copulas to jointly model
the observed data presents additional computational challenges compared to the
case of independent marginals or multivariate Gaussians, our second main
contribution is to develop a new optimization procedure in which we integrate
the _inference functions for margins (IFM)_ method of Joe and Xu [1996] within
the _ES algorithm_ of Elashoff and Ryan [2004].

In the next section we introduce the model along with numerical evidence of
the gains in decoding precision when using the copula-based joint modelling of
the observed data. The computational methods used for estimation are
introduced and theoretically justified in Sections 3 and 4. Section 5 contains
numerical experiments based on simulations as well as analysis of house
occupancy data. The paper ends with a summary and a discussion of future work.

##  2 A Copula-Based HMM

###  2.1 The General Model

Let $X_{1},X_{2},\ldots$ be an unobserved discrete-time first-order Markov
process taking values in a finite state space $\mathcal{X}=\\{1,2,\ldots,K\\}$
with initial distribution $\bm{\pi}:\mathcal{Q}\to[0,1]$ and transition
probability matrix $\bm{\Gamma}=[\gamma_{i,j}]_{i,j\in\mathcal{X}}$. Let
$\mathbf{Y}_{1},\mathbf{Y}_{2},\ldots\mathbf{Y}_{T}\in\mathbb{R}^{d}$
represent our observed data, assumed to satisfy the conditional independence
structure $(\mathbf{Y}_{i}\mid X_{i})\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}(\mathbf{Y}_{j}\mid X_{j})$ for $i\neq j$.

$\cdots$$X_{t-2}$$X_{t-1}$$X_{t}$$X_{t+1}$$X_{t+2}$$\cdots$$\mathbf{Y}_{t-2}$$\mathbf{Y}_{t-1}$$\mathbf{Y}_{t}$$\mathbf{Y}_{t+1}$$\mathbf{Y}_{t+2}$

Figure 1: Standard HMM dependence structure

We assume that for each state $k\in\mathcal{X}$, there exists a multivariate
distribution $H_{k}(\cdot):\mathbb{R}^{d}\to[0,1]$ defined by $d$ univariate
$k$’th state-dependent marginal distributions,
$F_{k,1}(\cdot\,;\lambda_{k,1}),\ldots,F_{k,d}(\cdot\,;\lambda_{k,d}):\mathbb{R}\to[0,1]$
and a $d$-dimensional parametric copula $C_{k}:[0,1]^{d}\to[0,1]$, such that

$\mathbf{Y}_{i}\mid(X_{i}=k)\sim
H_{k}(\cdot)=C_{k}\Big{(}F_{k,1}(\cdot\,;\lambda_{k,1}),\ldots,F_{k,d}(\cdot\,;\lambda_{k,d})\,\Big{|}\,\theta_{k}\Big{)}.$
(1)

This is the formulation of a generic finite-state time-homogeneous HMM (see
Figure 1). We further assume that each $C_{k}$ and $F_{k,h}$ admit densities
$c_{k}$ and $f_{k,h}$ with respect to Lebesgue measure (on $\mathbb{R}^{d}$
and $\mathbb{R}$, respectively), implying the existence of the $k$’th state-
dependent joint density

$h_{k}(\mathbf{y})=c_{k}(F_{k,1}(y_{1};\lambda_{k,1}),\ldots,F_{k,d}(y_{d};\lambda_{k,1})\mid\theta_{k})\cdot\prod_{h=1}^{d}f_{k,h}(y_{h};\lambda_{k,h}).$

Each $\lambda_{k,h}\in\mathbb{R}$ is a state- and margin-specific parameter,
whose underlying parameter space depends on both the state $k$ and the
component $h$. Thus, each state $k\in\mathcal{X}$ generates a $d$-dimensional
observation whose joint dependence structure is governed by a particular
state-specific copula $C_{k}$ with associated parameter $\theta_{k}$. Note
that while we assume that the $\lambda_{k,h}$’s and $\theta_{k}$’s are scalars
for ease of presentation, our theory easily generalizes to multivariate
parameters.

We further assume that each state-dependent copula $C_{k}$ is a member of a
parametric family whose log-densities (and derivatives thereof) can be
computed. The Gaussian copula is typically used to impose the same kind of
dependence structure as the multivariate Gaussian distribution, while the $t$
copula is more suitable for capturing dependent extreme values [Demarta and
McNeil, 2005]. _Archimedean_ copulas [Nelsen, 2007] are also popularly used in
statistical modelling. Each Archimedean copula can be explicitly represented
as
$C_{\theta}(\mathbf{u})=\psi_{\theta}[\psi_{\theta}^{-1}(u_{1})+\cdots+\psi_{\theta}^{-1}(u_{d})]$
for $\mathbf{u}\in[0,1]^{d}$, where the _generator_
$\psi_{\theta}:[0,\infty]\to[0,1]$ is a continuous decreasing function that
satisfies $\psi_{\theta}(0)=1$ and $\lim_{t\to\infty}\psi_{\theta}(t)=0$,
which is strictly decreasing on its support [Hofert et al., 2012]. The
explicit functional representation of Archimedean copulas does not necessarily
imply the analytic tractability of the usual quantities used for likelihood
inference, but Hofert et al. [2012] have derived explicit functional forms for
the densities, generator derivatives, and score functions in the five
“classical” one-parameter Archimedean copula families.

###  2.2 Coupling Benefits: an Illustration

HMMs are generally used for two different purposes: to model a data-generating
mechanism in which the Markov process may serve as a proxy for some process of
interest, and to establish a classification of the hidden states. In
applications where multivariate data regularly switch between different
complex dependence structures and hierarchies, we can model the data-
generating process as a finite-state HMM in which the state-dependent marginal
distributions are linked together via copulas. The copula-within-HMM framework
is highly flexible, because one has the freedom to vary both the copula itself
as well as the marginal distributions between states (although as with HMMs in
general, problems may arise if the state-dependent marginal distributions
within each dimensions do not remain within the same parametric family). Such
models have already been applied to financial data under the label of _regime-
switching copulas_ ; for example, Nasri et al. [2020] have used them for
option pricing and have developed a goodness-of-fit approach to selecting the
number of hidden states. In these applications, where the focus is on the data
generating mechanism, the HMM is assessed not by its capacity to predict the
states accurately, but rather via metrics that assess goodness-of-fit and the
model’s capacity to replicate key features of the data.

A copula-within-HMM framework can also be used for the second aim, in which
case the model would be assessed by its capacity to predict the hidden states.
In this situation, a simple example shows that an HMM performs poorly when
assuming contemporaneous conditional independence of the observation process,
even when the marginal state-dependent distributions are correctly specified.
Consider bivariate data arising from a 2-state finite mixture model of length
$T=100$ with equal persistence between states (which may be viewed as a
special case of an HMM), where the state-dependent distributions are Frank
copulas [Genest, 1987] with standard normal margins, and the state-dependent
copula parameters are $\theta_{1}=-\theta$ and $\theta_{2}=\theta$ for some
fixed value of $\theta\in(0,100)$, representing extreme negative dependence
and extreme positive dependence, respectively, when $|\theta|$ is large. In
Figure 2, we have plotted the zero-one loss for this model (i.e.,
$\ell_{01}=\frac{1}{T}\sum_{t=1}^{T}\mathbbm{1}_{\hat{X}_{t}\neq X_{t}}$,
where $\hat{X}_{t}$ is the predicted state at time $t$) based on local state
decoding as a function of $|\theta|$, as well as the zero-one loss based on an
incorrectly-specified model with independent marginals. The zero-one loss for
the independence model is, of course, constant; on the other hand, the zero-
one loss for the true model appears to decay exponentially in $|\theta|$.

Figure 2: Empirical zero-one loss (black dots) for local decoding of a 2-state
HMM with standard normal margins and Frank copulas with parameters $\theta$
and $-\theta$, for $|\theta|\in\\{1,2,\ldots,100\\}$. The red line is the
zero-one loss when predictions are made assuming independent margins (which is
constant in $\theta$); the blue line is the true expected zero-one loss at
$\theta$.

In fact, the expected zero-one loss as a function of $\theta$ is given by
$\ell_{01}(\theta)=\frac{1}{2}-\frac{2}{\theta}\hskip
1.0pt\mathrm{log}\left(\cosh{\frac{\theta}{4}}\right)$. Such a measure can
also be derived for other radially symmetric bivariate copula families within
the same finite mixture model (see A.1 within Subsection A.2); for example,
under the same setup but with Gauss copulas in place of Frank copulas with
state-dependent parameters $\rho_{1}=-\rho$ and $\rho_{2}=\rho\in(0,1)$, it
can be shown that $\ell_{0,1}(\rho)=\mathrm{cos}^{-1}(\rho)/{\pi}$. These are,
in fact, special cases of a far more general result:

######  Theorem 2.1.

Let $\nu_{t,k}=\mathbb{P}\left(X_{t}=k\right)=[\bm{\pi}\bm{\Gamma}^{t}]_{k}$.
The expected zero-one loss of the classifications made by local decoding is
given by

$\ell_{01}(\bm{\eta})=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,k}\cdot\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1\right|X_{t}=k\right)$
(2)

where $\kappa(x_{-t},\mathbf{y}_{-t})=\pi_{x_{1}}\prod_{s\neq
t}h_{x_{s}}(\mathbf{y}_{s})\cdot\prod_{s\neq t,t+1}\gamma_{x_{s-1},x_{s}}$.
Moreover, if the HMM constitutes a finite mixture model, then (2) reduces to

$\ell_{01}(\bm{\eta})=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,k}\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{\nu_{t,k}\cdot
h_{k}(\mathbf{y}_{t})}{\max_{j\neq k}\nu_{t,j}\cdot
h_{j}(\mathbf{y}_{t})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}_{t}).$
(3)

The complicated form of the HMM posterior state probabilities makes it
challenging to derive useful bounds for the zero-one loss of the most general
models; however, under certain conditions on the underlying Markov chain, one
can deduce upper bounds for individual components of the loss function. For
example, if the Markov chain spends most of its time travelling through state
$k$, then the contributions to the loss function provided by the $k$’th state
can be upper bounded by a more elegant expression:

######  Proposition 2.1.

Suppose that $\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid
X_{t-1},X_{t+1}\right)\geq\mathbb{P}_{\bm{\eta}}\left(X_{t}=j\mid
X_{t-1},X_{t+1}\right)$ for all $j\in\mathcal{X}$. Then

$\displaystyle\phantom{{}={}}\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1\right|X_{t}=k\right)$
$\displaystyle\leq\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{h_{k}(\mathbf{y}_{t})}{\max_{j\neq
k}h_{j}(\mathbf{y}_{t})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}_{t}).$

The effect of the copulas in the HMM can be analyzed more carefully through
the lens of the integral in 2.1; we have the following result for the
bivariate case that explicitly relates the integral to the strength of the
dependence within any particular state. Below, we write $\tau_{k}$ and
$\rho_{k}$ for Kendall’s $\tau$ and Spearman’s $\rho$ of the copula in the
$k$’th state, respectively.

######  Proposition 2.2.

Let $d=2$, and fix $k\in\mathcal{X}$. Then, as either $|\tau_{k}|\to 1$ or
$|\rho_{k}|\to 1$ while the other state-dependent copulas stay fixed, we have

$\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{\omega_{k}\cdot
h_{k}(\mathbf{y})}{\max_{j\neq k}\omega_{j}\cdot
h_{j}(\mathbf{y})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})\to 0$

for any $(\omega_{1},\ldots,\omega_{K})\in\mathcal{S}^{K-1}$ with
$\omega_{k}>0$.

In effect, Theorem 2.1, in combination with 2.1 and 2.2, says that as the
copula in a particularly “common” state approaches either of the Fréchet-
Hoeffding bounds, the observations produced by that state will become
sufficiently distinguished for the local decoding algorithm to detect the
state with complete accuracy, thereby eliminating its contribution to the
zero-one loss. Proofs of these results are provided in Subsection A.2.

##  3 Parameter Estimation

###  3.1 Pitfalls of the Vanilla EM Algorithm

We consider a parametric approach for the copula specification, in which the
parameter vector characterizes the state-dependent copulas and the marginal
distributions. Parameter estimation for HMMs usually relies on the Baum-Welch
algorithm [Baum et al., 1970], which is a particular instance of the
expectation-maximization (EM) algorithm [Dempster et al., 1977], with the
sequence of hidden states playing the role of missing data. For the model
proposed here, the E-step calculations are straightforward when all state-
dependent densities can be evaluated at the current parameter estimates; alas,
the maximization required to complete the M-step is considerably more
difficult because an exact solution is usually unavailable in closed form, and
one must then resort to numerical optimization. The task requires the
maximization of a nonlinear function of all state-dependent parameters. Copula
densities are particularly challenging to optimize directly due to the
algebraic complexities of their associated score functions. In the case of
Archimedean copulas, for example, a naive symbolic differentiation of the log-
likelihood function produces a number of terms that increase rapidly with the
size of each trajectory [Hofert et al., 2012]; this renders the maximization
of all parameters of the objective function computationally expensive. To
illustrate, we present here a succinct version of the relevant mathematical
equations; complete derivations may be found in Zucchini et al. [2017], for
example.

The set of parameters in our model consists of the initial probability mass
distribution $\\{\pi_{i}\\}_{i\in\mathcal{X}}$ and the transition
probabilities $\\{\gamma_{i,j}\\}_{i,j\in\mathcal{X}}$ of the latent Markov
chain, the parameters of the state-dependent marginal densities
$\\{\lambda_{i,h}\\}_{i\in\mathcal{X},h\in[d]}$, and the copula parameters
$\\{\theta_{i}\\}_{i\in\mathcal{X}}$. We write $\bm{\eta}\in\mathbb{R}^{p}$
for the complete collection of these unknown parameters arranged in a vector,
and $\bm{\eta}_{k}$ for the sub-vector of parameters associated with state
$k\in\mathcal{X}$ specifically. Our formulation here of the basic Baum-Welch
algorithm follows Zucchini et al. [2017]. The observed data at time $t$ is
denoted $\mathbf{y}_{t}$. Note that while our presentation deals with a single
trajectory $\\{(\mathbf{y}_{t},X_{t})\\}_{t}$, the algorithm can easily be
extended to handle multiple trajectories
$\\{\\{(\mathbf{y}_{t,1},X_{t,1})\\}_{t=1}^{T_{1}},\ldots,\\{(\mathbf{y}_{t,n},X_{t,n})\\}_{t=1}^{T_{n}}\\}$,
provided that they are independent and identically distributed (iid).

The complete-data log-likelihood for one trajectory of the copula HMM is given
by

$\displaystyle\ell_{\text{com}}\left(\bm{\eta}\mid\mathbf{y}_{1:T},X_{1:T}\right)$
$\displaystyle=\pi_{X_{1}}+\sum_{t=2}^{T}\log\gamma_{X_{t-1},X_{t}}$
$\displaystyle+\sum_{h=1}^{d}\log f_{X_{t},h}(y_{t,h};\lambda_{X_{t},h})$
$\displaystyle+\sum_{t=1}^{T}\log
c_{X_{t}}\left(F_{X_{t},1}(y_{t,1};\lambda_{X_{t},1}),\ldots,F_{X_{t},1}(y_{t,d};\lambda_{X_{t},d})\mid\theta_{X_{t}}\right).$
(4)

We define the latent indicators

$\displaystyle U_{k,t}$ $\displaystyle=\mathbbm{1}_{X_{t}=k},\quad
k\in\mathcal{X},\,t=1,\ldots,T$ $\displaystyle V_{j,k,t}$
$\displaystyle=\mathbbm{1}_{X_{t-1}=j,X_{t}=k},\quad
j,k\in\mathcal{X},\,t=2,\ldots,T,$

so that the complete data log-likelihood in (4) can be rewritten as

$\displaystyle\ell_{\text{com}}\left(\bm{\eta}\mid\mathbf{y}_{1:T},X_{1:T}\right)$
$\displaystyle=\sum_{k=1}^{K}U_{k,1}\cdot\log{\pi_{k}}+\sum_{j=1}^{K}\sum_{k=1}^{K}\left(\sum_{t=2}^{T}V_{j,k,t}\right)\log{\gamma_{j,k}}$
$\displaystyle+\sum_{k=1}^{K}\sum_{t=1}^{T}U_{k,t}\cdot\left(\sum_{h=1}^{d}\log{f_{k,h}(y_{t,h};\lambda_{k,h})}\right)$
$\displaystyle+\log{c_{k}\Big{(}F_{k,1}(y_{t,1};\lambda_{k,1}),\ldots,F_{k,d}(y_{t,d};\lambda_{k,d})\,\Big{|}\,\theta_{k}\Big{)}}.$
(5)

Now consider the E- and M-steps at the $(s+1)$th iteration of the algorithm,
assuming that the algorithm is initialized at $\bm{\eta}^{(0)}$. The E-step
requires the computation of the conditional expectation

$Q(\bm{\eta}\mid\bm{\eta}^{(s)})=\mathbb{E}_{\bm{\eta}^{(s)}}\left[\ell_{\text{com}}\left(\bm{\eta}\mid\mathbf{y}_{1:T},X_{1:T}\right)\right],$
(6)

and to do so requires computing the conditional expectations of all
$U_{k,t}$’s and $V_{j,k,t}$’s given the observed data $\mathbf{y}_{1:T}$ and
current parameter estimates $\bm{\eta}^{(s)}$, which are given by

$\hat{u}_{j,t}^{(s+1)}=\mathbb{P}_{\bm{\eta}^{(s)}}\left(X_{t}=j\mid\mathbf{y}_{1:T}\right)=\frac{\alpha_{j,t}(\mathbf{y}_{1:t};\bm{\eta}^{(s)})\cdot\beta_{j,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}{\sum_{l=1}^{K}\alpha_{l,t}(\mathbf{y}_{1:t};\bm{\eta}^{(s)})\cdot\beta_{l,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}$
(7)

and

$\hat{v}_{j,k,t}^{(s+1)}=\mathbb{P}_{\bm{\eta}^{(s)}}\left(X_{t-1}=j,X_{t}=k\mid\mathbf{y}_{1:T}\right)=\frac{\alpha_{j,t-1}(\mathbf{y}_{1:(t-1)};\bm{\eta}^{(s)})\cdot{\gamma}_{j,k}^{(s)}\cdot
h_{k}(\mathbf{y}_{t};\bm{\eta}^{(s)})\cdot\beta_{k,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}{\sum_{l=1}^{K}\alpha_{l,t}(\mathbf{y}_{1:t};\bm{\eta}^{(s)})\cdot\beta_{l,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})},$
(8)

where

$\alpha_{j,t}(\mathbf{y}_{1:t};\bm{\eta})=\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:t}=\mathbf{y}_{1:t},X_{t}=j\right)\quad\text{and}\quad\beta_{j,t}(\mathbf{y}_{(t+1):T};\bm{\eta})=\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{(t+1):T}=\mathbf{y}_{(t+1):T}\mid
X_{t}=j\right)$

are known as the _forward probabilities_ and _backward probabilities_ ,
respectively. These ingredients are essential for HMM classification
algorithms, and are well-known to admit recursive structures which allow for
their efficient computation via dynamic programming. The latter fact is not
altered by the addition of the copula density in (5).

The M-step requires finding the maximizer $\bm{\eta}^{(s+1)}$ of

$\displaystyle Q(\bm{\eta}|\bm{\eta}^{(s)})$
$\displaystyle=\sum_{k=1}^{K}\hat{u}_{k,1}^{(s)}\cdot\log{\pi_{k}}$
$\displaystyle+\sum_{j=1}^{K}\sum_{k=1}^{K}\left(\sum_{t=2}^{T}\hat{v}_{j,k,t}^{(s)}\right)\cdot\log{\gamma_{j,k}}$
$\displaystyle+\sum_{k=1}^{K}\sum_{t=1}^{T}\hat{u}_{k,t}^{(s)}\cdot\left(\log{c_{k}\Big{(}F_{k,1}(y_{t,1}\mid\lambda_{k,1}),\ldots,F_{k,d}(y_{t,d}\mid\lambda_{k,d})\,\Big{|}\,\theta_{k}\Big{)}}+\sum_{h=1}^{d}\log{f_{k,h}(y_{t,h}\mid\lambda_{k,h})}\right).$
(9)

Since the parameters in the three sums of (9) are functionally independent,
they can be maximized independently. Using Lagrange multipliers, one can
easily show that the maximizer for the initial distribution (as a vector in
the standard simplex $\mathcal{S}^{K-1}$) is

$\bm{\pi}^{(s+1)}=\underset{\bm{\pi}\in\mathcal{S}^{K-1}}{\mathrm{argmax}}\,\left(\sum_{k=1}^{K}\hat{u}_{k,1}^{(s)}\cdot\log{\pi_{k}}\right)=(\hat{u}_{1,1}^{(s)},\ldots,\hat{u}_{K,1}^{(s)}),$

while that for the vector of transition probabilities from state
$j\in\mathcal{X}$ is

$\bm{\gamma}_{j,\cdot}^{(s+1)}=\underset{\bm{\gamma}\in\mathcal{S}^{K-1}}{\mathrm{argmax}}\,\left(\sum_{k=1}^{K}\left(\sum_{t=2}^{T}\hat{v}_{j,k,t}^{(s)}\right)\cdot\log{\gamma_{k}}\right)=\left(\frac{\sum_{t=2}^{T}\hat{v}_{j,1,t}^{(s)}}{\sum_{k=1}^{K}\sum_{t=2}^{T}\hat{v}_{j,k,t}^{(s)}},\ldots,\frac{\sum_{t=2}^{T}\hat{v}_{j,K,t}^{(s)}}{\sum_{k=1}^{K}\sum_{t=2}^{T}\hat{v}_{j,k,t}^{(s)}}\right).$

The maximizer for the vector of parameters involved in the state-dependent
distributions consists of the marginal parameters
$\\{\lambda_{i,h}\\}_{i\in\mathcal{X},h\in[d]}$ and copula parameters
$\\{\theta_{i}\\}_{i\in\mathcal{X}}$ that jointly maximize

$\sum_{k=1}^{K}\sum_{t=1}^{T}\hat{u}_{k,t}^{(s-1)}\cdot\left[\log
c_{k}\Big{(}F_{k,1}(y_{t,1};\lambda_{k,1}),\ldots,F_{k,d}(y_{t,d};\lambda_{k,d})\,\Big{|}\,\theta_{k}\Big{)}+\sum_{h=1}^{d}\log
f_{k,h}(y_{t,h};\lambda_{k,h})\right].$ (10)

At this point, the M-step falters because the optimization of (10) is
difficult; when the dimension $d$ is even moderately high, naively applying a
numerical optimizer to (10) directly is likely to fail, for several reasons.
For one, initialization can be challenging; since any copula is grounded, even
a _correct_ initialization $\lambda_{k,h}^{(0)}=\lambda_{k,h}$ can led us
astray, for we can easily have $F_{k,h}(y_{t,h};\lambda_{k,h})\approx 0$ when
$X_{t}\neq k$, and hence
$c_{k}(\cdots,F_{k,h}(y_{t,h};\lambda_{k,h}),\cdots)\approx 0$. Unless
$\hat{u}_{k,t}^{(0)}\approx 0$ as well — which would be unusual when we have
no _a priori_ information about the true underlying states — evaluation of
(10) will immediately cause a numerical overflow. We emphasize that this can
occur even if all parameters are initialized to their true values.

A second challenge comes from the fact that the log-likelihood surface can
have wide regions of flatness; when Hessian-based methods (or quasi-Newton
methods) are used for optimization, this can cause instability. For example,
we consider a simple version of the model (1) in which only the copula
parameter contains information about the hidden state
$\mathbf{Y}_{t}\mid(X_{t}=k)\sim
C_{\mathrm{Frank}}\left(\mathcal{N}(\mu,1),\mathcal{N}(\mu,1)\mid\theta_{k}\right)$
for $k=1,2$. Figure 3 shows three log-likelihood surfaces, anchored by three
values of $\mu$. One can see that for moderate values of $\theta_{1}$ and
$\theta_{2}$, the log-likelihood surface is quite flat.

Figure 3: Log-likelihood surfaces based on the HMM
$\mathbf{Y}_{t}\mid(X_{t}=k)\sim
C_{\mathrm{Frank}}\left(\mathcal{N}(\mu,1),\mathcal{N}(\mu,1)\mid\theta_{k}\right)$
for $k=1,2$, with true parameters $\mu=0$, and $\theta_{1}=-\theta_{2}=10$ and
all initial probabilities and transition probabilities equal to $1/2$.

###  3.2 An Improved Algorithm

In contrast to jointly maximizing all parameters in (10), it is considerably
less challenging (if not necessarily easy) to estimate the copula parameters
$\\{\theta_{i}\\}_{i\in\mathcal{X}}$ and the marginal parameters
$\\{\lambda_{i,h}\\}_{i\in\mathcal{X},h\in[d]}$ _separately_. For instance,
likelihood inference for the five classical one-parameter Archimedean families
under _known_ marginals was studied by Hofert et al. [2012], who derived
concise functional representations for the copula densities and their score
functions. Altogether, these allow for maximum likelihood estimation using
gradient-based optimization methods. Meanwhile, the difficulty of performing
likelihood inference for the marginal distributions themselves depends highly
on their specific forms (as for the standard case of HMMs with univariate
state-dependent distributions). While closed-form solutions — or efficient
estimation procedures — for MLEs in the presence of iid data have been long-
established for most commonly encountered distributions, the situation is
different for the models we introduce here, where the function to be optimized
is essentially a weighted sum of log-densities, and closed-form expressions
for the maximizing values are rarely encountered.

If the HMM structure were removed from the model so that the data were assumed
to consist of iid observations from
$C(F_{1}(\cdot\,;\lambda_{1}),\ldots,F_{d}(\cdot\,;\lambda_{d})\big{|}\theta)$,
then (10) would reduce to

$\sum_{t=1}^{T}\left\\{\log{c\Big{(}F_{1}(y_{t,1}\,;\lambda_{1}),\ldots,F_{d}(y_{t,d}\,;\lambda_{d})\,\Big{|}\,\theta\Big{)}}+\log{f_{h}(y_{t,h};\lambda_{h})}\right\\}.$
(11)

Moreover, if the copula were known, an estimate for the marginal parameter
$\lambda_{h}$ could be inferred from the marginal data $y_{1:T,h}$ alone as
the standard MLE

$\tilde{\lambda}_{h}=\operatorname*{arg\,sup}_{\lambda}\,\sum_{t=1}^{T}\log{f_{h}(y_{t,h};\lambda)}.$
(12)

This is exactly what is done in the first step of the _inference functions for
margins (IFM)_ method of Joe and Xu [1996], which is itself grounded on the
theory of inference functions [McLeish and Small, 1988]. The IFM method starts
by finding the $\tilde{\lambda}_{h}$’s that separately maximize the marginal
log-likelihoods (12) for $h=1,\ldots,d$, and then proceeds by finding the
$\tilde{\theta}$ that maximizes the joint log-likelihood (11) in which the
marginal parameters are set to their previous estimates:

$\tilde{\theta}=\operatorname*{arg\,sup}_{\theta}\,\sum_{t=1}^{T}\log{c\Big{(}F_{1}(y_{t,1}\,;\tilde{\lambda}_{1}),\ldots,F_{d}(y_{t,d}\,;\tilde{\lambda}_{d})\,\Big{|}\,\theta\Big{)}}.$
(13)

Section 10.1 of Joe [1997] shows that under standard regularity conditions,
$(\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{d},\tilde{\theta})$ is a
consistent and asymptotically normal estimator of
$(\lambda_{1},\ldots,\lambda_{d},\theta)$, although it is not as efficient as
the usual maximum likelihood estimator.

We propose here a new algorithm, which we call the _expectation-IFM (EIFM)
algorithm_ , that alleviates the computational difficulties associated with
the estimation of $\bm{\eta}$. The main idea is to embed the IFM method within
the maximization step of the basic EM algorithm. Specifically, while the
E-step remains unchanged, in the M-step we replace the single-state log-
likelihoods in (12) and (13) with weighted averages that are computed over all
$K$ states. The log-likelihood for the $k$’th state is weighted by the
corresponding state membership probability estimates $\\{\hat{u}_{k,t}\\}$.
Starting with an initial guess $\bm{\eta}^{(0)}$, we perform the following
steps, iterating over $s\geq 1$ until we reach convergence:

  1. Step 1:

(E-Step) Calculate the conditional expectations:

    1. Step 1a:

Calculate

$\hat{u}^{(s+1)}_{j,t}=\frac{\alpha_{j,t}(\mathbf{y}_{1:t};\bm{\eta}^{(s)})\cdot\beta_{j,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}{\sum_{l=1}^{K}\alpha_{l,t}(\mathbf{y}_{1:t};\bm{\eta}^{(s)})\cdot\beta_{l,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}$

for $j\in\mathcal{X}$ and $t\in\\{1,\ldots,T\\}$.

    2. Step 1b:

Calculate

$\hat{v}^{(s+1)}_{j,k,t}=\frac{\alpha_{j,t-1}(\mathbf{y}_{1:(t-1)};\bm{\eta}^{(s)})\cdot{\gamma}_{j,k}^{(s)}\cdot
h_{k}(\mathbf{y}_{t};\bm{\eta}^{(s)})\cdot\beta_{k,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}{\sum_{l=1}^{K}\alpha_{l,t}(\mathbf{y}_{1:t};\bm{\eta}^{(s)})\cdot\beta_{l,t}(\mathbf{y}_{(t+1):T};\bm{\eta}^{(s)})}$

for $j,k\in\mathcal{X}$ and $t\in\\{2,\ldots,T\\}$.

  2. Step 2:

(IFM-Step) Estimate model parameters using an IFM approach:

    1. Step 2a:

Estimate the initial distribution and transition probabilities:

$\bm{\pi}^{(s+1)}=(\hat{u}_{1,1}^{(s)},\ldots,\hat{u}_{K,1}^{(s)})$

and

$\bm{\gamma}_{j,\cdot}^{(s+1)}=\left(\frac{\sum_{t=2}^{T}\hat{v}_{j,1,t}^{(s)}}{\sum_{k=1}^{K}\sum_{t=2}^{T}\hat{v}_{j,k,t}^{(s)}},\ldots,\frac{\sum_{t=2}^{T}\hat{v}_{j,K,t}^{(s)}}{\sum_{k=1}^{K}\sum_{t=2}^{T}\hat{v}_{j,k,t}^{(s)}}\right)$

for $j\in\mathcal{X}$.

    2. Step 2b:

Estimate the marginal parameters

$\lambda_{k,h}^{(s+1)}=\operatorname*{arg\,sup}_{\lambda}\sum_{t=1}^{T}\hat{u}_{k,t}^{(s+1)}\cdot\log{f_{k,h}(y_{t,h};\lambda)}$

for $k\in\\{1,\ldots,K\\}$ and $h\in\\{1,\ldots,d\\}$.

    3. Step 2c:

Estimate the copula parameters

$\tilde{\theta}_{k}^{(s+1)}=\operatorname*{arg\,sup}_{\theta}\sum_{t=1}^{T}\hat{u}_{k,t}^{(s+1)}\cdot\log{c_{k}\Big{(}F_{k,1}(y_{t,1};\lambda_{k,1}^{(s+1)}),\ldots,F_{k,d}(y_{t,d};\lambda_{k,d}^{(s+1)})\,\Big{|}\,\theta\Big{)}}$
(14)

for $k\in\\{1,\ldots,K\\}$.

Step 2b is an importance-weighted maximum likelihood estimation; it is
typically not more difficult than maximum likelihood estimation in the
unweighted case, and the optimal parameter values are sometimes available in
closed form. For example, if the marginal distribution is in an exponential
family with vector parameter $\bm{\lambda}\in\Theta\subseteq\mathbb{R}^{q}$,
that is,

$f_{k,h}(y_{t,h};\bm{\lambda})=g(y_{t,h})\cdot\hskip
1.0pt\mathrm{exp}\left(\bm{\eta}(\bm{\lambda})^{\top}\mathbf{T}(y_{t,h})-A(\bm{\eta})\right),$

it is straightforward to show that Step 2b amounts to solving the system

$\mathbf{J}_{\bm{\eta}}(\bm{\lambda})^{\top}\nabla
A(\bm{\lambda})=\frac{1}{\sum_{t=1}^{T}\hat{u}^{(s)}_{k,t}}\sum_{t=1}^{T}\hat{u}^{(s)}_{k,t}\mathbf{J}_{\bm{\eta}}(\bm{\lambda})^{\top}\mathbf{T}(y_{t,h})=\mathbf{0}$
(15)

in $\bm{\lambda}$. Here $\mathbf{J}_{\bm{\eta}}(\bm{\lambda})$ is the Jacobian
of $\bm{\eta}(\bm{\lambda})$, which is simply the identity matrix when the
exponential family is in canonical form.

The difficulty of Step 2c depends on the ease of evaluating the copula density
$c_{k}$ and its derivative with respect to $\theta$. For the five “classical”
one-parameter Archimedean copula families investigated by Hofert et al.
[2012], the functional forms for the densities, generator derivatives, and
score functions can be computed using the copula package [Hofert et al., 2020]
within R. These functions can be used in conjunction with the built-in optim
function or an appropriate numerical optimization library to approximate (14)
in a straightforward manner.

##  4 Theoretical Results

In this section, we demonstrate that the EIFM algorithm does not belong to the
class of generalized EM algorithms [Dempster et al., 1977], and therefore a
theoretical analysis is needed to justify its use. We start off with a purely
algorithmic perspective in Subsection 4.1, describing the circumstances under
which the sequence of estimates produced by the EIFM algorithm will converge.
In Subsection 4.2, we show that under relatively weak assumptions, the
estimator produced by the algorithm upon convergence will have desirable
statistical properties.

###  4.1 Analysis of the Algorithm

In order to justify the convergence of the EIFM algorithm, we temporarily set
aside its statistical content and view it as a deterministic nonlinear
optimization procedure. The algorithm is not an EM algorithm in the
traditional sense, because the expectation of the complete-data log-likelihood
is not maximized at each iteration. This is mainly a consequence of the fact
that

$\sum_{t=1}^{T}\hat{u}_{t}\cdot\log{f_{h}(y_{t,h};\lambda_{h}^{(s)})}\leq\sum_{t=1}^{T}\hat{u}_{t}\cdot\log{f_{h}(y_{t,h};\lambda_{h}^{(s+1)})},\qquad
h\in\\{1,\ldots,d\\}$ (16)

does not imply

$\sum_{t=1}^{T}\hat{u}_{t}\cdot\log{c\Big{(}F_{1}(y_{t,1};\lambda^{(s)}_{1}),\ldots,F_{d}(y_{t,d};\lambda^{(s)}_{d})\,\Big{|}\,\theta^{(s)}\Big{)}}\leq\sum_{t=1}^{T}\hat{u}_{t}\cdot\log{c\Big{(}F_{1}(y_{t,1};\lambda^{(s+1)}_{1}),\ldots,F_{d}(y_{t,d};\lambda^{(s+1)}_{d})\,\Big{|}\,\theta^{(s)}\Big{)}}.$

For a simple counterexample, one can take a single observation (i.e., $T=1$)
with $Y_{1}\sim\mathrm{Exp}(\lambda_{1})$ and
$Y_{2}\sim\mathrm{Exp}(\lambda_{2})$ such that $(F_{1}(Y_{1}),F_{2}(Y_{2}))$
has a Farlie–Gumbel–Morgenstern copula with density
$c(u,v\mid\theta)=1+\theta(2u-1)(2v-1)$. Suppose the algorithm were
initialized at $\lambda_{1}^{(0)}=\lambda_{2}^{(0)}=1$ and any
$\theta^{(0)}\in(0,1)$. Step 2b of the algorithm produces the standard
univariate MLEs $\lambda_{1}^{(1)}=y_{1}^{-1}$ and
$\lambda_{2}^{(1)}=y_{2}^{-1}$, which satisfy (16) by construction; however,
the inequality

$\log{c\left(F_{1}(y_{1};\lambda_{1}^{(0)}),F(y_{2};\lambda_{2}^{(0)})\mid\theta^{(0)}\right)}\leq\log{c\left(F_{1}(y_{1};\lambda_{1}^{(1)}),F(y_{2};\lambda_{2}^{(1)})\mid\theta^{(0)}\right)}$

is equivalent to

$(2e^{-y_{1}}-1)(2e^{-y_{2}}-1)\leq(2e^{-1}-1)^{2},$

which immediately fails when $y_{1},y_{2}>1$.

The above counterexample shows that despite sharing a similar strategy with
the ECM algorithm [Meng and Rubin, 1993], the EIFM algorithm cannot be
justified in the same way as any _generalized EM algorithm_ [Dempster et al.,
1977], as these all rely on increasing the objective function at each
iteration. Clearly, one needs to identify the conditions under which the EIFM
algorithm will converge, since its sequential updating rule does not offer
convergence guarantees without stronger assumptions. Our analysis here follows
roughly the one given for the _ES algorithm_ by Elashoff and Ryan [2004], as
the aim of that procedure is also to iteratively solve unbiased estimating
equations in the presence of missing data. However, the dependence inherent in
the HMM data precludes a direct application of the ES algorithm.

We regard the $\hat{u}_{j,t}$’s and $\hat{v}_{j,k,t}$’s as parameters
themselves (in a non-statistical context), and collect them and the parameter
of interest $\bm{\eta}$ into a larger parameter vector
$\bm{\xi}=(\hat{\mathbf{u}},\hat{\mathbf{v}},\bm{\eta})$ of length
$N:=TK+(T-1)K^{2}+K+K^{2}+dK+K$. For accounting purposes, let

$\displaystyle N_{1}$ $\displaystyle=0$ $\displaystyle N_{2}$
$\displaystyle=N_{1}+TK$ $\displaystyle N_{3}$
$\displaystyle=N_{2}+(T-1)K^{2}$ $\displaystyle N_{4}$ $\displaystyle=N_{3}+K$
$\displaystyle N_{5}$ $\displaystyle=N_{4}+K^{2}$ $\displaystyle N_{6}$
$\displaystyle=N_{5}+dK.$

We construct a vector-valued function
$\mathbf{g}(\bm{\xi})=\left(g_{1}(\bm{\xi}),\ldots,g_{N}(\bm{\xi})\right)^{\top}$
with each $g_{i}:\mathbb{R}^{N}\to\mathbb{R}$ defined as follows:

$\displaystyle g_{N_{1}+jt}(\bm{\xi})$
$\displaystyle=\hat{u}_{j,t}-\frac{\alpha_{j,t}(\mathbf{y}_{1:t};\bm{\eta})\cdot\beta_{j,t}(\mathbf{y}_{(t+1):T};\bm{\eta})}{\sum_{l=1}^{K}\alpha_{l,t}(\mathbf{y}_{1:t};\bm{\eta})\cdot\beta_{l,t}(\mathbf{y}_{(t+1):T};\bm{\eta})},$
$\displaystyle j\in\mathcal{X},\,\,t\in\\{1,\ldots,T\\}$ $\displaystyle
g_{N_{2}+jkt}(\bm{\xi})$
$\displaystyle=\hat{v}_{j,k,t}-\frac{\alpha_{j,t-1}(\mathbf{y}_{1:(t-1)};\bm{\eta})\cdot{\gamma}_{j,k}\cdot
h_{k}(\mathbf{y}_{t};\bm{\eta})\cdot\beta_{k,t}(\mathbf{y}_{(t+1):T};\bm{\eta})}{\sum_{l=1}^{K}\alpha_{l,t}(\mathbf{y}_{1:t};\bm{\eta})\cdot\beta_{l,t}(\mathbf{y}_{(t+1):T};\bm{\eta})},$
$\displaystyle j,k\in\mathcal{X},\,\,t\in\\{2,\ldots,T\\}$ $\displaystyle
g_{N_{3}+j}(\bm{\xi})$ $\displaystyle=\pi_{j}-\hat{u}_{j,1},$ $\displaystyle
j\in\mathcal{X}$ $\displaystyle g_{N_{4}+jk}(\bm{\xi})$
$\displaystyle=\gamma_{j,k}\cdot\sum_{l=1}^{K}\sum_{t=1}^{T}\hat{v}_{j,l,t}-\sum_{t=1}^{T}\hat{v}_{j,k,t},$
$\displaystyle j,k\in\mathcal{X}$ $\displaystyle g_{N_{5}+jh}(\bm{\xi})$
$\displaystyle=\sum_{t=1}^{T}\hat{u}_{j,t}\cdot\frac{\partial}{\partial\lambda}\log{f_{j,h}(y_{t,h};\lambda)},$
$\displaystyle j\in\mathcal{X},\,h\in\\{1,\ldots,d\\}$ $\displaystyle
g_{N_{6}+j}(\bm{\xi})$
$\displaystyle=\sum_{t=1}^{T}\hat{u}_{j,t}\cdot\frac{\partial}{\partial\theta}\log{c_{j}\left(F_{j,1}(y_{t,1};\lambda_{j,1}),\ldots,F_{j,d}(y_{t,d};\lambda_{j,d})\mid\theta\right)},$
$\displaystyle j\in\mathcal{X}.$

From the description of the EIFM algorithm in Subsection 3.2, it can be seen
that the entire $(s+1)$’th iteration of the algorithm exactly corresponds to
updating $\bm{\xi}^{(s)}\to\bm{\xi}^{(s+1)}$ by setting $\xi_{i}^{(s+1)}$ as
the solution to the univariate problem

$g_{i}(\xi^{(s+1)}_{1},\ldots,\xi_{i-1}^{(s+1)},\xi,\xi^{(s)}_{i+1},\ldots,\xi^{(s)}_{N})=0,$
(17)

for each $i=1,\ldots,N$ (in practice, many of these sub-updates are performed
in parallel — for example, the sub-updates $N_{4}+1$ to $N_{5}$ altogether
correspond to Step 2a of the algorithm). Supposing the sequence
$\\{\bm{\xi}^{(s)}\\}_{s\geq 1}$ converges to some
$\bm{\xi}^{*}\in\mathbb{R}^{N}$, the limiting vector will satisfy
$\mathbf{g}(\bm{\xi}^{*})=\mathbf{0}$, and the sub-vector $\bm{\eta}^{*}$ will
be taken as our estimator of $\bm{\eta}$, whose statistical properties are
studied in Subsection 4.2. The EIFM algorithm is thus an example of a
_nonlinear Gauss-Seidel method_ , or more generally, a _nonlinear successive
over-relaxation (SOR) method_ [Ortega and Rheinboldt, 2000]. The local
convergence of such algorithms depend on the behaviour of the Jacobian
$\mathbf{J}_{\mathbf{g}}(\bm{\xi})$ of $\mathbf{g}$ in a neighborhood of some
solution to $\mathbf{g}(\bm{\xi})=\mathbf{0}$ such that a certain
transformation $\tilde{\mathbf{g}}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ of
$\mathbf{g}$ is (locally) required to be a contraction mapping with
$\bm{\xi}^{*}$ as a fixed point. The critical theorem is the following:

######  Theorem 4.1 (Ortega and Rheinboldt [2000], Theorem 10.3.5).

Let $\mathbf{g}:D\subset\mathbb{R}^{N}\to\mathbb{R}^{N}$ be continuously
differentiable in an open neighborhood $S_{0}$ of $\bm{\xi}^{*}$ such that
$\mathbf{g}(\bm{\xi}^{*})=\mathbf{0}$. Decompose the Jacobian into diagonal,
strictly lower, and strictly upper matrices as
$\mathbf{J}_{\mathbf{g}}(\bm{\xi})=\mathbf{D}(\bm{\xi})-\mathbf{L}(\bm{\xi})-\mathbf{U}(\bm{\xi})$,
and suppose that $\mathbf{D}(\bm{\xi}^{*})$ is nonsingular. If
$\rho([\mathbf{D}(\bm{\xi}^{*})-\mathbf{L}(\bm{\xi}^{*})]^{-1}\mathbf{U}(\bm{\xi}^{*}))<1$,
then there exists an open ball $B(\bm{\xi}^{*},\delta)\subset S_{0}$ such that
for any $\bm{\xi}^{(0)}\in B(\bm{\xi}^{*},\delta)$, there exists a unique
sequence $\\{\bm{\xi}^{(s)}\\}\subset B(\bm{\xi}^{*},\delta)$ satisfying the
nonlinear SOR prescription, such that
$\lim_{s\to\infty}\bm{\xi}^{(s)}=\bm{\xi}^{*}$.

See Ortega and Rheinboldt [2000] for a proof. Note here that for a general
matrix $\mathbf{A}\in\mathbb{C}^{N\times N}$ with eigenvalues
$\nu_{1},\ldots,\nu_{N}$, the function $\rho(\mathbf{A})$ refers to the
_spectral radius_ of $\mathbf{A}$, defined as
$\max\\{|\nu_{1}|,\ldots,|\nu_{N}|\\}$. This theorem implies that the EIFM
algorithm, once sufficiently close to a local solution $\bm{\xi}^{*}$ of
$\mathbf{g}(\bm{\xi})=\mathbf{0}$, will produce a sequence converging
superlinearly to $\bm{\xi}^{*}$, provided that the spectral radius of the
matrix
$[\mathbf{D}(\bm{\xi}^{(s)})-\mathbf{L}(\bm{\xi}^{(s)})]^{-1}\mathbf{U}(\bm{\xi}^{(s)}))$
is less than 1. Verifying the nonsingularity of $\mathbf{D}(\bm{\xi}^{(s)})$
is straightforward; it is easy to see that the first $N_{5}$ diagonal entries
are equal to 1, while the remaining entries are of the form

$\left.\sum_{t=1}^{T}\hat{u}^{(s)}_{j,t}\cdot\frac{\partial^{2}}{\partial\lambda^{2}}\log{f_{j,h}(y_{t,h};\lambda)}\right|_{\lambda=\lambda^{(s)}_{j,h}}\quad\text{or}\quad\left.\sum_{t=1}^{T}\hat{u}^{(s)}_{j,t}\cdot\frac{\partial^{2}}{\partial\theta^{2}}\log{c_{j}\left(F_{j,1}(y_{t,1};\lambda^{(s)}_{j,1}),\ldots,F_{j,d}(y_{t,d};\lambda^{(s)}_{j,d})\mid\theta\right)}\right|_{\theta=\theta^{(s)}_{j}},$

which can be calculated directly. Some of these calculations may be
unnecessary with prior knowledge of the involved densities; for example, the
term on the left is certain to be negative if the mapping $\lambda\mapsto
f_{j,h}(y;\lambda)$ is known to be strictly concave. The remaining elements of
$\mathbf{J}_{\mathbf{g}}$ at $\bm{\xi}^{(s)}$ and the corresponding spectral
radius
$\rho([\mathbf{D}(\bm{\xi}^{(s)})-\mathbf{L}(\bm{\xi}^{(s)})]^{-1}\mathbf{U}(\bm{\xi}^{(s)}))$
can, in principle, be computed exactly following the $s$’th iteration of the
EIFM algorithm, although in practice approximating these quantities
numerically would be considerably easier than explicitly calculating all of
the required derivatives.

###  4.2 Analysis of the Estimator

We now discuss conditions under which the subvector $\bm{\eta}^{*}$ of the
solution $\bm{\xi}^{*}$ ultimately produced by the EIFM algorithm — now
regarded as a statistical estimator $\bm{\eta}^{*}(\mathbf{Y}_{1:T})$ — is
consistent for the true parameter $\bm{\eta}_{0}$ for the model described in
Subsection 2.1 which is assumed to have generated the data. While our analysis
of the algorithm’s convergence in Subsection 4.1 uses the same underlying
technique as Elashoff and Ryan [2004] for their ES algorithm, we cannot
emulate their proof of consistency. Unlike their data, which are assumed to be
iid, thus allowing them to use standard techniques (see Godambe [1960], for
example) to prove the consistency and asymptotic normality of their estimator,
the latent variables in our setup — namely, the conditional state membership
indicators $U_{k,t}$’s and the $V_{j,k,t}$’s — are not independent. Their
dependency is due to their strong connection to the underlying Markov
structure of the HMM. We instead follow the approach of Jensen [2011], who
studied the asymptotic properties of M-estimators used for HMMs.

Specifically, we assume that the EIFM algorithm has produced a vector
$\bm{\xi}^{*}=(\mathbf{u}^{*},\mathbf{v}^{*},\bm{\eta}^{*})$ which satisfies
$\mathbf{g}(\bm{\xi}^{*})=\mathbf{0}$, where $\mathbf{g}$ is as defined in
Subsection 4.1. To show that this estimator has desirable asymptotic
properties under mild regularity conditions, it suffices to show that the
lower $p=K+K^{2}+dK+K$ components of $\mathbf{g}$ — those corresponding to the
FM-step of the algorithm — constitute a system of unbiased estimating
equations for $\bm{\eta}$ conditional on $\mathbf{Y}_{1:T}$ after the latent
variables $X_{1:T}$ have been marginalized out. Specifically, we construct a
new vector-valued function $\bm{\psi}_{T}(\bm{\eta};\mathbf{Y}_{1:T})$ as
follows. First, we define the random $\mathbb{R}$-valued functions

$\displaystyle\phi_{N_{3}+j,t}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)$
$\displaystyle=\begin{cases}\pi_{j}-U_{j,1},&t=1\\\ 0,&t>1\end{cases},$
$\displaystyle j\in\mathcal{X}$ (18)
$\displaystyle\phi_{N_{4}+jk}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)$
$\displaystyle=\gamma_{j,k}\cdot\sum_{l=1}^{K}V_{j,l,t}-V_{j,k,t},$
$\displaystyle j,k\in\mathcal{X}$ (19)
$\displaystyle\phi_{N_{5}+jh}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)$
$\displaystyle=U_{j,t}\cdot\frac{\partial}{\partial\lambda}\log{f_{j,h}(Y_{t,h};\lambda)},$
$\displaystyle j\in\mathcal{X},\,h\in\\{1,\ldots,d\\}$ (20)
$\displaystyle\phi_{N_{6}+j}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)$
$\displaystyle=U_{j,t}\cdot\frac{\partial}{\partial\theta}\log{c_{j}\left(F_{j,1}(Y_{t,1};\lambda_{j,1}),\ldots,F_{j,d}(Y_{t,d};\lambda_{j,d})\mid\theta\right)},$
$\displaystyle j\in\mathcal{X}.$ (21)

Note that as functions, only (18) depends on $t$; this is required in order to
capture the initial distribution parameters $\bm{\pi}$. The remaining
functions are free of $t$ because of the time-homogeneity of the underlying
Markov chain. Now, let

$\bm{\phi}_{t}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)=\left(\phi_{N_{3}+1,t}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right),\ldots,\phi_{N_{3}+p}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\right)^{\top}$

and

$\bm{\psi}_{T}(\bm{\eta};\mathbf{Y}_{1:T})=\sum_{t=1}^{T}\mathbb{E}_{\bm{\eta}}\left[\bm{\phi}_{t}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\mid\mathbf{Y}_{1:T}\right].$

From (7) and (8) and the fact that the algorithm has converged to the fixed
point $\bm{\xi}^{*}$, it is clear that

$u^{*}_{j,t}=\mathbb{P}_{\bm{\eta}^{*}}\left(X_{t}=j\mid\mathbf{Y}_{1:T}\right)\quad\text{and}\quad
v^{*}_{j,k,t}=\mathbb{P}_{\bm{\eta}^{*}}\left(X_{t-1}=j,X_{t}=k\mid\mathbf{Y}_{1:T}\right),$

and from these it is easy to verify that $\bm{\eta}^{*}$ solves
$\bm{\psi}_{T}(\bm{\eta};\mathbf{Y}_{1:T})=\mathbf{0}$. Its use as an
estimator of $\bm{\eta}_{0}$ is justified by the following fundamental fact:

######  Theorem 4.2.

Under standard regularity conditions,
$\bm{\psi}_{T}(\bm{\eta};\mathbf{Y}_{1:T})$ defines an unbiased estimating
equation.

###### Proof.

It suffices to check that each of (18)–(21) have expectations (with respect to
the joint distribution of $\mathbf{Y}_{1:T}$ and $X_{1:T}$) of $0$. Below, all
expectations are taken with respect to a generic $\bm{\eta}\in\Theta$; we
suppress the subscript for ease of notation. To begin with,

$\mathbb{E}\left[\mathbb{E}\left[\phi_{N_{3}+j,1}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\mid\mathbf{Y}_{1:T}\right]\right]=\pi_{j}-\mathbb{P}\left(X_{1}=k\right)=0$

and, similarly,

$\displaystyle\mathbb{E}\left[\mathbb{E}\left[\phi_{N_{4}+jk}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\mid\mathbf{Y}_{1:T}\right]\right]$
$\displaystyle=\gamma_{j,k}\cdot\sum_{l=1}^{K}\mathbb{P}\left(X_{t-1}=j,X_{t}=l\right)-\mathbb{P}\left(X_{t-1}=j,X_{t}=k\right)$
$\displaystyle=\gamma_{j,k}\cdot\mathbb{P}\left(X_{t-1}=j\right)-\mathbb{P}\left(X_{t}=k\mid
X_{t-1}=j\right)\cdot\mathbb{P}\left(X_{t-1}=j\right)$ $\displaystyle=0.$

Moreover, under standard regularity conditions required for the Bartlett
identity

$\mathbb{E}\left[\frac{\partial}{\partial\lambda}\log{f_{k,h}(Y_{t,h};\lambda)}\right]=0$
(22)

to hold, we have

$\mathbb{E}\left[\mathbb{E}\left[\phi_{N_{5}+jh}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\mid
X_{1:T}\right]\right]=\mathbb{E}\left[U_{j,t}\cdot\mathbb{E}\left[\frac{\partial}{\partial\lambda}\log{f_{j,h}(Y_{t,h};\lambda)}\right]\right]=0.$

Likewise, assuming the same sorts of regularity conditions that permit the
interchange of integration and differentiation, we have

$\displaystyle\begin{split}\mathbb{E}\left[\phi_{N_{6}+j}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\right]&=\mathbb{E}\left[\frac{\partial}{\partial\theta}\log{c_{j}\left(\mathbf{F}_{j}(\mathbf{Y}_{t};\bm{\lambda}_{j})\mid\theta\right)}\right]\\\
&=\int_{\mathbb{R}^{d}}\frac{\partial}{\partial\theta}c_{j}\left(\mathbf{F}_{j}(\mathbf{y}_{t};\bm{\lambda}_{j})\mid\theta\right)\prod_{h=1}^{d}f_{j,h}(y_{t,h};\lambda_{j,h})\mathop{}\\!\mathrm{d}\mathbf{y}_{t}\\\
&=\int_{[0,1]^{d}}\frac{\partial}{\partial\theta}c_{j}\left(\mathbf{u}\mid\theta\right)\mathop{}\\!\mathrm{d}\mathbf{u}\\\
&=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}\theta}\int_{[0,1]^{d}}c_{j}\left(\mathbf{u}\mid\theta\right)\mathop{}\\!\mathrm{d}\mathbf{u}\\\
&=0.\end{split}$ (23)

where
$\mathbf{F}_{j}(\mathbf{y}_{t};\bm{\lambda}_{j}):=\left(F_{j,1}(y_{t,1};\lambda_{j,1}),\ldots,F_{j,d}(y_{t,d};\lambda_{j,d})\right)$
and we have made the usual change of variables
$u_{h}=F_{j,h}(y_{t,h};\lambda_{j,h})$ for $h\in\\{1,\ldots,d\\}$. This gives

$\mathbb{E}\left[\mathbb{E}\left[\phi_{N_{6}+j}\left(\bm{\eta};X_{(t-1):t},\mathbf{Y}_{t}\right)\mid
X_{1:T}\right]\right]=\mathbb{E}\left[U_{j,t}\cdot\mathbb{E}\left[\frac{\partial}{\partial\theta}\log{c_{j}\left(\mathbf{F}_{j}(\mathbf{y}_{t};\bm{\lambda}_{j})\mid\theta\right)}\right]\right]=0,$

which finishes the proof. ∎

The IFM method is justified by (22) and (23) [see Subsection 4.1 of Joe and
Xu, 1996]. The solution
$(\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{d},\tilde{\theta})$ of these
equations — at least in the presence of iid data — is shown to be a consistent
and asymptotically normal estimator of the true data-generating parameter as
$T\to\infty$. The asymptotic variance of the estimators is given by the
inverse Godambe information matrix. However, in order to account for the
underlying dependence structure of the HMM, stronger regularity conditions are
necessary for the same conclusion to hold. These conditions are encapsulated
in the findings of Jensen [2011], which establishes a central limit theorem
for the sequence of estimators produced by the EIFM algorithm. The relevant
assumptions are restated in Subsection A.1 using our notation; we present the
main theorem here:

######  Theorem 4.3 (Jensen [2011], Theorem 1).

Let
$\mathbf{G}_{t}=\mathrm{Var}_{\bm{\eta}_{0}}\left(\bm{\psi}_{t}(\bm{\eta}_{0};\mathbf{Y}_{1:t})\right)$
and
$\mathbf{H}_{t}=\mathbb{E}_{\bm{\eta}_{0}}\left[-\mathbf{J}_{\bm{\psi}_{t}(\cdot;\mathbf{Y}_{1:t})}(\bm{\eta}_{0})\right]$.
Under Assumptions 1-4 of Subsection A.1, there exists a consistent sequence
$\\{\hat{\bm{\eta}}_{T}\\}$ solving
$\bm{\psi}_{T}(\bm{\eta};\mathbf{Y}_{1:T})=\mathbf{0}$ such that

$\sqrt{T}\cdot\mathbf{G}_{T}^{-1/2}\mathbf{H}_{T}(\hat{\bm{\eta}}_{T}-\bm{\eta}_{0})\stackrel{{\scriptstyle
d}}{{\longrightarrow}}\mathcal{N}(\mathbf{0},I)$

as $T\to\infty$ under $\bm{\eta}_{0}$.

See Jensen [2011] for a proof.

In all but the most trivial situations, the variances and expectations needed
to calculate $\mathbf{G}_{T}^{-1}$ and, respectively, $\mathbf{H}_{T}$ are
analytically intractable. Therefore, we suggest two alternative approaches for
estimating standard errors and confidence intervals for $\bm{\eta}^{*}$. In
both cases, we start by sampling iid copies
$(\mathbf{Y}_{1:T}^{(1)},X_{1:T}^{(1)}),\ldots,(\mathbf{Y}_{1:T}^{(n)},X_{1:T}^{(n)})$
from the model under $\bm{\eta}^{*}$. The first approach involves estimating
$\mathbf{G}_{T}$ and $\mathbf{H}_{T}$ via Monte Carlo, using $\bm{\eta}^{*}$
as a plug-in estimator for $\bm{\eta}_{0}$. Using the generated samples, we
compute

$\hat{\mathbf{G}}_{T}=\frac{1}{n-1}\sum_{i=1}^{n}\left(\bm{\psi}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T}^{(i)})-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T})\right)\left(\bm{\psi}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T}^{(i)})-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T})\right)^{\top}$

where $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T})=\frac{1}{n}\sum_{i=1}^{n}\bm{\psi}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T}^{(i)})$,
and

$\hat{\mathbf{H}}_{T}=\frac{1}{n}\sum_{i=1}^{n}\left[-\mathbf{J}_{\bm{\psi}_{T}(\cdot;\mathbf{Y}_{1:T}^{(i)})}(\bm{\eta}^{*})\right].$

Note that the entries of each
$\bm{\psi}_{T}(\bm{\eta}^{*};\mathbf{Y}_{1:T}^{(i)})$ can be computed by
replacing the $U_{j,t}$’s and $V_{j,k,t}$’s in (18)–(21) with their
conditional expectations given $\mathbf{Y}_{1:T}$ (explicit forms for these
are provided by A.1). The Jacobian that is required to compute
$\hat{\mathbf{H}}_{T}$ can be challenging to derive using symbolic computation
software.

The second approach relies on the parametric bootstrap estimator of Efron and
Tibshirani [1994], as described by Zucchini et al. [2017] for HMMs.
Specifically, for each $i\in\\{1,\ldots,n\\}$, we re-run the EIFM algorithm on
$(\mathbf{Y}_{1:T}^{(i)},X_{1:T}^{(i)})$ to produce a new estimator
$\hat{\bm{\eta}}^{(i)}$, and then estimate the covariance matrix of
$\bm{\eta}^{*}$ via

$\hat{\mathbf{V}}_{\bm{\eta}^{*}}=\frac{1}{n-1}\sum_{i=1}^{n}\left(\hat{\bm{\eta}}^{(i)}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}\right)\left(\hat{\bm{\eta}}^{(i)}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}\right)^{\top},$

where $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}=\frac{1}{n}\sum_{i=1}^{n}\hat{\bm{\eta}}^{(i)}$. Computation time is
reasonable for both methods, because the required tasks are parallelizable.
Once all the first and second order derivatives are computed, however, the
Monte Carlo method is faster to implement (since it does not require re-
fitting the model), but may be prone to numerical errors when the true
parameter $\bm{\eta}_{0}$ lies near the boundary of $\bm{\Theta}$.

##  5 Simulations and Applications

###  5.1 A Simulation Study

The central aim of our simulation study was to demonstrate that introducing
the copula into the HMM formulation is advantageous when the dependence
structure between observed variables is informative about the latent process,
and that it is not harmful otherwise. Secondly, it provides a proof of concept
for the EIFM algorithm. Finally, we demonstrate the model’s performance under
copula misspecification.

All generative models have 3-state HMMs with the same transition probability
matrix

$\bm{\Gamma}=\begin{pmatrix}0.5&0.25&0.25\\\ 0.25&0.5&0.25\\\ 0.25&0.25&0.5\\\
\end{pmatrix}.$

and initial distribution $\bm{\pi}=(0,1,0)$. The marginal distributions are
Gaussian and the dependence structure is provided by a Frank copula, whose
copula parameter varies between models. Specifically, for $i\in\\{1,2,3,4\\}$,
the $i$’th model has the stochastic representation

$\mathbf{Y}_{t}^{(i)}\mid(X_{t}^{(i)}=k)\sim
C_{\text{Frank}}\left(\mathcal{N}(\mu_{k,1},\sigma=0.5),\mathcal{N}(\mu_{k,2},\sigma=0.5)\mid\theta_{i,k}\right),\quad
k\in\\{1,2,3\\}$

where $\mu_{k,h}=k+3\cdot\mathbbm{1}_{h=2}$ and

$\displaystyle\bm{\theta}_{1}$ $\displaystyle=(30,30,30)$
$\displaystyle\bm{\theta}_{2}$ $\displaystyle=(5,30,30)$
$\displaystyle\bm{\theta}_{3}$ $\displaystyle=(5,5,30)$
$\displaystyle\bm{\theta}_{4}$ $\displaystyle=(5,5,5).$

We refer to each $i$ as a dependence “scenario”; thus, the first scenario
features strong dependence in all states ($\theta=30$, or $\tau\approx
0.874$), while the second features weaker dependence in the first state
($\theta=5$ or $\tau\approx 0.457$), and so on. Note that in the first and
fourth scenarios, the copula parameters within each state-dependent
distribution are identical, and therefore only the marginals are informative
of the underlying states; this contrasts with the numerical example of
Subsection 2.2, in which _only_ the copula informed the hidden states. For the
$i$th scenario, we generated a time series of length $T=500$ from the
corresponding model with true parameter
$\bm{\eta}=(\mu_{1,1},\ldots,\mu_{3,2},\bm{\theta}_{i})$.

First, we assessed the EIFM algorithm’s ability to perform both parameter
estimation and state classification. Note that optimal values for the
$\mu_{k,h}$’s are available in closed form at each iteration of the algorithm
via (15), thereby limiting any potential numerical obstacles in the algorithm
(initializing values aside) to the estimation of the $\theta_{i,k}$’s in Step
2c. To provide the EIFM algorithm with plausible initial parameter values,
$\bm{\eta}_{i}^{(0)}$, for each scenario $i$, we first ran an initial EM
algorithm assuming independent Gaussian state-dependent distributions with
$\mu^{(0)}_{k,h}=(-1)^{h}\cdot|Z_{k,h}|$ and $Z_{k,h}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{N}(0,1)$. Using the parameter estimates from this
initial run, we used local decoding to make the classifications
$\hat{X}^{(i)}_{1},\ldots,\hat{X}^{(i)}_{T}$ and then set

$\hat{\tau}_{i,k}^{(0)}=\frac{2}{|C_{i,k}|\cdot(|C_{i,k}|-1)}\sum_{\begin{subarray}{c}s,t\,\in\,C_{i,k}\\\
s<t\end{subarray}}\mathrm{sgn}\left(y^{(i)}_{1,s}-y^{(i)}_{1,t}\right)\cdot\mathrm{sgn}\left(y^{(i)}_{2,s}-y^{(i)}_{2,t}\right),$

where $C_{i,k}=\\{t:\hat{X}^{(i)}_{t}=k\\}$; we then set
$\theta_{i,k}^{(0)}=\theta_{\text{Frank}}(\hat{\tau}_{i,k}^{(0)})$. In words,
our initial estimate for $\theta_{i,k}$ was based on the empirical Kendall’s
$\tau$ calculated on the subset of observations which was predicted (by the
independent Gaussians model) to have arisen from state $k$. The remaining
parameters, common to both the independent Gaussians model and the true data-
generating model, were initialized in the second stage algorithm as their
corresponding estimates from the initial run.

The entire process was independently replicated 20 times in parallel; we show
boxplots of the resulting parameter estimates in Figure 4. The estimates are
observed to be unbiased, and the copula parameters are generally estimated
well. Figure 5 shows the overall accuracy of the state classifications (via
local decoding) across the four dependence scenarios; label-switching was
addressed by selecting the best permutation of the predicted labels for each
scenario and each simulation, and then averaging the results. Figure 5 also
shows classification accuracy for the state predictions produced by two
misspecified models: one in which the Frank copula is replaced by the
independence copula (which corresponds to the independent Gaussians model used
to generate initial values for the EIFM algorithm), and another in which it is
replaced by a Gauss copula. It is clear that across virtually all states and
scenarios, the accuracy of the true model was equal or better than that of
both misspecified models.

Figure 4: Parameter estimates based on 20 independent simulations and EIFM
algorithm runs for the bivariate 3-state HMMs described in Subsection 5.1. The
$i$’th row corresponds to the $i$’th “scenario”. Colored horizontal lines
indicate true parameter values. Figure 5: Average percent of correct state
classifications for the bivariate 3-state HMMs under the true model (solid
bars) and the misspecified models featuring the independence copula (striped
bars) and the Gauss copula (hatched bars) described in Subsection 5.1.

###  5.2 Occupancy Detection

The ability to detect whether a room is occupied using sensor data (such as
temperature and $CO_{2}$ levels) can potentially reduce unnecessary energy
consumption by automatically controlling HVAC and lighting systems, without
the need for motion detectors or other methods that could constitute an
invasion of privacy [Candanedo and Feldheim, 2016]. To assess the performance
of the copula-within-HMM model and the EIFM algorithm on real-world data, we
consider three publicly-available labelled datasets presented by Candanedo and
Feldheim [2016] consisting of multivariate time series of four environmental
measurements (light, temperature, humidity, $CO_{2}$) and one derived metric
(the humidity ratio $W$, whose construction is described in Candanedo and
Feldheim [2016]) captured in an office room in Belgium, along with binary
indicators for whether the room was occupied at the time of measurement.

In Candanedo et al. [2017] various HMMs were fit to this data via the EM
algorithm using the depmixS4 package [Visser and Speekenbrink, 2010]; the
models featured either one or two (lagged) sensor measurements as the observed
data, and the occupancy indicators as the hidden states. The univariate data
was taken to be normally distributed, and the bivariate data was assumed to
feature independent normal margins. The performance of the models was assessed
by comparing predictions made by local decoding to the ground-truth labels via
the zero-one loss (as well as several other measures commonly used in binary
classification, such as sensitivity and specificity).

We also illustrate our model using bivariate data. We processed the data in
the same manner as Candanedo et al. [2017]: we averaged the data across five-
minute periods, and then calculated the lagged differences to obtain time
series of lengths 1629, 533, and 1951 respectively. We selected as our
observed data the series $\mathbf{Y}_{t}=((CO_{2})_{t},W_{t})$ and chose the
marginal distributions and copulas for each state using several heuristics.
First, we selected several common families of parametric copulas (the Frank,
Clayton, Gumbel, Joe, and Gauss families), and for each we carried out a
goodness-of-fit test based on the pseudo-observations using the multiplier
bootstrap method [Kojadinovic et al., 2011]; we then chose the parametric
family based on the lowest corresponding Cramér-von Mises test statistic. This
process yielded a Clayton copula for State 1 and a Frank copula for State 2
(see Table 1), choices which were corroborated by scatterplots of the pseudo-
observations (see Figure 6). We used a similar criterion to select the
marginal distributions, limiting our choice of distributions to those in
exponential families (with appropriate supports); the Gaussian distribution
gave the best results for both marginal components within both states.

Labelling the unoccupied state as ‘1’ and the occupied state as ‘2’, our model
thus has the stochastic representation

$\displaystyle\mathbf{Y}_{t}\mid(X_{t}=1)$ $\displaystyle\sim
C_{\text{Clayton}}\left(\mathcal{N}(\mu_{1,1},\sigma^{2}_{1,1}),\mathcal{N}(\mu_{1,2},\sigma^{2}_{1,2})\mid\theta_{1}\right)$
$\displaystyle\mathbf{Y}_{t}\mid(X_{t}=2)$ $\displaystyle\sim
C_{\text{Frank}}\left(\mathcal{N}(\mu_{2,1},\sigma^{2}_{2,1}),\mathcal{N}(\mu_{2,2},\sigma^{2}_{2,2})\mid\theta_{2}\right).$

State | Frank | Clayton | Gumbel | Joe | Gauss
---|---|---|---|---|---
1 | 0.356 | 0.255 | 0.423 | 0.770 | 0.345
2 | 0.018 | 0.433 | 0.038 | 0.206 | 0.045
Table 1: Cramér-von Mises test statistics based on pseudo-observations
computed from unoccupied (Row 1) and occupied (Row 2) subsets. Figure 6:
Pseudo-observations computed from unoccupied (Panel 1) and occupied (Panel 2)
subsets.

The parameter estimates produced by the EIFM algorithm (including those for
the initial distribution of the Markov chain $\bm{\pi}$ and the chain’s
transition probability matrix $\bm{\Gamma}$) are given in Table 2, along with
95% bootstrap confidence intervals produced using the method described in
Subsection 4.2. Local decoding was used to classify all three datasets, and
the accuracy of these classifications was assessed via the zero-one loss. To
compare the predictive performance of the copula-within-HMM model with that of
a more basic model, we repeated the experiment using an independence copula.
The performance of both models is shown in Table 3. The copula model performs
well, correctly classifying almost 90% of the states in the training dataset.
While the independence model also performs well on the training dataset
(similar to results in Candanedo et al. [2017]), it performs noticeably worse
for classifying the states in the test dataset.

Parameter | Lower | Estimate | Upper
---|---|---|---
$\pi_{1}$ | -0.336 | $2.243\times 10^{-8}$ | 0.336
$\pi_{2}$ | 0.664 | 1 | 1.336
$\gamma_{1,1}$ | 0.9882 | 0.9935 | 0.9988
$\gamma_{1,2}$ | 0.002489 | 0.01562 | 0.02876
$\gamma_{2,1}$ | 0.001206 | 0.006509 | 0.01181
$\gamma_{2,2}$ | 0.9712 | 0.9844 | 0.9975
$\mu_{1,1}$ | -0.4472 | -0.2685 | -0.08979
$\mu_{1,2}$ | $-8.785\times 10^{-7}$ | $-3.353\times 10^{-7}$ | $2.078\times 10^{-7}$
$\mu_{2,1}$ | -1.876 | 0.873 | 3.622
$\mu_{2,2}$ | $-6.644\times 10^{-7}$ | $2.562\times 10^{-6}$ | $5.788\times 10^{-6}$
$\sigma_{1,1}$ | 2.732 | 2.854 | 2.976
$\sigma_{1,2}$ | $9.482\times 10^{-6}$ | $9.894\times 10^{-6}$ | $1.031\times 10^{-5}$
$\sigma_{2,1}$ | 30.01 | 31.89 | 33.77
$\sigma_{2,2}$ | $3.418\times 10^{-5}$ | $3.613\times 10^{-5}$ | $3.808\times 10^{-5}$
$\theta_{1}$ | -0.2687 | 0.1209 | 0.5106
$\theta_{2}$ | 8.918 | 9.776 | 10.63
Table 2: Parameter estimates and 95% bootstrap confidence intervals for the copula-within-HMM model fit to the training data. Copula Model | Train | Test 1 | Test 2
---|---|---|---
Independence | 0.8950 | 0.8459 | 0.6795
Clayton/Frank | 0.8986 | 0.8515 | 0.6959
Table 3: Overall state classification accuracy for the training dataset and
the two test datasets using either the independence or Clayton/Frank copula.

##  6 Discussion

This paper proposes a copula-based approach for integration of information
from multivariate observations in a HMM setting. The dependence between the
components of the vector is captured using a copula, and we demonstrate that
ignoring the dependence can deteriorate the estimation capacity of the
statistical model. The computational challenges encountered when fitting the
copula HMM model are met using a new iterative algorithm, the EIFM algorithm,
which is designed broadly along the ideas of the EM algorithm but differs in
essential ways; most importantly, it makes computation possible when other
traditional methods fail.

The EIFM algorithm opens the door to developing other common HMM extensions
when considering copula-based state-dependent distributions, most of which
will impact the IFM step in the algorithm. One of the most common extensions
is to assume a time-inhomogeneous state process such that the dynamics of the
transition probability matrix $\bm{\Gamma}$ depend on time-varying covariates
of interest. A related extension is to allow for the state-dependent
parameters to vary over time. However, once the copula is introduced,
variations of the HMM will usually require specific customization, depending
on the application itself; for example, one might want to allow for the
marginal distributions or the copula parameters (or both) to vary over time,
or for situations in which components of the vector of observations are
observed at different times. One might also consider nonparametric estimation
of marginals.

Overall, our aim is to challenge the assumption of conditional contemporaneous
independence often made in HMM settings when working with multivariate data,
and when assuming longitudinal independence, to challenge the default
distributional assumption in the choice of the multivariate Gaussian. We
demonstrate the utility of both the Frank and Clayton copulas in an HMM
setting and demonstrate that with the aid of labelled data, we can even choose
different copula families across states. However, our work here is limited to
continuous-valued observation processes and to a bivariate setting. In cases
where the dimension of the observation process is greater than two, we will
investigate the use of vine copulas [Aas et al., 2009, Min and Czado, 2010].
Further, the support of the observation process may not always be the same
across dimensions. For example, having a multivariate process where one
dimension is discrete-valued, another continuous-valued and a third circular-
valued is not uncommon in ecological and environmental applications of HMMs
[Hodel and Fieberg, 2022]. Such situations, will require customizing the
copula techniques that were developed for mixed variables [Craiu and Sabeti,
2012] or discrete ones [Genest et al., 2014].

Lastly, in this work we limited our HMM presentation to a discrete-time
finite-state framework. However, the EIFM algorithm can be adapted to support
inference in other common HMM frameworks such as the continuous-time HMM and
hidden semi-Markov model (HSMM). In both of these scenarios, the time at which
the state evolves differs in distribution from that of the HMM presented here.
Continuous-time HMMs are commonly applied in medical applications where the
data collected over time at highly irregular time intervals. With the advent
of devices and tests that can collect multivariate observations, the copula-
within-HMM framework has the potential to better capture the correlation
between various aspects of someone’s health profile. The HSMM is a discrete-
time framework but with the added choice of selecting appropriate state-dwell
distributions, (i.e., the amount of time spent in a state before switching).
In a basic HMM, the state-dwell time distribution is necessarily geometric,
while in a HSMM it is common to select Poisson or negative binomial
distributions (although semi-parametric forms can also be specified [Pohle et
al., 2022]). In all cases, modifications to the EIFM algorithm will need to be
made to the forms presented in both the E-step and the IFM-step.

## Acknowledgments

RZ was supported by an Ontario Graduate Scholarship. RVC’s and VLB’s work have
been supported by the Natural Sciences and Engineering Research Council of
Canada. RZ thanks Michäel Lalancette and Yanbo Tang for helpful discussions.

## References

  * Hand [2020] David J Hand.  _Dark data: Why what you don’t know matters_.  Princeton University Press, 2020.
  * Rabiner and Juang [1986] Lawrence Rabiner and Biinghwang Juang.  An introduction to hidden markov models.  _ieee assp magazine_ , 3(1):4–16, 1986.
  * Zucchini et al. [2017] Walter Zucchini, Iain L MacDonald, and Roland Langrock.  _Hidden Markov models for time series: an introduction using R_.  CRC press, 2017.
  * Du et al. [2020] Kai Du, Steven Huddart, Lingzhou Xue, and Yifan Zhang.  Using a hidden markov model to measure earnings quality.  _Journal of Accounting and Economics_ , 69(2-3):101281, 2020.
  * Williams et al. [2020] Jonathan P Williams, Curtis B Storlie, Terry M Therneau, Clifford R Jack Jr, and Jan Hannig.  A bayesian approach to multistate hidden markov models: application to dementia progression.  _Journal of the American Statistical Association_ , 115(529):16–31, 2020.
  * Ailliot et al. [2009] Pierre Ailliot, Craig Thompson, and Peter Thomson.  Space–time modelling of precipitation by using a hidden markov model and censored gaussian distributions.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 58(3):405–426, 2009.
  * Ötting et al. [2021] Marius Ötting, Roland Langrock, and Antonello Maruotti.  A copula-based multivariate hidden markov model for modelling momentum in football.  _AStA Advances in Statistical Analysis_ , pages 1–19, 2021.
  * Cappé et al. [2005] Olivier Cappé, Eric Moulines, and Tobias Ryden.  _Inference in Hidden Markov Models (Springer Series in Statistics)_.  Springer-Verlag, Berlin, Heidelberg, 2005.  ISBN 0387402640.
  * Embrechts et al. [2001] Paul Embrechts, Filip Lindskog, and Alexander McNeil.  Modelling dependence with copulas.  _Rapport technique, Département de mathématiques, Institut Fédéral de Technologie de Zurich, Zurich_ , 14, 2001.
  * Sklar [1959] A. Sklar.  Fonctions de répartition à $n$ dimensions et leurs marges.  _Publications de l’Institut de Statistique de l’Université de Paris_ , 8:229–231, 1959.
  * Goetghebeur and Ryan [1995] E. Goetghebeur and L Ryan.  Analysis of competing risks survival data when some failure types are missing.  _Biometrika_ , 82:821–833, 1995.
  * Hougaard [2000] Philip Hougaard.  _Analysis of Multivariate Survival Data_.  Statistics for Biology and Health. Springer-Verlag, New York, 2000.  ISBN 0-387-98873-4.
  * Frees and Valdez [1998] E. W. Frees and E. A. Valdez.  Understanding relationships using copulas.  _North American Actuarial Journal_ , 2(1):1–25, 1998.
  * Genest et al. [2007] C. Genest, A.-C. Favre, J. Beliveau, and C. Jacques.  Metaelliptical copulas and their use in frequency analysis of multivariate hydrological data.  _Water Resources Research_ , 43(9), 9 2007.  doi: DOI10.1029/2006WR005275.
  * Dupuis [2007] D. J. Dupuis.  Using copulas in hydrology: Benefits, cautions, and issues.  _Journal of Hydrologic Engineering_ , 12(4):381–393, Jul-Aug 2007.  doi: DOI10.1061/(ASCE)1084-0699(2007)12:4(381).
  * Cherubini et al. [2004] Umberto Cherubini, Elisa Luciano, and Walter Vecchiato.  _Copula methods in finance_.  John Wiley & Sons, 2004.
  * da Silva Filho et al. [2012] Osvaldo Candido da Silva Filho, Flavio Augusto Ziegelmann, and Michael J Dueker.  Modeling dependence dynamics through copulas with regime switching.  _Insurance: Mathematics and Economics_ , 50(3):346–356, 2012.
  * Acar et al. [2011] Elif F Acar, Radu V Craiu, and Fang Yao.  Dependence calibration in conditional copulas: A nonparametric approach.  _Biometrics_ , 67(2):445–453, 2011.
  * Veraverbeke et al. [2011] N Veraverbeke, M. Omelka, and Irene Gijbels.  Estimation of a conditional copula and association measures.  _Scand. J. Statist._ , 38:766–780, 2011.
  * Chavez-Demoulin and Vatter [2015] V. Chavez-Demoulin and T. Vatter.  Generalized additive models for conditional copulas.  _J. Multivariate Anal._ , 141:147–167, 2015.
  * Craiu and Sabeti [2012] Radu V. Craiu and A. Sabeti.  In mixed company: Bayesian inference for bivariate conditional copula models with discrete and continuous outcomes.  _J. Multivariate Anal._ , 110:106–120, 2012.
  * Levi and Craiu [2018] Evgeny Levi and Radu V Craiu.  Bayesian inference for conditional copulas using gaussian process single index models.  _Computational Statistics & Data Analysis_, 122:115–134, 2018.
  * Aas et al. [2009] Kjersti Aas, Claudia Czado, Arnoldo Frigessi, and Henrik Bakken.  Pair-copula constructions of multiple dependence.  _Insurance Mathematics & Economics_, 44(2):182–198, 4 2009.  doi: DOI10.1016/j.insmatheco.2007.02.001.
  * Min and Czado [2010] Aleksey Min and Claudia Czado.  Bayesian inference for multivariate copulas using pair-copula constructions.  _Journal of Financial Econometrics_ , 8(4):511–546, 2010.
  * Joe and Xu [1996] Harry Joe and James Jianmeng Xu.  The estimation method of inference functions for margins for multivariate models.  Technical Report 166, Department of Statistics, University of British Columbia, 1996.
  * Elashoff and Ryan [2004] Michael Elashoff and Louise Ryan.  An em algorithm for estimating equations.  _Journal of Computational and Graphical Statistics_ , 13(1):48–65, 2004.
  * Demarta and McNeil [2005] Stefano Demarta and Alexander J McNeil.  The t copula and related copulas.  _International statistical review_ , 73(1):111–129, 2005.
  * Nelsen [2007] Roger B Nelsen.  _An introduction to copulas_.  Springer Science & Business Media, 2007.
  * Hofert et al. [2012] Marius Hofert, Martin Mächler, and Alexander J McNeil.  Likelihood inference for archimedean copulas in high dimensions under known margins.  _Journal of Multivariate Analysis_ , 110:133–150, 2012\.
  * Nasri et al. [2020] Bouchra R Nasri, Bruno N Remillard, and Mamadou Y Thioub.  Goodness-of-fit for regime-switching copula models with application to option pricing.  _Canadian Journal of Statistics_ , 48(1):79–96, 2020.
  * Genest [1987] Christian Genest.  Frank’s family of bivariate distributions.  _Biometrika_ , 74(3):549–555, 1987.
  * Baum et al. [1970] Leonard E Baum, Ted Petrie, George Soules, and Norman Weiss.  A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains.  _The annals of mathematical statistics_ , 41(1):164–171, 1970.
  * Dempster et al. [1977] Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * McLeish and Small [1988] Don L McLeish and Christopher G Small.  _The theory and applications of statistical interference functions_ , volume 44.  Springer Science & Business Media, 1988.
  * Joe [1997] Harry Joe.  _Multivariate models and multivariate dependence concepts_.  CRC Press, 1997.
  * Hofert et al. [2020] Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan.  _copula: Multivariate Dependence with Copulas_ , 2020.  URL https://CRAN.R-project.org/package=copula.  R package version 1.0-1.
  * Meng and Rubin [1993] Xiao-Li Meng and Donald B Rubin.  Maximum likelihood estimation via the ecm algorithm: A general framework.  _Biometrika_ , 80(2):267–278, 1993.
  * Ortega and Rheinboldt [2000] James M Ortega and Werner C Rheinboldt.  _Iterative solution of nonlinear equations in several variables_.  SIAM, 2000.
  * Godambe [1960] Vidyadhar P Godambe.  An optimum property of regular maximum likelihood estimation.  _The Annals of Mathematical Statistics_ , 31(4):1208–1211, 1960.
  * Jensen [2011] Jens Ledet Jensen.  Asymptotic normality of m-estimators in nonhomogeneous hidden markov models.  _Journal of Applied Probability_ , 48(A):295–306, 2011.
  * Efron and Tibshirani [1994] Bradley Efron and Robert J Tibshirani.  _An introduction to the bootstrap_.  CRC press, 1994.
  * Candanedo and Feldheim [2016] Luis M Candanedo and Véronique Feldheim.  Accurate occupancy detection of an office room from light, temperature, humidity and co2 measurements using statistical learning models.  _Energy and Buildings_ , 112:28–39, 2016.
  * Candanedo et al. [2017] Luis M Candanedo, Véronique Feldheim, and Dominique Deramaix.  A methodology based on hidden markov models for occupancy detection and a case study in a low energy residential building.  _Energy and Buildings_ , 148:327–341, 2017.
  * Visser and Speekenbrink [2010] Ingmar Visser and Maarten Speekenbrink.  depmixS4: An R package for hidden markov models.  _Journal of Statistical Software_ , 36(7):1–21, 2010.  URL https://www.jstatsoft.org/v36/i07/.
  * Kojadinovic et al. [2011] Ivan Kojadinovic, Jun Yan, and Mark Holmes.  Fast large-sample goodness-of-fit tests for copulas.  _Statistica Sinica_ , pages 841–871, 2011.
  * Hodel and Fieberg [2022] Florian H. Hodel and John R. Fieberg.  Circular–linear copulae for animal movement data.  _Methods in Ecology and Evolution_ , 13:1001–1013, 2022\.
  * Genest et al. [2014] Christian Genest, Johanna G Nešlehová, and Bruno Rémillard.  On the empirical multilinear copula process for count data.  _Bernoulli_ , 20(3):1344–1371, 2014.
  * Pohle et al. [2022] Jennifer Pohle, Timo Adam, and Larissa T. Beumer.  Flexible estimation of the state dwell-time distribution in hidden semi-markov models.  _Computational Statistics and Data Analysis_ , 172:107497, 2022.
  * Abramowitz et al. [1988] Milton Abramowitz, Irene A Stegun, and Robert H Romer.  Handbook of mathematical functions with formulas, graphs, and mathematical tables, 1988.

##  Appendix A Appendix

###  A.1 Regularity Conditions for Theorem 4.3

The following definition and assumptions are extracted from Jensen [2011].

######  Definition.

A set of functions
$\mathcal{A}:=\\{\mathbf{a}_{t}:\Theta\times\mathcal{Y}\times\mathcal{X}^{2}\to\mathbb{R}^{n}\\}_{t\geq
1}$ belongs to the class $C_{k}$ if there exist functions
$\\{a_{t}^{0}:\mathcal{Y}\to\mathbb{R}\\}_{t\geq 1}$, a constant
$\delta_{0}>0$, and a constant $K$ such that for all $t\geq 1$,

$\sup_{x_{1:2}\in\mathcal{X}^{2},\,\bm{\eta}\in
B(\bm{\eta}_{0},\delta_{0})}\left|\mathbf{a}_{t}(\bm{\eta},\mathbf{y},x_{1:2})\right|\leq
a^{0}_{t}(\mathbf{y})\qquad\text{and}\qquad\mathbb{E}_{\bm{\eta}_{0}}\left[a_{t}^{0}(\mathbf{Y}_{t})^{k}\right]\leq
K.$

Moreover, $\mathcal{A}$ belongs to the class $C_{k,m}$ if it belongs to
$C_{k}$, and in addition, there exist functions $\\{\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{a}_{t}:\mathcal{Y}\to\mathbb{R}\\}_{t\geq 1}$ a constant $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{0}>0$, and a constant $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{K}$
such that for all $\bm{\eta}\in B(\bm{\eta}_{0},\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{0})$,

$\left|\mathbf{a}_{t}(\bm{\eta},\mathbf{y},x_{1:2})-\mathbf{a}_{t}(\bm{\eta}_{0},\mathbf{y},x_{1:2})\right|\leq|\bm{\eta}-\bm{\eta}_{0}|\cdot\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{a}_{t}(\mathbf{y})\text{ for all
}x_{1:2}\in\mathcal{X}^{2}\qquad\text{and}\qquad\mathbb{E}_{\bm{\eta}_{0}}\left[\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{a}_{t}(\mathbf{Y}_{t})^{m}\right]\leq\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{K}.$

Below, we write $\omega_{t}(\bm{\eta};x_{(t-1):t},\mathbf{y}_{t})=\hskip
1.0pt\mathrm{log}\left(\gamma_{x_{t-1},x_{t}}\cdot
h_{x_{t}}(\mathbf{y}_{t})\right)$. Then the following assumptions are to be
made for Theorem 4.3 to hold:

  1. Assumption 1:

The entries of $\bm{\Gamma}_{0}$ are all strictly positive and
$0<\sum_{j=1}^{K}\nu_{j,t}\cdot h_{j}(\mathbf{y}_{t})<\infty$ for all
$\mathbf{y}_{t}\in\mathbb{R}^{d}$ and all $\bm{\eta}\in
B(\bm{\eta}_{0};\delta_{0})$.

  2. Assumption 2:

$\\{\bm{\phi}_{t}(\bm{\eta};x_{(t-1):t},\mathbf{y}_{t})\\}_{t\geq 1}$ is of
class $C_{4}$, and
$\\{\frac{\partial}{\partial\eta_{r}}\omega_{t}(\bm{\eta};x_{(t-1):t},\mathbf{y}_{t})\\}_{t\geq
1}$ and
$\\{\frac{\partial}{\partial\eta_{r}}\bm{\phi}_{t}(\bm{\eta};x_{(t-1):t},\mathbf{y}_{t})\\}_{t\geq
1}$ are of class $C_{4}$ and $C_{3,1}$, respectively, for each $r=1,\ldots,p$.

  3. Assumption 3:

There exist constants $K_{0}>0$ and $t_{0}$ such that for $t>t_{0}$,

$\mathbf{a}^{T}\mathbf{G}_{t}\mathbf{a}\geq
tK_{0}\left\lVert\mathbf{a}\right\rVert^{2}\qquad\text{for all
}\mathbf{a}\in\mathbb{R}^{p}.$

  4. Assumption 4:

There exist constants $c_{0}>0$ and $t_{0}$ such that for $t>t_{0}$, the
eigenvalues of $\mathbf{H}_{t}^{T}\mathbf{H}_{t}$ are bounded below by
$c_{0}$.

###  A.2 Proofs

In the following, we write
$X_{-t}=(X_{1},\ldots,X_{t-1},X_{t+1},\ldots,X_{T})$ and
$X_{-\\{t-1,t\\}}=(X_{1},\ldots,X_{t-2},X_{t+1},\ldots,X_{T})$.

######  Lemma A.1.

Under the model described by (1), we have

$\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid\mathbf{Y}_{1:T}=\mathbf{y}_{1:T}\right)=\frac{h_{k}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\sum_{j=1}^{K}h_{j}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}}.$
(24)

where

$\kappa(x_{-t},\mathbf{y}_{-t})=\pi_{x_{1}}\prod_{s\neq
t}h_{x_{s}}(\mathbf{y}_{s})\cdot\prod_{s\neq t,t+1}\gamma_{x_{s-1},x_{s}}.$

In particular, if the system constitutes a finite mixture model with
$\mathbb{P}_{\bm{\eta}}\left(X_{t}=j\right)=\nu_{j}$, we have

$\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid\mathbf{Y}_{1:T}=\mathbf{y}_{1:T}\right)=\frac{\nu_{k}\cdot
h_{k}(\mathbf{y}_{t})}{\sum_{j=1}^{K}\nu_{j}\cdot h_{j}(\mathbf{y}_{t})}.$
(25)

###### Proof.

Employing a standard abuse of notation in which densities are replaced by
probability measures where appropriate, Bayes’ rule gives us

$\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid\mathbf{Y}_{1:T}=\mathbf{y}_{1:T}\right)=\frac{\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:T}=\mathbf{y}_{1:T},X_{t}=k\right)}{\sum_{j=1}^{K}\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:T}=\mathbf{y}_{1:T},X_{t}=j\right)}.$
(26)

Now, for any $j\in\mathcal{X}$, the law of total probability and the
conditional independence of the observations $\mathbf{Y}_{1:T}$ respectively
give

$\displaystyle\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:T}=\mathbf{y}_{1:T},X_{t}=j\right)$
$\displaystyle=\sum_{x_{-t}\in\mathcal{X}^{T-1}}\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:T}=\mathbf{y}_{1:T},X_{t}=j,\mathbf{X}_{-t}=x_{-t}\right)$
$\displaystyle=\sum_{x_{-t}\in\mathcal{X}^{T-1}}\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:T}=\mathbf{y}_{1:T}\mid
X_{t}=j,X_{-t}=x_{-t}\right)\cdot\mathbb{P}_{\bm{\eta}}\left(X_{t}=j,X_{-t}=x_{-t}\right)$
$\displaystyle=\sum_{x_{-t}\in\mathcal{X}^{T-1}}h_{j}(\mathbf{y}_{t})\left(\prod_{s\neq
t}h_{x_{s}}(\mathbf{y}_{s})\right)\cdot\mathbb{P}_{\bm{\eta}}\left(X_{t}=j,X_{-t}=x_{-t}\right)$
$\displaystyle=\sum_{x_{-t}\in\mathcal{X}^{T-1}}h_{j}(\mathbf{y}_{t})\left(\prod_{s\neq
t}h_{x_{s}}(\mathbf{y}_{s})\right)\cdot\pi_{x_{1}}\prod_{s=2}^{t-1}\gamma_{x_{s-1},x_{s}}\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\cdot\prod_{s=t+2}^{T}\gamma_{x_{s-1},x_{s}}$
$\displaystyle=h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}.$
(27)

Inserting (27) into both numerator and denominator of (26) yields (24).

If the system constitutes a finite mixture model, then the Markov structure of
the underlying state sequence reduces to independence; that is,
$X_{t+1}\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X_{t}$ for all $t\in\\{1,\ldots,T\\}$, and in
particular, the rows of $\bm{\Gamma}$ are identical, with
$\gamma_{j,k}=\nu_{k}$ for all $j,k\in\mathcal{X}$. Then from (27) we have
that

$\displaystyle\mathbb{P}_{\bm{\eta}}\left(\mathbf{Y}_{1:T}=\mathbf{y}_{1:T},X_{t}=j\right)$
$\displaystyle=h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}$
$\displaystyle=h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\nu_{j}\cdot\nu_{x_{t+1}}$
$\displaystyle=C\cdot\nu_{j}\cdot h_{j}(\mathbf{y}_{t}),$ (28)

where $C=\sum_{x_{-t}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\nu_{x_{t+1}}$ is
constant with respect to both $j$ and $\mathbf{y}_{t}$. Inserting (28) into
both numerator and denominator of (26) yields (25) upon factoring out $C$. ∎

######  Theorem A.1.

Let $\nu_{t,k}=\mathbb{P}\left(X_{t}=k\right)=[\bm{\pi}\bm{\Gamma}^{t}]_{k}$.
The expected zero-one loss of the classifications made by local decoding is
given by

$\ell_{01}(\bm{\eta})=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,k}\cdot\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1\right|X_{t}=k\right)$
(29)

where $\kappa(x_{-t},\mathbf{y}_{-t})=\pi_{x_{1}}\prod_{s\neq
t}h_{x_{s}}(\mathbf{y}_{s})\cdot\prod_{s\neq t,t+1}\gamma_{x_{s-1},x_{s}}$.
Moreover, if the HMM constitutes a finite mixture model, then (2) reduces to

$\ell_{01}(\bm{\eta})=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,k}\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{\nu_{t,k}\cdot
h_{k}(\mathbf{y}_{t})}{\max_{j\neq k}\nu_{t,j}\cdot
h_{j}(\mathbf{y}_{t})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}_{t}).$
(30)

###### Proof.

Write $\hat{X}_{t}$ for the prediction made at time $t$. The local decoding
algorithm selects for $\hat{X}_{t}$ the state $k$ which maximizes the
conditional probability
$\mathbb{P}\left(X_{t}=k\mid\mathbf{Y}_{1:T}=\mathbf{y}_{1:T}\right)$, which
by A.1 is equivalent to

$\displaystyle\hat{X}_{t}$
$\displaystyle=\underset{k\in\mathcal{X}}{\mathrm{argmax}}\,\left(\frac{h_{k}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\sum_{j=1}^{K}h_{j}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}}\right)$
$\displaystyle=\underset{k\in\mathcal{X}}{\mathrm{argmax}}\,\left(h_{k}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}\right).$

Thus, foregoing the possibility of ties (which almost surely do not occur), we
have that

$\displaystyle\phantom{{}\iff{}}\hat{X}_{t}\neq k$ $\displaystyle\iff
h_{k}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}\neq\underset{j\neq
k}{\mathrm{argmax}}\,\left\\{h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}$
$\displaystyle\iff
h_{k}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}<\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}$
$\displaystyle\iff\frac{h_{k}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1.$
(31)

Taking expectations with respect to $(\mathbf{Y}_{1:T},X_{1:T})$, the expected
zero-one loss satisfies

$\displaystyle\ell_{0,1}(\bm{\eta})$
$\displaystyle=\mathbb{E}\left[\frac{1}{T}\sum_{t=1}^{T}\mathbbm{1}\left({\hat{X}_{t}\neq
X_{t}}\right)\right]$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\mathbb{P}_{\bm{\eta}}\left(\hat{X}_{t}\neq
X_{t}\right)$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,k}\cdot\mathbb{P}_{\bm{\eta}}\left(\hat{X}_{t}\neq
k\mid X_{t}=k\right)$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,j}\cdot\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1\right|X_{t}=k\right).$

If the system constitutes a finite mixture model, then (31) reduces to

$\frac{\nu_{t,k}\cdot h_{k}(\mathbf{y}_{t})}{\max_{j\neq k}\nu_{t,j}\cdot
h_{k}(\mathbf{y}_{t})}<1$

by (30), and the same derivation leads to

$\ell_{01}(\bm{\eta})=\frac{1}{T}\sum_{t=1}^{T}\sum_{k=1}^{K}\nu_{t,k}\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{\nu_{t,k}\cdot
h_{k}(\mathbf{y}_{t})}{\max_{j\neq k}\nu_{t,j}\cdot
h_{j}(\mathbf{y}_{t})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}_{t})$

as desired. ∎

######  Proposition A.1.

Suppose that $\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid
X_{t-1},X_{t+1}\right)\geq\mathbb{P}_{\bm{\eta}}\left(X_{t}=j\mid
X_{t-1},X_{t+1}\right)$ for all $j\in\mathcal{X}$. Then

$\displaystyle\phantom{{}={}}\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{X_{t-1},k}\cdot\gamma_{k,X_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{X_{t-1},j}\cdot\gamma_{j,X_{t+1}}\right\\}}<1\right|X_{t}=k\right)$
$\displaystyle\leq\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{h_{k}(\mathbf{y}_{t})}{\max_{j\neq
k}h_{j}(\mathbf{y}_{t})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}_{t}).$

###### Proof.

To begin with, note that for any $x_{t-1},x_{t+1}\in\mathcal{X}$,

$\displaystyle\phantom{{}={}}\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid
X_{t-1}=x_{t-1},X_{t+1}=x_{t+1}\right)$
$\displaystyle=\frac{\mathbb{P}_{\bm{\eta}}\left(X_{t-1}=x_{t-1},X_{t}=k,X_{t+1}=x_{t+1}\right)}{\mathbb{P}_{\bm{\eta}}\left(X_{t-1}=x_{t-1},X_{t+1}=x_{t+1}\right)}$
(32)
$\displaystyle=\frac{\mathbb{P}_{\bm{\eta}}\left(X_{t}=k,X_{t+1}=x_{t+1}\mid
X_{t-1}=x_{t-1}\right)\cdot\mathbb{P}_{\bm{\eta}}\left(X_{t-1}=x_{t-1}\right)}{\mathbb{P}_{\bm{\eta}}\left(X_{t-1}=x_{t-1},X_{t+1}=x_{t+1}\right)}$
$\displaystyle=\frac{\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid
X_{t-1}=x_{t-1}\right)\cdot\mathbb{P}_{\bm{\eta}}\left(X_{t+1}=x_{t+1}\mid
X_{t}=k\right)\cdot\mathbb{P}_{\bm{\eta}}\left(X_{t-1}=x_{t-1}\right)}{\mathbb{P}_{\bm{\eta}}\left(X_{t-1}=x_{t-1},X_{t+1}=x_{t+1}\right)}$
$\displaystyle=\frac{\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\mathbb{P}_{\bm{\eta}}\left(X_{t+1}=x_{t+1}\mid
X_{t-1}=x_{t-1}\right)},$

and therefore the Markov chain condition
$\mathbb{P}_{\bm{\eta}}\left(X_{t}=k\mid
X_{t-1},X_{t+1}\right)\geq\mathbb{P}_{\bm{\eta}}\left(X_{t}=j\mid
X_{t-1},X_{t+1}\right)$ for all $j\in\mathcal{X}$ is equivalent to

$\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}\geq\max_{j\neq
k}\left\\{\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}$ (33)

for all $x_{t-1},x_{t+1}\in\mathcal{X}$. Thus, for any
$\mathbf{y}_{1:T}\in\mathbb{R}^{d\times T}$, we have that

$\displaystyle\phantom{{}\iff{}}h_{k}(\mathbf{y}_{t})\geq\max_{j\neq
k}h_{j}(\mathbf{y}_{t})$ $\displaystyle\iff
h_{k}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\geq\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\right\\}$
$\displaystyle\implies
h_{k}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}\geq\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\right\\}\cdot\max_{j\neq
k}\left\\{\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}$
$\displaystyle\implies
h_{k}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}\geq\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\},$

where we have used (33) and the standard fact that
$\sup_{x}f(x)\cdot\sup_{x}g(x)\geq\sup_{x}\left(f(x)\cdot g(x)\right)$ when
$f,g>0$. It follows that as subsets of $\mathbb{R}^{d\times T}$,

$\left\\{\mathbf{y}_{1:T}:\frac{h_{k}(\mathbf{y}_{t})}{\max_{j\neq
k}h_{j}(\mathbf{y}_{t})}\geq
1\right\\}\subseteq\left\\{\mathbf{y}_{1:T}:\frac{h_{k}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}\geq
1\right\\},$

or equivalently,

$\left\\{\mathbf{y}_{1:T}:\frac{h_{k}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{y}_{-t})\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1\right\\}\subseteq\left\\{\mathbf{y}_{1:T}:\frac{h_{k}(\mathbf{y}_{t})}{\max_{j\neq
k}h_{j}(\mathbf{y}_{t})}<1\right\\}.$

Thus,

$\displaystyle\phantom{{}={}}\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},k}\cdot\gamma_{k,x_{t+1}}}{\max_{j\neq
k}\left\\{h_{j}(\mathbf{Y}_{t})\cdot\sum_{x_{-t}\in\mathcal{X}^{T-1}}\kappa(x_{-t},\mathbf{Y}_{-t})\cdot\gamma_{x_{t-1},j}\cdot\gamma_{j,x_{t+1}}\right\\}}<1\right|X_{t}=k\right)$
$\displaystyle\leq\mathbb{P}_{\bm{\eta}}\left(\left.\frac{h_{k}(\mathbf{Y}_{t})}{\max_{j\neq
k}h_{j}(\mathbf{Y}_{t})}<1\right|X_{t}=k\right)$
$\displaystyle=\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{h_{k}(\mathbf{y}_{t})}{\max_{j\neq
k}h_{j}(\mathbf{y}_{t})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}_{t}).$

∎

######  Proposition A.2.

Let $d=2$ and fix $k\in\mathcal{X}$. As either $|\tau_{k}|\to 1$ or
$|\rho_{k}|\to 1$ while the other state-dependent parameters stay fixed, we
have

$\int_{\mathbb{R}^{d}}\mathbbm{1}\left\\{\frac{\omega_{k}\cdot
h_{k}(\mathbf{y})}{\max_{j\neq k}\omega_{j}\cdot
h_{j}(\mathbf{y})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})\to 0$

for any $(\omega_{1},\ldots,\omega_{K})\in\mathcal{S}^{K-1}$ with
$\omega_{k}>0$.

###### Proof.

For ease of notation, we omit the subscript $t$ in what follows. First, note
that by Theorem 3.6 of Embrechts et al. [2001], $\rho_{k}=1$ if and only if
$\tau_{k}=1$, either of which are equivalent to $C_{k}$ being the
comonotonicity copula. Similarly, $\rho_{k}=-1$ if and only if $\tau_{k}=-1$,
corresponding to the countermonotonicity copula.

Let $\epsilon>0$, and start by assuming $\rho_{k}\to 1$. Then

$\displaystyle\int_{\mathbb{R}^{d}}\mathbbm{1}\left({\frac{\omega_{k}\cdot
h_{k}(\mathbf{y})}{\max_{j\neq k}\omega_{j}\cdot
h_{j}(\mathbf{y})}<1}\right)\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})$
$\displaystyle=\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\max_{j\neq
k}\omega_{j}\cdot h_{j}(\mathbf{y})>\omega_{k}\cdot
h_{k}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})$
$\displaystyle=\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\bigvee_{j\neq
k}\omega_{j}\cdot h_{j}(\mathbf{y})>\omega_{k}\cdot
h_{k}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})$
$\displaystyle\leq\sum_{j\neq
k}\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\omega_{j}\cdot
h_{j}(\mathbf{y})>\omega_{k}\cdot
h_{k}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y}),$ (34)

where we have used Boole’s inequality in the last line. Consider now the
integral

$\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\omega_{j}\cdot
h_{j}(\mathbf{y})>\omega_{k}\cdot
h_{k}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})$ (35)

inside the summation of (34). If $\omega_{j}=0$ then (35) is zero and
contributes nothing to the sum, so we may assume $\omega_{j}>0$ and let
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}=\omega_{j}/\omega_{k}$. Further, let
$d_{j}(\mathbf{u})=c_{j}(F_{j,1}(F_{k,1}^{-1}(u_{1})),F_{j,2}(F_{k,2}^{-1}(u_{2}))\cdot\prod_{h=1}^{2}f_{j,h}(F_{k,h}^{-1}(u_{h}))/f_{k,h}(F_{k,h}^{-1}(u_{h}))$
and $\mathbf{u}=\left(F_{k,1}(y_{1}),F_{k,2}(y_{2})\right)$ so that by a
change of variables, (35) becomes

$\displaystyle\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}h_{j}(\mathbf{y})>h_{k}(\mathbf{y})}\right)h_{k}(\mathbf{y})\mathop{}\\!\mathrm{d}\mathbf{y}$
$\displaystyle=\int_{[0,1]^{2}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}$
$\displaystyle=\int_{B_{\delta}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}+\int_{B^{c}_{\delta}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u},$

where $B_{\delta}=\\{\mathbf{u}\in[0,1]^{2}:|u_{1}-u_{2}|\leq\delta\\}$, for
some $\delta>0$ to be specified later. Now for $(U_{1},U_{2})\sim C_{k}$, we
have, by Markov’s inequality,

$\int_{B^{c}_{\delta}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}\leq\int_{B^{c}_{\delta}}c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}=\mathbb{P}\left(|U_{1}-U_{2}|>\delta\right)\leq\frac{\mathbb{E}\left[(U_{1}-U_{2})^{2}\right]}{\delta^{2}}=\frac{1-\text{Corr}(U_{1},U_{2})}{6\delta^{2}},$
(36)

On the other hand,

$\displaystyle\int_{B_{\delta}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}$
$\displaystyle<\int_{B_{\delta}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}$
$\displaystyle\leq\int_{B_{\delta}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}$
$\displaystyle=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}\int_{B_{\delta}}c_{j}(F_{j,1}(F_{k,1}^{-1}(u_{1})),F_{j,2}(F_{k,2}^{-1}(u_{2}))\cdot\prod_{h=1}^{2}f_{j,h}(F_{k,h}^{-1}(u_{h}))/f_{k,h}(F_{k,h}^{-1}(u_{h}))\mathop{}\\!\mathrm{d}\mathbf{u}.$
(37)

Let
$B^{\prime}_{\delta}=\\{\mathbf{w}\in\mathbb{R}^{2}:|F_{k,1}(w_{1})-F_{k,2}(w_{2})|\leq\delta\\}$.
Making the substitution
$\mathbf{w}=(F_{k,1}^{-1}(u_{1}),F_{k,2}^{-1}(u_{2}))$, (37) becomes

$\displaystyle\phantom{{}={}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}\int_{B^{\prime}_{\delta}}c_{j}(F_{j,1}(w_{1}),F_{j,2}(w_{2}))\cdot\prod_{h=1}^{2}\left(f_{j,h}(w_{h})/f_{k,h}(w_{h})\right)\cdot\prod_{h=1}^{2}f_{k,h}(w_{h})\mathop{}\\!\mathrm{d}\mathbf{w}$
$\displaystyle=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}\int_{B^{\prime}_{\delta}}c_{j}(F_{j,1}(w_{1}),F_{j,2}(w_{2}))\cdot\prod_{h=1}^{2}f_{j,h}(w_{h})\mathop{}\\!\mathrm{d}\mathbf{w}$
$\displaystyle=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}\int_{B^{\prime}_{\delta}}h_{j}(\mathbf{w})\mathop{}\\!\mathrm{d}\mathbf{w}$
$\displaystyle=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}\cdot\mu_{j}(B^{\prime}_{\delta}),$

where $\mu_{j}$ is the measure on $\mathbb{R}^{2}$ induced by $H_{j}$. Since
$\mu_{j}$ is absolutely continuous with respect to Lebesgue measure
$\Lambda(\cdot)$ on $\mathbb{R}^{2}$, there exists some $\delta_{\epsilon}$
such that whenever $\Lambda(B^{\prime}_{\delta})<\delta_{\epsilon}$, we have
$\mu_{j}(B^{\prime}_{\delta})<\epsilon/((K-1)\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j})$. Since we have $\Lambda(B^{\prime}_{\delta})\to 0$ as $\delta\to
0$, there must exist some $\delta^{*}$ small enough that
$\Lambda(B^{\prime}_{\delta^{*}})<\delta_{\epsilon}$. Making the choice
$\delta=\delta^{*}$ thus gives

$\int_{B_{\delta}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}d_{j}(\mathbf{u})>c_{k}(\mathbf{u})}\right)c_{k}(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}<\frac{\epsilon}{K-1}.$
(38)

Now, $\rho_{k}=\text{Corr}(U_{1},U_{2})$ because $H_{k}$ has continuous
margins [Nelsen, 2007], and so upon combining (36) and (38), we obtain

$\limsup_{\rho_{k}\to
1}\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{}_{j}h_{j}(\mathbf{y})>h_{k}(\mathbf{y})}\right)h_{k}(\mathbf{y})\mathop{}\\!\mathrm{d}\mathbf{y}<\limsup_{\rho_{k}\to
1}\left(\frac{1-\rho_{k}}{6\delta^{2}}+\frac{\epsilon}{K-1}\right)\leq\frac{\epsilon}{K-1}.$
(39)

Inserting (39) into (34) gives

$\limsup_{\rho_{k}\to
1}\int_{\mathbb{R}^{2}}\mathbbm{1}\left\\{\frac{\omega_{k}\cdot
h_{k}(\mathbf{y})}{\max_{j\neq k}\omega_{j}\cdot
h_{j}(\mathbf{y})}<1\right\\}\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})\leq\epsilon.$
(40)

Finally, if instead we have $\rho_{k}\to-1$, we require the numerator of (36)
to be $1+\text{Corr}(U_{1},U_{2})$, and this only needs the modification
$B_{\delta}=\\{\mathbf{u}\in[0,1]^{2}:|u_{1}+u_{2}-1|\leq\delta\\}$ and
similarly for $B^{\prime}_{\delta}$, aside from which the proof is identical.
In either case, (40) can be made arbitrarily close to 0, which completes the
proof. ∎

######  Corollary A.1.

Let $d=2$ and $K=2$. Suppose
$\mathcal{C}=\\{C_{\theta}(\cdot,\cdot):\theta\in\Theta\\}$ is a parametric
family of radially symmetric bivariate copulas satisfying
$C_{-\theta}(u,v)=u-C_{\theta}(u,1-v)=v-C_{\theta}(1-u,v)$. If our HMM is a
finite mixture model with $C_{1}=C_{\theta}$ and $C_{2}=C_{-\theta}$ with
identical marginals and
$\mathbb{P}\left(X_{t}=1\right)=\mathbb{P}\left(X_{t}=2\right)=1/2$ for all
$t\geq 1$, then

$\ell_{01}(\theta)=1-2\cdot C_{\theta}\left(\frac{1}{2},\frac{1}{2}\right).$
(41)

In particular, for $\theta>0$ we have
$\ell_{01}^{\mathrm{Frank}}(\theta)=1-2\hskip
1.0pt\mathrm{log}\left(\cosh{\frac{\theta}{4}}\right)/\theta$,
$\ell_{01}^{\mathrm{Gauss}}(\rho)=\cos^{-1}(\rho)/\pi$, and
$\ell_{01}^{\mathrm{FGM}}(\theta)=1/2-\theta/8$.

######  Remark.

In all three cases, $\ell_{01}(0)=1/2$, which is in line with intuition. While
$\ell_{01}^{\mathrm{Frank}}(\pm\infty)=\ell_{01}^{\mathrm{Gauss}}(\pm 1)=0$,
we only have $\ell_{01}^{\mathrm{FGM}}(1)=3/8$. That
$\ell_{01}^{\mathrm{FGM}}(1)\neq 0$ arises from the fact that the FGM copula
can model only weak dependence [Nelsen, 2007], and thus is not comprehensive.
The Plackett family of copulas is radially symmetric, but with parameter space
$(0,\infty)$, it does not satisfy the relationship in the statement of the
corollary. It does, however, satisfy $C_{1/\theta}(u,v)=u-C_{\theta}(u,v)$
[Nelsen, 2007], for which a modification of (41) can be developed.

###### Proof.

Since $\nu_{t,k}=1/2$ for all $t\geq 1$ and $k\in\mathcal{X}$, using Theorem
A.1, the zero-one loss in this situation reduces to

$\displaystyle\ell_{01}(\bm{\eta})$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{2}\int_{\mathbb{R}^{2}}\mathbbm{1}\left({\frac{h_{k}(\mathbf{y})}{\max_{j\neq
k}h_{j}(\mathbf{y})}<1}\right)\mathop{}\\!\mathrm{d}H_{k}(\mathbf{y})$
$\displaystyle=\frac{1}{2}\left(\int_{\mathbb{R}^{2}}\mathbbm{1}\left({h_{1}(\mathbf{y})<h_{2}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{1}(\mathbf{y})+\int_{\mathbb{R}^{2}}\mathbbm{1}\left({h_{2}(\mathbf{y})<h_{1}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{2}(\mathbf{y})\right)$
$\displaystyle=1-\frac{1}{2}\left(\int_{\mathbb{R}^{2}}\mathbbm{1}\left({h_{1}(\mathbf{y})>h_{2}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{1}(\mathbf{y})+\int_{\mathbb{R}^{2}}\mathbbm{1}\left({h_{2}(\mathbf{y})>h_{1}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{2}(\mathbf{y})\right)$
$\displaystyle=1-\frac{1}{2}\left(\int_{\mathbb{R}^{2}}\mathbbm{1}\left({h_{\theta}(\mathbf{y})>h_{-\theta}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{\theta}(\mathbf{y})+\int_{\mathbb{R}^{2}}\mathbbm{1}\left({h_{-\theta}(\mathbf{y})>h_{\theta}(\mathbf{y})}\right)\mathop{}\\!\mathrm{d}H_{-\theta}(\mathbf{y})\right),$
(42)

where we have made the change of notation $h_{\theta}:=h_{1}$ and
$h_{-\theta}:=h_{2}$ for convenience. Now,
$h_{\theta}(\mathbf{y})=c_{\theta}\left(F_{1}(y_{1}),F_{2}(y_{2})\right)f_{1}(y_{1})f_{2}(y_{2})$,
so by cancelling out marginal densities and writing
$u_{h}=F_{h}(y_{h})\in[0,1]$ for $h=1,2$, we have

$h_{\theta}(\mathbf{y})>h_{-\theta}(\mathbf{y})\iff
c_{\theta}(u_{1},u_{2})>c_{-\theta}(u_{1},u_{2}),$

and therefore after applying the change of variables
$\mathbf{u}=(F_{1}(y_{1}),F_{2}(y_{2}))$, (42) becomes

$1-\frac{1}{2}\left(\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{\theta}(u_{1},u_{2})>c_{-\theta}(u_{1},u_{2})}\right)c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}\mathbf{u}+\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{-\theta}(u_{1},u_{2})>c_{\theta}(u_{1},u_{2})}\right)c_{-\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}\mathbf{u}\right).$
(43)

Now, the two integrals in (43) are _equal_ ; to show this, we make the
temporary change of variables $w=1-u_{2}$ and repeatedly use the implication
$c_{-\theta}(u_{1},u_{2})=c_{\theta}(u_{1},1-u_{2})=c_{\theta}(1-u_{1},u_{2})$
to obtain

$\displaystyle\phantom{{}={}}\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{-\theta}(u_{1},u_{2})>c_{\theta}(u_{1},u_{2})}\right)c_{-\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}\mathbf{u}$
$\displaystyle=-\int_{0}^{1}\int_{1}^{0}\mathbbm{1}\left({c_{-\theta}(u_{1},1-w)>c_{\theta}(u_{1},1-w)}\right)c_{-\theta}(u_{1},1-w)\mathop{}\\!\mathrm{d}w\mathop{}\\!\mathrm{d}u_{1}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{\theta}(u_{1},w)>c_{-\theta}(u_{1},w)}\right)c_{\theta}(u_{1},w)\mathop{}\\!\mathrm{d}w\mathop{}\\!\mathrm{d}u_{1}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{\theta}(u_{1},u_{2})>c_{-\theta}(u_{1},u_{2})}\right)c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}\mathbf{u}$

Moreover, the common value of the two integrals is equal to

$\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{\theta}(u_{1},u_{2})>c_{\theta}(u_{1},1-u_{2})}\right)c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}\mathbf{u},$
(44)

so

$\ell_{01}(\bm{\eta})=1-\int_{0}^{1}\int_{0}^{1}\mathbbm{1}\left({c_{\theta}(u_{1},u_{2})>c_{\theta}(u_{1},1-u_{2})}\right)c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}\mathbf{u}.$

It remains to evaluate (44). Let
$A=\\{\mathbf{u}:c_{\theta}(u_{1},u_{2})>c_{\theta}(u_{1},1-u_{2})\\}$, and
partition it as $A=A_{1}\cup A_{2}\cup A_{3}\cup A_{4}$ with

$\displaystyle A_{1}$ $\displaystyle=A\cap(0,1/2)^{2}$ $\displaystyle A_{2}$
$\displaystyle=A\cap\left((0,1/2)\times(1/2,1)\right)$ $\displaystyle A_{3}$
$\displaystyle=A\cap\left((1/2,1)\times(0,1/2)\right)$ $\displaystyle A_{4}$
$\displaystyle=A\cap(1/2,1)^{2}.$

Then

$\displaystyle A_{2}$
$\displaystyle=\\{\mathbf{u}:c_{\theta}(u_{1},u_{2})>c_{\theta}(u_{1},1-u_{2}),u_{1}\in(0,1/2),u_{2}\in(1/2,1)\\}$
$\displaystyle=\\{\mathbf{u}:c_{\theta}(1-u_{1},1-u_{2})>c_{\theta}(u_{1},1-u_{2}),u_{1}\in(0,1/2),u_{2}\in(1/2,1)\\}$
$\displaystyle=\\{\mathbf{u}:c_{\theta}(1-u_{1},u_{2})>c_{\theta}(u_{1},u_{2}),u_{1}\in(0,1/2),u_{2}\in(0,1/2)\\}$
$\displaystyle=\\{\mathbf{u}:c_{\theta}(u_{1},u_{2})\leq
c_{\theta}(u_{1},1-u_{2}),u_{1}\in(0,1/2),u_{2}\in(0,1/2)\\}$
$\displaystyle=A^{c}\cap(0,1/2)^{2}$

and so $A_{1}\cup A_{2}=(0,1/2)^{2}$. With similar manipulations, we can
obtain $A_{3}\cup A_{4}=(1/2,1)^{2}$ and hence $A=(0,1/2)^{2}\cup(1/2,1)^{2}$.
Thus, (44) reduces to

$\displaystyle\phantom{{}={}}\int_{0}^{1/2}\int_{0}^{1/2}c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}u_{2}\mathop{}\\!\mathrm{d}u_{1}+\int_{1/2}^{1}\int_{1/2}^{1}c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}u_{2}\mathop{}\\!\mathrm{d}u_{1}$
$\displaystyle=\int_{0}^{1/2}\int_{0}^{1/2}c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}u_{2}\mathop{}\\!\mathrm{d}u_{1}+\int_{1/2}^{1}\int_{1/2}^{1}c_{\theta}(1-u_{1},1-u_{2})\mathop{}\\!\mathrm{d}u_{2}\mathop{}\\!\mathrm{d}u_{1}$
$\displaystyle=\int_{0}^{1/2}\int_{0}^{1/2}c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}u_{2}\mathop{}\\!\mathrm{d}u_{1}+\int_{0}^{1/2}\int_{0}^{1/2}c_{\theta}(u_{1},u_{2})\mathop{}\\!\mathrm{d}u_{2}\mathop{}\\!\mathrm{d}u_{1}$
$\displaystyle=2\cdot C_{\theta}\left(\frac{1}{2},\frac{1}{2}\right),$

which proves the general result. The specific cases for the Frank and FGM
families follow from direct calculation, and that for the Gauss copula follows
from a well-known result about quadrant probabilities for the bivariate normal
distribution (Formula 23.6.19 of Abramowitz et al. [1988]). ∎

######  Remark.

This expression is very close to _Blomqvist’s $\beta$_. See Section 5.1.4 of
[Nelsen, 2007].