# Modelling multi-scale state-switching functional data with hidden Markov
models

Evan Sidrow
Department of Statistics
University of British Columbia
Vancouver, Canada
<EMAIL_ADDRESS>
&Nancy Heckman
Department of Statistics
University of British Columbia
Vancouver, Canada
&Sarah M.E. Fortune
Marine Mammal Research Unit
University of British Columbia
Vancouver, Canada
&Andrew W. Trites
Institute for the Oceans and Fisheries
University of British Columbia
Vancouver, Canada
&Ian Murphy
Department of Statistics
University of British Columbia
Vancouver, Canada
&Marie Auger-Méthé
Department of Statistics
University of British Columbia
Vancouver, Canada

###### Abstract

Data sets comprised of sequences of curves sampled at high frequencies in time
are increasingly common in practice, but they can exhibit complicated
dependence structures that cannot be modelled using common methods of
Functional Data Analysis (FDA). We detail a hierarchical approach which treats
the curves as observations from a hidden Markov model (HMM). The distribution
of each curve is then defined by another fine-scale model which may involve
auto-regression and require data transformations using moving-window summary
statistics or Fourier analysis. This approach is broadly applicable to
sequences of curves exhibiting intricate dependence structures. As a case
study, we use this framework to model the fine-scale kinematic movement of a
northern resident killer whale (Orcinus orca) off the coast of British
Columbia, Canada. Through simulations, we show that our model produces more
interpretable state estimation and more accurate parameter estimates compared
to existing methods.

_K_ eywords Accelerometer data $\cdot$ animal movement $\cdot$ biologging
$\cdot$ diving behaviour $\cdot$ hierarchical modelling $\cdot$ killer whales
$\cdot$ state-switching $\cdot$ statistical ecology $\cdot$ time series

##  1 INTRODUCTION

Biologging technology now provides researchers with kinematic data collected
almost continuously in time (Hooten et al.,, 2017). The collection and
analysis of data from devices such as accelerometers have brought new insights
to areas ranging from monitoring machine health (Getman et al.,, 2009) to
understanding physical activity levels in children (Morris et al.,, 2006). The
study of animal movement in particular has been transformed by tracking
devices that record kinematic information in a variety of environments (Börger
et al.,, 2020; Jeanniard-du Dot et al.,, 2016). Tags can record over 50
observations per second, resulting in time series containing millions of
observations over the course of several hours. These data contain a wealth of
information about human and animal behaviour, but parsing these large data
sets poses a challenge for statisticians and biologists.

Biologging data are frequently modelled as a set of curves analyzed by methods
of Functional Data Analysis, or FDA – see, e.g., Ramsay and Silverman, (2005).
For example, Morris et al., (2006) view a child’s activity level as a set of
daily curves where metabolic activity is a function of time. Similarly, Fu and
Heckman, (2019) view the dive profile of a southern elephant seal (Mirounga
leonina) as a set of dive curves whereby the amplitude and phase variation of
the dive is used to classify dive types.

FDA was originally developed to process curves assumed to be independent
replicates (i.e., there is no between-curve dependence), and within-curve
fine-scale structure is not usually incorporated in FDA models. However, sets
of curves often exhibit complex sequential dependencies both between curves
and within curves, especially in the case of biologging data (Leos-Barajas et
al.,, 2017). For example, dive profiles of marine animals show distinct dive
“types” that cluster together in time (Tennessen et al., 2019a, ) while
simultaneously displaying bouts of short-term periodicity within dives (Adam
et al.,, 2019). Fine-scale periodicity within a larger process is also common
in fields ranging from machine health (Xin et al.,, 2018; Lucero et al.,,
2019) to speech recognition (Juang and Rabiner,, 1991).

Some FDA models account for between-curve dependence occurring when multiple
curves arise from separate groups of individuals, but these models are
inadequate when modelling certain types of dependence in time. For example,
previous studies have used multilevel models with random effects to model
variation between and within individuals in daily activity levels of children
(Morris et al.,, 2006) or in menstrual cycles of adults (Brumback and Rice,,
1998). More recent work involving multilevel models includes that of Di et
al., (2009), Crainiceanu et al., (2009), and Chen and Müller, (2012). However,
these multilevel models are not appropriate in many biologging applications
since they do not account for different curve types. Further, the first two
papers do not account for temporal between-curve dependence and the third
paper’s model of temporal dependence is not appropriate in most biologging
applications. In addition to multi-level models, FDA researchers have also
used functional time series to model dependence in a sequence of curves.
Functional time series extend the ideas of classic time series to model the
evolution of one curve into the next (Kokoszka and Reimherr,, 2018), but they
do not account for sequences of time-series curves whose distributions are
determined by well-defined hidden states.

Traditional FDA techniques similarly fall short when modelling complicated
within-curve data. In particular, within-curve structure is usually modelled
by a generic smooth mean function and a covariance function (Yao et al.,,
2005) or with random regression (Rice and Wu,, 2001). However, time-series
data exhibiting both sharp behavioural changes and periodic fine-scale
structure are difficult to model with these classic FDA techniques.

To accommodate both temporal dependence and changes in curve type, we turn to
the field of animal movement modelling (Hooten et al.,, 2017), where one of
the most prevalent techniques of late is the hidden Markov model, or HMM
(Patterson et al.,, 2017; McClintock et al.,, 2020). HMMs interpret animal
movement data as arising from a sequence of unobserved behavioural states,
allowing biologists to infer the underlying behaviour of an animal from
sequential observations of its position. While ubiquitous in ecology
literature, HMMs have seen little use in non-parametric functional modelling
with a few notable exceptions. In particular, Langrock et al., (2018) take a
non-parametric approach to model the distributions of HMM observations with
B-splines, but this approach does not account for certain types of temporal
correlation.

While useful, HMMs alone are also not sufficient to model intricate fine-scale
time-series data for three primary reasons. Firstly, HMMs assume that
subsequent observations are independent given an underlying hidden state
process, but this is often not the case when observations are taken at
extremely high frequencies. Several solutions have been proposed in the
ecology literature such as the hidden movement Markov model (Whoriskey et
al.,, 2016) and the conditionally auto-regressive hidden Markov model, or
CarHMM (Lawler et al.,, 2019). Secondly, classic HMMs fail to model
simultaneous behavioural processes that occur at different time scales (i.e.,
both between and within curves). To address this issue, statistical ecologists
have employed hierarchical hidden Markov models (HHMMs) (Leos-Barajas et al.,,
2017; Adam et al.,, 2019), which model both scales with conditionally
dependent HMMs. Thirdly, traditional HMMs, CarHMMs, and HHMMs cannot easily
capture complicated dependence structures at short time scales. For example,
Adam et al., (2019) fail to capture fine-scale periodic swimming pattern of
horn sharks (Heterodontus francisci) using a traditional HHMM. Heerah et al.,
(2017) successfully use Fourier analysis within an HMM to account for daily
behavioural cycles in marine mammals, and Fourier analysis has previously been
used with accelerometer data to explain animal behaviour (Fehlmann et al.,,
2017; Alex Shorter et al.,, 2017). Thus, incorporating Fourier analysis within
the structure of an HMM appears to be a promising approach to account for
fine-scale periodic structures.

In this paper, we combine existing methods from statistical ecology literature
in novel ways to account for complex, temporally dependent functional data.
The resulting suite of methods makes up a “tool box” that can be used to build
arbitrarily complex hierarchical models to explain multi-scale functional and
time-series data with intricate dependence structure. We begin in Section 2 by
describing HMMs as well as two variants, CarHMMs and HHMMs, and discuss how
Fourier analysis can handle fine-scale dependence structures. We also show how
these methods can be combined to analyze increasingly complex data. In Section
3 we fit several candidate models to data from a killer whale (Orcinus orca)
from the threatened northern resident population off the coast of British
Columbia, Canada. Section 4 details a simulation study based on these
candidate models, and in Section 5 we discuss our results.

##  2 MODELS AND PARAMETER ESTIMATION

Consider a sequence of $T$ curves, where any particular curve $t$ is
characterized by a curve-level (or coarse-scale) observation $Y_{t}$ as well
as a sequence of $T^{*}_{t}$ within-curve (or fine-scale) observations
$Y^{*}_{t}$. Namely,
$Y^{*}_{t}\equiv\big{\\{}Y^{*}_{t,1},\ldots,Y^{*}_{t,T^{*}_{t}}\big{\\}}$ is
made up fine-scale quantities derived from curve $t$ indexed using $t^{*}$.
Both $Y_{t}$ and $Y^{*}_{t,t^{*}}$ can be either vectors or scalars. We call
the sequence of coarse-scale observations
$Y\equiv\big{\\{}Y_{1},\ldots,Y_{T}\big{\\}}$ and the collection of all fine-
scale observations $Y^{*}\equiv\big{\\{}Y^{*}_{1},\ldots,Y^{*}_{T}\big{\\}}$.
To develop our model for this data, we detail the structure of a traditional
HMM followed by three variations which generalize its base structure. We then
show how each of these generalized HMMs can be synthesized to form a wide
variety of more complicated models.

###  2.1 HMMs as a base structure

Hidden Markov models describe state-switching Markovian processes in discrete
time and are the core structure we use to model both $Y$ and $Y^{*}$. For
simplicity we focus on $Y$ to introduce the model. An HMM is comprised of a
sequence of unobserved states $X\equiv\big{\\{}X_{1},\ldots,X_{T}\big{\\}}$
together with an observation sequence
$Y\equiv\big{\\{}Y_{1},\ldots,Y_{T}\big{\\}}$, where $X_{t}$ is associated
with the observation $Y_{t}$. The $Y_{t}$’s are often referred to as
“emissions” and the index $t$ typically refers to time. The $X_{t}$’s form a
Markov chain and can take integer values between $1$ and $N$. Their
distribution is governed by the distribution of the initial state $X_{1}$ and
the $N\times N$ transition probability matrix $\Gamma$, where
$\Gamma_{ij}=\Pr(X_{t+1}=j|X_{t}=i)$. We assume that $X_{1}$ follows the
chain’s stationary distribution, which is denoted by an $N$-dimensional row
vector $\delta$, where $\delta_{i}=\Pr(X_{1}=i).$ A Markov chain’s stationary
distribution is determined by its probability transition matrix via
$\delta=\delta\Gamma$ and $\sum_{i=1}^{N}\delta_{i}=1$. The distribution of an
emission $Y_{t}$ conditioned on the corresponding hidden state $X_{t}$ does
not depend upon any other observation or hidden state. If $X_{t}=i$ then we
denote the conditional density or probability mass function of $Y_{t}$ as
$f^{(i)}(\cdot;\theta^{(i)})$ or simply $f^{(i)}(\cdot)$, where $\theta^{(i)}$
is a state-dependent parameter describing the emission distribution.

Using observation emissions, denoted here as
$y\equiv\\{y_{1},\ldots,y_{T}\\}$, we can find the maximum likelihood
estimates of the parameters $\Gamma$ and
$\theta\equiv\\{\theta^{(1)},\ldots,\theta^{(N)}\\}$. The likelihood
$\mathcal{L}_{\text{HMM}}$ can be evaluated using the well-known forward
algorithm (Zucchini et al.,, 2016):

$\mathcal{L}_{\text{HMM}}(\theta,\Gamma;y)=\delta
P(y_{1};\theta)\prod_{t=2}^{T}\Gamma P(y_{t};\theta)\mathbf{1}_{N},$

where $\mathbf{1}_{N}$ is an $N$-dimensional column vector of ones and
$P(y_{t};\theta)$ is an $N\times N$ diagonal matrix with $(i,i)^{th}$ entry
$f^{(i)}(y_{t};\theta^{(i)})$.

Following Leos-Barajas et al., (2017), we reparameterize the $N\times N$
transition probability matrix $\Gamma$ such that the entries of the matrix are
forced to be non-negative and the rows sum to 1:

$\Gamma_{ij}=\frac{\exp(\eta_{ij})}{\sum_{k=1}^{N}\exp(\eta_{ik})},$

where $i,j=1,\ldots,N$ and $\eta_{ii}$ is set to zero for identifiability.
This simplifies likelihood maximization by removing constraints in the
optimization problem. For simplicity, we will continue to use $\Gamma$ in our
notation, suppressing the reparameterization in terms of $\eta$. Figure 1a
represents the dependence structure of an HMM.

(a) Hidden Markov Model (HMM)

(b) Conditionally Autoregressive HMM (CarHMM)

(c) Hierarchical HMM (HHMM)

(d) HMM with Discrete Fourier Transform (HMM-DFT)

Figure 1: Dependence structures of a standard HMM (a) and three HMM variants
(b, c, and d). Hidden state sequences are denoted as $X$ on the coarse-scale
and $X^{*}$ on the fine scale. Likewise, observations (or emissions) are
denoted by $Y$ on the coarse scale and $Y^{*}$ on the fine scale. The HHMM
shown here has a traditional HMM as both the coarse- and fine-scale model. In
(d), observations are transformed using a moving window and denoted as
$\tilde{Y}^{*}$ with corresponding hidden states $\tilde{X}^{*}$.

