# Continuous-time multivariate analysis

Biplab Paul, Philip T. Reiss
Department of Statistics, University of Haifa;
Erjia Cui
Division of Biostatistics and Health Data Science, University of Minnesota;
and
Noemi Foà
Department of Statistics, University of Haifa  The authors thank Hadar Fisher,
Ciprian Crainiceanu, Prince Osei, Laura Sangalli and Dror Arbiv for helpful
discussions, and the Editor, Associate Editor and referees for valuable
feedback. The work of Paul, Reiss and Foà was supported by Israel Science
Foundation under Grant 1076/19.

###### Abstract

The starting point for much of multivariate analysis (MVA) is an $n\times p$
data matrix whose $n$ rows represent observations and whose $p$ columns
represent variables. Some multivariate data sets, however, may be best
conceptualized not as $n$ discrete $p$-variate observations, but as $p$ curves
or functions defined on a common time interval. Here we introduce a framework
for extending techniques of multivariate analysis to such settings. The
proposed continuous-time multivariate analysis (CTMVA) framework rests on the
assumption that the curves can be represented as linear combinations of basis
functions such as $B$-splines, as in the Ramsay-Silverman representation of
functional data; but whereas functional data analysis extends MVA to the case
of observations that are curves rather than vectors — heuristically, $n\times
p$ data with $p$ infinite — we are instead concerned with what happens when
$n$ is infinite. We present continuous-time extensions of the classical MVA
methods of covariance and correlation estimation, principal component
analysis, Fisher’s linear discriminant analysis, and $k$-means clustering. We
show that CTMVA can improve on the performance of classical MVA, in particular
for correlation estimation and clustering, and can be applied in some settings
where classical MVA cannot, including variables observed at disparate time
points. CTMVA is illustrated with a novel perspective on a well-known Canadian
weather data set, and with applications to data sets involving international
development, brain signals, and air quality. The proposed methods are
implemented in the publicly available R package ctmva.

Keywords: $B$-splines; correlation matrix; Fisher’s linear discriminant
analysis; functional data; $k$-means clustering; principal component analysis

##  1 Introduction

The generic mathematical object that underpins much of multivariate analysis
(MVA) is an $n\times p$ data matrix $X$ whose $n$ rows represent observations
and whose $p$ columns represent variables. While classical MVA typically
assumes the $n$ observations to be independent, there is a long history of
applications of MVA to multivariate time series, notably in econometrics and
atmospheric science (Jolliffe, 2002). In this paper we argue that some time-
indexed multivariate data sets may be best conceptualized not as a time series
of $n$ discrete $p$-variate observations, but as $p$ curves or functions
defined on a common time interval; and we introduce a framework for
straightforwardly extending several techniques of MVA to this setting.

Heuristically, our proposed framework, henceforth _continuous-time
multivariate analysis_ (CTMVA), is the “transpose” of functional data analysis
(FDA; Ramsay and Silverman, 2005), in the following sense. In a programmatic
address that introduced the term “functional data” and foreshadowed the FDA
paradigm, Ramsay (1982) argued that some modern data sets call for replacing
the traditional $p$ variables with infinitely many points on a continuum, as
in the upper right picture in Figure 1, and suggested that spline methodology
could play a key role in handling such infinite-dimensional data. He briefly
noted that instead of taking $p=\infty$ as in FDA, one might consider
situations in which $n=\infty$, as in the lower left picture in Figure 1; but
he did not develop the $n=\infty$ case. The CTMVA framework is an attempt to
do so, while maintaining an emphasis on the role of spline and other basis-
function representations as in Ramsay (1982) and Ramsay and Silverman (2005).

Figure 1: Schematic diagram presented by Ramsay (1982) to explain “the concept
of a functional datum.” The diagram appeared as Figure 2 of that paper above
the caption: “Possible domains for statistical observations. These domains
depend on whether the number of replications or cases ($n$) is finite or
infinite and whether the number of variables or points of observation ($p$) is
finite or corresponds to the points on a continuum.” The upper right
represents functional data, whereas the lower left scenario, mentioned but not
developed by Ramsay (1982), is the continuous-time multivariate analysis setup
of the present paper. Figure reproduced with permission from Springer Nature.

To illustrate the contrast between CTMVA and FDA it may be instructive to
reconsider what is perhaps the most famous example of functional data: the
Canadian temperature data of Ramsay and Silverman (2005), consisting of mean
temperatures on each day of the year at 35 weather stations in Canada. In FDA,
these are viewed not as $n=35$ observations of dimension $p=365$, but as
smooth curves $x_{1}(t),\ldots,x_{35}(t)$ on the time interval
$\mathcal{I}=[0,365]$, so that in effect $p=\infty$ (uncountably infinite).
Figure 2 displays both the raw data (upper panel) and the curves following
roughness-penalty smoothing with respect to a 45-element Fourier basis (middle
panel), with the smoothing parameter chosen by restricted maximum likelihood
(Ruppert et al., 2003; Reiss and Ogden, 2009; Wood, 2011) . In CTMVA we again
view the data as smooth curves $x_{1}(t),\ldots,x_{35}(t)$, but these are now
understood as $n=\infty$ observations of $p=35$ variables.

In FDA, given a sample of curves arising from a stochastic process on a finite
interval $\mathcal{I}\subset\mathbb{R}$, much interest focuses on
$\Gamma(s,t)$, the covariance function of the underlying process. In
particular, in functional principal component analysis one estimates the
eigenfunctions of the associated covariance operator. In CTMVA, by contrast,
the covariances or correlations of interest are those among the
curves—representing, in this example, relationships among the various weather
stations. In §4 below we discuss the matrices shown at the bottom of Figure 2,
representing correlations among the 35 stations’ temperature curves, and in
Appendix A we consider correlations among precipitation curves from the same
data set.

Figure 2: Top: Mean daily temperatures at 35 Canadian weather stations (Ramsay
and Silverman, 2005). Middle: The same data after smoothing with respect to a
Fourier basis. Below: Continuous-time correlations among daily mean
temperatures in the 35 weather stations, visualized by the corrgram method
(Friendly, 2002) implemented in Wei and Simko (2021), before (left) and after
(right) removing the common seasonal trend.

The practical impetus for CTMVA is twofold. First, the strategy we present
below, of representing a multivariate time series as a set of smooth curves
rather than a discrete set of MV observations, can sometimes improve
performance of MVA techniques. We demonstrate this below for correlation
estimation, in the simulation study of §8, and for clustering, with the
environmetric data analyzed in §9. Second, CTMVA can sometimes be applied when
classical MVA cannot—for example, if the $p$ variables are not measured at a
uniform set of time points (as in §3.2 below).

After presenting the basic setup and notation in §2, in §3 we define
continuous-time versions of the sample mean vector and the sample covariance
and correlation matrices. As part of this development, §4 clarifies the basic
but possibly confusing role of centering, which can have an important effect
on correlation estimation. In §5–§7 we develop continuous-time (CT) versions
of three key MVA techniques, namely principal component analysis, Fisher’s
linear discriminant analysis and $k$-means clustering. A simulation study and
an application to air pollution data are presented in §8 and §9, respectively.
The paper concludes with a brief discussion in §10.

##  2 Setup

In what follows we assume
$\mbox{\boldmath$x$}(t)\equiv[x_{1}(t),\ldots,x_{p}(t)]^{T}$ to be a
realization of a stochastic vector process
$\mathcal{X}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{I}\rightarrow\mathbb{R}^{p}$. Depending on the application,
the process $\mathcal{X}(t)=[X_{1}(t),\ldots,X_{p}(t)]^{T}$ may or may not be
stationary. If it is, the mean vector
$\mbox{\boldmath$\mu$}(t)=E[\mathcal{X}(t)]\in\mathbb{R}^{p}$ and covariance
matrix
$\mbox{\boldmath$\Sigma$}(t)=E[\\{\mathcal{X}(t)-\mbox{\boldmath$\mu$}(t)\\}\\{\mathcal{X}(t)-\mbox{\boldmath$\mu$}(t)\\}^{T}]\in\mathbb{R}^{p\times
p}$ do not depend on $t$ and thus may be denoted by
$\mbox{\boldmath$\mu$},\mbox{\boldmath$\Sigma$}$. In some applications the $p$
functions exhibit a common trend, i.e. a time-varying smooth function
$m(t)=E[\frac{1}{p}\sum_{u=1}^{p}X_{u}(t)]$, an estimate of which should be
subtracted from each component of $\mbox{\boldmath$x$}(t)$ before proceeding
with data analysis. We explain this with reference to an example in §4.

Following Ramsay and Silverman (2005), we make the fundamental simplifying
assumption that, for $u=1,\ldots,p$, the $u$th function can be represented by
means of an expansion
$x_{u}(t)=\mbox{\boldmath$c$}_{u}^{T}\mbox{\boldmath$\phi$}(t)$ where
$\mbox{\boldmath$\phi$}(t)=[\phi_{1}(t),\ldots,\phi_{K}(t)]^{T}$ is a vector
of _a priori_ smooth basis functions on $\mathcal{I}$. Thus, if $C$ is the
$K\times p$ coefficient matrix with $u$th column $\mbox{\boldmath$c$}_{u}$,
then

$\mbox{\boldmath$x$}(t)=\mbox{\boldmath$C$}^{T}\mbox{\boldmath$\phi$}(t).$ (1)

Advantages of basis functions over other smoothing techniques include
flexibility, compact representation of even large MV data sets, and not
requiring equally spaced time points or a uniform set of times for the
component functions.

The observed data from which this basis expansion is derived may be a
traditional $n\times p$ data matrix, representing $p$ variables observed at a
common set of $n$ time points, with smoothing applied to each of the $p$
columns (see §3.3 below). But in some applications (e.g., in §3.2), the
variables may be observed at different time points. The only requirement here
is that the same basis is employed for all $p$ curves. For curves that are
cyclic it is natural for $\phi$ to be a Fourier basis: for the weather data,
for example, this ensures fits that transition smoothly from December 31 to
January 1. For other applications, $B$-spline bases, most often cubic
$B$-splines, are a very popular choice due to their computational efficiency
and favorable approximation properties. In practice, the representation is
typically obtained by penalized smoothing (see §3.3), and this renders the
results fairly insensitive to the choice of $K$ (e.g., Ruppert et al., 2003).

Our implementation of CTMVA in the R environment (R Core Team, 2023) depends
on the $B$-spline and Fourier basis functions from the fda package for
functional data analysis (Ramsay et al., 2009, 2023). In some applications it
might be desirable to employ data-driven basis functions such as functional
principal components (Silverman, 1996).

If $\mathcal{S}$ is a subinterval, or finite union of subintervals, of
$\mathcal{I}$, we define

$\bar{\mbox{\boldmath$\phi$}}_{\mathcal{S}}=|\mathcal{S}|^{-1}\int_{\mathcal{S}}\mbox{\boldmath$\phi$}(t)dt\in\mathbb{R}^{K}\quad\mbox{and}\quad\mbox{\boldmath$Q$}_{\mathcal{S}}=|\mathcal{S}|^{-1}\int_{\mathcal{S}}[\mbox{\boldmath$\phi$}(t)-\bar{\mbox{\boldmath$\phi$}}_{\mathcal{S}}][\mbox{\boldmath$\phi$}(t)-\bar{\mbox{\boldmath$\phi$}}_{\mathcal{S}}]^{T}dt\in\mathbb{R}^{K\times
K}.$ (2)

Let $\bar{\mbox{\boldmath$\phi$}}=\bar{\mbox{\boldmath$\phi$}}_{\mathcal{I}}$
and $\mbox{\boldmath$Q$}=\mbox{\boldmath$Q$}_{\mathcal{I}}$. These expressions
will be used repeatedly in what follows. Computation of
$\bar{\mbox{\boldmath$\phi$}}$ and $Q$ is discussed in Appendix B.

##  3 Continuous-time sample mean vector and covariance and correlation
matrices

###  3.1 Definitions

In what follows we motivate the continuous-time sample mean vector and sample
covariance matrix as analogues, or limits, of their classical counterparts.
These will be seen, moreover, to be natural estimates of corresponding
stochastic process parameters (see §3.3 below).

In classical multivariate analysis, given a data matrix
$\mbox{\boldmath$X$}=(x_{ij})_{1\leq i\leq n,1\leq j\leq p}$ representing $n$
observations of a $p$-dimensional random vector, the $p$-dimensional sample
mean vector and $p\times p$ covariance matrix are given, respectively, by
$\bar{\mbox{\boldmath$x$}}=(\bar{x}_{1},\ldots,\bar{x}_{p})^{T}$ with
$\bar{x}_{u}=n^{-1}\sum_{i=1}^{n}x_{iu}$ and
$\mbox{\boldmath$S$}=\left[\frac{1}{n}\sum_{i=1}^{n}(x_{iu}-\bar{x}_{u})(x_{iv}-\bar{x}_{v})\right]_{1\leq
u,v\leq p}$ [here we follow Mardia et al. (1979) in dividing by $n$ rather
than by $n-1$].