###  2.2 Relaxing conditional independence with the CarHMM

A conditionally auto-regressive hidden Markov model, or CarHMM (Lawler et
al.,, 2019), is a generalization of an HMM which explicitly models auto-
correlation in the observation sequence beyond the correlation induced by the
hidden state process. Like a traditional HMM, a CarHMM is made up of a Markov
chain of unobserved states $X_{1},\ldots,X_{T}$ that can take integer values
between $1$ and $N$. A CarHMM also has a transition probability matrix
$\Gamma$ and initial distribution $\delta$ equal to the stationary
distribution of $\Gamma$. Unlike a traditional HMM, the CarHMM assumes that
the distribution of $Y_{t}$ conditioned on $X_{1},\ldots,X_{T}$ and
$Y_{1},\ldots,Y_{t-1}$ depends on both $X_{t}$ and $Y_{t-1}$ rather than only
$X_{t}$. The first emission $Y_{1}$ is treated as a fixed initial value which
does not depend upon $X_{1}$. We denote the conditional density or probability
mass function of $Y_{t}$ given $Y_{t-1}=y_{t-1}$ and $X_{t}=i$ as
$f^{(i)}(\cdot|y_{t-1};\theta^{(i)})$ or simply $f^{(i)}(\cdot|y_{t-1})$. This
model is highly general, as $f^{(i)}(\cdot|y_{t-1})$ can be any valid density
or probability mass function that depends upon the parameters $\theta^{(i)}$
and the previous observation $Y_{t-1}$. As a concrete example, if $Y_{t}$ is a
scalar, then one may assume that $Y_{t}$ given $X_{t}=i$ is Normally
distributed with parameters
$\theta^{(i)}=\\{\mu^{(i)},\sigma^{(i)},\phi^{(i)}\\}$, where:

$\displaystyle\begin{split}\mathbb{E}(Y_{t}|Y_{t-1}=y_{t-1},X_{t}=i)&=\phi^{(i)}~{}y_{t-1}~{}+~{}(1-\phi^{(i)})~{}\mu^{(i)},\\\
\mathbb{V}(Y_{t}|Y_{t-1}=y_{t-1},X_{t}=i)&=(\sigma^{(i)})^{2}.\end{split}$ (1)

A CarHMM which follows Equation (1) can be viewed as a discrete time version
of a state-switching Ornstein-Uhlenbeck process (Michelot and Blackwell,,
2019). This follows in the same way that an AR(1) process is the discrete-time
version of a traditional Ornstein-Uhlenbeck process.

As previously, the likelihood corresponding to a general CarHMM can be easily
calculated using the forward algorithm. If $y$ is the sequence of observed
emissions, then

$\mathcal{L}_{\text{CarHMM}}(\theta,\Gamma;y)=\delta\prod_{t=2}^{T}\Gamma
P(y_{t}|y_{t-1};\theta)\mathbf{1}_{N},$

where $P(y_{t}|y_{t-1};\theta)$ is an $N\times N$ diagonal matrix with
$(i,i)^{th}$ entry equal to $f^{(i)}(y_{t}|y_{t-1};\theta^{(i)})$. Figure 1b
shows a graphical representation of the dependence structure of a CarHMM.

###  2.3 Incorporating multiple scales with the HHMM

A hierarchical hidden Markov model, or HHMM, accounts for processes occurring
simultaneously at different scales by modelling both the coarse-scale process
and fine-scale process with either HMMs (Leos-Barajas et al.,, 2017; Adam et
al.,, 2019) or CarHMMs. The coarse-scale model is either an HMM and CarHMM as
defined in Sections 2.1 and 2.2, where $X_{1},\ldots,X_{T}$ make up an
unobserved Markov chain with $N$ possible states and $Y_{1},\ldots,Y_{T}$ are
the corresponding observations with state-dependent parameters $\theta^{(i)}$
for $i=1,\ldots,N$. In the hierarchical setting, however, each state $X_{t}$
also emits another sequence of fine-scale unobserved states,
$X_{t}^{*}\equiv\\{X_{t,1}^{*},\ldots,X_{t,T_{t}^{*}}\\}$, which in turn emits
a sequence of fine-scale observations
$Y_{t}^{*}\equiv\\{Y_{t,1}^{*},\ldots,Y_{t,T_{t}^{*}}\\}$. For each curve $t$,
the fine-scale process $\\{X_{t}^{*},Y_{t}^{*}\\}$ then follows another HMM
(or CarHMM) whose parameters depend on the value of $X_{t}$. If $X_{t}=i$,
then the components of $X_{t}^{*}$ make up a Markov chain with $N^{*(i)}$
possible states, an $N^{*(i)}\times N^{*(i)}$ transition probability matrix
$\Gamma^{*(i)}$, and an initial distribution $\delta^{*(i)}$ which we assume
is equal to the stationary distribution of the chain. The distribution of
$Y^{*}_{t,t^{*}}$ given $Y^{*}_{t,t^{*}-1}=y^{*}_{t,t^{*}-1}$,
$X^{*}_{t,t^{*}}=i^{*}$, and $X_{t}=i$ is governed by the parameter
$\theta^{*(i,i^{*})}$ and has density or probability mass function denoted
$f^{*(i,i^{*})}\left(\cdot|y^{*}_{t,t^{*}-1};\theta^{*(i,i^{*})}\right)$ or
simply $f^{*(i,i^{*})}(\cdot|y^{*}_{t,t^{*}-1})$. We denote the set of fine-
scale emission parameters corresponding to $X_{t}=i$ as
$\theta^{*(i)}=\big{\\{}\theta^{*(i,1)},\ldots,\theta^{*\left(i,N^{*(i)}\right)}\big{\\}}$.
In summary:

$\displaystyle\\{Y,X\\}\text{ follows a (Car)HMM with
}\Gamma\in\mathbb{R}^{N\times N},$
$\displaystyle(Y_{t}|Y_{t-1}=y_{t-1},X_{t}=i)\text{ has density
}f^{(i)}(\cdot|y_{t-1};\theta^{(i)}),$
$\displaystyle\\{Y^{*}_{t},X^{*}_{t}|X_{t}=i\\}\text{ follows a (Car)HMM with
}\Gamma^{*(i)}\in\mathbb{R}^{N^{*(i)}\times N^{*(i)}},$
$\displaystyle(Y^{*}_{t,t^{*}}|Y^{*}_{t,t^{*}-1}=y^{*}_{t,t^{*}-1},X^{*}_{t,t^{*}}=i^{*},X_{t}=i)\text{
has density }f^{*(i,i^{*})}(\cdot|y^{*}_{t,t^{*}-1};\theta^{(i,i^{*})}).$

Given the coarse-scale hidden state sequence $X$, the $T+1$ sets
$\\{X_{1}^{*},Y_{1}^{*}\\},\ldots,\\{X_{T}^{*},Y_{T}^{*}\\}$, and
$\\{Y_{1},\ldots,Y_{T}\\}$ are assumed to be independent of one another.

Forcing certain parameters to be shared can reduce complexity and increase
interpretability of the HHMM. For example, in our killer whale case study (see
Section 3), we take $N^{*(i)}=N^{*}$ for all $i$. We also share the fine-scale
emission parameters across the $N$ coarse-scale hidden states
$\left(\text{i.e.,
}\theta^{*(1,i^{*})}=\cdots=\theta^{*(N,i^{*})}=\theta^{*(\cdot,i^{*})}\text{
for all }i^{*}=1,\ldots,N^{*}\right)$. Coarse-scale hidden states therefore
differ only in their coarse-scale emission parameters $\theta^{(i)}$ and fine-
scale probability transition matrices $\Gamma^{*(i)}$.

Due to the nested structure of the HHMM, the likelihood is easily calculated
using the forward algorithm. Let $y$ be the sequence of observed coarse-scale
emissions and $y^{*}\equiv\\{y^{*}_{1},\ldots,y^{*}_{T}\\}$ be the collection
of $T$ observed fine-scale emission vectors. In addition, let
$\theta^{*}\equiv\\{\theta^{*(1)},\ldots,\theta^{*(N)}\\}$ denote the
collection of all fine-scale emission parameters and
$\Gamma^{*}\equiv\\{\Gamma^{*(1)},\ldots,\Gamma^{*(N)}\\}$ denote the
collection of all fine-scale transition probability matrices. The likelihood
of the observed data is then

$\mathcal{L}_{\text{HHMM}}(\theta,\theta^{*},\Gamma,\Gamma^{*};y,y^{*})=\delta
P_{m}(y_{1},y_{1}^{*};\theta,\theta^{*},\Gamma^{*})\prod_{t=2}^{T}\Gamma
P_{m}(y_{t},y_{t}^{*}|y_{t-1};\theta,\theta^{*},\Gamma^{*})\mathbf{1}_{N}$ (2)

where $P_{m}(y_{t},y_{t}^{*}|y_{t-1};\theta,\theta^{*},\Gamma^{*})$ is an
$N\times N$ diagonal matrix whose exact structure depends upon the coarse- and
fine-scale models. If the coarse-scale model is an HMM,
$P_{m}(y_{1},y_{1}^{*};\theta,\theta^{*},\Gamma^{*})$ and
$P_{m}(y_{t},y_{t}^{*}|y_{t-1};\theta,\theta^{*},\Gamma^{*})$ for $t\geq 2$
both have $(i,i)^{th}$ entries equal to
$f^{(i)}(y_{t})\mathcal{L}_{\text{fine}}\left(\theta^{*(i)},\Gamma^{*(i)};y_{t}^{*}\right)$.
If the coarse-scale model is a CarHMM,
$P_{m}(y_{1},y_{1}^{*};\theta,\theta^{*},\Gamma^{*})$ has $(i,i)^{th}$ entry
equal to
$\mathcal{L}_{\text{fine}}\left(\theta^{*(i)},\Gamma^{*(i)};y_{1}^{*}\right)$
and $P_{m}(y_{t},y_{t}^{*}|y_{t-1};\theta,\theta^{*},\Gamma^{*})$ for $t\geq
2$ has $(i,i)^{th}$ entry equal to
$f^{(i)}(y_{t}|y_{t-1})\mathcal{L}_{\text{fine}}\left(\theta^{*(i)},\Gamma^{*(i)};y_{t}^{*}\right)$.
The fine-scale likelihood $\mathcal{L}_{\text{fine}}$ corresponds to the
likelihood of the fine-scale model, which can be either a CarHMM or an HMM.
Figure 1c graphically displays the dependence structure for an HHMM.

###  2.4 Transforming fine-scale observations with the HMM-DFT

In many applications where data are collected at high frequencies, intricate
dependency structures arise within the fine-scale process that cannot be
adequately modelled with the HMM-based models described thus far. To handle
these additional fine-scale structures, we recommend replacing
$Y_{t}^{*}=\\{Y^{*}_{t,1},\ldots,Y^{*}_{t,T^{*}_{t}}\\}$ with relevant
statistics that summarize any non-Markovian behaviour. To maintain the
temporal structure of the fine-scale process, local summary statistics can be
calculated from a moving window over the elements of $Y_{t}^{*}$. Subject
matter experts are often required to determine the specific summary statistics
employed as well as the optimal window size and stride length of the moving
window. Stride length refers to the distance between the first element of
consecutive windows, so a stride length of $h$ indicates that the first window
starts at $Y^{*}_{t,1}$, the second at $Y^{*}_{t,h+1}$, and so on. Larger
stride lengths result in a loss of information but also reduce the dimension
of the fine-scale process, which allows for faster model fitting. In addition,
setting the stride length equal to the window size avoids artificial residual
correlation arising from overlapping windows. For our case study, we use the
discrete Fourier transform (DFT) of a moving forward window of width $h$ and
stride $h$ across $Y^{*}_{t}$. Namely:

$\displaystyle
DFT\\{Y^{*}_{t,t^{*}},\ldots,Y^{*}_{t,t^{*}+h-1}\\}(k)=\sum_{n=0}^{h-1}Y^{*}_{t,t^{*}+n}\exp\left(-\frac{i2\pi}{h}kn\right)$
(3)

for $t^{*}=1,h+1,2h+1,\ldots$ and $k=0,1,\ldots,h-1$ with $i=\sqrt{-1}$. If
$Y^{*}_{t,t^{*}}$ is a vector then the DFT is taken component-wise. We omit
the final window if $t^{*}+h-1$ exceeds $T^{*}_{t}$, denote the total number
of windows as $\tilde{T}^{*}_{t}$, and index each window with
$\tilde{t}^{*}=1,\ldots,\tilde{T}^{*}_{t}$. Next, we calculate transformed
observations
$\tilde{Y}^{*}_{t,\tilde{t}^{*}}\equiv\\{\tilde{A}^{*}_{t,\tilde{t}^{*}},\tilde{W}^{*}_{t,\tilde{t}^{*}}\\}$:

$\tilde{A}^{*}_{t,\tilde{t}^{*}}\equiv\frac{1}{h}\sum_{n=1}^{h}Y^{*}_{t,h(\tilde{t}^{*}-1)+n},\quad\tilde{W}^{*}_{t,\tilde{t}^{*}}\equiv\sum_{k=1}^{\tilde{\omega}}\Big{|}\Big{|}DFT\\{Y^{*}_{t,h(\tilde{t}^{*}-1)+1},\ldots,Y^{*}_{t,h\tilde{t}^{*}}\\}(k)\Big{|}\Big{|}^{2},$
(4)

where $\tilde{\omega}\leq h-1$ is a problem-specific tuning parameter
corresponding to the maximum recorded frequency within each window. In words,
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ is the average value of $Y^{*}_{t}$ within
window $\tilde{t}^{*}$ and $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ is the squared
two-norm of the component of the window that can be attributed to frequencies
between one and $\tilde{\omega}$ periods per window. More intuitively,
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$ corresponds to the “wiggliness” of the fine-
scale data within curve $t$ and window $\tilde{t}^{*}$.

When performing this transformation, the fine-scale HMM (or CarHMM) must be
redefined since $\tilde{Y}^{*}_{t}$ exists on a coarser scale than $Y^{*}_{t}$
itself. As a result, there are only $\tilde{T}_{t}^{*}$ hidden states
associated with
$\tilde{Y}^{*}_{t}=\big{\\{}\tilde{Y}^{*}_{t,1},\ldots,\tilde{Y}^{*}_{t,\tilde{T}^{*}_{t}}\big{\\}}$,
which we denote as
$\tilde{X}^{*}_{t}=\big{\\{}\tilde{X}^{*}_{t,1},\ldots,\tilde{X}^{*}_{t,\tilde{T}^{*}_{t}}\big{\\}}$.
The fine-scale probability transition matrices $\Gamma^{*(i)}$ and probability
density functions $f^{*(i,i^{*})}$ are then applied directly to
$\tilde{X}^{*}_{t}$ and $\tilde{Y}^{*}_{t}$ instead of $X^{*}_{t}$ and
$Y^{*}_{t}$. The likelihood of this model is therefore identical to that of
the original HMM or CarHMM as defined in Sections 2.1 and 2.2, but $Y^{*}$ is
replaced with $\tilde{Y}^{*}$ and $X^{*}$ is replaced with $\tilde{X^{*}}$. To
clearly differentiate models, we refer to an HMM with $\tilde{Y}^{*}$ as
observations and $\tilde{X^{*}}$ as hidden states as an HMM-DFT. Figure 1d
displays the dependence structure of an HMM-DFT.

###  2.5 Generalized hierarchical Markov models

Traditional HHMMs treat both the coarse-scale and the fine-scale processes as
realizations of an HMM or CarHMM. However, the fine-scale observations for a
particular dive $Y^{*}_{t}$ can be modelled using a large variety of
parametric models which admit easy-to-compute likelihoods or penalized
likelihoods. As such, the fine-scale HMM likelihood term
$\mathcal{L}_{\text{fine}}$ in Equation (2) can be replaced by the likelihood
of a general fine-scale model whose parameters depend upon the coarse-scale
hidden state. Possible candidates for the fine-scale model include any of the
models described in the previous subsections in addition to many others not
described here. For example, Bebbington, (2007) and Borchers et al., (2013)
investigate data sets with count onsets as observations, so they use
variations of a Poisson process as their fine-scale model. If the fine-scale
model is a simple Poisson process, then this approach is equivalent to a
Markov-modulated Poisson process (Fischer and Meier-Hellstern,, 1993). The
fine-scale process can also be modelled as in (Langrock et al.,, 2018), who
use B-splines to model the emission distribution of an HMM. This non-
parametric approach uses a penalized likelihood term which can easily replace
the usual fine-scale likelihood term in Equation (2). Another class of fine-
scale models is the set of continuous time methods such as a continuous-time
HMM (CTHMM) (Liu et al.,, 2015) or a state-switching Ornstein-Uhlenbeck
process (Michelot and Blackwell,, 2019). A continuous-time HMM may be
appropriate if observations are not equi-spaced in time (Liu et al.,, 2015).
Xu et al., (2020) model high-frequency biologging accelerometer data of
individuals by incorporating a CTHMM into a hierarchical model similar to
ours. However, they assume that individuals are partitioned into subgroups a
priori whereas we use an HMM to infer the coarse-scale hidden states.

These examples include a few of many fine-scale models that can act as initial
building blocks in a practitioner’s toolbox to construct increasingly complex
hierarchical models. A myriad of possible models can be built using this
framework, but these models can quickly become complicated and computationally
expensive to fit. Therefore, models should be constructed with care to achieve
an adequate fit of the data while avoiding over-fitting and high computational
costs.

##  3 KILLER WHALE CASE STUDY

To illustrate the process of constructing a model using these building blocks,
we analyze the dive behaviour of a northern resident killer whale in Queen
Charlotte Sound off the coast of British Columbia, Canada, and construct
several candidate models to categorize and describe its diving behaviour.

Understanding animal behaviour is important for conservation efforts, as
environmental changes caused by anthropogenic activity can directly impact
animal behaviour (Sutherland,, 1998). HMMs have been used to understand how
diving behaviours of various species are affected by disturbances (e.g.,
DeRuiter et al., (2017) and Isojunno et al., (2017)). For killer whales, we
are interested in categorizing different diving behaviours to identify
potential foraging dives. Northern resident killer whales feed almost
exclusively on calorie-rich Chinook salmon (Oncorhynchus tshawytscha) (Ford
and Ellis,, 2006), which typically occur deeper and are less numerous than
smaller types of salmon (Ford et al.,, 2009). Northern resident killer whales
therefore must expend significant amounts of energy to capture Chinook
(Williams and Noren,, 2009; Noren,, 2011; Wright et al.,, 2017). Acceleration
data can be used to estimate an animal’s energy expenditure (Green et al.,,
2009; Wilson et al.,, 2019), but the animal’s behavioural state must be
accounted for in order to obtain accurate estimates (Jeanniard du Dot et al.,,
2016). Therefore, understanding both the behavioural state of the killer whale
and the distribution of acceleration within each behavioural state is needed
to determine the true energetic requirements of the animal.

###  3.1 Data collection and preprocessing

The data we use were collected on September 2, 2019 from 12:49 pm to 6:06 pm
PDT and consist of depth and acceleration over time. Observations were
collected at a rate of 50 Hz using a CATs biologger (Customizable Animal
Tracking Solutions, www.cats.is). Acceleration was measured in three
dimensions, which together represent the complete range of movement of an
animal (forward/backward, upward/downward, and right/left). Tri-axial
acceleration readings are common in these types of tags and are often used to
infer animal behaviour such as foraging (Cade et al.,, 2018; Fehlmann et al.,,
2017; Wright et al.,, 2017). The act of attaching and detaching the tag caused
anomalous behaviour before 1:20 pm and after 6:00 pm, so observations taken
during these time periods are ignored. There were also periods of time when
the tag failed to record observations, resulting in data gaps between 2:25 pm
and 2:37 pm and between 4:07 pm and 5:07 pm. To preprocess the data, we smooth
the depth and acceleration curves by taking a moving average within a window
of $1/10^{th}$ of a second. We then define a killer whale “dive” as any
continuous interval of data that occurs below 0.5 meters in depth and lasts
for at least 10 seconds. Data are preprocessed in part with the divebomb
package in Python (Nunes,, 2019). The preprocessed data contain a total of 267
dives, all of which are displayed in Figure 2. Each dive is treated as one
curve, and the sequence of dives makes up the coarse-scale process.
Specifically, the sequence of observed coarse-scale observations
$y=\big{\\{}y_{1},\ldots,y_{267}\big{\\}}$ is a time series of dive durations
in seconds. For dive $t$, the fine-scale observations are contained in
$y^{*}_{t}\equiv\big{\\{}y^{*}_{t,1},\ldots,y^{*}_{t,T^{*}_{t}}\big{\\}}$,
which is a sequence containing the within-dive acceleration data in units of
meters per second squared. The collection of all acceleration data is denoted
as $y^{*}=\big{\\{}y^{*}_{1},\ldots,y^{*}_{267}\big{\\}}$.

Figure 2: Dive depth (top panel) and three-dimensional acceleration (bottom
three panels) from a killer whale over approximately 5 hours. An exact
physical interpretation of each of component acceleration is difficult due to
variations in tag orientation. There are data gaps occurring from around $1.5$
to $1.8$ hours and from around $3.2$ to $4.5$ hours. Both data gaps were
excluded from the analysis.

###  3.2 Model definition and selection

Defining a suitable model to describe this killer whale kinematic data
involves selecting an appropriate number of hidden states, model structure,
and emission distributions for both the coarse- and fine-scale observations.

We do not use information criteria to select the number of dive types $N$
since these metrics tend to overestimate the number of behavioural states in
biological processes (Pohle et al.,, 2017). We instead plot the duration of
each dive versus the duration of the dive preceding it ($y_{t}$ versus
$y_{t-1}$ for $t\geq 2$). This type of visualization is known as a lag plot.
If the emission distributions of the hidden states are well-separated, a lag
plot should reveal $N$ distinct patterns, where each pattern corresponds to
one dive type (Lawler et al.,, 2019). This is unfortunately not the case for
our killer whale data, as there is one cluster of data centred at
approximately $y_{t}=y_{t-1}=30$ seconds. However, longer dives appear to be
characterized by bouts of less “wiggly” behaviour in the acceleration data
compared to shorter dives, so we choose $N=2$ to differentiate these dive
types. The absence of a more principled method to select $N$ highlights the
importance of model validation techniques in lieu of information criteria (see
Section 3.4). Prior to fitting the model, lag plots reveal no significant
auto-correlation between dive duration observations (see Figure 1 of the
supplementary material), and visual inspection shows no obvious complicated
dependence. Therefore, we select a simple HMM to model the coarse-scale
process since neither a CarHMM nor a moving-window transformation is called
for. Given that dive $t$ is of type $i$, we assume that the dive duration
$Y_{t}$ follows a Gamma distribution with unknown parameters $\mu^{(i)}$ and
$\sigma^{(i)}$:

$\mathbb{E}(Y_{t}|X_{t}=i)=\mu^{(i)},\qquad\mathbb{V}(Y_{t}|X_{t}=i)=\left(\sigma^{(i)}\right)^{2}.$

This is consistent with previous studies, including that of Leos-Barajas et
al., (2017).

We then select a model corresponding to the fine-scale observations of
acceleration. Similarly to the coarse model, we rely on lag plots and visual
inspection to select $N^{*}=3$ subdive states. Although $N^{*}$ is selected
heuristically, we test the validity of this model in Section 3.4. In contrast
to the coarse-scale observations, the fine-scale acceleration data exhibit
significant sinusoidal behaviour. Thus, we transform each fine-scale
observation sequence $y_{t}^{*}$ into $\tilde{y}^{*}_{t}$ using Equation (4)
with a window size of $h=100$ (two seconds) and a maximum frequency of
$\tilde{\omega}=10$ (5 Hz). We then have that
$\tilde{y}^{*}_{t,\tilde{t}^{*}}=\\{\tilde{a}^{*}_{t,\tilde{t}^{*}},\tilde{w}^{*}_{t,\tilde{t}^{*}}\\}$,
where $\tilde{a}^{*}_{t,\tilde{t}^{*}}$ is a three dimensional vector of
component-wise average acceleration and $\tilde{w}^{*}_{t,\tilde{t}^{*}}$ is a
scalar describing the “wiggliness” of a particular window. Even after
transforming the raw acceleration data, there is still strong auto-correlation
within each component of $\tilde{a}^{*}_{t,\tilde{t}^{*}}$ prior to fitting
the model (see Figure 1 of the supplementary material). Therefore, we choose a
CarHMM as defined in Section 2.2 as the fine-scale model.

We then select the specific emission distribution of
$\tilde{Y}^{*}_{t,\tilde{t}^{*}}$. First, we assume that
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$ and all three components of
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ are independent of one another when
conditioned on the dive types and subdive states. To reduce model complexity,
we also assume that the three sets of fine-scale emission parameters are
shared across the two dive types
$\left(\theta^{*(1,i^{*})}=\theta^{*(2,i^{*})}\equiv\theta^{*(\cdot,i^{*})}\text{
for }i^{*}=1,2,3\right)$. This implies that the subdive states within dive
type 1 have the same interpretation as those within dive type 2. To specify
the emission distribution of $\tilde{A}^{*}_{t,\tilde{t}^{*}}$, consider the
sequence
$\big{\\{}\tilde{A}^{*}_{t,1},\ldots,\tilde{A}^{*}_{t,T^{*}_{t}}\big{\\}}$ for
a particular dive $t$. We assume that each of the three components of this
sequence are Normally distributed as in Equation (1) and all components are
independent of one another when conditioned on the subdive states
$\big{\\{}\tilde{X}_{t,1},\ldots,\tilde{X}_{t,\tilde{T}^{*}_{t}}\big{\\}}$.
Each component is assumed to have its own mean and variance parameters, but
all components share the same auto-correlation parameter. Thus the
distribution of $\tilde{A}^{*}_{t,\tilde{t}^{*}}$ given
$\tilde{A}^{*}_{t,\tilde{t}^{*}-1}$ and $X^{*}_{t,\tilde{t}^{*}}=i^{*}$ has
parameters $\mathbf{\mu}_{A}^{*(\cdot,i^{*})}\in\mathbb{R}^{3}$,
$\mathbf{\sigma}_{A}^{*(\cdot,i^{*})}\in\mathbb{R}^{3}$, and
$\phi_{A}^{*(\cdot,i^{*})}\in[0,1]$. To specify the emission distribution of
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$, we assume that given
$\tilde{X}^{*}_{t,\tilde{t}^{*}}=i^{*}$, $\tilde{W}^{*}_{t,\tilde{t}^{*}}$
follows a Gamma distribution parameterized by its mean
$\mu_{W}^{*(\cdot,i^{*})}$ and standard deviation
$\sigma_{W}^{*(\cdot,i^{*})}$. In addition,
$\tilde{W}^{*}_{t,1},\ldots,\tilde{W}^{*}_{t,T^{*}_{t}}$ are assumed to be
independent of one another given the subdive state sequence
$\big{\\{}\tilde{X}_{t,1},\ldots,\tilde{X}_{t,T_{t}^{*}}\big{\\}}$. We do not
include $\tilde{W}^{*}_{t,\tilde{t}^{*}-1}$ in the distribution of
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$ because the auto-correlation evident from
the lag plot is not severe and may be explained by subsequent observations
occurring within the same subdive state.