The CT (sample) mean vector $\bar{\mbox{\boldmath$x$}}^{*}$ and covariance
matrix $\mbox{\boldmath$S$}^{*}$ (here and below we use ∗ to signify CT
analogues of classical MVA estimators) are simply the limits, as
$n\rightarrow\infty$, of their classical counterparts
$\bar{\mbox{\boldmath$x$}},\mbox{\boldmath$S$}$, when the rows of $X$ are
$\mbox{\boldmath$x$}(t_{1})^{T},\ldots,\mbox{\boldmath$x$}(t_{n})^{T}$, where
$t_{1},\ldots,t_{n}$ form a grid of points spanning $\mathcal{I}$ at equal
intervals of $|\mathcal{I}|/n$ (e.g., $t_{i}=i/n$ for each $i$, for
$\mathcal{I}=[0,1]$). In that case
$|\mathcal{I}|\bar{\mbox{\boldmath$x$}}=\frac{|\mathcal{I}|}{n}\sum_{i=1}^{n}\mbox{\boldmath$x$}(t_{i})$
is a Riemann sum converging to $\int_{\mathcal{I}}\mbox{\boldmath$x$}(t)dt$ as
$n\rightarrow\infty$; dividing by $|\mathcal{I}|$ yields
$\bar{\mbox{\boldmath$x$}}\rightarrow\bar{\mbox{\boldmath$x$}}^{*}$, where

$\bar{\mbox{\boldmath$x$}}^{*}=(\bar{x}^{*}_{1},\ldots,\bar{x}^{*}_{p})^{T}\quad\mbox{with}\quad\bar{x}^{*}_{u}=|\mathcal{I}|^{-1}\int_{\mathcal{I}}x_{u}(t)dt.$
(3)

Consequently,
$|\mathcal{I}|\mbox{\boldmath$S$}=\frac{|\mathcal{I}|}{n}\sum_{i=1}^{n}[\mbox{\boldmath$x$}(t_{i})-\bar{\mbox{\boldmath$x$}}][\mbox{\boldmath$x$}(t_{i})-\bar{\mbox{\boldmath$x$}}]^{T}$
is asymptotically equivalent to
$\frac{|\mathcal{I}|}{n}\sum_{i=1}^{n}[\mbox{\boldmath$x$}(t_{i})-\bar{\mbox{\boldmath$x$}}^{*}][\mbox{\boldmath$x$}(t_{i})-\bar{\mbox{\boldmath$x$}}^{*}]^{T}$,
which is a Riemann sum converging to
$\int_{\mathcal{I}}[\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}][\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}]^{T}dt$.
Dividing by $|\mathcal{I}|$ leads to
$\mbox{\boldmath$S$}\rightarrow\mbox{\boldmath$S$}^{*}$, where

$\mbox{\boldmath$S$}^{*}=\left[|\mathcal{I}|^{-1}\int_{\mathcal{I}}[x_{u}(t)-\bar{x}^{*}_{u}]][x_{v}(t)-\bar{x}^{*}_{v}]dt\right]_{1\leq
u,v\leq
p}=|\mathcal{I}|^{-1}\int_{\mathcal{I}}[\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}][\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}]^{T}dt.$
(4)

The basis representation (1) makes the integrals in (3), (4) computable even
though $\mbox{\boldmath$x$}(t)$ is observed at only finitely many times $t$,
and (2) yields convenient and succinct expressions for the CT mean and
covariance matrix:
$\bar{\mbox{\boldmath$x$}}^{*}=\mbox{\boldmath$C$}^{T}\bar{\mbox{\boldmath$\phi$}}$
and
$\mbox{\boldmath$S$}^{*}=\mbox{\boldmath$C$}^{T}\mbox{\boldmath$Q$}\mbox{\boldmath$C$}$.

The CT (sample) correlation between the $u$th and $v$th functions is similarly
defined as

$r^{*}_{uv}=\frac{\int_{\mathcal{I}}x^{c}_{u}(t)x^{c}_{v}(t)dt}{\left[\int_{\mathcal{I}}x^{c}_{u}(t)^{2}dt\int_{\mathcal{I}}x^{c}_{v}(t)^{2}dt\right]^{1/2}},$
(5)

where $x^{c}_{u}(t)=x_{u}(t)-\bar{x}^{*}_{u}$. As in classical MVA, the matrix
of CT correlations equals
$\mbox{\boldmath$R$}^{*}=\mbox{\boldmath$D$}^{*-1/2}\mbox{\boldmath$S$}^{*}\mbox{\boldmath$D$}^{*-1/2}$,
where $\mbox{\boldmath$D$}^{*}$ is the diagonal matrix with the same main
diagonal as $\mbox{\boldmath$S$}^{*}$.

###  3.2 An example: World Development Indicators

To illustrate the use of CT correlation for time series with irregularly
spaced times and substantial missingness, we consider several of the World
Development Indicators (WDI; https://datacatalog.worldbank.org), variables
that are estimated annually for each nation in the world. Quantitative
political scientists use such national time series to investigate
relationships among different indicators. Missing values are common, and may
be addressed with ad hoc fixes (e.g., Baum and Lake, 2003) or complex
imputation schemes (e.g., Honaker and King, 2010). CT correlation offers a
simple and fast way to assess relationships between a potentially large number
of indicators, while allowing for substantial missing data. Such applications
of CT correlation may prove highly useful for exploratory data analysis and
hypothesis generation, which are increasingly important goals in light of the
popularity of online resources such as Gapminder and Our World in Data, and
the large number of available candidate variables.

Figure 3: Left: Domestic fixed investment (above) and female secondary school
enrollment (below) for DR Congo, with spline smooths $\pm$ 2 standard errors.
CT correlation is computed for 1991–2021, the overlap of the two variables’
time ranges, resulting in an estimate of 0.977. Right: Histograms of country-
level CT correlations for each pair of variables among the five WDIs
considered.

We focus here on five variables taken from Baum and Lake (2003): fertility,
gross domestic product (GDP) per capita, gross domestic fixed investment,
female life expectancy, and female secondary school enrollment. At left in
Figure 3 are plots of investment and female enrollment in the Democratic
Republic of the Congo, which illustrate the difficulty of assessing
correlation between highly incomplete time series. In only 12 years are both
variables estimated, so that ordinary correlation can be computed with only 12
observations. But the CT correlation is obtained by using all observations for
1991–2021, the overlap of the two variables’ time ranges, to estimate the
trend in each variable.

The resulting correlation estimate of 0.977 is not very precise, with a 95%
confidence interval from 0.500 to 0.997 based on 999 bootstrap replicates. (In
§10 we mention a weighting approach that might help to improve the precision
of CT correlation.) But while individual correlation estimates from this data
set lack precision, it is possible to discern interesting patterns from a
collection of such estimates. We computed CT correlation for each pair of
variables, for each of the 216 countries in the data set, subject to the
requirements that (i) each variable have at least 8 observations and (ii) the
overlap of the two time ranges be at least 12 years. These requirements were
met for 85% of the $216{5\choose 2}$ possible correlations, resulting in 1840
CT correlations, which were computed in 27 seconds on a MacBook Pro with 32 GB
of RAM. At right in Figure 3 are histograms of the country-specific CT
correlations for each pair of variables. Notwithstanding the above results for
DR Congo, there is no clear international pattern of positive or negative
correlation between investment and any of the other variables. On the other
hand, the four histograms in the top right corner suggest that in most
countries, historical trends in two measures of women’s development, female
life expectancy and female school enrollment, are negatively correlated with
fertility trends and positively correlated with GDP trends.

###  3.3 Estimators or estimands?

In classical MVA, the sample mean $\bar{\mbox{\boldmath$x$}}$ is an estimate,
based on observations $\mbox{\boldmath$x$}_{1},\ldots,\mbox{\boldmath$x$}_{n}$
from a common multivariate distribution, of the population mean
$\mbox{\boldmath$\mu$}=E(\mathcal{X})$ of a random vector $\mathcal{X}$
arising from that distribution, while the sample covariance matrix $S$ is an
estimate of the population covariance matrix
$\mbox{\boldmath$\Sigma$}=E[(\mathcal{X}-\mbox{\boldmath$\mu$})(\mathcal{X}-\mbox{\boldmath$\mu$})^{T}]$.
Likewise, in CTMVA, $\bar{\mbox{\boldmath$x$}}^{*},\mbox{\boldmath$S$}^{*}$
are natural estimators of $\mbox{\boldmath$\mu$},\mbox{\boldmath$\Sigma$}$ as
defined in §2, assuming stationarity. If the process is not stationary,
$\bar{\mbox{\boldmath$x$}}^{*},\mbox{\boldmath$S$}^{*}$ can be viewed as
estimators of $|\mathcal{I}|^{-1}\int_{\mathcal{I}}\mbox{\boldmath$\mu$}(t)dt$
and $|\mathcal{I}|^{-1}\int_{\mathcal{I}}\mbox{\boldmath$\Sigma$}(t)dt$,
respectively.

But in many applications, interest centers not on the underlying process but
on the specific instance $\mbox{\boldmath$x$}(\cdot)$. It is then more
appropriate to regard $\bar{\mbox{\boldmath$x$}}^{*},\mbox{\boldmath$S$}^{*}$,
as well as $\mbox{\boldmath$R$}^{*}$, not as _estimators_ but as _estimands_ ,
since in practice the function $\mbox{\boldmath$x$}(\cdot)$ on which they are
based is not observed perfectly, for three reasons: (i) the $p$ component
functions are observed only at finitely many time points, (ii) they are
generally observed with noise and/or measurement error, and (iii) the basis
representation (1) introduces approximation error.

Formally, for $u=1,\ldots,p$, at finitely many times $t_{iu}$ we observe not
the true $x_{u}(t_{iu})$ but

$z_{iu}=x_{u}(t_{iu})+\varepsilon_{iu},$ (6)

where $\varepsilon_{iu}$ represents zero-mean error. We apply penalized
smoothing (Ruppert et al., 2003; Wood, 2017) to the $z_{iu}$’s for each $u$,
to obtain estimates
$\hat{\mbox{\boldmath$x$}}(t)=\widehat{\mbox{\boldmath$C$}}^{T}\mbox{\boldmath$\phi$}(t)\in\mathbb{R}^{p}$
for the $p$ functions. We note that for given $u$, the $\varepsilon_{iu}$’s
may not be independent. The penalized spline literature includes methods to
take residual dependence into account, as well as some evidence that
likelihood-based smoothness selection is robust to such dependence
(Krivobokova and Kauermann, 2007). In this setup,
$\bar{\mbox{\boldmath$x$}}^{*},\mbox{\boldmath$S$}^{*},\mbox{\boldmath$R$}^{*}$
are unobservable but are estimated, respectively, by

$\hat{\bar{\mbox{\boldmath$x$}}}^{*}=|\mathcal{I}|^{-1}\int_{\mathcal{I}}\hat{\mbox{\boldmath$x$}}(t)dt=\widehat{\mbox{\boldmath$C$}}^{T}\bar{\mbox{\boldmath$\phi$}},\quad\quad\widehat{\mbox{\boldmath$S$}}^{*}=|\mathcal{I}|^{-1}\int_{\mathcal{I}}[\hat{\mbox{\boldmath$x$}}(t)-\hat{\bar{\mbox{\boldmath$x$}}}^{*}][\hat{\mbox{\boldmath$x$}}(t)-\hat{\bar{\mbox{\boldmath$x$}}}^{*}]^{T}dt=\widehat{\mbox{\boldmath$C$}}^{T}\mbox{\boldmath$Q$}\widehat{\mbox{\boldmath$C$}},$

and
$\widehat{\mbox{\boldmath$R$}}^{*}=\widehat{\mbox{\boldmath$D$}}^{*-1/2}\widehat{\mbox{\boldmath$S$}}^{*}\widehat{\mbox{\boldmath$D$}}^{*-1/2}$,
where $\widehat{\mbox{\boldmath$D$}}^{*}$ is the diagonal matrix with the same
main diagonal as $\widehat{\mbox{\boldmath$S$}}^{*}$. For simplicity, in what
follows we omit the “hats” and treat $x$ and consequently
$\bar{\mbox{\boldmath$x$}}^{*},\mbox{\boldmath$S$}^{*},\mbox{\boldmath$R$}^{*}$
as given, except when presenting the simulation study of §8, which shows the
advantage of CT over ordinary correlation for stochastic processes observed
with noise.

To restate the above argument more succinctly and concretely, let us consider
a bivariate process that is stationary, so that the covariance and hence the
correlation $\rho$ between the two components are time-independent. Classical
MVA distinguishes the population parameter $\rho$ from the sample quantity
$r$. But in CTMVA there are not two but three quantities: the process
correlation $\rho$ gives rise to the CT correlation $r^{*}$ between a pair of
functions drawn from the process, which is estimated from data by
$\hat{r}^{*}$—or schematically,

$\rho\xrightarrow{(i)}r^{*}\xrightarrow{(ii)}\hat{r}^{*}$ (7)

(where the arrows informally denote “gives rise to” rather than denoting
convergence). Depending on the context we may be more interested in step (i)
or in step (ii) of (7).

##  4 Two forms of centering

The mean-centered variable vector
$\mbox{\boldmath$x$}^{c}(t)=\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}$
appeared in §3.1 and has a number of other uses, such as for CT principal
component analysis (see §5 below). By (1),
$\mbox{\boldmath$x$}^{c}(t)=\mbox{\boldmath$C$}^{T}[\mbox{\boldmath$\phi$}(t)-\bar{\mbox{\boldmath$\phi$}}]$,
i.e. the centered variables are linear combinations of
$\phi_{1}(t)-\bar{\phi}_{1},\ldots,\phi_{K}(t)-\bar{\phi}_{K}$, but it would
be more convenient if we could represent them as linear combinations of the
original (not centered) basis functions $\phi_{1}(t),\ldots,\phi_{K}(t)$. It
turns out that we can do so, under the simple condition that