In total, the parameters to estimate are

$\displaystyle\Gamma,\qquad\Gamma^{*}=\\{\Gamma^{*(1)},\Gamma^{*(2)}\\}\qquad\text{(probability
transition matrices)},$
$\displaystyle\theta=\\{\mu^{(1)},\sigma^{(1)},\mu^{(2)},\sigma^{(2)}\\}\qquad\text{($Y$
emission parameters), and}$
$\displaystyle\theta^{*}=\\{\theta^{*(\cdot,1)},\theta^{*(\cdot,2)},\theta^{*(\cdot,3)}\\}\qquad\text{($\tilde{Y}^{*}$
emission parameters), where}$
$\displaystyle\theta^{*(\cdot,i^{*})}=\\{\mu_{A}^{*(\cdot,i^{*})},\sigma_{A}^{*(\cdot,i^{*})},\phi_{A}^{*(\cdot,i^{*})},\mu_{W}^{*(\cdot,i^{*})},\sigma_{W}^{*(\cdot,i^{*})}\\}.$

Recall that $\theta^{*(\cdot,i^{*})}$ is the set of parameters describing the
distribution of $\tilde{Y}^{*}_{t,\tilde{t}^{*}}$ conditioned on
$\tilde{X}^{*}_{t,\tilde{t}^{*}}=i^{*}$. We refer to this final model as the
CarHHMM-DFT since it includes a CarHMM, hierarchical HMM, and DFT-based
transformation. The likelihood of this model is easily calculated using the
forward algorithm and can be maximized with respect to the parameters above
(see the appendix for details). Figure 3 shows the dependence structure of the
full CarHHMM-DFT.

Figure 3: Graphical representation of the conditionally auto-regressive
hierarchical hidden Markov model with discrete Fourier transform (CarHHMM-DFT)
used in the simulation and case study. The type of dive $t$ is denoted by
$X_{t}$ and $Y_{t}$ represents the associated dive duration. The raw
acceleration vector associated with dive $t$ and time stamp $t^{*}$ is denoted
by $Y^{*}_{t,t^{*}}$. The subdive state of the killer whale during dive $t$
and window $\tilde{t}^{*}$ is denoted as $\tilde{X}^{*}_{t,\tilde{t}^{*}}$,
and the corresponding transformed observation is denoted by
$\tilde{Y}^{*}_{t,\tilde{t}^{*}}$.

In addition to the CarHHMM-DFT, we consider three variations for comparison.
As in the full model, each of the following models assume that all components
of $\tilde{Y}^{*}_{t,\tilde{t}^{*}}$ are conditionally independent of one
another given the dive types and subdive states:

  1. 1.

An HHMM-DFT, which models the coarse-scale observations with an HMM and
transforms the fine-scale observations using Equation (4), but models
$\tilde{Y}^{*}_{t,\tilde{t}^{*}}$ as emissions from a simple HMM rather than a
CarHMM.

  2. 2.

A CarHHMM, which models the coarse-scale observations with an HMM, transforms
the fine-scale observations using Equation (4), and models
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ as emissions of a CarHMM. However, the
“wiggliness” $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ is omitted from this model
altogether.

  3. 3.

A CarHMM-DFT, which models the coarse-scale observations as an independent and
identically distributed sequence of dives, transforms the fine-scale
observations using Equation (4), and models $\tilde{Y}^{*}_{t,\tilde{t}^{*}}$
as emissions of a CarHMM. This model assumes that there is only one dive type.

Each of the three candidate models above leaves out one important aspect of
the full CarHHMM-DFT: the HHMM-DFT assumes there is no auto-correlation
between fine-scale observations, the CarHHMM does not incorporate “wiggliness”
$\big{(}\tilde{W}^{*}_{t,\tilde{t}^{*}}\big{)}$, and the CarHMM-DFT lacks a
hierarchical structure and thus does not distinguish between dive types.

###  3.3 Case study results

To illustrate an application of this method and compare the candidate models,
we fit all four models to the data shown in Figure 2. We first report the
results from the full CarHHMM-DFT in detail and assess the quality of the fit.
We then compare these results with those from the other candidate models.

(a) Estimated Gamma probability density functions of a killer whale’s dive
duration ($Y_{t}$) corresponding to dive types 1 and 2.

(b) Estimated Normal conditional densities of the three components of a killer
whale’s acceleration
$\left(\tilde{A}^{*}_{t,\tilde{t}^{*}}|\tilde{A}^{*}_{t,\tilde{t}^{*}-1}=\mu_{A}^{*(\cdot,i^{*})}\right)$
plotted on a linear scale and estimated Gamma probability density of the
killer whale’s wiggliness $\left(\tilde{W}^{*}_{t,\tilde{t}^{*}}\right)$
plotted on a log-log scale. Density functions correspond to subdive states
$i^{*}=1,2,$ and $3$. Density functions corresponding to acceleration are
conditioned on $\tilde{A}^{*}_{t,\tilde{t}^{*}-1}=\mu_{A}^{*(\cdot,i^{*})}$
because the density function of $\tilde{A}^{*}_{t,\tilde{t}^{*}}$ itself
depends on $\tilde{A}^{*}_{t,\tilde{t}^{*}-1}$.

Figure 4: Estimated probability density functions for coarse-scale and fine-scale emissions corresponding to observations of killer whale behaviour. Densities are estimated by fitting the CarHHMM-DFT to the case study data (see Table 1). Table 1: Estimates with standard errors for the parameters of the distributions of dive duration $\big{(}Y_{t}\big{)}$, acceleration $\big{(}\tilde{A}^{*}_{t,\tilde{t}^{*}}\big{)}$, and wiggliness $\big{(}\tilde{W}^{*}_{t,\tilde{t}^{*}}\big{)}$ of the killer whale kinematic data using the full CarHHMM-DFT. The parentheses refer to standard errors estimated using the observed information matrix. Feature | Dive / Subdive Type | Parameter Estimate (Standard Error)
---|---|---
$\hat{\mu}$ | $\hat{\sigma}$ | $\hat{\phi}$
Dive Duration $(s)$ – $Y_{t}$ | 1 | $25.68\ (0.60)$ | $9.57\ (0.51)$ | —
2 | $104.6\ (9.4)$ | $64.7\ (7.5)$ | —
$x$-Acc. $(m/s^{2})$ – $\left(\tilde{A}^{*}_{t,\tilde{t}^{*}}\right)_{x}$ | 1 | $0.020\ (0.042)$ | $0.034\ (0.001)$ | $0.976\ (0.007)$
2 | $0.244\ (0.013)$ | $0.079\ (0.001)$ | $0.886\ (0.005)$
3 | $0.218\ (0.028)$ | $0.265\ (0.007)$ | $0.626\ (0.029)$
$y$-Acc. $(m/s^{2})$ – $\left(\tilde{A}^{*}_{t,\tilde{t}^{*}}\right)_{y}$ | 1 | $0.469\ (0.052)$ | $0.044\ (0.001)$ | $0.976\ (0.009)$
2 | $0.436\ (0.014)$ | $0.082\ (0.001)$ | $0.886\ (0.012)$
3 | $0.384\ (0.033)$ | $0.321\ (0.009)$ | $0.626\ (0.034)$
$z$-Acc. $(m/s^{2})$ – $\left(\tilde{A}^{*}_{t,\tilde{t}^{*}}\right)_{z}$ | 1 | $-0.683\ (0.061)$ | $0.052\ (0.001)$ | $0.976\ (0.005)$
2 | $-0.593\ (0.016)$ | $0.096\ (0.001)$ | $0.886\ (0.009)$
3 | $-0.366\ (0.033)$ | $0.317\ (0.009)$ | $0.626\ (0.033)$
Wiggliness - $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ | 1 | $23.34\ (0.29)$ | $12.95\ (0.27)$ | —
2 | $301.2\ (3.2)$ | $330.1\ (4.2)$ | —
3 | $10200\ (210)$ | $15300\ (350)$ | —

The coarse-scale parameter estimates suggest that the killer whale has at
least two distinct dive behaviours (see Table 1 and Figure 4(a)). Dive type 1
corresponds to shorter, shallower dives which likely reflect resting,
travelling, and to a lesser extent searching for prey. Dive type 2 is longer,
deeper, and may be associated with behaviours such as hunting (Tennessen et
al., 2019a, ), but it is unclear whether any of the dives in this study are
successful foraging dives. No dive in this data set has a maximum depth deeper
than 30 meters and a study by Wright et al., (2017) of killer whales in
Johnstone Strait found that most prey captures occur at depths deeper than 100
meters. However, the killer whale studied here was tagged north of Johnstone
Strait in Queen Charlotte Sound, and significant numbers of fish have been
observed to be caught near the surface (Fortune and Trites, unpublished data).
Dive type 2 could also be associated with behaviours such as socializing that
can take place several meters below the surface (Tennessen et al., 2019a, ).

The means of “wiggliness” $\big{(}\tilde{W}^{*}_{t,\tilde{t}^{*}}\big{)}$
associated with each subdive state are separated by an order of magnitude (see
Table 1 and Figure 4(b)). Subdive state 1 has the smallest mean corresponding
to $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ and the smallest variance corresponding
to $\tilde{A}^{*}_{t,\tilde{t}^{*}}$. It also has the highest auto-correlation
in $\tilde{A}^{*}_{t,\tilde{t}^{*}}$. This implies less overall activity and
more consistent acceleration compared to the other subdive states. Subdive
state 2 has a mean “wiggliness”
$\big{(}\tilde{W}^{*}_{t,\tilde{t}^{*}}\big{)}$ one order of magnitude higher
than subdive state 1 and its acceleration has about twice the variance
compared to subdive state 1. The auto-correlation of acceleration is also
slightly lower than subdive state 1. We therefore hypothesize that subdive
state 2 corresponds to fluking (active swimming), as strong sinusoidal
behaviour in acceleration is characteristic of this behaviour in marine
mammals (Simon et al.,, 2012). Finally, the mean of
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$ and variance of
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ in subdive state 3 are both much higher than
in the other two states, and the auto-correlation of
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ is also much lower. This corresponds to
vigorous swimming activity, especially as the killer whale begins or ends a
dive (see Figure 5).

Figure 5: The $x$–component of acceleration $\left(y^{*}_{t,t^{*}}\right)_{x}$
(top two panels) and dive depth (bottom two panels) of a northern resident
killer whale for a sequence of six selected dives. Each panel is partitioned
into dives by vertical black lines. The curve colours in the first and third
panels correspond to the estimated dive types while the curve colours of the
second and fourth panels correspond the estimated subdive states. Both the
dive types and subdive states are estimated by fitting the CarHHMM-DFT to the
data and performing the forward-backward algorithm to determine the hidden
state with the highest probability.

The estimated probability transition matrices and associated stationary
distributions on the coarse scale are

$\hat{\Gamma}=\begin{pmatrix}0.788&0.212\\\ 0.809&0.191\end{pmatrix}\text{ and
}$ $\hat{\delta}=\begin{pmatrix}0.792&0.208\end{pmatrix}$

for the transitions between dives. The estimated probability transition
matrices and stationary distributions on the fine scale are

$\hat{\Gamma}^{*(1)}=\begin{pmatrix}0.679&0.321&0.000\\\ 0.038&0.904&0.058\\\
0.000&0.232&0.768\end{pmatrix},\qquad\hat{\Gamma}^{*(2)}=\begin{pmatrix}0.859&0.141&0.000\\\
0.114&0.841&0.045\\\ 0.000&0.216&0.784\end{pmatrix},$
$\hat{\delta}^{*(1)}=\begin{pmatrix}0.087&0.731&0.182\end{pmatrix},\enspace\text{and}\enspace\hat{\delta}^{*(2)}=\begin{pmatrix}0.401&0.496&0.103\end{pmatrix}$