$\mbox{there exists }\mbox{\boldmath$w$}\in\mathbb{R}^{K}\mbox{ such that
}\mbox{\boldmath$w$}^{T}\mbox{\boldmath$\phi$}(t)=1\mbox{ for all }t,$ (8)

since then
$\mbox{\boldmath$x$}^{c}(t)=\mbox{\boldmath$C$}^{T}[\mbox{\boldmath$\phi$}(t)-\bar{\mbox{\boldmath$\phi$}}]=\mbox{\boldmath$C$}^{T}[\mbox{\boldmath$I$}_{K}-\bar{\mbox{\boldmath$\phi$}}\mbox{\boldmath$w$}^{T}]\mbox{\boldmath$\phi$}(t)$,
i.e., the centered curves have a basis representation as in (1), but with
coefficient matrix
$[\mbox{\boldmath$I$}_{K}-\mbox{\boldmath$w$}\bar{\mbox{\boldmath$\phi$}}^{T}]\mbox{\boldmath$C$}$
in place of $C$. Condition (8) holds for $B$-spline bases, with
$\mbox{\boldmath$w$}=(1,\ldots,1)^{T}$, and for Fourier bases, with
$\mbox{\boldmath$w$}=(|\mathcal{I}|^{1/2},0,\ldots,0)^{T}$.

Centering in the above sense is the continuous-time version of column-
centering an $n\times p$ data matrix. As alluded to in §2, in some
applications, a CT analogue of _row_ -centering is warranted, i.e., an
estimated common trend $\widehat{m}(t)=\frac{1}{p}\sum_{u=1}^{p}x_{u}(t)$
should be subtracted from each of the $p$ curves before computing the
covariance or correlation. Consider, for example, the Canadian temperature
curves in Figure 2. Without removal of the common trend, the CT correlation
matrix $\mbox{\boldmath$R$}^{*}$ for this example has all elements above 0.9,
since the temperatures at all stations are relatively high in the summer and
low in the winter, yielding the uninformative correlation matrix in the lower
left panel of Figure 2. But with detrending applied, the result is the more
interesting correlation matrix in the lower right panel. The 35 stations are
listed geographically, roughly from east to west, except for the last six
stations, which are located in the Canadian North. The top and middle panels
of Figure 2 suggest that in many cases, stations with nearly the same
longitude have similar temperature curves. This explains why the CT
correlation matrix has several blocks of very high positive correlations along
the main diagonal: for instance, the first six stations are those located in
the Atlantic provinces, and the correlations among them are all above 0.78 and
mostly above 0.93. In Appendix C we propose a simple diagnostic to help
determine, for a given data set, whether detrending is advisable.

Since the Canadian weather curves are spatially indexed, Scheipl et al. (2015)
modeled the spatial correlation to improve precision in an FDA context (see
also Liu et al., 2017). In the above analysis, however, the correlations among
stations are the estimands of interest, rather than a means to a further
inferential end as in spatial FDA.

##  5 CT principal component analysis

In classical MVA, the (sample) principal components are linear combinations
$\mbox{\boldmath$v$}_{1}^{T}\mbox{\boldmath$x$},\ldots,\mbox{\boldmath$v$}_{p}^{T}\mbox{\boldmath$x$}$
that have maximal variance, subject to mutual orthonormality of the
$\mbox{\boldmath$v$}_{j}$’s. If the sample covariance matrix $S$ has
orthonormal eigenvectors
$\mbox{\boldmath$e$}_{1},\ldots,\mbox{\boldmath$e$}_{p}$ with corresponding
eigenvalues $\lambda_{1}\geq\ldots\geq\lambda_{p}$, then by a standard result
on Rayleigh quotient maximization, the solution to this iterative optimization
problem is $\mbox{\boldmath$v$}_{m}=\mbox{\boldmath$e$}_{m}$ for
$m=1,\ldots,p$.

Continuous-time PCA seeks linear combination _functions_
$\mbox{\boldmath$v$}_{m}^{T}\mbox{\boldmath$x$}(t)$ that achieve maximal
variance in the above sense, but with the variance of
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$x$}(t)$ defined as in (4) with respect
to time, as
$|\mathcal{I}|^{-1}\int_{\mathcal{I}}[\mbox{\boldmath$v$}^{T}\mbox{\boldmath$x$}^{c}(t)]^{2}dt=\mbox{\boldmath$v$}^{T}\mbox{\boldmath$S$}^{*}\mbox{\boldmath$v$}$.
By the same argument as in the classical case, but with
$\mbox{\boldmath$S$}^{*}$ in place of $S$, the solution is given by the
eigenvectors of $\mbox{\boldmath$S$}^{*}$: that is,
$\mbox{\boldmath$v$}_{m}=\mbox{\boldmath$e$}^{*}_{m}$ for $m=1,\ldots,p$,
where $\mbox{\boldmath$e$}^{*}_{1},\ldots,\mbox{\boldmath$e$}^{*}_{p}$ are
orthonormal eigenvectors of $\mbox{\boldmath$S$}^{*}$ in descending order of
the corresponding eigenvalues $\lambda^{*}_{1}\geq\ldots\geq\lambda^{*}_{p}$.

Similarly to the classical setting, we then have the expansion
$\mbox{\boldmath$x$}^{c}(t)=\sum_{j=1}^{p}\mbox{\boldmath$e$}^{*}_{j}s_{j}(t)$,
where $s_{j}(t)$ is the $j$th score function, defined by
$s_{j}(t)=\mbox{\boldmath$e$}^{*T}_{j}\mbox{\boldmath$x$}^{c}(t)$, and it
follows straightforwardly from this definition that the CT covariance matrix
(4) of $[s_{1}(t),\ldots,s_{p}(t)]^{T}$ is
$\mbox{Diag}(\lambda^{*}_{1},\ldots,\lambda^{*}_{p})$. Assuming (8), since
each component of $\mbox{\boldmath$x$}^{c}(t)$ is a linear combination of the
basis functions $\phi_{1}(t),\ldots,\phi_{K}(t)$, the same is true for each of
the scores $s_{1}(t),\ldots,s_{p}(t)$. Appendix D presents an application of
CT PCA to electroencephalography data.

Canonical correlation analysis (CCA) is another classical MVA method that is
somewhat related to PCA. Here one has two sets of variables and seeks pairs of
linear combinations, one for each set, having maximal correlation in an
iterative sense. A continuous-time version of CCA is presented in Appendix E.

##  6 CT Fisher’s linear discriminant analysis

In its classical form, Fisher’s linear discriminant analysis (LDA; see
Appendix F) seeks linear combinations
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$x$}$ that optimally separate $G\geq 2$
sets of observations $\mbox{\boldmath$X$}_{1},\ldots,\mbox{\boldmath$X$}_{G}$
into which the data set is divided _a priori_. Like PCA and CCA, Fisher’s LDA
finds the desired linear combinations of the variables by solving an iterative
optimization problem, in this case based on a partition of the total sums of
squares and cross-products matrix $\mbox{\boldmath$T$}=n\mbox{\boldmath$S$}$
into _within-groups_ and _between-groups_ components,
$\mbox{\boldmath$T$}=\mbox{\boldmath$W$}+\mbox{\boldmath$B$}$.

We propose a continuous-time variant of Fisher’s LDA that optimally
discriminates among $G\geq 2$ _subintervals_
$\mathcal{I}_{1},\ldots,\mathcal{I}_{G}$ into which $\mathcal{I}$ is
partitioned. For $g\in\\{1,\ldots,G\\}$, let
$\bar{\mbox{\boldmath$x$}}^{*}_{g}=\mbox{\boldmath$C$}^{T}\bar{\mbox{\boldmath$\phi$}}_{\mathcal{I}_{g}}=|\mathcal{I}_{g}|^{-1}\int_{\mathcal{I}_{g}}\mbox{\boldmath$x$}(t)dt$.
Analogously to the classical setting, we have the decomposition
$\mbox{\boldmath$T$}^{*}=\mbox{\boldmath$W$}^{*}+\mbox{\boldmath$B$}^{*}$,
where

$\displaystyle\mbox{\boldmath$T$}^{*}$ $\displaystyle=$
$\displaystyle\int_{\mathcal{I}}[\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}][\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}]^{T}dt,$
$\displaystyle\mbox{\boldmath$W$}^{*}$ $\displaystyle=$
$\displaystyle\sum_{g=1}^{G}\int_{\mathcal{I}_{g}}[\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}_{g}][\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}_{g}]^{T}dt,$
$\displaystyle\mbox{\boldmath$B$}^{*}$ $\displaystyle=$
$\displaystyle\sum_{g=1}^{G}|\mathcal{I}_{g}|(\bar{\mbox{\boldmath$x$}}^{*}_{g}-\bar{\mbox{\boldmath$x$}}^{*})(\bar{\mbox{\boldmath$x$}}^{*}_{g}-\bar{\mbox{\boldmath$x$}}^{*})^{T}.$

(see Appendix F for a derivation). The CT Fisher’s linear discriminants are
then functions
$\mbox{\boldmath$v$}_{1}^{T}\mbox{\boldmath$x$}(t),\ldots,\mbox{\boldmath$v$}_{s}^{T}\mbox{\boldmath$x$}(t)$
that optimally separate the subintervals in the sense of iteratively
maximizing
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$B$}^{*}\mbox{\boldmath$v$}/\mbox{\boldmath$v$}^{T}\mbox{\boldmath$W$}^{*}\mbox{\boldmath$v$}$,
subject to
$\mbox{\boldmath$v$}_{i}^{T}\mbox{\boldmath$W$}^{*}\mbox{\boldmath$v$}_{j}=0$
for $i\neq j$. By the same argument as in the classical setting, the optimally
separating vectors $\mbox{\boldmath$v$}_{1},\ldots,\mbox{\boldmath$v$}_{s}$
are eigenvectors of $\mbox{\boldmath$W$}^{*-1}\mbox{\boldmath$B$}^{*}$
corresponding to the positive eigenvalues in descending order. Once again, our
implementation relies on the basis representation (1) along with (2): we have

$\displaystyle\mbox{\boldmath$T$}^{*}$ $\displaystyle=$
$\displaystyle|\mathcal{I}|\mbox{
}\mbox{\boldmath$C$}^{T}\mbox{\boldmath$Q$}\mbox{\boldmath$C$},$
$\displaystyle\mbox{\boldmath$W$}^{*}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$C$}^{T}\left(\sum_{g=1}^{G}|\mathcal{I}_{g}|\mbox{\boldmath$Q$}_{\mathcal{I}_{g}}\right)\mbox{\boldmath$C$},$
(9) $\displaystyle\mbox{\boldmath$B$}^{*}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$C$}^{T}\left[\sum_{g=1}^{G}|\mathcal{I}_{g}|(\bar{\mbox{\boldmath$\phi$}}_{\mathcal{I}_{g}}-\bar{\mbox{\boldmath$\phi$}})(\bar{\mbox{\boldmath$\phi$}}_{\mathcal{I}_{g}}-\bar{\mbox{\boldmath$\phi$}})^{T}\right]\mbox{\boldmath$C$}.$

##  7 CT $k$-means clustering

Continuous-time $k$-means clustering partitions $\mathcal{I}$ into _temporal
clusters_ , each of which is a subinterval, or a finite union of subintervals,
of $\mathcal{I}$. Our algorithm is a CT extension of the basic $k$-means
procedure introduced in a 1957 manuscript later published as Lloyd (1982):

  1. 1.

Randomly select centers
$\mbox{\boldmath$m$}_{1},\ldots,\mbox{\boldmath$m$}_{k}\in\mathbb{R}^{p}$, by
evaluating $\mbox{\boldmath$x$}(t)$ at $k$ randomly chosen time points.

  2. 2.

Until a convergence criterion is met:

    1. (a)