for dive types 1 and 2. In summary, about 79% of dives are short dives of type
1, and the whale performs an average of 4.72 short type 1 dives before
switching to dive type 2 and an average of 1.24 longer type 2 dives before
switching back to dive type 1. This finding is consistent with those of
Tennessen et al., 2019a  and Williams and Noren, (2009), both of whom describe
common bouts of short resting dives before a killer whale performs a longer,
more energy-intensive deep dive. Further, this killer whale is in the less
active subdive state 1 40% of the time during a dive of type 1 compared to
only only 9% of the time during a dive of type 2. Less active swimming
behaviour is consistent with the need for marine mammals to conserve energy
when diving to depth and holding their breath for long periods (Williams et
al.,, 1999; Hastie et al.,, 2006). Figure 5 shows the decoded dive behaviour
of six selected dives, and Section 1.3 of the supplementary material also
shows the probability of each dive type and subdive state given the data and
the fitted model.

###  3.4 Model validation

We use two visual tools to evaluate the CarHHMM-DFT: pseudoresidual plots and
empirical histograms. The pseudoresidual of a coarse-scale observation $y_{t}$
is equal to
$\Phi^{-1}\left(\Pr(Y_{t}<y_{t}|\\{Y_{1},\ldots,Y_{T},\tilde{Y}^{*}_{1},\ldots,\tilde{Y}^{*}_{T}\\}/\\{Y_{t}\\})\right)$
and the pseudoresidual of a fine-scale observation
$\tilde{y}^{*}_{t,\tilde{t}^{*}}$ is
$\Phi^{-1}\left(\Pr(\tilde{Y}^{*}_{t,t^{*}}<\tilde{y}^{*}_{t,\tilde{t}^{*}}|\\{Y_{1},\ldots,Y_{T},\tilde{Y}^{*}_{1},\ldots,\tilde{Y}^{*}_{T}\\}/\\{\tilde{Y}^{*}_{t,\tilde{t}}\\})\right)$,
where $\Phi$ is the cumulative distribution function of a standard Normal
distribution. If the model is correct, then all pseudoresiduals are
independent and follow a standard Normal distribution. Histograms of the
pseudoresiduals mostly support that the CarHHMM-DFT is well-specified. One
exception is $\tilde{W}^{*}_{t,\tilde{t}^{*}}$, whose pseudoresiduals are
noticeably right-skewed (see Figure 6(a)). This implies that the true
distribution of $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ may follow a heavier-tailed
distribution than the Gamma distribution used in the case study. See Sections
1.4 through 1.6 of the supplementary material for pseudoresidual plots
corresponding to all observations and models.

We also plot histograms of dive duration corresponding to each dive type in
Figure 6(b). Each observation of dive duration is weighted by the estimated
probability that it corresponds to a particular dive type as decoded by the
forward-backward algorithm (Zucchini et al.,, 2016). This procedure results in
two histograms – one corresponding to dive type 1 and another corresponding to
dive type 2. Each histogram is then plotted together with the corresponding
emission distribution estimated by the CarHHMM-DFT. Analogous histograms
corresponding to the fine-scale observations are contained in Section 1 of the
supplementary material. Our results mostly show that the CarHHMM-DFT explains
the data well, but there are some exceptions. In histograms corresponding to
subdive state 3, $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ is right-skewed and
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ has heavier tails compared to a normal
distribution, indicating the existence of rare events corresponding to
exceptionally sudden changes in acceleration of the killer whale. These
outliers are potential subjects for future study and may indicate biologically
relevant phenomena such as prey capture (Tennessen et al., 2019b, ).

(a) Histogram of pseudoresiduals of $\tilde{W}^{*}_{t,\tilde{t}^{*}}$

(b) Empirical distribution of $Y_{t}$

Figure 6: Pseudoresiduals of wiggliness
($\Phi^{-1}\big{(}\Pr(\tilde{W}^{*}_{t,\tilde{t}^{*}}<\tilde{w}^{*}_{t,\tilde{t}^{*}}|Y,\tilde{Y}^{*}/\\{\tilde{Y}^{*}_{t,\tilde{t}^{*}}\\})\big{)}$,
left) plotted over a standard normal density as well as a weighted empirical
distribution of dive duration ($Y_{t}$, right) plotted over the corresponding
fitted Gamma distributions. Both plots are generated by fitting the CarHHMM-
DFT to the killer whale case study data and performing the forward-backward
algorithm.

###  3.5 Comparison with candidate models

The HHMM-DFT, which ignores auto-correlation in acceleration, decodes dive
types and subdive states similarly to the CarHHMM-DFT, but it is less likely
to categorize the behaviour at the beginning and end of dives as subdive state
3 (see Figures 4 and 5 of the supplementary material). In addition, for all
three components of each of $\sigma_{A}^{*(\cdot,1)}$,
$\sigma_{A}^{*(\cdot,2)}$, and $\sigma_{A}^{*(\cdot,3)}$, the HHMM-DFT
produces estimates which are approximately $50$ to $100$ percent larger than
those of the CarHHMM-DFT. The estimated uncertainties of the three components
of each of $\hat{\mu}_{A}^{*(\cdot,1)}$, $\hat{\mu}_{A}^{*(\cdot,2)}$, and
$\hat{\mu}_{A}^{*(\cdot,3)}$ are also less than half of those for the CarHHMM-
DFT (see Tables 1 and 2 of the supplementary material). This suggests that
including auto-correlation in the model significantly affects parameter
estimates. Further, the pseudoresiduals of the HHMM-DFT are noticeably light-
tailed and do not follow a standard normal distribution (see Figure 15 of the
supplementary material). These findings suggest that the HHMM-DFT is a
significantly worse fit to these data than the full CarHHMM-DFT.

The CarHHMM does not model the “wiggliness” of the acceleration data, so it
regularly fails to pick up obvious behavioural changes corresponding to the
periodicity shown in Figure 5 (see Figure 6 of the supplementary material).
These results essentially disqualify the CarHHMM as a viable model for this
data set. The pseudoresiduals of acceleration are also light-tailed relative
to a Normal distribution (see Figure 16 of the supplementary material).

Finally, the CarHMM-DFT, which lacks a hierarchical structure, produces fine-
scale parameter estimates and subdive state estimates similar to those of the
CarHHMM-DFT. However, its lack of hierarchical structure means that it fails
to differentiate between short and long dives. This model therefore does not
infer the dive-level Markov chain or the relationship between the dive and
subdive levels. In particular, the CarHMM-DFT does not indicate that the whale
is more likely to be in subdive state 1 when engaged in longer dives compared
to shorter dives.

A more complete set of results for each of the candidate models is presented
in Section 1 of the supplementary material.

##  4 SIMULATION STUDY

We perform a simulation study based on data generated from the full CarHHMM-
DFT as defined in Section 3.2 to evaluate each candidate model when the
ground-truth is known. The parameters used to generate the data are based on
those estimated in the case study (see Table 1), with slight modifications
made for simplicity. In particular, we set the number of subdive states to
$N^{*}=2$ and $\tilde{A}^{*}_{t,\tilde{t}^{*}}$ to a scalar instead of a three
dimensional vector. We then fit all four models to the simulated data. Metrics
used to evaluate each model include decoding accuracy of hidden states, bias
in parameter estimates, empirical standard errors of parameter estimates, and
fitting times. To assess the accuracy of uncertainty estimates, we also
compare the empirical standard errors of a given model’s parameter estimates
with the standard errors estimated using the inverse of the observed Fisher
information.

###  4.1 Simulation procedure

We generate 500 independent training data sets using the CarHHMM-DFT as a
generative model. Each training data set consists of a sequence of 100 curves
which we call a sequence of killer whale dives. Each dive can be one of $N=2$
dive types based on a Markov chain with probability transition matrix

$\Gamma=\begin{pmatrix}0.79&0.21\\\ 0.81&0.19\end{pmatrix}.$

Dive duration is Gamma distributed and the coarse-scale emission parameters
are

$\mu^{(1)}=25.7s,\enspace\sigma^{(1)}=9.6s,\enspace\mu^{(2)}=104.6s,\enspace\sigma^{(2)}=64.7s.$

After generating the dive durations for all 100 dives in a data set, dive $t$
is broken into a sequence of $\tilde{T}^{*}_{t}=\lfloor Y_{t}/2\rfloor$ two-
second windows, where the last $Y_{t}-2\tilde{T}^{*}_{t}$ seconds of each
simulated dive are ignored. Each two-second segment is assigned one of
$N^{*}=2$ behaviours according to a fine-scale Markov chain
$\tilde{X}^{*}_{t}\equiv\big{\\{}\tilde{X}^{*}_{t,1},\ldots,\tilde{X}^{*}_{t,\tilde{T}^{*}_{t}}\big{\\}}$
with probability transition matrices

$\Gamma^{*(1)}=\begin{pmatrix}0.68&0.32\\\ 0.05&0.95\end{pmatrix}\quad\text{
and }\quad\Gamma^{*(2)}=\begin{pmatrix}0.86&0.14\\\ 0.15&0.85\end{pmatrix}$

for dive types 1 and 2, respectively. Instead of generating the raw
observations $Y^{*}_{t,t^{*}}$, we directly simulate the fine-scale
transformed observations
$\tilde{Y}^{*}_{t,\tilde{t}^{*}}=\big{\\{}\tilde{A}^{*}_{t,\tilde{t}^{*}},\tilde{W}^{*}_{t,\tilde{t}^{*}}\big{\\}}$.
Recall from Section 3.2 that we must specify the mean, standard deviation, and
auto-correlation parameters corresponding to
$\big{\\{}\tilde{A}^{*}_{t,1},\ldots,\tilde{A}^{*}_{t,\tilde{T}_{t}^{*}}\big{\\}}$
as well as the mean and standard deviation parameters corresponding to
$\big{\\{}\tilde{W}^{*}_{t,1},\ldots,\tilde{W}^{*}_{t,\tilde{T}_{t}^{*}}\big{\\}}$.
We select the following parameters in line with the results from the case
study:

$\displaystyle\mu_{A}^{*(\cdot,1)}=0.0s,\enspace\sigma_{A}^{*(\cdot,1)}=0.034s,\enspace\phi_{A}^{*(\cdot,1)}=0.98,$
$\displaystyle\mu_{A}^{*(\cdot,2)}=0.0s,\enspace\sigma_{A}^{*(\cdot,2)}=0.079s,\enspace\phi_{A}^{*(\cdot,2)}=0.87,$
$\displaystyle\mu_{W}^{*(\cdot,1)}=23.3,\quad\sigma_{W}^{*(\cdot,1)}=13.0,$
$\displaystyle\mu_{W}^{*(\cdot,2)}=301.2,\quad\sigma_{W}^{*(\cdot,2)}=330.1.$

It is not possible to uniquely reconstruct the raw accelerometer data $Y^{*}$
from $\tilde{Y}^{*}$ alone, but we describe one possible mapping from
$\tilde{Y}^{*}$ to $Y^{*}$ in the appendix. Figure 21 of the supplementary
material shows one realization of $\tilde{Y}^{*}$ for five dives of one
simulated data set along with the corresponding reconstructed realization of
$Y^{*}$.

The two simulated dive types differ in that dives of type 1 are much shorter
on average (26 seconds) than dives of type 2 (105 seconds). The two simulated
subdive states differ primarily due to $\mu_{W}^{*}$ and $\sigma_{W}^{*}$
since both are much higher for subdive state 2 than for subdive state 1. These
larger parameter values correspond to much more vigorous and variable periodic
behaviour in the acceleration data.

Table 2: Average decoding accuracies and training times for all models used to categorize dive type and subdive state in the simulation study. Each of the four models was fit to 500 training data sets comprised of 100 simulated dives and tested on test data sets also comprised of 100 simulated dives. Reported values are averages, and $\pm$ refers to the sample standard deviation across the 500 data sets. Rows labelled Both/Both correspond to overall average decoding accuracy. Model | Train Time (min) | Dive Type | Subdive State | Dive Accuracy | Subdive Accuracy
---|---|---|---|---|---
CarHMM-DFT | $70\pm 11$ | Both | Both | ————- | $0.93\pm 0.01$
1 | 1 | ————- | $0.79\pm 0.04$
1 | 2 | $0.96\pm 0.01$
2 | 1 | ————- | $0.92\pm 0.01$
2 | 2 | $0.94\pm 0.01$
HHMM-DFT | $209\pm 51$ | Both | Both | $0.94\pm 0.04$ | $0.88\pm 0.04$
1 | 1 | $0.97\pm 0.03$ | $0.63\pm 0.12$
1 | 2 | $0.96\pm 0.01$
2 | 1 | $0.85\pm 0.10$ | $0.79\pm 0.16$
2 | 2 | $0.92\pm 0.03$
CarHHMM | $236\pm 52$ | Both | Both | $0.87\pm 0.17$ | $0.74\pm 0.05$
1 | 1 | $0.87\pm 0.21$ | $0.77\pm 0.21$
1 | 2 | $0.71\pm 0.09$
2 | 1 | $0.82\pm 0.14$ | $0.87\pm 0.23$
2 | 2 | $0.66\pm 0.10$
CarHHMM-DFT | $132\pm 40$ | Both | Both | $0.94\pm 0.04$ | $0.93\pm 0.01$
1 | 1 | $0.96\pm 0.03$ | $0.76\pm 0.04$
1 | 2 | $0.96\pm 0.01$
2 | 1 | $0.87\pm 0.09$ | $0.93\pm 0.01$
2 | 2 | $0.93\pm 0.01$