Partition $\mathcal{I}$ into clusters $\mathcal{C}_{1},\ldots,\mathcal{C}_{k}$
defined by $\mathcal{C}_{i}=\\{t\in\mathcal{I}\mathrel{\mathop{\mathchar
58\relax}}d_{i}(t)<d_{j}(t)\mbox{ for all }j\neq i\\}$ for each $i$, where
$d_{i}(t)=\|\mbox{\boldmath$x$}(t)-\mbox{\boldmath$m$}_{i}\|$.

    2. (b)

For $i=1,\ldots,k$, update
$\mbox{\boldmath$m$}_{i}=|\mathcal{C}_{i}|^{-1}\int_{\mathcal{C}_{i}}\mbox{\boldmath$x$}(t)dt=\mbox{\boldmath$C$}^{T}\bar{\mbox{\boldmath$\phi$}}_{\mathcal{C}_{i}}$.

The nontrivial part is step 2(a). Since
$d_{i}(t)^{2}=\|\mbox{\boldmath$C$}^{T}\mbox{\boldmath$\phi$}(t)-\mbox{\boldmath$m$}_{i}\|^{2}=\mbox{\boldmath$\phi$}(t)^{T}\mbox{\boldmath$C$}\mbox{\boldmath$C$}^{T}\mbox{\boldmath$\phi$}(t)-2\mbox{\boldmath$\phi$}(t)^{T}\mbox{\boldmath$C$}\mbox{\boldmath$m$}_{i}+\mbox{\boldmath$m$}_{i}^{T}\mbox{\boldmath$m$}_{i}$,
whose first term does not depend on $i$, it follows that finding the cluster
center nearest to $\mbox{\boldmath$x$}(t)$ reduces to finding the
$i\in\\{1,\ldots,k\\}$ that minimizes

$A_{i}(t)\equiv-2\mbox{\boldmath$\phi$}(t)^{T}\mbox{\boldmath$C$}\mbox{\boldmath$m$}_{i}+\mbox{\boldmath$m$}_{i}^{T}\mbox{\boldmath$m$}_{i}.$

Suppose for simplicity that $k=2$. Then each of
$\mathcal{C}_{1}=\\{t\in\mathcal{I}\mathrel{\mathop{\mathchar
58\relax}}A_{1}(t)<A_{2}(t)\\}$ and
$\mathcal{C}_{2}=\\{t\in\mathcal{I}\mathrel{\mathop{\mathchar
58\relax}}A_{1}(t)>A_{2}(t)\\}$ is a union of subintervals of $\mathcal{I}$,
where the boundaries between these subintervals are the points $t$ at which
$A_{1}(t)$ and $A_{2}(t)$ cross, i.e., the zero crossings of

$A_{1}(t)-A_{2}(t)=2\mbox{\boldmath$\phi$}(t)^{T}\mbox{\boldmath$C$}(\mbox{\boldmath$m$}_{2}-\mbox{\boldmath$m$}_{1})+\mbox{\boldmath$m$}_{1}^{T}\mbox{\boldmath$m$}_{1}-\mbox{\boldmath$m$}_{2}^{T}\mbox{\boldmath$m$}_{2}.$
(10)

If the basis functions are cubic splines then this expression is a cubic
polynomial in $t$ on each of the intervals between knots. Thus finding the
transition points between the two clusters reduces to finding, for the cubic
polynomial on each of these intervals, any zero crossings that occur within
that interval. In this way the problem of partitioning infinitely many time
points among $k$ clusters is reduced to identifying the zeros of a finite set
of polynomials. Appendix G provides a formula for the polynomial coefficients,
and discusses how this CT approach to $k$-means clustering may improve
computational efficiency. The $k>2$ case is similar to the above, but not
every crossing point of any pair of $A_{i}(t)$’s is then a transition point.

Analogously to ordinary $k$-means, a CT $k$-means clustering decomposes the
total integrated sum of squared distances from $\bar{\mbox{\boldmath$\phi$}}$
into between-cluster and within-cluster components. The latter component, more
precisely the total within-cluster integrated squared distance from the
cluster mean, is the objective function of CT $k$-means and is readily shown
to equal
$\mbox{tr}\left[\mbox{\boldmath$C$}^{T}\left(\sum_{i=1}^{k}|\mathcal{C}_{i}|\mbox{\boldmath$Q$}_{\mathcal{C}_{i}}\right)\mbox{\boldmath$C$}\right]$
[cf. (9)]. Like ordinary $k$-means, the above CT algorithm is not guaranteed
to converge to the global minimum of this objective function. It is thus
recommended to run the algorithm repeatedly with different starting values,
and choose the solution that minimizes the objective.

In Appendix H we present a continuous-time variant of the silhouette width
(Rousseeuw, 1987), which indicates how clearly a given time point belongs to
its assigned cluster. The silhouette width can be plotted for each of a grid
of time points, and its average value provides a simple index of how tight the
clusters are. As with ordinary $k$-means clustering, the average silhouette
width may be used to guide the choice of $k$.

##  8 A simulation study

We conducted a simulation study to assess the performance of CT correlation
estimation relative to classical (discrete-observation) correlation between
$p=2$ variables. True curves $\mbox{\boldmath$x$}(t)$ were generated from a
bivariate Gaussian process on $[0,1]$ (see Appendix I) with between-curve
covariance matrix $\mbox{\boldmath$\Sigma$}=\left(\begin{array}[]{cc}1&0.5\\\
0.5&1\end{array}\right)$, and thus correlation $\rho=0.5$, and within-curve
covariance function $\Gamma$ having the squared exponential (radial basis
function) form

$\Gamma(s,t)=\exp[-(s-t)^{2}/(2\ell^{2})]$ (11)

where the _length parameter_ $\ell>0$ determines the wiggliness of the curves,
with higher $\ell$ implying greater smoothness. For each pair of curves we
obtain noisy observations (6) at each of 500 common observation times, with
noise standard deviation $\sigma=0.5$. For a given pair of curves
$\mbox{\boldmath$x$}(\cdot)$, the true CT correlation between them is
$r^{*}=r^{*}_{12}$ [see (5)]; the CT correlation _estimate_ , based on an
estimate $\hat{\mbox{\boldmath$x$}}(\cdot)$ of the pair of curves, is denoted
by $\hat{r}^{*}$ (see §3.3). We compare the ordinary correlation $\hat{r}$,
based on $n=500$ bivariate observations, versus CT correlation $\hat{r}^{*}$,
based on treating the number of observations $n$ as infinite and smoothing
both curves with respect to a $B$-spline basis of dimension $K=40$.

Figure 4: Top row: Examples of bivariate Gaussian processes on $[0,1]$ with
underlying between-curve correlation 0.5 and squared exponential within-curve
correlation (11), with length parameter $\ell=0.02,0.1,0.3$. Middle row:
Noisily observed versions of the same example pairs of curves. Bottom row:
Realized correlations $r^{*}$, and estimates $\hat{r},\hat{r}^{*}$ by ordinary
and continuous-time correlation, for simulated bivariate processes, with each
panel displaying results for 50 simulated noisy curve pairs akin to that shown
in the panel directly above it. Solid line is the line of identity; dashed
vertical line indicates the process correlation 0.5.

Figure 4 demonstrates how the CT correlation estimate $\hat{r}^{*}$
outperforms the ordinary sample correlation $\hat{r}$. The top row displays
example curve pairs $x_{1}(t),x_{2}(t)$ generated from a bivariate Gaussian
process on $[0,1]$ with underlying between-curve correlation 0.5 as above, and
within-curve (auto)correlation (11) for length parameter $\ell=0.02,0.1,0.3$.
For each of the above three values of $\ell$, we generated 50 pairs of curves
$x_{1},x_{2}$, like those shown in the top row of Figure 4, and for each pair
we sampled noisy observations (as described in the previous paragraph), as
displayed in the middle row. In the bottom row we have plotted the realized
correlations $r^{*}$ of the noiseless curves $x_{1},x_{2}$, against the
corresponding correlation estimates, $\hat{r}$ (ordinary) and $\hat{r}^{*}$
(continuous-time); each of the 50 points in each panel represents a pair of
curves like those in the top row.

The horizontal spread of the points in the bottom row represents variation of
the realized $r^{*}$ about the underlying process correlation 0.5 [step (i) of
(7)], which is seen to increase as the length parameter $\ell$ increases.
Deviation of the points from the line of identity represents error in the
estimates with respect to the realized $r^{*}$ [step (ii) of (7)]. This error
also increases with $\ell$, in particular for ordinary correlation, which is
markedly biased toward zero. This bias is easy to explain: the noise term in
(6) inflates variance without affecting between-curve covariance, resulting in
diminished correlation. Such attenuation of correlation due to noise is a
long-established phenomenon (e.g., Spearman, 1910) and is sometimes corrected
by multiplying the estimated correlation by the geometric mean of the two
variables’ reliabilities (Cohen et al., 2013). For time series as opposed to
independent observations, such a correction is not straightforward to
implement, but the CT correlation approach proceeds by instead denoising the
observed curves, resulting here in consistently accurate estimation of the
between-curve correlation $r^{*}$. This simulation study, then, offers clear
and concrete evidence for the benefits of CTMVA, even when ordinary MVA is
feasible. Appendix J displays the results of (i) simulations as above but with
underlying correlation 0.2 and 0.8, and (ii) simulations in which we vary the
noise standard deviation $\sigma$ as well as $n$, the number of points at
which the curves are observed. The overall conclusion is that unless $n$ is
small and/or the curves are very bumpy, the CT correlation outperforms
ordinary correlation, in some cases dramatically.

##  9 Example: Temperature and air pollution in Chicago

Figure 5: Top: The four standardized and smoothed variables from the Chicago
air pollution data. Vertical lines indicate transitions among the clusters
found by CT 3-means clustering, with cluster numbers shown along the top of
the plot. Middle: The same cluster transitions displayed separately by year.
Bottom: Trajectory of the first two PCs of the Chicago air pollution data.
Each label ($1987,\ldots,2000$) marks the end of the year indicated.

Wood (2017) performed several analyses relating air pollution to death rates,
using a data set of daily indicators for Chicago for the years 1987 to 2000
made available by Peng and Welty (2004). Here we focus on the four predictors
in these analyses, namely PM10 (median 2.5- to 10-$\mu$m particles per m3),
ozone (in parts per billion), SO2 (median sulfur dioxide measurement), and
temperature. Our goal here is exploratory analyses of these four metrics using
the CTMVA methods of §5–§7. All four variables were standardized and were
smoothed using a 200-dimensional $B$-spline basis and residuals with an AR(1)
residual structure (Wood, 2017).

Figure 5 displays results from continuous-time 3-means clustering. The top
panel shows the four standardized variables’ time series with vertical lines
indicating transitions among the three clusters. As can be discerned more
clearly in the middle panel, which shows the clustering separately by year,
there is a fairly consistent seasonal pattern: cluster 1 corresponds roughly
to the winter, cluster 2 to the summer, and cluster 3 to the spring and
autumn. Thus each calendar year is divided into portions belonging to clusters
1, 3, 2, 3, 1, in that order, with one exception: the winter of 1991–92, in
which cluster 1 does not appear. The clear seasonal pattern in the clusters
seems to be driven largely by temperature variations, along with the high
positive correlations of PM10 and ozone with temperature. SO2, on the other
hand, is negatively correlated with temperature. In the winter of 1991–92,
however, SO2 uncharacteristically dropped along with the other three
variables, and this may help explain the non-occurrence of cluster 1. Appendix
H provides further details of this analysis, in particular the silhouette
width and the stability of the 3-means solution. We note that when ordinary
3-means clustering is applied to the raw time series, there are over 1100
transitions among the resulting three clusters, making it difficult at best to
detect the cyclic seasonal pattern that emerges clearly from the continuous-
time analysis.

The above results are complemented by CT principal component analysis, which
found that most of the variation is explained by two PCs: 77.3% by the first
PC, with loadings 0.524, 0.505, -0.417, 0.544 for PM10, ozone, SO2 and
temperature, respectively, and 15.0% by the second PC, with loadings -0.351,
-0.347, -0.870, -0.006. In view of the correlations noted above, the first PC
is expected to be highest in the summer and lowest in the winter. This is
confirmed by the bottom panel of Figure 5, which displays the trajectory of
the first two PCs for the 14-year period, with labels ($1987,\ldots,2000$)
marking the end of each year. The colors are lightest at the beginning and end
of each year and darkest in the middle of the year, and this makes clear the
annual oscillation of the first PC. Both PCs’ year-end scores are highest for
1991, suggesting again that the winter of 1991–92 was atypical.

Figure 6: Left: The first two linear discriminants for the Chicago air
pollution data, displayed as functions of date. Center: Mean values of the
first two discriminants for each of the 14 years. Right: Mean PM10 and mean
SO2 for each year.