We calculate maximum likelihood estimates
$\\{\hat{\theta},\hat{\Gamma},\hat{\theta}^{*},\hat{\Gamma}^{*}\\}$ for all
four candidate models for each of the 500 data sets using the Cedar Compute
Canada cluster with 1 CPU and 4 GB of dedicated memory per data set. For each
of the 500 training data sets, we simulate a test data set to assess how well
each model predicts the hidden states, as follows. Each test data set consists
of a sequence of 100 dives and is created from the generative model with true
parameters $\\{\theta,\Gamma,\theta^{*},\Gamma^{*}\\}$. To assess the coarse-
scale hidden state prediction, we estimate
$p_{t}(i|y,\tilde{y}^{*})\equiv\Pr(X_{t}=i|Y=y,\tilde{Y}^{*}=\tilde{y}^{*})$,
$i=1,2$, $t=1,\ldots,100$ using the test-set observations $(y,\tilde{y}^{*})$
and training-set maximum likelihood estimates. These estimates are found using
the forward-backward algorithm (Zucchini et al.,, 2016). We compare these
estimated conditional probabilities to $\\{x_{1},\ldots,x_{100}\\}$, the true
coarse-scale state realizations in the test data, by calculating the average
dive decoding accuracy for a single training/test data set pair,
$\sum_{t=1}^{100}\hat{p}_{t}(x_{t}|y,\tilde{y}^{*})/100$. We then report the
average of these over the 500 training/test data set pairs. Analogously, to
assess prediction of the fine-scale states, we estimate
$p^{*}_{t,\tilde{t}^{*}}(i^{*}|y,\tilde{y}^{*})\equiv\Pr(\tilde{X}^{*}_{t,\tilde{t}^{*}}=i|Y=y,\tilde{Y}^{*}=\tilde{y}^{*})$,
$i^{*}=1,2$, $\tilde{t}^{*}=1,\ldots,\tilde{T}^{*}_{t}$, $t=1,\ldots,100$,
using the test-set observations, the training-set maximum likelihood
estimates, and the forward-backward algorithm. Denoting the true fine-scale
state realizations from the test data set and dive $t$ as
$\\{\tilde{x}^{*}_{t,1},\ldots,\tilde{x}^{*}_{t,\tilde{T}^{*}_{t}}\\}$, we
define the overall average subdive decoding accuracy as the average value of
$\hat{p}^{*}_{t,\tilde{t}^{*}}(\tilde{x}^{*}_{t,\tilde{t}^{*}}|y,\tilde{y}^{*})$
across all simulated test data sets, dives, and windows. The conditional
probabilities are estimated according to one of the four models under study,
using the maximum likelihood estimates from the training data set in
conjunction with the forward-backward algorithm.

###  4.2 Simulation results

Figure 7: Estimated probabilities that each dive is of type 2 for five
selected dives of a simulated data set of killer whale dive behaviour. Each
panel is partitioned into dives by vertical black lines. The colour of the
curve corresponds to the true dive type while the colour of the background
corresponds to the true subdive state. The CarHMM-DFT is omitted because it
assumes that there is only one dive type. Figure 8: Estimated probabilities
that each window corresponds to subdive state 2 for five selected dives of a
simulated data set of killer whale dive behaviour. Each panel is partitioned
into dives by vertical black lines. The colour of the curve corresponds to the
true subdive state while the colour of the background corresponds to the true
dive type.

The full CarHHMM-DFT is the best performing model of the four candidates since
it is the generating model. Its average dive decoding accuracy and average
subdive decoding accuracy are both greater than 0.9. All parameter estimates
of fine-scale mean values and probability transition matrices,
$\hat{\mu}^{*}$, $\hat{\Gamma}^{*}$ and $\hat{\Gamma}$, respectively, have
statistically insignificant biases. In addition, the biases of $\hat{\sigma}$,
$\hat{\phi}$, and $\hat{\mu}$ are either as small or smaller than the other
three candidate models. The empirical standard errors of all parameter
estimates ($\hat{\theta}$, $\hat{\Gamma}$, $\hat{\theta}^{*}$,
$\hat{\Gamma}^{*}$) are well-approximated by the inverse of the observed
Fisher information matrix, although the estimated standard errors tend to be
slightly smaller than the empirical standard errors. This underestimation is
especially noticeable for parameters associated with the wiggliness
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$, where the empirical standard error can be
up to twice as large as the estimated standard error. See Tables 5 through 9
of the supplementary material for detailed results.

The HHMM-DFT performs similarly to the CarHHMM-DFT in most respects. Its
average dive decoding accuracy is comparable to the CarHHMM-DFT while its
average subdive decoding accuracy is worse by approximately 5 percentage
points. Its parameter estimates are comparable to the CarHHMM-DFT with the
notable exception that it greatly overestimates $\sigma_{A}^{*(\cdot,1)}$ and
$\sigma_{A}^{*(\cdot,2)}$. In addition, the estimated standard errors of
$\hat{\mu}_{A}^{*(\cdot,1)}$, $\hat{\mu}_{A}^{*(\cdot,2)}$,
$\hat{\sigma}_{A}^{*(\cdot,1)}$, and $\hat{\sigma}_{A}^{*(\cdot,2)}$ are much
smaller than the associated empirical standard errors (see Table 7 of the
supplementary material). These results suggest that the estimates of standard
deviation can be too large and estimates of standard errors can be too small
when auto-correlation is ignored. This finding is consistent with the results
of the case study, where the HHMM-DFT produced larger estimates of
$\sigma_{A}^{*}$ and smaller estimates of standard error compared to the
CarHHMM-DFT. When standard errors are underestimated, the associated
confidence intervals are too narrow, implying that researchers may be
overconfident in their parameter estimates.

The CarHHMM is the worst-performing model in terms of accuracy, as its average
dive decoding accuracy is below $0.9$ and its average subdive decoding
accuracy is below $0.8$. This result is consistent with expectations because
the CarHHMM does not model the “wiggliness” of the fine-scale process, which
is the most distinct difference between the subdive states. In addition to its
poor average decoding accuracy, the CarHHMM is also the worst of the four
candidate models at estimating parameters. Parameter estimates associated with
subdive type 2 ($\theta^{*(\cdot,2)}$) are especially poor. See Section 2 of
the supplementary material for more detailed results.

Finally, the CarHMM-DFT is nearly identical to the CarHHMM-DFT in terms of
average subdive decoding accuracy, fine-scale parameter biases, and both
estimated and empirical standard error for the fine-scale parameter estimates.
In addition, the time required to fit the CarHMM-DFT is less than half of that
of the other models (see Table 2). However, this model cannot estimate dive
type as it lacks any hierarchical structure. The CarHMM-DFT nonetheless fits a
(misspecified) single Gamma distribution over the dive duration of all dives.
The resulting parameter estimates ($\hat{\mu}$ and $\hat{\sigma}$) are highly
correlated (See Figure 25 of the supplementary material).

Figures 7 and 8 display five dives of one simulated data set as well as the
decoded dive types and subdive states associated with each model. The CarHHMM-
DFT and CarHMM-DFT produce similar estimates of subdive state while the HHMM-
DFT is slightly more likely to predict that a given window corresponds to
subdive state 2. As expected, the CarHHMM is the least accurate model when
predicting subdive state. The CarHHMM-DFT and HHMM-DFT yield similar estimates
of dive type while the CarHHMM is less accurate and even misclassifies the
dive type of the second dive. The CarHMM-DFT does not estimate dive type.

##  5 DISCUSSION

Current Functional Data Analysis literature addresses dependence between
curves either with multilevel models (Chen and Müller,, 2012; Di et al.,,
2009), which lack a time component, or with functional time series, which
overlook the possibility that curves have several distinct “types” (Kokoszka
and Reimherr,, 2018). Our work addresses these issues and introduces a
flexible framework to model functional time-series data using HMMs. We suggest
handling temporal dependence between curves by using either an HMM or a CarHMM
to model the curve sequence. We then suggest viewing each individual curve as
an HMM emission whose distribution is described by a fine-scale model. Here we
use a CarHMM as the fine-scale model, but there are a wide range of possible
fine-scale models, including a Poisson process or continuous time approach
similar to that of Michelot and Blackwell, (2019). We also incorporate a
moving-window transformation at the fine scale to capture intricate dependence
structures on short time scales. Together, the coarse- and fine-scale models
make up a hierarchical structure which can account for simultaneous processes
taking place at different time scales. Provided the construction is not overly
complex, a hierarchical model created using our method can be both flexible
and easy to fit using maximum likelihood estimation.

We demonstrate the usefulness of this framework using a
biomechanical/ecological example, where we use HMMs to classify the coarse-
and fine-scale diving behaviour of a northern resident killer whale in Queen
Charlotte Sound off the coast of British Columbia, Canada. Our analysis gives
a deeper understanding of a killer whale’s tri-axial movement and thus its
behaviour and energy expenditure (Gleiss et al.,, 2011; Qasem et al.,, 2012),
both of which are important for understanding the foraging ecology and
nutritional status of northern resident killer whales (Noren,, 2011). Our
model is also applicable to many diving animals such as sharks (Adam et al.,,
2019), seals (Jeanniard du Dot et al.,, 2016), and porpoises (Leos-Barajas et
al.,, 2017). In addition, since complicated state-switching processes with
temporal dependence are common in settings ranging from speech recognition
(Juang and Rabiner,, 1991) and neuroscience (Langrock et al.,, 2013) to
oceanography (Bulla et al.,, 2012) and ecology (Adam et al.,, 2019), we
believe that researchers can adapt our methodology for the analysis of a wide
range of time series data in a variety of fields.

##  6 ACKNOWLEDGEMENTS

The killer whale data was collected under University of British Columbia
Animal Care Permit no. A19-0053 and Fisheries and Oceans Canada Marine Mammal
Scientific License for Whale Research no. XMMS 6 2019. This research was
enabled in part by support provided by WestGrid (www.westgrid.ca) and Compute
Canada (www.computecanada.ca). We acknowledge the support of the Natural
Sciences and Engineering Research Council of Canada (NSERC) as well as the
support of Fisheries and Oceans Canada (DFO). Marie Auger-Méthé and Nancy
Heckman thank the NSERC Discovery program, and Marie Auger-Méthé additionally
thanks the Canadian Research Chair program. Evan Sidrow thanks the University
of British Columbia and the Four-Year Doctoral Fellowship program. We are
grateful to Dr. Joe Watson for his constructive suggestions.