Continuous-time Fisher’s linear discriminant analysis, with the $G=14$
calendar years as subintervals among which we wish to discriminate, offers a
way to compare and contrast the different years. The coefficients of PM10,
ozone, SO2 and temperature, respectively, are given by
$\mbox{\boldmath$v$}_{1}=(0.0302,-0.0037,-0.0186,-0.0380)^{T}$ for the first
discriminant and $\mbox{\boldmath$v$}_{2}=(-0.0235,0.0134,-0.0105,0.0034)^{T}$
for the second. The plot of the time series of the first two discriminants in
Figure 6 suggests that, as in the cluster analysis, the year 1992 stands out:
the value of LD1 tends to be very high during that year. This impression is
confirmed by the middle panel, which shows the average values of the first two
discriminants in each year and reveals that 1992 is a clear outlier for LD1.
This finding is explained by the very high mean PM10 and extremely low mean
SO2 for 1992, as seen in the right panel.

##  10 Discussion

Our simulation study in §8 and 3-means cluster analysis in §9 suggest that the
denoised, smooth representation in CTMVA can, at least in some cases, lead to
improved performance over standard MVA with the raw data. In addition, as
illustrated by the WDI example of §3.2, CTMVA can be applied in certain
settings for which classical MVA is unavailable, such as variables observed at
disparate times. Another such setting is correlated point processes, such as
inhomogeneous Poisson processes. Ramsay and Silverman (2005) present a smooth
basis function approach to estimating a time-varying intensity function
$\iota(t)$. Given $p$ inhomogeneous Poisson processes with intensities
$\iota_{1}(t),\ldots,\iota_{p}(t)$ on a common time interval, one can compute
CT correlations among estimates of these intensities. However, more work is
required on reliable smoothing parameter estimation for intensity functions.

While the applicability of CT covariance and correlation even for pairs of
curves whose observation times have little or no overlap (as in the WDI
example of §3.2) is a clear advantage of the CT approach, there is an
accompanying pitfall: as seen in the left panel of Figure 3, curves are
estimated less precisely within time ranges with fewer observations, and
intuitively these times should be weighted less when estimating covariance.
Formally, the covariance between the $u$th and $v$th curves would then be
defined as
$\int_{\mathcal{I}}w(t)[x_{u}(t)-\bar{x}^{*}_{u}]][x_{v}(t)-\bar{x}^{*}_{v}]dt\mbox{
}/\int_{\mathcal{I}}w(t)dt$ for some function $w(t)$ [generalizing (4)].
Inverse variance is a standard choice for such a weight function, but it is
not obvious how to incorporate the inverse variances of two function estimates
into a single function $w(t)$. Minimizing the coefficient of variation, as in
Chen et al. (2014), is a possible way forward.

A key objective of ongoing work is interval estimation for CT correlation, via
analytic and resampling-based approaches. Optimal weighting, as described in
the previous paragraph, may sometimes yield narrower confidence intervals.

Another aim of future work is to develop _multilevel CTMVA_ for multiple
instances of continuous-time multivariate data, such as the full WDI and EEG
data sets of §3.2 and Appendix D, respectively. This may be contrasted with
two complex data structures extending the $n\times\infty$ setup of traditional
FDA: _multilevel_ and _multivariate_ FDA. As a crude oversimplification, we
may say that multilevel FDA (e.g. Crainiceanu et al., 2009) considers
$(n\times m)\times\infty$ data, in which a given function is observed $m$
times in each of $n$ individuals, whereas multivariate FDA (e.g. Happ and
Greven, 2018) concerns $n\times(v\times\infty)$ data, with $n$ instances of a
collection of $v$ functions. The multilevel CTMVA setup, of multiple instances
of $\infty\times p$ data, differs from the either multilevel or multivariate
FDA and may be more similar to the setup of Dubin and Müller (2005). But
whereas these authors estimate population dynamical correlation from a sample
of $\infty\times p$-type observations, the aim of multilevel CT correlation
would be to decompose the correlations among $p$ variables into “within” and
“between” components. Multilevel data will also allow for assessing the
reliability of CTMVA results, via extensions of intraclass correlation (Xu et
al., 2021).

Most of the methods of this paper are CT extensions of MVA methods involving
$p\times p$ matrices. Certain other MVA methods, such as multidimensional
scaling, begin with an $n\times n$ matrix, often representing dissimilarities
or distances among $n$ observations, or alternatively similarities among them.
Continuous-time ($n=\infty$) extensions are then more challenging than for
methods based on $p\times p$ matrices, since instead of $n\times n$ matrices
we must work with their infinite-dimensional counterparts, namely bounded
linear operators on $\mathcal{I}$. We are pursuing CT methods of this type in
ongoing work.

Our methods are implemented in the R package ctmva (Paul and Reiss, 2024), and
code for our real-data analyses can be found at
https://github.com/reissphil/ctmva.

## APPENDICES

##  Appendix A Daily average precipitation in Canada

Figure A7: Base-10 log mean daily precipitation in mm for the 35 Canadian
weather stations of Ramsay and Silverman (2005), in raw form (above) and after
smoothing (middle). Below: Correlations among the 35 weather stations,
obtained by treating the data as $n=365$ discrete observations (left) and by
the continuous-time approach (right).

The simulation study of §8 demonstrated the contrast between ordinary versus
CT correlations. As a real-data illustration of this phenomenon, we consider
the Canadian precipitation data of Ramsay and Silverman (2005). As in Figure
2, the data consist of daily values for the 35 Canadian weather stations, but
here the variable of interest is base-10 log of the average precipitation in
millimeters (see Ramsay et al., 2023). Figure A7 displays the data in raw form
(upper panel) and after smoothing (middle panel), using the same strategy as
for the temperature data in Figure 2. The lower left panel of Figure A7
presents the ordinary Pearson correlation among the 35 stations’ precipitation
levels, when the data are treated as $n=365$ discrete observations of $p=35$
variables.

The lower right panel shows the CT correlations based on the smooth curves.
The CT correlations are much stronger than raw correlations, and the
simulation results of §8 support viewing the former as more accurate than the
latter. For each correlation matrix, the stations are ordered by the angles
formed by the leading two eigenvectors; this ordering helps to reveal patterns
and groupings among the correlations (Friendly and Kwan, 2003). In both
correlation matrices the reordering puts the same nine stations at the end:
the first six of the original ordering, which are located on the Atlantic
coast, along with Vancouver, Victoria and Prince Rupert, on the Pacific coast.
The CT correlation matrix suggests that these nine stations form a tight
cluster, with positive correlations among them and negative correlations
between them and the inland stations. This pattern is much less clear in the
ordinary correlation matrix, due perhaps to the attenuation phenomenon
discussed in §8. As can be seen in the middle panel of Figure A7, the Atlantic
and Pacific stations tend to have the least precipitation in the summer, while
the inland stations have an opposite pattern.

The CT correlation matrix of the log mean precipitation values for the 35
Canadian weather stations, displayed in Figure A7, reveals that the nine
coastal stations (six on the Atlantic and three on the Pacific coast) tend to
be highly correlated with each other and anticorrelated with the 26 inland
stations. Figure A8 provides a map of this division of the 35 weather stations
into three subsets, using the same color-coding as in Figure A7. We remark
that Ramsay and Silverman (2005), as well as R package fda (Ramsay et al.,
2023) and a number of subsequent analyses of the Canadian weather data, divide
the stations into four regions, including 15 “Atlantic” and 5 “Pacific”
stations. The nine stations referred to in the previous paragraph, and
indicated in Figure A8, meet a stricter criterion of proximity to the Atlantic
or Pacific.

Figure A8: The 35 Canadian weather stations, color-coded by location on the
Atlantic or Pacific coast or inland.

##  Appendix B Computing (2)

For simplicity we consider only $\mathcal{S}=\mathcal{I}$. It is readily shown
that the matrix
$\mbox{\boldmath$Q$}=|\mathcal{I}|^{-1}\int_{\mathcal{I}}[\mbox{\boldmath$\phi$}(t)-\bar{\mbox{\boldmath$\phi$}}][\mbox{\boldmath$\phi$}(t)-\bar{\mbox{\boldmath$\phi$}}]^{T}dt$
of (2) can be written as

$\mbox{\boldmath$Q$}=|\mathcal{I}|^{-1}\mbox{\boldmath$G$}-\bar{\mbox{\boldmath$\phi$}}\bar{\mbox{\boldmath$\phi$}}^{T},$

where
$\mbox{\boldmath$G$}=\int_{\mathcal{I}}\mbox{\boldmath$\phi$}(t)\mbox{\boldmath$\phi$}(t)^{T}dt$,
i.e., the $K\times K$ matrix with $(k,\ell)$ element
$g_{k\ell}=\int_{\mathcal{I}}\phi_{k}(t)\phi_{\ell}(t)dt$.

If $\phi$ is an orthonormal Fourier basis on $\mathcal{I}$ then
$\mbox{\boldmath$G$}=\mbox{\boldmath$I$}_{K}$ and
$\bar{\mbox{\boldmath$\phi$}}=(|\mathcal{I}|^{-1/2},0,\ldots,0)^{T}$, so that
$Q$ is simply the diagonal matrix with main diagonal
$(0,|\mathcal{I}|^{-1},\ldots,|\mathcal{I}|^{-1})^{T}$.

For a $B$-spline basis, $\bar{\mbox{\boldmath$\phi$}}$ and $G$ are obtained by
numerical integration, with quadrature points $t_{1},\ldots,t_{Q}$ and
corresponding weights $w_{1},\ldots,w_{Q}$. Let $\Phi$ be the $Q\times K$
matrix of evaluation of the basis functions at the quadrature points, i.e.,
the $(q,k)$ entry of $\Phi$ is $\phi_{k}(t_{q})$. For basis objects in the R
package fda (Ramsay et al., 2023), this matrix can be generated with a call to
function eval.basis from that package. Let $W$ be the $Q\times Q$ diagonal
matrix whose diagonal elements are the quadrature weights. Then for
$k=1,\ldots,K$, $\mbox{\boldmath$\Phi$}^{T}\mbox{\boldmath$W$}\mbox{\bf
1}_{Q}$ has $k$th element
$\sum_{q=1}^{Q}w_{q}\phi_{k}(t_{q})\approx\int_{\mathcal{I}}\phi_{k}(t)dt$ and
thus

$\bar{\mbox{\boldmath$\phi$}}\approx|\mathcal{I}|^{-1}\mbox{\boldmath$\Phi$}^{T}\mbox{\boldmath$W$}\mbox{\bf
1}_{Q}.$

Similarly,
$\mbox{\boldmath$\Phi$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$\Phi$}$ has
$(k,\ell)$ element
$\sum_{q=1}^{Q}w_{q}\phi_{k}(t_{q})\phi_{\ell}(t_{q})\approx\int_{\mathcal{I}}\phi_{k}(t)\phi_{\ell}(t)dt$,
and thus

$\mbox{\boldmath$G$}\approx\mbox{\boldmath$\Phi$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$\Phi$}.$

Our implementation uses the the Newton-Cotes 7-point rule for numerical
quadrature, as suggested by Wand and Ormerod (2008). For cubic or lower-order
$B$-splines, the Newton-Cotes 7-point rule makes the above integrals exact
rather than approximate, so that $Q$ is obtained exactly.

##  Appendix C Diagnosing the need for detrending

We present here a simple diagnostic to indicate, for a given data set, whether
removing a common trend before computing the CT covariance or correlation (see
§4) would be beneficial. If we adopt an FDA viewpoint toward the $p$ functions
$x_{1}(t),\ldots,x_{p}(t)$, a strong common trend implies that all $p$
functions are quite similar to an underlying mean function $\mu_{x}(t)$. This
similarity can be quantified by using the function pffr (Ivanescu et al.,
2015; Scheipl et al., 2015) from the R package refund (Goldsmith et al., 2023)
to fit the simplest possible functional-response model, namely
$x_{i}(t)=\mu_{x}(t)+\varepsilon_{i}(t)$ for $i=1,\ldots,p$. Here $\mu_{x}(t)$
is a “functional intercept,” in line with the syntax for this model, pffr(X
$\sim$ 1) where X is a $p$-row matrix containing the functions. The (adjusted)
$R^{2}$ outputted by pffr for this model (based on the underlying model fitted
by the R package mgcv of Wood, 2017) gives the variance explained, for the
collection of all points $x_{i}(t)$ in all $p$ functions, by the fitted values
$\hat{x}_{i}(t)=\hat{\mu}_{x}(t)$. This “pointwise” $R^{2}$ (as opposed to the
_functional_ $R^{2}$ of Müller and Yao, 2008, which equals zero for this
model) provides a measure of the strength of the common trend. For the
Canadian weather data set, $R^{2}$ is found to be 71% for the temperature
curves, but only 4.9% for the log precipitation curves; hence our decision to
detrend the former (see §4) but not the latter (see Appendix A).

##  Appendix D Electroencephalography data

Figure A9: Loadings of the first three continuous-time principal components
for an example EEG trial. panel titles give the percent variance explained by
each PC.