## References

  * Adam et al., (2019) Adam, T., Griffiths, C., Leos Barajas, V., Meese, E., Lowe, C., Blackwell, P., Righton, D., and Langrock, R. (2019).  Joint modelling of multi-scale animal movement data using hierarchical hidden Markov models.  Methods in Ecology and Evolution, 10(9):1536–1550.
  * Alex Shorter et al., (2017) Alex Shorter, K., Shao, Y., Ojeda, L., Barton, K., Rocho-Levine, J., van der Hoop, J., and Moore, M. (2017).  A day in the life of a dolphin: Using bio-logging tags for improved animal health and well-being.  Marine Mammal Science, 33(3):785–802.
  * Bebbington, (2007) Bebbington, M. S. (2007).  Identifying volcanic regimes using hidden Markov models.  Geophysical Journal International, 171(2):921–942.
  * Borchers et al., (2013) Borchers, D. L., Zucchini, W., Heide-Jørgensen, M. P., Cañadas, A., and Langrock, R. (2013).  Using hidden Markov models to deal with availability bias on line transect surveys.  Biometrics, 69(3):703–713.
  * Brumback and Rice, (1998) Brumback, B. A. and Rice, J. A. (1998).  Smoothing spline models for the analysis of nested and crossed samples of curves.  Journal of the American Statistical Association, 93(443):961–976.
  * Bulla et al., (2012) Bulla, J., Lagona, F., Maruotti, A., and Picone, M. (2012).  A multivariate hidden Markov model for the identification of sea regimes from incomplete skewed and circular time series.  Journal of Agricultural, Biological and Environmental Statistics, 17:544–567.
  * Börger et al., (2020) Börger, L., Bijleveld, A., Fayet, A., Machovsky-Capuska, G., Patrick, S., Street, G., and Wal, E. (2020).  Biologging special feature.  Journal of Animal Ecology, 89(1):6–15.
  * Cade et al., (2018) Cade, D. E., Barr, K. R., Calambokidis, J., Friedlaender, A. S., and Goldbogen, J. A. (2018).  Determining forward speed from accelerometer jiggle in aquatic environments.  Journal of Experimental Biology, 221(2):jeb170449.
  * Chen and Müller, (2012) Chen, K. and Müller, H.-G. (2012).  Modeling repeated functional observations.  Journal of the American Statistical Association, 107(500):1599–1609.
  * Crainiceanu et al., (2009) Crainiceanu, C. M., Staicu, A.-M., and Di, C.-Z. (2009).  Generalized multilevel functional regression.  Journal of the American Statistical Association, 104(488):1550–1561.
  * DeRuiter et al., (2017) DeRuiter, S. L., Langrock, R., Skirbutas, T., Goldbogen, J. A., Calambokidis, J., Friedlaender, A. S., and Southall, B. L. (2017).  A multivariate mixed hidden Markov model for blue whale behaviour and responses to sound exposure.  The Annals of Applied Statistics, 11(1):362–392.
  * Di et al., (2009) Di, C.-Z., Crainiceanu, C. M., Caffo, B. S., and Punjabi, N. M. (2009).  Multilevel functional principal component analysis.  The Annals of Applied Statistics, 3(1):458–488.
  * Fehlmann et al., (2017) Fehlmann, G., O’Riain, J., Hopkins, P. W., O’Sullivan, J., Holton, M. D., Shepard, E. L. C., and King, A. J. (2017).  Identification of behaviours from accelerometer data in a wild social primate.  Animal Biotelemetry, 5:6.
  * Fischer and Meier-Hellstern, (1993) Fischer, W. and Meier-Hellstern, K. S. (1993).  The Markov-modulated Poisson process (MMPP) cookbook.  Perform. Evaluation, 18(2):149–171.
  * Ford and Ellis, (2006) Ford, J. K. B. and Ellis, G. M. (2006).  Selective foraging by fish-eating killer whales _Orcinus orca_ in British Columbia.  Marine Ecology-Progress Series, 316:185–199.
  * Ford et al., (2009) Ford, J. K. B., Ellis, G. M., Olesiuk, P. F., and Balcomb, K. C. (2009).  Linking killer whale survival and prey abundance: Food limitation in the oceans’ apex predator?  Biology Letters, 6(1):139–42.
  * Fu and Heckman, (2019) Fu, E. and Heckman, N. (2019).  Model-based curve registration via stochastic approximation EM algorithm.  Computational Statistics and Data Analysis, 131:159–175.
  * Getman et al., (2009) Getman, A., Cooper, C. D., Key, G., Zhou, H., and Frankle, N. (2009).  Detection of mobile machine damage using accelerometer data and prognostic health monitoring techniques.  In 2009 IEEE Workshop on Computational Intelligence in Vehicles and Vehicular Systems, pages 101–104.
  * Gleiss et al., (2011) Gleiss, A. C., Wilson, R. P., and Shepard, E. L. C. (2011).  Making overall dynamic body acceleration work: On the theory of acceleration as a proxy for energy expenditure.  Methods in Ecology and Evolution, 2(1):23–33.
  * Green et al., (2009) Green, J. A., Halsey, L. G., Wilson, R. P., and Frappell, P. B. (2009).  Estimating energy expenditure of animals using the accelerometry technique: Activity, inactivity and comparison with the heart-rate technique.  Journal of Experimental Biology, 212(5):745–746.
  * Hastie et al., (2006) Hastie, G. D., Rosen, D. A. S., and Trites, A. W. (2006).  The influence of depth on a breath-hold diver: Predicting the diving metabolism of Steller sea lions (_Eumetopias jubatus_).  Journal of Experimental Marine Biology and Ecology, 3:163–170.
  * Heerah et al., (2017) Heerah, K., Woillez, M., Fablet, R., Garren, F., Martin, S., and De Pontual, H. (2017).  Coupling spectral analysis and hidden Markov models for the segmentation of behavioural patterns.  Movement Ecology, 5:20.
  * Hooten et al., (2017) Hooten, M. B., King, R., and Langrock, R. (2017).  Guest editor’s introduction to the special issue on “animal movement modeling”.  Journal of Agricultural, Biological and Environmental Statistics, 22(3):224–231.
  * Isojunno et al., (2017) Isojunno, S., Sadykova, D., DeRuiter, S., Curé, C., Visser, F., Thomas, L., Miller, P. J. O., and Harris, C. M. (2017).  Individual, ecological, and anthropogenic influences on activity budgets of long‐finned pilot whales.  Ecosphere, 8(12):e02044.
  * Jeanniard-du Dot et al., (2016) Jeanniard-du Dot, T., Guinet, C., Arnould, J. P., Speakman, J. R., and Trites, A. W. (2016).  Accelerometers can measure total and activity-specific energy expenditures in free-ranging marine mammals only if linked to time-activity budgets.  Functional Ecology, 31(2):377–386.
  * Jeanniard du Dot et al., (2016) Jeanniard du Dot, T., Trites, A. W., Arnould, J. P. Y., Speakman, J. R., and Guinet, C. (2016).  Activity-specific metabolic rates for diving, transiting, and resting at sea can be estimated from time-activity budgets in free-ranging marine mammals.  Ecology and Evolution, 7(9):2969–2976.
  * Juang and Rabiner, (1991) Juang, B. H. and Rabiner, L. R. (1991).  Hidden Markov models for speech recognition.  Technometrics, 33(3):251–272.
  * Kokoszka and Reimherr, (2018) Kokoszka, P. and Reimherr, M. (2018).  Introduction to Functional Data Analysis, chapter 8.  Chapman and Hall/CRC.
  * Langrock et al., (2018) Langrock, R., Adam, T., Leos-Barajas, V., Mews, S., Miller, D. L., and Papastamatiou, Y. P. (2018).  Spline‐based nonparametric inference in general state‐switching models.  Statistica Neerlandica, 72(3):179–200.
  * Langrock et al., (2013) Langrock, R., Swihart, B. J., Caffo, B. S., Punjabi, N. M., and Crainiceanu, C. M. (2013).  Combining hidden Markov models for comparing the dynamics of multiple sleep electroencephalograms.  Statistics in Medicine, 32(19):3342–3356.
  * Lawler et al., (2019) Lawler, E., Whoriskey, K., Aeberhard, W. H., Field, C., and Mills Flemming, J. (2019).  The conditionally autoregressive hidden Markov model (CarHMM): Inferring behavioural states from animal tracking data exhibiting conditional autocorrelation.  Journal of Agricultural, Biological and Environmental Statistics, 24(4):651–668.
  * Leos-Barajas et al., (2017) Leos-Barajas, V., Gangloff, E. J., Adam, T., Langrock, R., van Beest, F. M., Nabe-Nielsen, J., and Morales, J. M. (2017).  Multi-scale modeling of animal movement and general behavior data using hidden Markov models with hierarchical structures.  Journal of Agricultural, Biological and Environmental Statistics, 22(3):232–248.
  * Liu et al., (2015) Liu, Y.-Y., Li, S., Li, F., Song, L., and Rehg, J. M. (2015).  Efficient learning of continuous-time hidden Markov models for disease progression.  Advances in Neural Information Processing Systems, 28:3599–3607.
  * Lucero et al., (2019) Lucero, P., Sánchez, R., Macancela, J., Cabrera, D., Cerrada, M., Li, C., and Alonso, H. R. (2019).  Accelerometer placement comparison for crack detection in railway axles using vibration signals and machine learning.  In 2019 Prognostics and System Health Management Conference (PHM-Paris), pages 291–296.
  * McClintock et al., (2020) McClintock, B. T., Langrock, R., Gimenez, O., Cam, E., Borchers, D. L., Glennie, R., and Patterson, T. A. (2020).  Uncovering ecological state dynamics with hidden Markov models.  Ecology Letters, 23(12):1878–1903.
  * Michelot and Blackwell, (2019) Michelot, T. and Blackwell, P. G. (2019).  State-switching continuous-time correlated random walks.  Methods in Ecology and Evolution, 10(5):637–649.
  * Morris et al., (2006) Morris, J. S., Arroyo, C., Coull, B. A., Ryan, L. M., Herrick, R., and Gortmaker, S. L. (2006).  Using wavelet-based functional mixed models to characterize population heterogeneity in accelerometer profiles.  Journal of the American Statistical Association, 101(476):1352–1364.
  * Noren, (2011) Noren, D. P. (2011).  Estimated field metabolic rates and prey requirements of resident killer whales.  Marine Mammal Science, 27(1):60 – 77.
  * Nunes, (2019) Nunes, A. (2019).  Divebomb version 1.0.7.  Ocean Tracking Network, gitlab.oceantrack.org/anunes/divebomb.
  * Patterson et al., (2017) Patterson, T. A., Parton, A., Langrock, R., Blackwell, P. G., Thomas, L., and King, R. (2017).  Statistical modelling of individual animal movement: an overview of key methods and a discussion of practical challenges.  Advances in Statistical Analysis, 101:399–438.
  * Pohle et al., (2017) Pohle, J., Langrock, R., van Beest, M., and Schmidt, N. M. (2017).  Selecting the number of states in hidden Markov models: Pragmatic solutions illustrated using animal movement.  Journal of Agricultural, Biological and Environmental Statistics, 22:1–24.
  * Qasem et al., (2012) Qasem, L., Cardew, A., Wilson, A., Griffiths, I., Halsey, L. G., Shepard, E. L. C., Gleiss, A. C., and Wilson, R. (2012).  Tri-axial dynamic acceleration as a proxy for animal energy expenditure; should we be summing values or calculating the vector?  PloS one, 7(2):e31187.
  * Ramsay and Silverman, (2005) Ramsay, J. O. and Silverman, B. W. (2005).  Functional Data Analysis.  Springer.
  * Rice and Wu, (2001) Rice, J. A. and Wu, C. O. (2001).  Nonparametric mixed effects models for unequally sampled noisy curves.  Biometrics, 57(1):253–9.
  * Simon et al., (2012) Simon, M., Johnson, M., and Madsen, P. T. (2012).  Keeping momentum with a mouthful of water: Behavior and kinematics of humpback whale lunge feeding.  Journal of Experimental Biology, 215(21):3786–3798.
  * Sutherland, (1998) Sutherland, W. J. (1998).  The importance of behavioural studies in conservation biology.  Animal Behaviour, 56(4):801–809.
  * (47) Tennessen, J., Holt, M. M., Ward, E. J., Hanson, M. B., Emmons, C. K., Giles, D. A., and Hogan, J. T. (2019a).  Hidden Markov models reveal temporal patterns and sex differences in killer whale behavior.  Scientific Reports, 9:14951.
  * (48) Tennessen, J. B., Holt, M. M., Hanson, M. B., Emmons, C. K., Giles, D. A., and Hogan, J. T. (2019b).  Kinematic signatures of prey capture from archival tags reveal sex differences in killer whale foraging activity.  Journal of Experimental Biology, 222(3).
  * Whoriskey et al., (2016) Whoriskey, K., Auger-Méthé, M., Albertsen, C. M., Whoriskey, F. G., Binder, T. R., Krueger, C. C., and Mills Flemming, J. (2016).  A hidden Markov movement model for rapidly identifying behavioral states from animal tracks.  Ecology and Evolution, 7(7):2112–2121.
  * Williams and Noren, (2009) Williams, R. and Noren, D. P. (2009).  Swimming speed, respiration rate, and estimated cost of transport in adult killer whales.  Marine Mammal Science, 25(2):327–350.
  * Williams et al., (1999) Williams, T. M., Haun, J. E., and Friedl, W. A. (1999).  The diving physiology of bottlenose dolphins (_Tursiops truncatus_): I. balancing the demands of exercise for energy conservation at depth.  Journal of Experimental Biology, 202(20):2739–2748.
  * Wilson et al., (2019) Wilson, R. P., Börger, L., Holton, M. D., Scantlebury, D. M., Gómez-Laich, A., Quintana, F., Rosell, F., Graf, P. M., Williams, H., Gunner, R., Hopkins, L., Marks, N., Geraldi, N. R., Duarte, C. M., Scott, R., Strano, M. S., Robotka, H., Eizaguirre, C., Fahlman, A., and Shepard, E. L. C. (2019).  Estimates for energy expenditure in free-living animals using acceleration proxies: A reappraisal.  Journal of Animal Ecology, 89(1):161–172.
  * Wright et al., (2017) Wright, B. M., Ford, J. K. B., Ellis, G. M., Deecke, V. B., Shapiro, A. D., Battaile, B. C., and Trites, A. W. (2017).  Fine-scale foraging movements by fish-eating killer whales (_Orcinus orca_) relate to the vertical distributions and escape responses of salmonid prey (_Oncorhynchus_ spp.).  Movement Ecology, 5:3.
  * Xin et al., (2018) Xin, G., Hamzaoui, N., and Antoni, J. (2018).  Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals.  Measurement, 127:141–166.
  * Xu et al., (2020) Xu, Z., Laber, E. B., and Staicu, A.-M. (2020).  Hierarchical Continuous Time Hidden Markov Model, with Application in Zero-Inflated Accelerometer Data, pages 125–142.  Springer International Publishing, Cham.
  * Yao et al., (2005) Yao, F., Müller, H.-G., and Wang, J.-L. (2005).  Functional data analysis for sparse longitudinal data.  Journal of the American Statistical Association, 100(470):577–590.
  * Zucchini et al., (2016) Zucchini, W., Macdonald, I. L., and Langrock, R. (2016).  Hidden Markov Models for Time Series - An Introduction Using R.  CRC Press.

##  Appendix A Detailed description of data simulation from Section 4.1

We easily simulate realizations of the coarse-scale HMM ($X$ and $Y$) given
the parameters $\Gamma$ and $\theta$. For each dive $t$, we also easily
generate the fine-scale hidden Markov chain
$\tilde{X}^{*}_{t}\equiv\left\\{\tilde{X}^{*}_{t,1},\ldots,\tilde{X}^{*}_{t,\tilde{T}_{t}^{*}}\right\\}$
according to one of the probability transition matrices $\Gamma^{*(1)}$ or
$\Gamma^{*(2)}$, depending upon the value of $X_{t}$. This determines the
sequence of fine-scale hidden states corresponding to each window. Recall that
the fine-scale model is based on a sequence of $\tilde{T}_{t}^{*}$ two-second
windows, each containing 100 observations, and our model is formulated in
terms of quantities derived from the raw data within each window (namely, the
average acceleration, $\tilde{A}^{*}_{t,\tilde{t}^{*}}$, and wiggliness,
$\tilde{W}^{*}_{t,\tilde{t}^{*}}$). Generating the raw acceleration data from
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ and $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ is not
straightforward. In our simulation study, we generate raw acceleration data so
that we can visualize our results in terms of the underlying data curves for
each dive. Here, we explain how we generate the acceleration curves so that
$\tilde{A}^{*}_{t,\tilde{t}^{*}}$ and $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ both
follow the specified model. A key component is the discrete Fourier
transformation of the 100 raw acceleration values in window $\tilde{t}^{*}$ of
dive $t$:

$\hat{Y}^{*(k)}_{t,\tilde{t}^{*}}\equiv
DFT\left\\{Y^{*}_{t,100(\tilde{t}^{*}-1)+1},\ldots,Y^{*}_{t,100\tilde{t}^{*}}\right\\}(k)$

for $k\geq 0$, as defined in Equation (3).

We simulate the raw acceleration data for dive $t$ in three steps: (1)
simulate the average acceleration within each window
$\left(\hat{Y}^{*(0)}_{t,\tilde{t}^{*}}\right)$, (2) simulate all other
Fourier coefficients within each window
$\left(\hat{Y}^{*(k)}_{t,\tilde{t}^{*}},k=1,\ldots,99\right)$, and (3) take
the inverse discrete Fourier transform of $\hat{Y}^{*}_{t,\tilde{t}^{*}}$:

$\\{Y^{*}_{t,100(\tilde{t}^{*}-1)+1},\ldots,Y^{*}_{t,100\tilde{t}^{*}}\\}\equiv
IDFT\left\\{\hat{Y}^{*(0)}_{t,\tilde{t}^{*}},\ldots,\hat{Y}^{*(99)}_{t,\tilde{t}^{*}}\right\\}{\rm{~{}for~{}}}\tilde{t}^{*}=1,\ldots,\tilde{T}^{*}_{t}.$

The details of steps (1) and (2) are described below.

For step (1), we generate
$\hat{Y}^{*(0)}_{t,1},\ldots,\hat{Y}^{*(0)}_{t,\tilde{T}_{t}^{*}}$ as a CarHMM
with underlying Markov state sequence
$\tilde{X}^{*}_{t,1},\ldots,\tilde{X}^{*}_{t,\tilde{T}_{t}^{*}}$ and random
first emission. Specifically, we let

$\hat{Y}^{*(0)}_{t,1}|\tilde{X}^{*}_{t,1}=i^{*}~{}~{}\sim~{}~{}\mathcal{N}\left(0,\left(100\sigma_{A}^{*(\cdot,i^{*})}\right)^{2}\right)~{}~{}{\rm{and}}$
$\displaystyle\hat{Y}^{*(0)}_{t,\tilde{t}^{*}}|\tilde{X}^{*}_{t,\tilde{t}^{*}}=i^{*},\hat{Y}^{*(0)}_{t,\tilde{t}^{*}-1}$
$\displaystyle~{}~{}\sim~{}~{}\mathcal{N}\left(\phi_{A}^{*(\cdot,i^{*})}\hat{Y}^{*(0)}_{t,\tilde{t}^{*}-1},\left(100\sigma_{A}^{*(\cdot,i^{*})}\right)^{2}\right),$
(5) $\displaystyle\tilde{t}^{*}=2,\ldots,\tilde{T}^{*}_{t}$

where $\sigma_{A}^{*(\cdot,1)}=0.034s$, $\sigma_{A}^{*(\cdot,2)}=0.079s$.
$\phi_{A}^{*(\cdot,1)}=0.98$ and $\phi_{A}^{*(\cdot,2)}=0.87$.

For step (2), we first construct $\hat{Y}^{*(k)}_{t,\tilde{t}^{*}}$,
$k=1,\ldots,49$, as

$\hat{Y}^{*(k)}_{t,\tilde{t}^{*}}=a_{t,\tilde{t}^{*}}^{(k)}i\sqrt{b^{(k)}_{t,\tilde{t}^{*}}}$
(6)

where the $a^{(k)}_{t,\tilde{t}^{*}}$’s are independent and equal to either 1
or -1 each with probability 1/2 and $i$ is the imaginary number. We include
$i$ in Equation (6) to force all variation within a window to take the form of
a sine wave, which reduces the variation between the endpoints of windows
compared to a cosine wave. Given the fine scale states, the
$b^{(k)}_{t,\tilde{t}^{*}}$’s are independent and independent of the
$a^{(k)}_{t,\tilde{t}^{*}}$’s. The distribution of $b^{(k)}_{t,\tilde{t}^{*}}$
is

$\displaystyle\begin{split}b^{(k)}_{t,\tilde{t}^{*}}|\tilde{X}^{*}_{t,\tilde{t}^{*}}=1&~{}~{}\sim~{}~{}{\rm{Gamma}}(16.38/k^{3},36.23)\\\
b^{(k)}_{t,\tilde{t}^{*}}|\tilde{X}^{*}_{t,\tilde{t}^{*}}=2&~{}~{}\sim~{}~{}{\rm{Gamma}}(4.20/k^{3},1825.53).\\\
\end{split}$ (7)

The first argument of ${\rm{Gamma}}\left(\cdot,\cdot\right)$ is the shape
parameter and the second is the scale parameter. The squared magnitude of the
$k^{th}$ Fourier coefficient is equal to $b^{(k)}_{t,\tilde{t}^{*}}$, which
decays like $1/k^{3}$ to “smooth out” the raw acceleration data.

We then define the remaining 50 Fourier coefficients:

$\hat{Y}^{*(50)}_{t,\tilde{t}^{*}}=0{\rm{~{}~{}and~{}~{}}}\hat{Y}^{*(k)}_{t,\tilde{t}^{*}}=-\hat{Y}^{*(100-k)}_{t,\tilde{t}^{*}},k=51,\ldots,99.$

This guarantees that the inverse discrete Fourier transform is real-valued.

We now show that this construction of the raw acceleration data results in the
distributions listed in Section 4.1. It suffices to show that the construction
of the discrete Fourier transformations, the
$\hat{Y}^{(k)}_{t,\tilde{t}^{*}}$’s, yields the desired distributions.

First, since
$\hat{Y}^{(0)}_{t,\tilde{t}^{*}}=\sum_{n=1}^{100}Y^{*}_{t,100(\tilde{t}^{*}-1)+n}$
$=100\tilde{A}^{*}_{t,\tilde{t}^{*}}$, Equation (5) implies that
$\tilde{A}^{*}_{t,1},\ldots,\tilde{A}^{*}_{t,\tilde{T}^{*}}$ follows a CarHMM
with Normal emissions distributions and parameters
$\mu_{A}^{*(\cdot,1)}=\mu_{A}^{*(\cdot,2)}=\mu_{A}^{*(\cdot,i^{*})}=0$,
$\sigma_{A}^{*(\cdot,1)}=0.034s$, $\phi_{A}^{*(\cdot,1)}=0.98$,
$\sigma_{A}^{*(\cdot,2)}=0.079s$, and $\phi_{A}^{*(\cdot,2)}=0.87$.

From Equations (4) and (6), the wiggliness within window $\tilde{t}^{*}$ of
dive $t$ is

$\tilde{W}^{*}_{t,\tilde{t}^{*}}=\sum_{k=1}^{\tilde{\omega}}\big{|}\big{|}\hat{Y}^{(k)}_{t,\tilde{t}^{*}}\big{|}\big{|}^{2}=\sum_{k=1}^{\tilde{\omega}}b^{(k)}_{t,\tilde{t}^{*}}.$

If $\tilde{\omega}<50$, then $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ is the sum of
independent Gamma-distributed random variables with identical scale
parameters, so the distribution of $\tilde{W}^{*}_{t,\tilde{t}^{*}}$ is also
Gamma. Thus, by Equation (7)

$\tilde{W}^{*}_{t,\tilde{t}^{*}}|\tilde{X}^{*}_{t,\tilde{t}^{*}}=1~{}~{}\sim~{}~{}{\rm{Gamma}}\left(\sum_{k=1}^{\tilde{\omega}}16.38/k^{3},36.23\right)\text{
and }$
$\tilde{W}^{*}_{t,\tilde{t}^{*}}|\tilde{X}^{*}_{t,\tilde{t}^{*}}=2~{}~{}\sim~{}~{}{\rm{Gamma}}\left(\sum_{k=1}^{\tilde{\omega}}4.20/k^{3},1825.53\right).$

Setting $\tilde{\omega}$ to 10 and carrying out simple calculation of the mean
and variance of a Gamma distribution yields $\mu_{W}^{*(\cdot,1)}=23.3$,
$\sigma_{W}^{*(\cdot,1)}=13.0$, $\mu_{W}^{*(\cdot,2)}=301.2$, and
$\sigma_{W}^{*(\cdot,2)}=330.1$.

##  Appendix B Likelihood of CarHHMM-DFT

The overall likelihood of the CarHHMM-DFT model is as follows:

$\mathcal{L}_{\text{CarHHMM-
DFT}}(\theta,\theta^{*},\Gamma,\Gamma^{*};y,\tilde{y}^{*})=\delta
P(y_{1},\tilde{y}^{*}_{1};\theta,\theta^{*},\Gamma^{*})\prod_{t=2}^{T}\Gamma
P(y_{t},\tilde{y}^{*}_{t};\theta,\theta^{*},\Gamma^{*})\mathbf{1}_{N}.$

In particular,

$\displaystyle
P(y_{t},\tilde{y}^{*}_{t};\theta,\theta^{*},\Gamma^{*})=\text{diag}\Big{[}$
$\displaystyle
f^{(1)}(y_{t};\theta^{(1)})\mathcal{L}_{\text{fine}}\left(\theta^{*},\Gamma^{*(1)};\tilde{y}^{*}_{t}\right),\ldots,$
$\displaystyle
f^{(N)}(y_{t};\theta^{(N)})\mathcal{L}_{\text{fine}}\left(\theta^{*},\Gamma^{*(N)};\tilde{y}^{*}_{t}\right)\Big{]},$

where $f^{(i)}(y_{t};\theta^{(i)})$ is the emission distribution of dive
duration given that $X_{t}=i$. The likelihood $\mathcal{L}_{\text{fine}}$
corresponds to the fine-scale model and is equal to the following:

$\mathcal{L}_{\text{fine}}\left(\theta^{*},\Gamma^{*(i)};\tilde{y}^{*}_{t}\right)=\delta^{*(i)}\prod_{\tilde{t}^{*}=2}^{\tilde{T}^{*}_{t}}\Gamma^{*(i)}P(\tilde{y}^{*}_{t,\tilde{t}^{*}}|\tilde{y}^{*}_{t,\tilde{t}^{*}-1};\theta^{*})\mathbf{1}_{N^{*}},$

where
$P(\tilde{y}^{*}_{t,\tilde{t}^{*}}|\tilde{y}^{*}_{t,\tilde{t}^{*}-1};\theta^{*})$
is an $N^{*}\times N^{*}$ diagonal matrix with $(i^{*},i^{*})^{th}$ entry
equal to
$f^{*(\cdot,i^{*})}(\tilde{y}^{*}_{t,\tilde{t}^{*}}|\tilde{y}^{*}_{t,\tilde{t}^{*}-1};\theta^{*(\cdot,i^{*})})$.
Recall that
$f^{*(\cdot,i^{*})}(\cdot|\tilde{y}^{*}_{t,t^{*}-1};\theta^{*(\cdot,i^{*})})$
is the probability density function of $\tilde{Y}^{*}_{t,\tilde{t}^{*}}$ when
$\tilde{X}^{*}_{t,\tilde{t}^{*}}=i^{*}$ and
$\tilde{Y}^{*}_{t,\tilde{t}^{*}-1}=\tilde{y}^{*}_{t,\tilde{t}^{*}-1}$.

See pages 1- of CarHHMM_ArXiv_sup.pdf