We applied continuous-time PCA to a portion of the Begleiter
electroencephalography (EEG) data set from the UCI Machine Learning
Repository, made available in R by Helwig (2018). The data consist of multiple
EEG trials for each of a set of individuals, with each trial consisting of
electrical signals at $p=64$ channels with 256 time points within 1 second.
Loadings for the first 3 PCs for a randomly selected trial are shown in Figure
A9, and these have reasonably clear interpretations: PC1 roughly captures the
balance between signal in brain regions that are more central versus more
peripheral (as viewed from above), PC2 represents anterior versus posterior
regions, and PC3 represents left versus right hemisphere.

As suggested in §10, a natural next step would be multilevel CT-PCA to analyze
multiple trials nested within individuals, as are available with this data
set. The aim would be similar to the multilevel, longitudinal and structured
variants of functional PCA (Di et al., 2009; Greven et al., 2010; Shou et al.,
2015; Cui et al., 2023), but rather than decomposing a covariance operator,
the structure being decomposed is a $p\times p$ covariance or correlation
matrix as in Timmerman and Kiers (2002). We leave this development to future
work.

##  Appendix E Continuous-time canonical correlation analysis

In classical canonical correlation analysis (CCA) we are given an $n\times p$
matrix $X$ and an $n\times q$ matrix $Y$, and we seek _pairs_ of unit vectors
$\mbox{\boldmath$a$}\in\mathbb{R}^{p},\mbox{\boldmath$b$}\in\mathbb{R}^{q}$
such that the sample correlation of
$\mbox{\boldmath$a$}^{T}\mbox{\boldmath$x$},\mbox{\boldmath$b$}^{T}\mbox{\boldmath$y$}$
is maximized in an iterative sense:

  1. 1.

$\mbox{\boldmath$a$}_{1}^{T}\mbox{\boldmath$x$},\mbox{\boldmath$b$}_{1}^{T}\mbox{\boldmath$y$}$
attain the maximal sample correlation among all pairs
$(\mbox{\boldmath$a$},\mbox{\boldmath$b$})$ of unit vectors
$\mbox{\boldmath$a$}\in\mathbb{R}^{p},\mbox{\boldmath$b$}\in\mathbb{R}^{q}$.

  2. 2.

For $m=2,\ldots,\min\\{p,q\\}$,
$\mbox{\boldmath$a$}_{m}^{T}\mbox{\boldmath$x$},\mbox{\boldmath$b$}_{m}^{T}\mbox{\boldmath$y$}$
attain the maximal sample correlation among all pairs
$(\mbox{\boldmath$a$},\mbox{\boldmath$b$})$ of unit vectors such that
$\mbox{\boldmath$a$}^{T}\mbox{\boldmath$x$}$ is uncorrelated with
$\mbox{\boldmath$a$}_{1}^{T}\mbox{\boldmath$x$},\ldots,\mbox{\boldmath$a$}_{m-1}^{T}\mbox{\boldmath$x$}$
and $\mbox{\boldmath$b$}^{T}\mbox{\boldmath$y$}$ is uncorrelated with
$\mbox{\boldmath$b$}_{1}^{T}\mbox{\boldmath$y$},\ldots,\mbox{\boldmath$b$}_{m-1}^{T}\mbox{\boldmath$y$}$.

Then
$(\mbox{\boldmath$a$}_{m}^{T}\mbox{\boldmath$x$},\mbox{\boldmath$b$}_{m}^{T}\mbox{\boldmath$y$})$
are the $m$th pair of canonical variables, and their correlation is the $m$th
canonical correlation.

The solution to the iterative maximization problem is given in terms of the
covariance matrices of the two sets of variables, as well as their cross-
covariance matrix:

$\mbox{\boldmath$S$}_{xx}=n^{-1}\mbox{\boldmath$X$}^{cT}\mbox{\boldmath$X$}^{c},\quad\mbox{\boldmath$S$}_{yy}=n^{-1}\mbox{\boldmath$Y$}^{cT}\mbox{\boldmath$Y$}^{c},\quad\mbox{\boldmath$S$}_{xy}=n^{-1}\mbox{\boldmath$X$}^{cT}\mbox{\boldmath$Y$}^{c},\quad\mbox{\boldmath$S$}_{yx}=\mbox{\boldmath$S$}_{xy}^{T},$

where $\mbox{\boldmath$X$}^{c},\mbox{\boldmath$Y$}^{c}$ are column-centered
versions of the two matrices and
$\mbox{\boldmath$S$}_{xx},\mbox{\boldmath$S$}_{xy}$ are assumed to be
nonsingular. It can be shown (e.g., Rencher and Christensen, 2012) that

$\mbox{\boldmath$M$}_{1}=\mbox{\boldmath$S$}^{-1}_{xx}\mbox{\boldmath$S$}_{xy}\mbox{\boldmath$S$}^{-1}_{yy}\mbox{\boldmath$S$}_{yx}\quad\mbox{and}\quad\mbox{\boldmath$M$}_{2}=\mbox{\boldmath$S$}^{-1}_{yy}\mbox{\boldmath$S$}_{yx}\mbox{\boldmath$S$}^{-1}_{xx}\mbox{\boldmath$S$}_{xy}$

have the same positive eigenvalues $r_{1}^{2}\geq r_{2}^{2}\geq\ldots\geq
r^{2}_{\min\\{p,q\\}}$, and that if
$\mbox{\boldmath$a$}_{i},\mbox{\boldmath$b$}_{i}$ are eigenvectors of
$\mbox{\boldmath$M$}_{1},\mbox{\boldmath$M$}_{2}$ respectively corresponding
to the $i$th largest eigenvector, then
$\mbox{\boldmath$a$}_{i}^{T}\mbox{\boldmath$x$}$,
$\mbox{\boldmath$b$}_{i}^{T}\mbox{\boldmath$y$}$ are the $i$th canonical
variates, with canonical correlation $r_{i}$, which can be forced to be
positive by multiplying either $\mbox{\boldmath$a$}_{i}$ or
$\mbox{\boldmath$b$}_{i}$ by -1 if needed.

To develop a continuous-time version of CCA, suppose that two sets of
functions, $\mbox{\boldmath$x$}(t)=[x_{1}(t),\ldots,x_{p}(t)]^{T}$ and
$\mbox{\boldmath$y$}(t)=[y_{1}(t),\ldots,y_{q}(t)]^{T}$, are defined on
$\mathcal{I}$. We now seek pairs $(\mbox{\boldmath$a$},\mbox{\boldmath$b$})$
of unit vectors defining linear combinations
$\mbox{\boldmath$a$}^{T}\mbox{\boldmath$x$}(t),\mbox{\boldmath$b$}^{T}\mbox{\boldmath$y$}(t)$
of the two sets of functions that have maximal correlation in the above
iterative sense, but here, generalizing (4), the covariance of
$\mbox{\boldmath$a$}^{T}\mbox{\boldmath$x$}(t),\mbox{\boldmath$b$}^{T}\mbox{\boldmath$y$}(t)$
is given by

$|\mathcal{I}|^{-1}\int_{\mathcal{I}}[\mbox{\boldmath$a$}^{T}\\{\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}\\}][\mbox{\boldmath$b$}^{T}\\{\mbox{\boldmath$y$}(t)-\bar{\mbox{\boldmath$y$}}\\}]dt,$

and their correlation is then given by

$\frac{\int_{\mathcal{I}}[\mbox{\boldmath$a$}^{T}\\{\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}\\}][\mbox{\boldmath$b$}^{T}\\{\mbox{\boldmath$y$}(t)-\bar{\mbox{\boldmath$y$}}\\}]dt}{\sqrt{\int_{\mathcal{I}}[\mbox{\boldmath$a$}^{T}\\{\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}\\}]^{2}dt\int_{\mathcal{I}}[\mbox{\boldmath$b$}^{T}\\{\mbox{\boldmath$y$}(t)-\bar{\mbox{\boldmath$y$}}\\}]^{2}dt}}.$

These pairs of linear combinations may be called the _canonical functions_.

We assume for simplicity that both sets of functions are represented with
respect to a common $K$-dimensional basis $\phi$, i.e.,

$\mbox{\boldmath$x$}(t)=\mbox{\boldmath$C$}^{T}_{x}\mbox{\boldmath$\phi$}(t),\quad\mbox{\boldmath$y$}(t)=\mbox{\boldmath$C$}^{T}_{y}\mbox{\boldmath$\phi$}(t)$

for some matrices $\mbox{\boldmath$C$}_{x},\mbox{\boldmath$C$}_{y}$ of
dimension $K\times p$, $K\times q$ respectively. As with PCA, the continuous-
time solution to the iterative maximization problem ensues upon “translating”
the relevant matrices to the continuous-time setting. We do so by defining

$\mbox{\boldmath$M$}^{*}_{1}=\mbox{\boldmath$S$}^{*-1}_{xx}\mbox{\boldmath$S$}^{*}_{xy}\mbox{\boldmath$S$}^{*-1}_{yy}\mbox{\boldmath$S$}^{*}_{yx}\quad\mbox{and}\quad\mbox{\boldmath$M$}^{*}_{2}=\mbox{\boldmath$S$}^{*-1}_{yy}\mbox{\boldmath$S$}^{*}_{yx}\mbox{\boldmath$S$}^{*-1}_{xx}\mbox{\boldmath$S$}^{*}_{xy},$

where

$\mbox{\boldmath$S$}^{*}_{xx}=\mbox{\boldmath$C$}^{T}_{x}\mbox{\boldmath$Q$}\mbox{\boldmath$C$}_{x},\quad\mbox{\boldmath$S$}^{*}_{yy}=\mbox{\boldmath$C$}^{T}_{y}\mbox{\boldmath$Q$}\mbox{\boldmath$C$}_{y},\quad\mbox{\boldmath$S$}^{*}_{xy}=\mbox{\boldmath$C$}^{T}_{x}\mbox{\boldmath$Q$}\mbox{\boldmath$C$}_{y}^{T},\quad\mbox{\boldmath$S$}^{*}_{yx}=\mbox{\boldmath$S$}_{xy}^{*T},$

with $Q$ defined as in (2). Then, as in the classical case,
$\mbox{\boldmath$M$}^{*}_{1},\mbox{\boldmath$M$}^{*}_{2}$ have the same
positive eigenvalues; and if $\mbox{\boldmath$a$}_{i},\mbox{\boldmath$b$}_{i}$
are eigenvectors of $\mbox{\boldmath$M$}^{*}_{1},\mbox{\boldmath$M$}^{*}_{2}$
respectively corresponding to the $i$th largest eigenvector, then
$\mbox{\boldmath$a$}_{i}^{T}\mbox{\boldmath$x$}(t)$,
$\mbox{\boldmath$b$}_{i}^{T}\mbox{\boldmath$y$}(t)$ are the $i$th canonical
functions.

##  Appendix F Classical Fisher’s linear discriminant analysis, and the CT
within-between decomposition

As noted in §6, classical Fisher’s linear discriminant analysis (LDA) seeks
linear combinations $\mbox{\boldmath$v$}^{T}\mbox{\boldmath$x$}$ that
optimally separate $G\geq 2$ _a priori_ subsets
$\mbox{\boldmath$X$}_{1},\ldots,\mbox{\boldmath$X$}_{G}$ of the data set.
Assume that, for $g=1,\ldots,G$, $\mbox{\boldmath$X$}_{g}$ is $n_{g}\times p$
with $i$th row $\mbox{\boldmath$x$}_{gi}^{T}$; let

$\mbox{\boldmath$X$}=\left[\begin{array}[]{c}\mbox{\boldmath$X$}_{1}\\\
\vdots\\\ \mbox{\boldmath$X$}_{G}\end{array}\right]$

be the $n\times p$ combined data set, where $n=n_{1}+\ldots+n_{G}$; and let
$\bar{\mbox{\boldmath$x$}}=n^{-1}\mbox{\boldmath$X$}^{T}\mbox{\bf 1}_{n}$ and
$\bar{\mbox{\boldmath$x$}}_{g}=n_{g}^{-1}\mbox{\boldmath$X$}_{g}^{T}\mbox{\bf
1}_{n_{g}}$ denote the overall and group-$g$ mean vectors. We proceed by
partitioning the total sums of squares and cross-products matrix
$\mbox{\boldmath$T$}=n\mbox{\boldmath$S$}=\sum_{g=1}^{G}\sum_{i=1}^{n_{g}}(\mbox{\boldmath$x$}_{gi}-\bar{\mbox{\boldmath$x$}})(\mbox{\boldmath$x$}_{gi}-\bar{\mbox{\boldmath$x$}})^{T}$
into _within-groups_ and _between-groups_ components,
$\mbox{\boldmath$T$}=\mbox{\boldmath$W$}+\mbox{\boldmath$B$}$, where

$\mbox{\boldmath$W$}=\sum_{g=1}^{G}\sum_{i=1}^{n_{g}}(\mbox{\boldmath$x$}_{gi}-\bar{\mbox{\boldmath$x$}}_{g})(\mbox{\boldmath$x$}_{gi}-\bar{\mbox{\boldmath$x$}}_{g})^{T},\quad\mbox{\boldmath$B$}=\sum_{g=1}^{G}n_{g}(\bar{\mbox{\boldmath$x$}}_{g}-\bar{\mbox{\boldmath$x$}})(\bar{\mbox{\boldmath$x$}}_{g}-\bar{\mbox{\boldmath$x$}})^{T}.$

With these definitions,
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$v$}$ is the sum of
squared deviations of the linear combinations
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$x$}_{gi}$ from their within-group
means, while $\mbox{\boldmath$v$}^{T}\mbox{\boldmath$B$}\mbox{\boldmath$v$}$
is a weighted sum of these within-group means from the overall mean, with
weight $n_{g}$ for group $g$. (Some authors formulate Fisher’s LDA
differently, with these weights omitted. Moreover, some authors divide $W$ and
$B$ by constants that allow them to be interpreted as within- and between-
group covariance matrices, but this does not change the solution.) The linear
combination $\mbox{\boldmath$v$}_{1}^{T}\mbox{\boldmath$x$}$ where
$\mbox{\boldmath$v$}_{1}=\operatornamewithlimits{arg\,max}_{v}\mbox{\boldmath$v$}^{T}\mbox{\boldmath$B$}\mbox{\boldmath$v$}/\mbox{\boldmath$v$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$v$}$
can be said to maximally separate the groups, relative to the within-group
variation. It can be shown that $\mbox{\boldmath$v$}_{1}$ is the eigenvector
of $\mbox{\boldmath$W$}^{-1}\mbox{\boldmath$B$}$ corresponding to its largest
eigenvalue. If $s$ $(\leq\min\\{G-1,p\\})$ is the number of positive
eigenvalues of $\mbox{\boldmath$W$}^{-1}\mbox{\boldmath$B$}$, then for
$m=2,\ldots,s$, $\mbox{\boldmath$v$}_{m}$, the maximizer of
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$B$}\mbox{\boldmath$v$}/\mbox{\boldmath$v$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$v$}$
subject to
$\mbox{\boldmath$v$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$v$}_{1}=\ldots=\mbox{\boldmath$v$}^{T}\mbox{\boldmath$W$}\mbox{\boldmath$v$}_{m-1}=0$,
is given by the eigenvector of $\mbox{\boldmath$W$}^{-1}\mbox{\boldmath$B$}$
corresponding to its $m$th-largest eigenvalue. The linear combinations
$\mbox{\boldmath$v$}_{1}^{T}\mbox{\boldmath$x$},\ldots,\mbox{\boldmath$v$}_{s}^{T}\mbox{\boldmath$x$}$,
which optimally separate the groups in the above iterative sense, are known as
Fisher’s linear discriminants.

Similarly to the derivations of
$\bar{\mbox{\boldmath$x$}}^{*},\mbox{\boldmath$S$}^{*}$ in §3.1, the matrices
$\mbox{\boldmath$T$}^{*},\mbox{\boldmath$W$}^{*},\mbox{\boldmath$B$}^{*}$ in
§6 can be derived as limits of their classical counterparts
$\mbox{\boldmath$T$},\mbox{\boldmath$W$},\mbox{\boldmath$B$}$ times
$|\mathcal{I}|/n$ for observations on a uniform grid of $n$ time points, as
$n\rightarrow\infty$. Thus the CT decomposition
$\mbox{\boldmath$T$}^{*}=\mbox{\boldmath$W$}^{*}+\mbox{\boldmath$B$}^{*}$ can
be inferred from the classical decomposition
$\mbox{\boldmath$T$}=\mbox{\boldmath$W$}+\mbox{\boldmath$B$}$. A more direct
proof of the CT decomposition is as follows:

$\displaystyle\mbox{\boldmath$T$}^{*}$ $\displaystyle=$
$\displaystyle\sum_{g=1}^{G}\int_{\mathcal{I}_{g}}[\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}][\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}]^{T}dt$
$\displaystyle=$
$\displaystyle\sum_{g=1}^{G}\int_{\mathcal{I}_{g}}[\\{\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}_{g}\\}+\\{\bar{\mbox{\boldmath$x$}}^{*}_{g}-\bar{\mbox{\boldmath$x$}}^{*}\\}][\\{\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}_{g}\\}+\\{\bar{\mbox{\boldmath$x$}}^{*}_{g}-\bar{\mbox{\boldmath$x$}}^{*}]\\}^{T}dt$
$\displaystyle=$
$\displaystyle\mbox{\boldmath$W$}^{*}+\mbox{\boldmath$B$}^{*}+2\sum_{g=1}^{G}[\bar{\mbox{\boldmath$x$}}^{*}_{g}-\bar{\mbox{\boldmath$x$}}^{*}]\int_{\mathcal{I}_{g}}[\mbox{\boldmath$x$}(t)-\bar{\mbox{\boldmath$x$}}^{*}_{g}]dt.$

Since the integral in the last line equals zero by the definition of
$\bar{\mbox{\boldmath$x$}}^{*}_{g}$, we conclude that
$\mbox{\boldmath$T$}^{*}=\mbox{\boldmath$W$}^{*}+\mbox{\boldmath$B$}^{*}$ as
claimed.

##  Appendix G Computational aspects of CT $k$-means clustering

If the basis $\mbox{\boldmath$\phi$}(t)$ consists of cubic splines then,
within the $j$th inter-knot interval (for $j=1,\ldots,K-3$; see Ramsay and
Silverman, 2005), $\mbox{\boldmath$\phi$}(t)^{T}=(1\quad t\quad t^{2}\quad
t^{3})\mbox{\boldmath$L$}_{j}$ for a $4\times K$ matrix
$\mbox{\boldmath$L$}_{j}$ of polynomial coefficients. Substituting this into
(10), we conclude that $A_{1}(t)-A_{2}(t)=(1\quad t\quad t^{2}\quad
t^{3})\mbox{\boldmath$v$}_{j}$ where

$\mbox{\boldmath$v$}_{j}=2\mbox{\boldmath$L$}_{j}\mbox{\boldmath$C$}(\mbox{\boldmath$m$}_{2}-\mbox{\boldmath$m$}_{1})+\left[\begin{array}[]{c}\mbox{\boldmath$m$}_{1}^{T}\mbox{\boldmath$m$}_{1}-\mbox{\boldmath$m$}_{2}^{T}\mbox{\boldmath$m$}_{2}\\\
0\\\ 0\\\ 0\end{array}\right].$

Consequently, the zeroes of the cubic polynomial (10) in the $j$th inter-knot
interval can be found by inputting this vector $\mbox{\boldmath$v$}_{j}$ to R
function polyroot, which implements the Jenkins-Traub zero-finding algorithm.

Once the data have been reduced to a spline representation, the computing time
for CT $k$-means is insensitive to the number $n$ of observations in the raw
data. In experiments with 20 spline basis functions and $p=4$ (using the same
laptop that was used for the WDI analysis in §3.2), it took about 5 seconds to
run CT $k$-means for $k=2,3,\ldots,15$. By contrast, for ordinary $k$-means
with the same range of values of $k$, as implemented by the standard R
function kmeans, we found the computing time for the default algorithm of
Hartigan and Wong (1979) to increase slightly faster than linearly with $n$,
reaching about 80 seconds when $n\approx 3\times 10^{6}$. Moreover, in many
instances, convergence was not attained at the “quick-transfer” stage (step 6
in Hartigan and Wong, 1979). Computing times were slightly lower with two less
refined algorithms implemented in kmeans, one due to Lloyd (1982) and Forgy
(1965) and the other due to MacQueen (1967), but these often struggled to
converge. It must be acknowledged that for large $n$ the time required to
smooth the data with respect to the basis may at least partly offset the
computation-time advantage of the CT method. This depends, however, on how one
chooses the smoothing parameter in penalized basis smoothing. The R package
fda (Ramsay et al., 2023) offers the function Data2fd with a simple default
choice of the smoothing parameter. Wood (2017) summarizes some techniques for
optimal smoothness selection with large $n$, while Li et al. (2024) propose a
faster technique suitable for even larger $n$.

##  Appendix H Silhouette width and stability of the 3-means analysis

Given a clustering of $n$ discrete observations into $k>1$ clusters, Rousseeuw
(1987) defines the silhouette, or silhouette width, of the $i$th observation
as $s(i)=\frac{b(i)-a(i)}{\max\\{a(i),b(i)\\}}$, where $a(i)$ is the mean
distance from the given observation to other elements of the same cluster, and
$b(i)$ is the minimum, over the other $k-1$ clusters, of the mean distance
from the $i$th observation to elements of the cluster. When $b(i)>a(i)$, as is
usually but not always the case, the silhouette reduces to $1-a(i)/b(i)$, an
intuitive measure of how clearly the $i$th observation belongs to its assigned
cluster. An obvious CT extension is to define, for $t\in\mathcal{I}$,
$s(t)=\frac{b(t)-a(t)}{\max\\{a(t),b(t)\\}}$, where the mean distances
$a(t),b(t)$ are defined as above but as integrals over clusters, divided by
the length of the relevant cluster. It is straightforward to approximate these
integrals by the trapezoidal rule, and thereby obtain $s(t)$, at each of a
grid of equally spaced time points spanning $\mathcal{I}$. The left panel of
Figure A10 displays $s(t)$ for the CT 3-means clustering of the Chicago air
pollution data of §9. Within most of the subintervals into which each cluster
is divided, the silhouette width is low near the left boundary (as one would
expect near a transition between clusters), rises to a peak, then falls again
near the right boundary. For the long 1991–92 subinterval of cluster 3
discussed in the main text, there are two peaks with a trough in between. This
trough occurs around December 1991 and represents times that are on the
borderline between the assigned cluster (3) and the cluster to which the
winter is ordinarily assigned (1).

Figure A10: At left, a plot of the silhouette $s(t)$ for the CT 3-means
clustering of the Chicago air pollution data presented in §9. At right,
average silhouette width of $k$-means clusterings of the data, for
$k=2,3,\ldots,9$.

The right panel of Figure A10 shows the average of $s(t)$, over 6000 grid
points, for CT $k$-means clustering of the air pollution data with $k$ ranging
from 2 to 9. Rousseeuw (1987) suggests that “one way to choose $k$
‘appropriately’ is to select that value of $k$ for which $\bar{s}(k)$ [mean
silhouette width] is as large as possible.” According to this widely adopted
proposal, the plot indicates that $k=3$, while preferable to $k-4,\ldots,9$,
gives a less successful clustering of this data set than $k=2$ (due,
apparently, to the relatively low $s(t)$ within the spring/autumn cluster). On
the other hand, Rousseeuw (1987) did not advocate a rigid policy of treating
the mean-silhouette-maximizing $k$ as the true value to the exclusion of all
others; he left open the possibility that different numbers of clusters might
provide alternative informative views of the data. This is implied both by the
scare quotes in the above quotation, and the fact that Rousseeuw (1987)
himself presented both 2- and 3-cluster partitions of an example data set from
political science. Adopting a similarly flexible approach, we felt that
notwithstanding the mean silhouette results, $k=3$ gives a more interesting
cluster partition than $k=2$ for the air pollution data, and we therefore
chose to present the 3-means clustering in the text.

The 3-means clustering of this data set was run 50 times. In 46 instances,
there were 54 cluster transitions, exactly as in the solution described in §9.
The timings of each of these transitions never differed by more than one day
between any two of these 46 clusterings. In the other 4 instances, there were
40 cluster transitions, whose timings were again very consistent, with
discrepancies that were uniformly less than one day. The 54-transition
solution that was usually chosen had lower objective function (total within-
cluster integrated squared distance) than the 40-transition solution; the
ratio of between to total integrated squared distance, a metric that is
analogous to a standard measure of clustering quality, is about 88% for the
former solution and about 80% for the latter. These results suggest that for
this application, 3-means clustering is quite stable, with about 90% of
starting values leading to an optimal partition. At the same time, the
suboptimal solution attained in several instances underscores the importance
of trying multiple starting values.

##  Appendix I Multivariate Gaussian processes

The true functions $\mbox{\boldmath$x$}(\cdot)$ in the simulation study of §8
arise from a multivariate Gaussian process, which can be defined as follows
(Chen et al., 2020). Given a mean function
$\mbox{\boldmath$\mu$}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{I}\rightarrow\mathbb{R}^{p}$, within-curve covariance
function $\Gamma\mathrel{\mathop{\mathchar
58\relax}}\mathcal{I}\times\mathcal{I}\rightarrow\mathbb{R}$ and $p\times p$
between-curve covariance matrix $\Sigma$, we say that
$\mbox{\boldmath$x$}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{I}\rightarrow\mathbb{R}^{p}$ arises from the multivariate
Gaussian process MGP$(\mbox{\boldmath$\mu$},\Gamma,\mbox{\boldmath$\Sigma$})$
if, for any $t_{1},\ldots,t_{n}\in\mathcal{I}$, the $n\times p$ matrix
$\mbox{\boldmath$X$}_{t_{1},\ldots,t_{n}}\equiv[x_{u}(t_{i})]_{1\leq i\leq
n,1\leq u\leq p}$ has the matrix-variate normal distribution (Dawid, 1981;
Gupta and Nagar, 1999) with $n\times p$ mean matrix having $i$th row
$\mbox{\boldmath$\mu$}(t_{i})^{T}$, between-row covariance matrix
$\mbox{\boldmath$\Gamma$}_{t_{1},\ldots,t_{n}}\equiv[\Gamma(t_{i},t_{j})]_{1\leq
i,j\leq n}$ and between-column covariance matrix $\Sigma$.

The above process is unidentifiable in the sense that it is equal to
MGP$(\mbox{\boldmath$\mu$},h\Gamma,h^{-1}\mbox{\boldmath$\Sigma$})$ for any
$h>0$. But if we assume stationarity, identifiability can be established by
letting $\Gamma$ be an autocorrelation function $\Gamma(s,t)=a(|s-t|)$ with
$a(0)=1$, e.g., the squared exponential (11).

##  Appendix J Additional simulations

We conducted further simulations that were identical to those in §8 except for
the between-curve correlation. More specifically, as in §8, bivariate Gaussian
processes as in (11), with $\ell=0.02,0.1,0.3$, were observed at $n=500$
points with noise having standard deviation $\sigma=0.5$; but instead of
between-curve correlation $\rho=0.5$, we generated 50 replicates with
$\rho=0.2$ and 50 replicates with $\rho=0.8$. The results, analogous to those
in the bottom row of Figure 4, are displayed in Figure A11. As in §8, the
ordinary correlation estimates are attenuated and the CT correlations are much
more accurate (points are closer to the line of identity).

Figure A11: Simulation results, analogous to those in the bottom row of Figure
4, but with underlying between-curve correlation 0.2 (above) and 0.8 (below).

These results, along with those of Figure 4, provide strong visual evidence of
the superiority of CT over ordinary correlation when $n=500$ and $\sigma=0.5$.
We performed a more extensive simulation study, with $\rho=0.5$ and
$\ell=0.02,0.1,0.3$ as in Figure 4, but now with 50 replicates for each
combination of $n=50,100,200,500,1000,2000$ and $\sigma=0.2,0.5,0.8$. The
upper part of Figure A12 shows the median of the absolute error in estimating
the correlation $r^{*}$ between the two functions, i.e., $|\hat{r}-r^{*}|$ for
ordinary correlation and $|\hat{r}^{*}-r^{*}|$ for CT correlation, along with
error bars from the 5th to the 95th percentiles, for each setting. The
performance of ordinary correlation worsens as $\sigma$ increases, and does
not depend noticeably on $n$. CT correlation, on the other hand, performs
poorly in most settings with $n=50$, but as $n$ increases its median error
drops well below that for ordinary correlation. The lower part of Figure A12
shows the ratio of the root mean square error (RMSE) for CT versus ordinary
correlation. For $\ell=0.02$, regardless of $\sigma$, CT correlation performs
much worse with $n=50$ and somewhat worse with $n=100$ in terms of RMSE, but
outperforms ordinary correlation for $n\geq 200$. For $\ell=0.1,0.3$, i.e. for
smoother curves, CT correlation has lower RMSE for all values of $n$. In
summary, CT correlation is more accurate than ordinary correlation except when
the function is too bumpy and/or is observed too sparsely.

Figure A12: Above: Median absolute error in estimating the between-curve
correlation $r^{*}$ by ordinary and CT correlation, with error bars for 5th
and 95th percentiles. Horizontal axis represents the number $n$ of points at
which the curves are observed with noise, and $\sigma$ is the standard error
of the noise. “Length” denotes the parameter $\ell$ in (11). Below: Ratio of
root mean square error in estimating $r^{*}$ by CT versus ordinary
correlation. Several values lie beyond the vertical axis range, and are
indicated in parentheses.

Computing time for the CT correlation estimate $\hat{r}^{*}$, including
$B$-spline smoothing of the two functions, depended primarily on $n$. On an HP
Pavilion x360 with 8 GB of RAM, with the other settings fixed at
$\ell=0.02,\sigma=0.5$, the computing time in seconds was 0.35, 0.53, 0.97,
1.42 for $n=100,500,1000,2000$, respectively.

## References

  * Baum and Lake (2003) Baum, M. A. and D. A. Lake (2003).  The political economy of growth: democracy and human capital.  American Journal of Political Science 47(2), 333–347.
  * Chen et al. (2014) Chen, H., P. T. Reiss, and T. Tarpey (2014).  Optimally weighted $L^{2}$ distance for functional data.  Biometrics 70(3), 516–525.
  * Chen et al. (2020) Chen, Z., B. Wang, and A. N. Gorban (2020).  Multivariate Gaussian and Student-t process regression for multi-output prediction.  Neural Computing and Applications 32(8), 3005–3028.
  * Cohen et al. (2013) Cohen, J., P. Cohen, S. G. West, and L. S. Aiken (2013).  Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences (3rd ed.).  Routledge.
  * Crainiceanu et al. (2009) Crainiceanu, C. M., A.-M. Staicu, and C.-Z. Di (2009).  Generalized multilevel functional regression.  Journal of the American Statistical Association 104, 1550–1561.
  * Cui et al. (2023) Cui, E., R. Li, C. M. Crainiceanu, and L. Xiao (2023).  Fast multilevel functional principal component analysis.  Journal of Computational and Graphical Statistics 32(2), 366–377.
  * Dawid (1981) Dawid, A. P. (1981).  Some matrix-variate distribution theory: Notational considerations and a Bayesian application.  Biometrika 68(1), 265–274.
  * Di et al. (2009) Di, C.-Z., C. M. Crainiceanu, B. S. Caffo, and N. M. Punjabi (2009).  Multilevel functional principal component analysis.  Annals of Applied Statistics 3(1), 458.
  * Dubin and Müller (2005) Dubin, J. A. and H.-G. Müller (2005).  Dynamical correlation for multivariate longitudinal data.  Journal of the American Statistical Association 100, 872–881.
  * Forgy (1965) Forgy, E. W. (1965).  Cluster analysis of multivariate data: efficiency versus interpretability of classifications.  Biometrics 21, 768–769.
  * Friendly (2002) Friendly, M. (2002).  Corrgrams: Exploratory displays for correlation matrices.  The American Statistician 56(4), 316–324.
  * Friendly and Kwan (2003) Friendly, M. and E. Kwan (2003).  Effect ordering for data displays.  Computational Statistics & Data Analysis 43(4), 509–539.
  * Goldsmith et al. (2023) Goldsmith, J., F. Scheipl, L. Huang, J. Wrobel, C. Di, J. Gellar, J. Harezlak, M. W. McLean, B. Swihart, L. Xiao, C. Crainiceanu, and P. T. Reiss (2023).  refund: Regression with Functional Data.  R package version 0.1-30.
  * Greven et al. (2010) Greven, S., C. Crainiceanu, B. Caffo, and D. Reich (2010).  Longitudinal functional principal component analysis.  Electronic Journal of Statistics 4, 1022–1054.
  * Gupta and Nagar (1999) Gupta, A. K. and D. K. Nagar (1999).  Matrix Variate Distributions.  Chapman and Hall/CRC.
  * Happ and Greven (2018) Happ, C. and S. Greven (2018).  Multivariate functional principal component analysis for data observed on different (dimensional) domains.  Journal of the American Statistical Association 113, 649–659.
  * Hartigan and Wong (1979) Hartigan, J. and M. Wong (1979).  A k-means clustering algorithm.  Applied Statistics 28, 100–108.
  * Helwig (2018) Helwig, N. E. (2018).  eegkit: Toolkit for Electroencephalography Data.  R package version 1.0-4.
  * Honaker and King (2010) Honaker, J. and G. King (2010).  What to do about missing values in time-series cross-section data.  American Journal of Political Science 54(2), 561–581.
  * Ivanescu et al. (2015) Ivanescu, A. E., A.-M. Staicu, F. Scheipl, and S. Greven (2015).  Penalized function-on-function regression.  Computational Statistics 30, 539–568.
  * Jolliffe (2002) Jolliffe, I. T. (2002).  Principal Component Analysis (2nd ed.).  Springer.
  * Krivobokova and Kauermann (2007) Krivobokova, T. and G. Kauermann (2007).  A note on penalized spline smoothing with correlated errors.  Journal of the American Statistical Association 102, 1328–1337.
  * Li et al. (2024) Li, M., J. Zhang, and C. Meng (2024).  Nonparametric additive models for billion observations.  _Journal of Computational and Graphical Statistics_ , in press.
  * Liu et al. (2017) Liu, C., S. Ray, and G. Hooker (2017).  Functional principal component analysis of spatially correlated data.  Statistics and Computing 27, 1639–1654.
  * Lloyd (1982) Lloyd, S. P. (1982).  Least squares quantization in PCM.  IEEE Transactions on Information Theory 28(2), 129–137.
  * MacQueen (1967) MacQueen, J. (1967).  Some methods for classification and analysis of multivariate observations.  In L. M. Le Cam and J. Neyman (Eds.), Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1, pp. 281–297. University of California Press.
  * Mardia et al. (1979) Mardia, K. V., J. T. Kent, and J. M. Bibby (1979).  Multivariate Analysis.  New York: Academic Press.
  * Müller and Yao (2008) Müller, H.-G. and F. Yao (2008).  Functional additive models.  Journal of the American Statistical Association 103, 1534–1544.
  * Paul and Reiss (2024) Paul, B. and P. T. Reiss (2024).  ctmva: Continuous-time multivariate analysis.  R package version 1.4.0.
  * Peng and Welty (2004) Peng, R. D. and L. J. Welty (2004).  The nmmapsdata package.  R News 4(2), 10–14.
  * R Core Team (2023) R Core Team (2023).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Ramsay (1982) Ramsay, J. O. (1982).  When the data are functions.  Psychometrika 47(4), 379–396.
  * Ramsay et al. (2009) Ramsay, J. O., G. Hooker, and S. Graves (2009).  Functional Data Analysis with R and MATLAB.  New York: Springer.
  * Ramsay et al. (2023) Ramsay, J. O., G. Hooker, and S. Graves (2023).  fda: Functional Data Analysis.  R package version 6.1.4.
  * Ramsay and Silverman (2005) Ramsay, J. O. and B. W. Silverman (2005).  Functional Data Analysis (2nd ed.).  New York: Springer.
  * Reiss and Ogden (2009) Reiss, P. T. and R. T. Ogden (2009).  Smoothing parameter selection for a class of semiparametric linear models.  Journal of the Royal Statistical Society: Series B 71(2), 505–523.
  * Rencher and Christensen (2012) Rencher, A. C. and W. F. Christensen (2012).  Methods of Multivariate Analysis.  John Wiley & Sons, Inc.
  * Rousseeuw (1987) Rousseeuw, P. J. (1987).  Silhouettes: A graphical aid to the interpretation and validation of cluster analysis.  Journal of Computational and Applied Mathematics 20, 53–65.
  * Ruppert et al. (2003) Ruppert, D., M. P. Wand, and R. J. Carroll (2003).  Semiparametric Regression.  Cambridge, UK: Cambridge University Press.
  * Scheipl et al. (2015) Scheipl, F., A.-M. Staicu, and S. Greven (2015).  Functional additive mixed models.  Journal of Computational and Graphical Statistics 24(2), 477–501.
  * Shou et al. (2015) Shou, H., V. Zipunnikov, C. M. Crainiceanu, and S. Greven (2015).  Structured functional principal component analysis.  Biometrics 71(1), 247–257.
  * Silverman (1996) Silverman, B. W. (1996).  Smoothed functional principal components analysis by choice of norm.  Annals of Statistics 24(1), 1–24.
  * Spearman (1910) Spearman, C. (1910).  Correlation calculated from faulty data.  British Journal of Psychology 3(3), 271–295.
  * Timmerman and Kiers (2002) Timmerman, M. E. and H. A. L. Kiers (2002).  Three-way component analysis with smoothness constraints.  Computational Statistics & Data Analysis 40(3), 447–470.
  * Wand and Ormerod (2008) Wand, M. P. and J. T. Ormerod (2008).  On semiparametric regression with O’Sullivan penalized splines.  Australian & New Zealand Journal of Statistics 50(2), 179–198.
  * Wei and Simko (2021) Wei, T. and V. Simko (2021).  R package ‘corrplot’: Visualization of a Correlation Matrix.  Version 0.92.
  * Wood (2011) Wood, S. N. (2011).  Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.  Journal of the Royal Statistical Society: Series B 73(1), 3–36.
  * Wood (2017) Wood, S. N. (2017).  Generalized Additive Models: An Introduction with R (2nd ed.).  Boca Raton, Florida, USA: CRC Press.
  * Xu et al. (2021) Xu, M., P. T. Reiss, and I. Cribben (2021).  Generalized reliability based on distances.  Biometrics 77(1), 258–270.