# Extracting Common Time Trends from Concurrent Time Series: Maximum
Autocorrelation Factors with Applications

Matz A. Haugen Electronic address<EMAIL_ADDRESS>Stanford University
Bala Rajaratnam  Stanford University  Paul Switzer  Stanford University

###### Abstract

Concurrent time series commonly arise in various applications, including when
monitoring the environment such as in air quality measurement networks,
weather stations, oceanographic buoys, or in paleo form such as lake
sediments, tree rings, ice cores, or coral isotopes, with each monitoring or
sampling site providing one of the time series. The goal in such applications
is to extract a common time trend or signal in the observed data. Other
examples where the goal is to extract a common time trend for multiple time
series are in stock price time series, neurological time series, and quality
control time series. For this purpose we develop properties of MAF [Maximum
Autocorrelation Factors] that linearly combines time series in order to
maximize the resulting SNR [signal-to-noise-ratio] where there are multiple
smooth signals present in the data. Equivalence is established in a regression
setting between MAF and CCA [Canonical Correlation Analysis] even though MAF
does not require specific signal knowledge as opposed to CCA. We proceed to
derive the theoretical properties of MAF and quantify the SNR advantages of
MAF in comparison with PCA [Principal Components Analysis], a commonly used
method for linearly combining time series, and compare their statistical
sample properties. MAF and PCA are then applied to real and simulated data
sets to illustrate MAFs efficacy.

##  1 Introduction and Preliminaries

A common goal in the analysis of a collection of $p$ concurrent time series
$Z_{j}(t)$, $j=1,\dots,p$, observed at times $t=1,\dots,n$, is to extract a
common time trend which we refer to as the signal. Specifically, we look at
optimizing a linear combination $\bm{Y}(t)=\bm{w}^{\prime}\bm{Z}(t)$,
$t=1,\dots,n$, where $\bm{w}$ is an optimized coefficient $p$-vector. For
example, if the goal is to maximize variance over time of the combined series,
$Y(t)$, then this is equivalent to finding the first principal component in a
PCA (Principal Component Analysis). Then the coefficient vector $\bm{w}_{PCA}$
is the principal eigenvector of the cross-covariance matrix, $\bm{S}$, where
$\bm{S}_{ij}$ is the covariance over time between the pair of time series
$Z_{i}(t)$ and $Z_{j}(t)$. The idea of PCA is to reduce dimensionality through
retaining linear combinations of the data which have the highest variability.
Some applications of PCA to multiple time series analysis are given in Li et
al. (2007); Briffa et al. (2008); McShane and Wyner (2011); Jansen and
Rajaratnam (2014) Find references outside earth sciences. However, maximizing
variance across time, as PCA seeks to do, will not necessarily be well suited
to revealing coherent underlying latent time trends because PCA does not make
use of the specific time order of the data or optimize any property dependent
on temporal coherence. If the time order of the time series were permuted,
say, then the covariance matrix $\bm{S}$ and the coefficient vector
$\bm{w}_{PCA}$ are unchanged.

Arguably, an optimization criterion for the coefficient vector $\bm{w}$ for
combining the $p$ concurrent time series should specifically maximize a
measure of temporal coherence of the transformed time series, rather than the
time variance used in PCA.

###  1.1 MAF - Maximum Autocorrelation Factors

An alternative to PCA is Maximum Autocorrelation Factors (MAF) (Switzer and
Green, 1984; Shapiro and Switzer, 1989) where variance maximization is
replaced by autocorrelation maximization, which explicitly does depend on the
time ordering of the $p$-variate observations. The motivation for MAF is that
smoothly evolving time trends contained in time series data will enhance
autocorrelation. We show in Appendix B that the MAF-optimized coefficient
vector $\bm{w}_{MAF}$ is obtained as the leading eigenvector of the matrix

$\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2},$ (1.1)

where $\bm{S}_{\Delta}$ is the $p\times p$ covariance matrix of the time-
differenced time series. Any rescaling of the original time series,
$\bm{Z}(t)$, will preserve the MAF time series. This invariance property for
MAF is also derived in Appendix B. On the other hand, PCA component time
series are not invariant to rescaling or recombining of the original data.

Some applications of MAF to multiple time series analysis are given in Switzer
and Green (1984); Shapiro and Switzer (1989); Gallagher et al. (2014). Our
interest in MAF derives from applications to the analysis of multiple time
series of climate proxy data from tree ring measurements, described in Section
6. A fuller discussion of the analysis of tree ring data will be presented in
a separate paper. In this paper we shall focus on the methodological
development of the MAF framework.

To intuitively appreciate the difference between MAF and PCA, suppose we have
$p=2$ time series, one that is pure white noise and the other that is a linear
time trend without noise, with both series having unit variance over time.
Since PCA looks for a combined time series with maximal variance, it is
indifferent between the noisy time series with zero autocorrelation and the
clean time series with unit autocorrelation. On the other hand, MAF will put
all its weight on the noiseless linear time trend. If the two original time
series contained each a mixture of time trend and noise, then the MAF time
series will amplify the time trend relative to the noise.

###  1.2 An Illustration

Figure 1 shows an example with four parallel time series, rescaled to have
zero mean and unit variance. These 150-year time series are extracted from the
database used in Mann et al. (2008) and represent tree-ring time series. To
measure temporal coherence we introduce an empirical signal-to-noise ratio
(SNR), which is obtained by taking the ratio of two standard deviations; that
of a smoothed version of the time series and that of the associated residuals
after the smooth has been subtracted from the original. Standard deviations
are calculated by summing over the time steps. The annotated empirical SNR
suggest that the first two time series exhibit more evident temporal structure
than the last two time series. The corresponding PCA and MAF time series are
shown in Figure 2, and these are also rescaled to have zero mean and unit
variance. The MAF time series appears to concentrate the temporal structure
whereas PCA seems to exhibit more temporal noise. The empirical SNR of the MAF
time series is $1.46$ while that of the PCA time series is $0.92$. PCA and MAF
coefficient matrices are shown in Table 1 and we see that the MAF time series
up-weights the first two data time series and down-weights the last two data
time series.

Figure 1: Four tree ring time series, each one scaled to have unit variance and zero mean. Autocorrelation is annotated above each figure. Figure 2: MAFs and PCs of the time series shown in Figure 1. | MAF | PCA
---|---|---
1 | 0.80 | 0.59
2 | 0.30 | 0.58
3 | 0.24 | 0.42
4 | -0.47 | 0.37
Table 1: MAF and PCA coefficients of 4 time series.

###  1.3 Summary of results

In Section 2, we introduce the signal-plus-noise model and show that under
general conditions, the MAF time series yields the highest signal-to-noise
ratio among all possible combined time series. Equivalently, MAF also
maximizes the correlation between the combined time series and the underlying
signal time series. The PCA time series, on the other hand, maximizes signal
_plus_ noise variance rather than the ratio. PCA does not generally share the
MAF “oracle property”, i.e. finding the linear combination of time series
which is maximally correlated with the underlying signal. We show that the SNR
of the MAF time series is equal or greater than that of the PCA time series in
all situations involving one or more signals. Only in the trivial setting when
the noise is _iid_ , i.e. with zero cross-correlation and equal variance are
MAF and PCA equivalent. Otherwise, MAF increases the SNR compared to PCA.

We then extend the model to having $q\leq p$ multiple signals, where we
establish that first $q$ MAFs and Canonical Correlation Factors (CCFs) span
the same subspace that contains any linear combination of the underying signal
time series, thus extending the “oracle property” to the case of multiple
signals. Consequently, in a regression setting with one response time series
and a set of predictor time series, where the latter contains multiple
signals, MAF regression with $q$ factors will be optimal in a ‘least squares’
sense. It is assumed that the response signal is a particular linear
combination of the underlying set of signals present in the predictors. On the
other hand, since the first $q$ Principal Components (PCs) do not span the
subspace of signals, their regression on the response will be suboptimal in
the least squares sense.

In Section 3, a specific illustration is given where two groups of time series
are considered, each with different signal strengths present in combination
with noise. Explicit expressions are given for both MAF and PCA where we
replace the sample covariance matrices by their expected values under ther
model. We then derive the explicit form of the MAF and PCA coefficients
vectors in other models. Doing so allows us to investigate how the
coefficients change as functions of the noise cross-correlation, relative
signals strengths contained in each time series, and total number of time
series. We find that the leading MAF SNR improves compared to PCA as noise
cross-correlation, number of time series and/or signal strength differences
increase(s).

Section 5 explores the statistical properties of MAF and shows that MAF
coefficient estimates are consistent as the number of time steps are increased
while keeping the number of time series constant. Illustrations are also given
to quantify the difference between MAF and PCA regarding their correlations
with the underlying time trend. To determine the presence of a signal in the
data, a hypothesis testing procedure is presented where the null hypothesis is
a pure noise time series. Using resampling, we illustrate the power of the
test at different sample sizes and significance levels.

Application to tree ring time series in western North America is shown in
Section 6. We illustrate MAF and PCA for these time series. Both MAF and PCA
suggest underlying common time trends, but MAF appears to show these trends
more clearly. A null hypothesis test is highly significant and suggests the
presence of time trends in the data. Concluding remarks are presented in
Section 7.

##  2 The signal-plus-noise model

###  2.1 Preliminaries

We now formally define the Maximum Autocorrelation Factor (MAF). For a given
set of $p$ observed concurrent time series,
$\bm{Z}(t):\mathbb{N}\rightarrow\mathbb{R}^{p}$, the leading MAF coefficient
vector, is defined as the linear combination of the time series in $\bm{Z}(t)$
such that

$\bm{w}_{MAF}(\bm{Z})\equiv\bm{w}_{MAF}=\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{\text{Cor}(\bm{w}^{\prime}\bm{Z}(t),\bm{w}^{\prime}\bm{Z}(t+1))\\}.$
(2.1)

Similary, the leading Principal Component coefficient vector are defined as

$\bm{w}_{PCA}(\bm{Z})\equiv\bm{w}_{PCA}=\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{\text{Var}(\bm{w}^{\prime}\bm{Z}(t))\\}.$
(2.2)

Note that the MAF yields the optimal linear combination such that the
autocorrelation is maximized while the PC yields the linear combination that
maximizes variance. Furthermore, the leading MAF factor is defined as follows,

$\bm{Y}_{MAF}(\bm{Z}(t))=\bm{w}^{\prime}_{MAF}\bm{Z}(t)\text{ for }t=1,...,n.$
(2.3)

This single time series is the linear combination of the original $p$ time
series with maximal autocorrelation. With these definitions, we now proceed to
derive various properties related to these two techniques.

###  2.2 The model

Suppose $f(t):\mathbb{N}\rightarrow\mathbb{R}$ is a fixed but unknown
normalized underlying signal time series with zero mean and Euclidian norm
equal to 1 over the observation period $t=1,...,n$. We have $p$ observed
concurrent time series, $\bm{Z}(t):\mathbb{N}\rightarrow\mathbb{R}^{p}$, that
are represented as

$\bm{Z}(t)=\bm{s}(t)+\bm{\varepsilon}(t)=f(t)\cdot\bm{b}+\bm{\varepsilon}(t),$
$\text{i.e. }\quad\sum_{t}f(t)=0\text{ and }\sum_{t}f^{2}(t)=1,$ (2.4)
$\text{with}\>E[\bm{\varepsilon}(t)]=0\>\text{and
}Var[\bm{\varepsilon}(t)]=\bm{\Sigma}_{\varepsilon},\quad\forall t$

where $\bm{\varepsilon}(t)$ is the random $p$-variate covariance-stationary
residual noise time series and $\bm{b}=(b_{1},b_{2},\dots,b_{p})$ is a
coefficient vector, fixed and unknown. The quantities $f(t)$,
$\bm{\varepsilon}(t)$ and $\bm{b}$ are all unobserved and unknown. We call
this the ‘S+N model’. A linear combination of the $p$ observed times series
$\bm{Z}(t)$ is another time series $\bm{w}^{\prime}\bm{Z}(t)$, with
$\bm{w}^{\prime}\in\mathbb{R}^{p}$. The signal-to-noise ratio for the combined
time series is denoted $SNR(\bm{w})$ and defined as

$\displaystyle\frac{\text{Average signal mean square}}{\text{Average noise
variance}}=$
$\displaystyle\text{SNR}(\bm{w}):=\frac{\frac{1}{n}\sum_{t=1}^{n}Var[\bm{w}^{\prime}\bm{s}(t)]}{\frac{1}{n}\sum_{t=1}^{n}Var[\bm{w}^{\prime}\bm{\varepsilon}(t)]}$
$\displaystyle=$
$\displaystyle\frac{\frac{1}{n}\sum_{t=1}^{n}f^{2}(t)[(\bm{w}^{\prime}\bm{b})^{2}]}{\frac{1}{n}\sum_{t=1}^{n}\bm{w}^{\prime}\Sigma_{\varepsilon}\bm{w}}$
$\displaystyle=$
$\displaystyle\frac{(\bm{w}^{\prime}\bm{b})^{2}}{\bm{w}^{\prime}\Sigma_{\varepsilon}\bm{w}}.$
(2.5)

The MAF and PCA time series are examples of such linearly combined time series
with particular choices for $\bm{w}$.

Now define

$\frac{1}{n-1}\sum_{t=1}^{n-1}\varepsilon(t)\varepsilon(t+1)=k_{\varepsilon},\;\;\text{and
}\frac{1}{n-1}\sum_{t=1}^{n-1}f(t)f(t+1)=k_{f},$ (2.6)

where $k_{f}$ may be regarded as a measure of signal coherence. We now show
conditions under which the MAF time series maximizes SNR$(\bm{w})$ over
$\bm{w}$.

######  Proposition 1.

Suppose that the stationary time series model for the residual noise is such
that the $p\times p$ residual autocovariance matrix has the proportional form,

$\displaystyle\text{Cov}(\bm{\varepsilon}(t),\bm{\varepsilon}(t+1))=k_{\varepsilon}\text{Cov}(\bm{\varepsilon}(t))=k_{\varepsilon}\bm{\Sigma}_{\varepsilon},\quad\text{for
some }k_{\varepsilon}\in\mathbb{R},$ (2.7)

such that

$\text{Var}(\bm{w}^{\prime}\bm{Z}(t),\bm{w}^{\prime}\bm{Z}(t+1))=(\bm{w}^{\prime}\bm{b})^{2}k_{f}+k_{\varepsilon}\bm{w}^{\prime}\Sigma_{\bm{\varepsilon}}\bm{w},\quad\text{with
}k_{f}>k_{\varepsilon},$ (2.8)

where $k_{f}$ is the lag-1 autocorrelation of a normalized signal $f(t)$, as
given in Equation 2.6. Then MAF maximizes S/N and PCA maximizes S+N, i.e.,

$\displaystyle\bm{w}_{MAF}=$
$\displaystyle\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\frac{(\bm{w}^{\prime}\bm{b})^{2}}{\bm{w}^{\prime}\bm{\Sigma}_{\varepsilon}\bm{w}},\;\;\text{and}$
$\displaystyle\bm{w}_{PCA}=$
$\displaystyle\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{(\bm{w}^{\prime}\bm{b})^{2}+\bm{w}^{\prime}\bm{\Sigma}_{\varepsilon}\bm{w}\\}.$
(2.9)

###### Proof.

We show that maximizing SNR$(\bm{w})$ over linear combinations if $\bm{w}$ is
equivalent to maximizing the lagged autocorrelation, denoted $r(\bm{w})$, of
the combined time series $\bm{w}^{\prime}\bm{Z}(t)$. Now define the following,

$\displaystyle r(\bm{w})=$
$\displaystyle\frac{Cov(\bm{w}^{\prime}\bm{Z}(t),\bm{w}^{\prime}\bm{Z}(t+1)}{Cov(\bm{w}^{\prime}\bm{Z}(t))}=1-\frac{\text{Var}(\bm{w}^{\prime}\Delta\bm{Z}(t))}{2\text{Var}(\bm{w}^{\prime}\bm{Z}(t))},$
$\displaystyle\Delta\bm{Z}(t)=$
$\displaystyle\bm{Z}(t)-\bm{Z}(t-1)=\bm{b}^{\prime}\Delta
f(t)+\Delta\bm{\varepsilon}(t)\quad\text{is the time differenced data
vector,}$ $\displaystyle\Delta f(t)=$ $\displaystyle f(t)-f(t-1)\quad\text{is
the time-differenced signal,}$ $\displaystyle\Delta\bm{\varepsilon}(t)=$
$\displaystyle\bm{\varepsilon}(t)-\bm{\varepsilon}(t-1)\quad\text{is the time
differenced noise vector.}$ (2.10)

We can write

$\displaystyle\text{Var}(\bm{w}’\bm{Z}(t))=$
$\displaystyle(\bm{w}^{\prime}\bm{b})^{2}+\bm{w}^{\prime}\Sigma_{\bm{\varepsilon}}\bm{w}.$
(2.11)

Using (LABEL:*)eq:R, (LABEL:*)eq2, (LABEL:*)eq3 we can express the model
autocorrelation, $r(\bm{w})$, of the combined time series
$\bm{w}^{\prime}\bm{Z}(t)$ as

$r(\bm{w})=\frac{\text{SNR}(\bm{w})\cdot
k_{f}+k_{\varepsilon}}{\text{SNR}(\bm{w})+1}$ (2.12)

which is a monotone function of SNR$(\bm{w})$, if $k_{f}>k_{\varepsilon}$.
Hence, maximizing $r(\bm{w})$ is equivalent to maximizing SNR$(\bm{w})$. Since
MAF maximizes autocorrelation, MAF will also maximize the signal-to-noise
variance ratio over combinations of $p$ observable cross-correlated time
series, where each observable time series is a sum of a signal contribution
and a random noise contribution.

PCA, one the other hand, is defined as

$\bm{w}_{PCA}:=\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{\text{Var}(\bm{w}^{\prime}\bm{Z}(t))\\}=\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{Var[\bm{w}^{\prime}(\bm{b}^{\prime}f(t)+\bm{\varepsilon}(t))]\\}=\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{(\bm{w}^{\prime}\bm{b})^{2}+\bm{w}^{\prime}\Sigma_{\varepsilon}\bm{w}\\}.$
(2.13)

∎

The above theorem has important consequences. In signal extraction, maximizing
SNR$(\bm{w})$ is arguably more desirable than maximizing overall variance of a
linear combination of the input time series as in PCA. The MAF optimization
criterion is clearly more suited to the goal of extracting a common signal
component from multiple time series. It is also important to note that the MAF
time series is invariant to any rescaling of the input time series, shown in
Appendix B, whereas the PCA time series is scale dependent.

We now proceed to state the theoretical properties of MAF time series in terms
of four lemmas. First, we show that the MAF time series is maximally
correlated with the underlying signal time series under the S+N model. This
property is fundamentally important and is henceforth referred to as the
“oracle property” of MAF.

######  Lemma 1.

Consider the model given in Proposition 1, then

$\bm{w}_{MAF}=\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\;\;Cor[f(t),\bm{w}^{\prime}\bm{Z}(t)]$
(2.14)

###### Proof.

First note that the squared cross-correlation is given by

$Cor[f(t),\bm{w}^{\prime}\bm{Z}(t)]^{2}=\frac{\text{SNR}(\bm{w})}{\text{SNR}(\bm{w})+1}.$
(2.15)

The proof now follows immediately by letting $k_{f}=1$ and $k_{\varepsilon}=0$
in Equation 2.12 in the proof of Proposition 1. ∎

Remark: Note that Lemma 1 above is incidentally the defining property of
Canonical Correlation Analysis (CCA) with one underlying signal. However,
there is a fundamental difference: CCA requires the knowledge of the signal,
$f(t)$, while MAF does not, hence the above lemma being called the “oracle
property” of MAF.

We now proceed to show that the leading MAF time series is invariant to any
rescaling of the input time series.

######  Lemma 2.

Consider a data matrix $\bm{Z}\in R^{n\times p}$ where each column of $\bm{Z}$
represents a single time series of length $n$ and $Y_{MAF}(\bm{Z})$ represents
the leading MAF factor of $\bm{Z}$. Now let $\bm{A}$ be an invertible matrix
such that $\tilde{\bm{Z}}=\bm{Z}\bm{A}$. Then,

$Y_{MAF}(\bm{Z})=Y_{MAF}(\tilde{\bm{Z}}).$ (2.16)

###### Proof.

We shall show,

$Y_{MAF}(\bm{Z}):=\bm{Z}\bm{w}_{MAF}(\bm{Z})=\tilde{\bm{Z}}\bm{w}_{MAF}(\tilde{\bm{Z}})=:Y_{MAF}(\tilde{\bm{Z}}).$
(2.17)

First from Equation 2.1,

$\displaystyle\bm{w}_{MAF}(\bm{Z})=$
$\displaystyle\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\\{\text{Cor}(\bm{w}^{\prime}\bm{Z}(t),\bm{w}^{\prime}\bm{Z}(t+1))\\}$
$\displaystyle=$
$\displaystyle\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\frac{\bm{w}^{\prime}\bm{S}_{\delta}\bm{w}}{\bm{w}^{\prime}\bm{S}\bm{w}}$
(2.18)

where $\bm{S}_{\delta}=Cov(\bm{Z}(t),\bm{Z}(t+1))$ and
$\bm{S}=Cov(\bm{Z}(t))$. Then note,

$\displaystyle\bm{w}_{MAF}(\bm{Z}\bm{A})=$
$\displaystyle\operatorname*{argmax}_{\bm{w}\in\mathbb{R}^{p}}\frac{\bm{w}^{\prime}\bm{A}^{\prime}\bm{S}_{\delta}\bm{A}\bm{w}}{\bm{w}^{\prime}\bm{A}^{\prime}\bm{S}\bm{A}\bm{w}}$
$\displaystyle=$
$\displaystyle\operatorname*{argmax}_{\bm{u}\in\mathbb{R}^{p}}\frac{\bm{u}^{\prime}\bm{S}_{\delta}\bm{u}}{\bm{u}^{\prime}\bm{S}\
\bm{u}},\text{with }\bm{u}=\bm{A}\bm{w},$ (2.19)

as $\bm{A}$ is invertible. The above then gives

$\bm{w}_{MAF}(\bm{Z}\bm{A})=\bm{A}^{-1}\bm{w}_{MAF}(\bm{Z}).$ (2.20)

Thus,

$Y_{MAF}(\bm{Z})=:\tilde{\bm{Z}}\bm{w}_{MAF}(\tilde{\bm{Z}})=\bm{Z}\bm{A}\bm{A}^{-1}\bm{w}_{MAF}(\bm{Z})=\bm{Z}\bm{w}_{MAF}(\bm{Z})=:Y_{MAF}(\tilde{\bm{Z}}).$
(2.21)

∎

We now proceed to give an analytic representation of the MAF coefficient
vector.

######  Lemma 3.

Consider the ‘S+N model’ in Equation 2.4. It follows that the MAF coefficient
vector can be expressed as

$\bm{w}_{MAF}=\bm{\Sigma}_{\varepsilon}^{-1}\bm{b},$ (2.22)

###### Proof.

See proof of Lemma 4. ∎

Lastly, we show that the SNR of the MAF time series under the S+N model is
proportional to the expected value of a likelihood ratio statistic for a
Gaussian noise specification.

######  Lemma 4.

Consider the ‘S+N model’ in Equation 2.4 and the following set of hypotheses,

$\displaystyle H_{0}:$ $\displaystyle\bm{Z}_{n}(t)=\bm{b}+\bm{\varepsilon}(t)$
$\displaystyle H_{A}:$
$\displaystyle\bm{Z}_{n}(t)=f(t)\bm{b}+\bm{\varepsilon}(t),\quad
f(t)\neq\text{constant},$ (2.23)

such that, $\sum_{t}f(t)=0$ and $\sum_{t}f^{2}(t)=1$ as defined in
(LABEL:*)eq:model, and $\bm{\varepsilon}(t)\sim N(0,\bm{\Sigma})$. Then,

$SNR(\bm{w}_{MAF})=E[\ln(L_{A})-\ln(L_{0})],$ (2.24)

where $L_{A}$ and $L_{0}$ are the likelihoods of the two hypotheses given the
data matrix $\bm{Z}\in\mathbb{R}^{n\times p}$.

###### Proof.

See Appendix B. ∎

Switzer and Green (1984) show that MAF and PCA are equivalent in the special
and restrictive case where the noise covariance matrix is given by

$Cov(\bm{\varepsilon}(t))=\Sigma_{\varepsilon}=\sigma^{2}\bm{I},$ (2.25)

i.e., the noise component of each input has the same variance and these
$p\times p$ noise components have no cross-correlation. However, this
equivalence between MAF and PCA does not hold when there is noise cross-
correlation or heterogeneous noise variance.

###  2.3 Multiple signals model

We can generalize the S+N model to allow for multiple underlying signal time
series. Each of the $p$ observed concurrent time series is made up of its own
unknown smooth signal time series and its own superposed noise time series
representing short term fluctuations. The specific structure of the problem
represents each of these $p$ signal time series in terms of $q<p$ underlying
orthogonal factor time series, representing the reduced dimensionality of the
signal structure. The goal is to find $q$ new time series which are linear
combinations of the observed time series. These $q$ new time series aim to
recover the underlying orthogonal factor time series, i.e. the signals. We
show conditions under which the MAF linear combinations of the observed time
series achieve this concentration of the underlying signal information.

Our strategy for showing that the $q$ MAF time series jointly capture the
available signal information contained in the observed $p$-variate time series
is to demonstrate that the $q$-space spanned by MAF is the same as the
$q$-space obtained from a canonical correlation analysis (CCA) of the $p$
observed times series one the one hand and the $q$ unobserved signal time
series on the other hand. Theorem 1 below shows this equivalence, under
specific conditions for the additive noise component of the observed time
series. The equivalence, using the modeled noise covariance structure, implies
that MAF, which is computed _without_ specifying the underlying signal, can
capture the same signal information as a canonical correlation analysis which
_requires_ the signal specification. In this sense MAF can be said to have an
oracle property under the specified conditions insofar as covariances and
lagged covariances computed from the observed data approximate their modeled
structure. Thus, MAF is able to achieve the same result as CCA by taking
advantage of time order.

Consider a $p$-variate set of time series, $\bm{Z}(t),\>t=1,\dots n>p$,
comprised of $q\leq p$ normalized underlying smooth orthogonal signals,
$\bm{F}(t)=(f_{1}(t),f_{2}(t),...,f_{q}(t))$, and zero-mean $p$-variate noise,
$\bm{\varepsilon}(t)\sim\bm{\Sigma}_{\varepsilon}$. For the signal, we assume

$\frac{1}{n-1}\sum_{t=1}^{n-1}f_{i}(t)f_{j}(t+1)=k_{i}\delta_{ij}$ (2.26)

where $1>k_{1}\geq k_{2}\geq\dots\geq k_{q}>k_{\varepsilon}$ and $\delta_{ij}$
us the familiar Kronecker delta function111We neglect any non-orthogonality
that might arise between lagged versions of the signal time series.. For the
noise, assume a proportional covariance model
$Cov[\bm{\varepsilon}(t),\bm{\varepsilon}(t+1)]=k_{\varepsilon}\bm{\Sigma}_{\varepsilon}$.
A $p$-length signal strength vector, $\bm{b}_{i}$, describes the amount of
signal $f_{i}(t)$ present in each of the $p$ original time series
$\bm{Z}_{j}(t)$. With $\bm{B}\in\mathbb{R}^{(p\times q)}$ with columns
$(\bm{b}_{1},\bm{b}_{2},\dots,\bm{b}_{q})$, the full model is

$\displaystyle\bm{Z}(t)=\bm{B}\bm{F}(t)+\bm{\varepsilon}(t).$ (2.27)

Letting $diag(\bm{k})$ be the matrix formed by the $q$-vector
$\bm{k}=(k_{1},k_{2},\dots,k_{q})$ in the diagonal and zeros in the off
diagonal, we can write

$\displaystyle\bm{\Sigma}_{Z}=$ $\displaystyle
Cov[\bm{Z}(t),\bm{Z}(t)]=\bm{B}\bm{B}^{\prime}+\bm{\Sigma}_{\varepsilon}$
$\displaystyle\bm{\Sigma}_{\delta\bm{Z}}=$ $\displaystyle
Cov[\bm{Z}(t),\bm{Z}(t+1)]$ $\displaystyle=$
$\displaystyle\bm{B}Cov[\bm{F}(t),\bm{F}(t+1)]\bm{B}^{\prime}+Cov[\bm{\varepsilon}(t),\bm{\varepsilon}(t+1)]$
$\displaystyle=$
$\displaystyle\bm{B}diag(\bm{k})\bm{B}^{\prime}+k_{\varepsilon}\bm{\Sigma}_{\varepsilon},$

both assumed to be positive definite.

Canonical Correlation Analysis (CCA) looks for linear combinations of the
columns of $\bm{Z}$ which maximize correlation between linear combinations of
the signals contained in $F(t)$, while being orthogonal to each other. We
shall refer to these combinations as Canonical Correlation Factors (CCFs).

In Appendix B, we show that the first $q$ MAF and CCA factor coefficients for
$\bm{Z}$ are both contained in the range of
$\bm{\Sigma}_{\varepsilon}^{-1}\bm{B}$, formalized in the following theorem.

######  Theorem 1.

If $\min_{i}(k_{i})>k_{\varepsilon}$, the first $q$ CCA and MAF coefficient
vectors span the same hyperplane of dimension $q$ in $\mathbb{R}^{p}$.

Consequently, the first $q$ MAFs are optimal as regressors in the following
sense. Since the first $q$ CCFs maximize correlations of $q$ different linear
combinations of $f_{1}(t),\dots,f_{q}(t)$, we can construct a maximizer of any
linear combination of $f_{1}(t),\dots,f_{q}(t)$. Moreover, for a response
variable $y(t)=\sum_{q}\alpha_{i}f_{i}(t)$, there exists a linear combination
of the CCFs, $\hat{y}(t)$, for which the corresponding correlation,
Cor$(y(t),\hat{y}(t))$, is maximized. Thus, $\hat{y}(t)$ is an optimal
predictor of $y(t)$ using $q$ time series in a least squares sense. And
because MAF spans the same $q$-subspace as the first $q$ CCFs, by trasitivity,
MAF is also optimal in this sense. The benefit of MAF is that no knowledge of
the underlying signal is needed for its computation as opposed to CCA.

For PCA in the multiple signal case we find the eigenvectors of

$\bm{\Sigma}_{Z}=\bm{B}\bm{B}^{\prime}+\Sigma_{\varepsilon}.$ (2.29)

If the data time series has been normalized by their respective variance,
$\bm{\sigma}_{Z}^{2}=\text{diag}(\bm{\Sigma}_{Z})$ the diagonal of the
covariance matrix, the corresponding normalized PCA would we the eigenvectors
of

$\text{diag}(\bm{\sigma}_{Z})^{-0.5}\left[\sum_{i}\bm{b}_{i}\bm{b}_{i}^{\prime}+\Sigma_{\varepsilon}\right]\text{diag}(\bm{\sigma}_{Z})^{-0.5},$
(2.30)

where $\text{diag}(\bm{\sigma}_{Z})^{-0.5}$ has the variance in the diagonal
and zeros in the off-diagonal. In both cases there is no closed form for the
PCA eigenvectors. Moreover, the space spanned by the first $q$ PC coefficient
vectors are not the same as the space spanned by CCA, and thus PCA is sub-
optimal in this setting. This can easily be seen by noting that
$\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}_{i}$ is not an eigenvector of either
matrix. An important aspect of this sub-optimality comes from PCAs lack of
invariance under linear transformations. Looking at Equation 2.29, the only
situation in which MAF and PCA are equivalent is if
$\Sigma_{\varepsilon}=\bm{I}$.

In a situation with no noise, MAF will recover a multivariate mixture of
orthogonal signals into their separate components without loss of information.
For example, if two time series are supplied, both with a combination of a
linear and a quadratic signal and both mutually orthogonal, then the MAFs will
decompose these two into their separate forms.

######  Property 1.

Let $\bm{z}(t)=\bm{B}\bm{F}(t)$ and $\bm{F}(t)$ represents q concurrent
unknown and uncorrelated time series at time $t$ such that $q\leq p$ sorted in
decreasing autocorrelations, $k_{1}\geq...\geq k_{q}$, $\bm{k}$ in vector
form. And let $\bm{B}$ be an unknown $p\times q$ matrix. Then the MAF will
recover $\bm{F}(t)$. If $q=p$, $\bm{B}$ will also be recovered.

###### Proof.

See Appendix B. ∎

##  3 Illustrations of MAF/PCA comparisons in S+N model

We now consider two models where it is possible to derive closed form
expressions for the SNRs of the leading MAF and PC. These expressions allow us
to quantify the improvement that MAF yields over PCA, and get a firm
understanding of how each model parameter affects the different SNRs. A closed
form expression is also derived for the leading MAF coefficient,
$\bm{w}_{MAF}$.

###  3.1 Model I: Two groups of time series

Consider a scenario with two groups of concurrent time series following the
signal-plus-noise model, with time-independent noise. Both groups contain $q$
time series each. The following lemma gives the relationship between the
coefficients of each group of time series both for the MAF and the PCA case.

######  Lemma 5.

Consider two groups of $q$ time series, with SNR equal to $b_{1}$ and $b_{2}$,
with $b_{2}<b_{1}$ respectively, and where noise has equal variance and a
common cross-correlation of $\rho>\frac{-1}{p-1}$ between each $2q$ time
series. Let the total number of time series be represented by the $2q$-vectors
$\bm{Z}(t)$, for $t=1,...,n$. Then consider a linear combination of these $2q$
time series $\bm{w}^{\prime}\bm{Z}(t)$. The associated SNR of this linear
combination is given by

$\displaystyle\text{SNR}(\bm{w})=\text{SNR}(w_{1},w_{2})=$
$\displaystyle\frac{b_{1}^{2}q[1+\nu\gamma]^{2}}{(1-\rho)(1+\nu^{2})+\rho
q(1+\nu)^{2}}$ $\displaystyle\nu=$ $\displaystyle
w_{2}/w_{1},\quad\gamma=b_{2}/b_{1},$ (3.1)

where $w_{1}$ and $w_{2}$ represent the coefficient for each group of time
series. The maximum SNR, and also the MAF SNR, occurs when

$\nu_{MAF}:=\frac{w_{2}}{w_{1}}=\frac{\gamma(1-\rho+\rho q)-\rho
q}{1-\rho+\rho q-\gamma\rho q},$ (3.2)

which we call the MAF coefficient ratio.

Similarly, the PCA coefficients are determined by maximizing total variance,

$S+N(w_{1},w_{2})=S(w_{1},w_{2})+N(w_{1},w_{2})=\max_{\nu}\left\\{(qw_{1}b_{1}+qw_{2}b_{2})^{2}+(1-\rho)+\rho(qw_{1}+qw_{2})^{2}\right\\},$
(3.3)

which is maximized when

$\nu_{PCA}=\frac{w_{2}}{w_{1}}=\sqrt{\alpha^{2}+1}-\alpha,\text{ with
}\alpha=\frac{b_{1}^{2}-b_{2}^{2}}{2(b_{1}b_{2}+\rho)}.$ (3.4)

###### Proof.

For the MAF result, substitute the specific parameter values of the above
model into the general expression for SNR$(\bm{w})$ in (LABEL:*)eq:R to obtain
(LABEL:*)eq:5. Thereafter, find the maximum of the quadratic expression in
(LABEL:*)eq:5. Note that the minimum is attained when $\nu=-1/\gamma$. For the
PCA result, find the values of $w_{1}$ and $w_{2}$ for which (LABEL:*)eq:7 is
maximized under the constraint that $w_{1}^{2}+w_{2}^{2}=1$. ∎

Note that the input parameters investigated here are the cross-correlation in
the noise, $\rho$, the relative differences in the two groups’ SNR, $\gamma$,
and the overall number of time series, $p=2q$.

In particular, consider what happens to MAF and PCA SNR when changing the
number of time series in each group, $q$. Note that $\nu_{PCA}$ does not
depend on $q$, while $\nu_{MAF}$ does. Taking limits in $q$,

$\text{SNR}(\bm{w}_{MAF})\sim q\frac{b_{1}^{2}(1-\gamma)^{2}}{2(1-\rho)}\text{
as }q\rightarrow\infty.$ (3.5)

Similarly,

$\text{SNR}(\bm{w}_{PCA})\sim\frac{b_{1}^{2}(1+\nu_{PCA}\gamma)}{\rho(1+\nu_{PCA})^{2}}\text{
as }q\rightarrow\infty.$ (3.6)

Thus, the associated SNR of PCA approaches a constant, while SNR of MAF will
grow linearly with $q$. This implies that the MAF SNR continues to improve as
the number of time series increases, unlike PCA which reaches a plateau.
Furthermore, $\lim_{q\to\infty}\nu_{MAF}=-1$, a result that intuitively
follows from the fact that the noise has equal variance across the groups,
unlike the signal. Thus, if $w_{1}=-w_{2}$ and $q$ is large enough the noise
component will cancel while the signal remains.

In Figure 3, MAF and PCA SNR values are compared as $\gamma$ and $\rho$ are
changed. The ratios of the SNRs are plotted in a contour plot. Each panel
shows a different $p$, the total number of time series. We see that increasing
the cross-correlation, $\rho$, increases the difference between MAF and PCA
SNR while increasing $\gamma$ has the opposite effect. Increasing the number
of time series will exacerbate the difference between the SNRs, as explained
in the asymptotic analysis above.

Figure 3: The SNR of MAF divided by the SNR of PCA when $q={1,5,25}$, with
$\gamma$ and $\rho$ changing.

###  3.2 Model II: A model with common cross-correlation and different
variances

To generalize Model I, we allow each time series to have a unique signal
strength and noise variance. The following lemma derives the form the MAF
coefficient vector, $\bm{w}_{MAF}$, takes in this model.

######  Lemma 6.

Consider the multivariate time series model,

$\bm{Z}(t)=f(t)\bm{b}+\bm{\varepsilon}(t),$ (3.7)

where $f(t)$ is the signal time series, $b$ is the vector of signal strengths
for each time series, and $Cor[\varepsilon_{i}(t),\varepsilon_{j}(t)]=\rho$
and $Var[\varepsilon_{i}(t)]=\sigma_{i}^{2}$.

Then, the MAF coefficient vector $\bm{w}=(w_{1},\dots,w_{p})$, is given by

$w_{i}\propto\;\frac{b_{i}}{\sigma_{i}^{2}}-\frac{\rho}{1+\rho(p-1)}\sum_{j=1}^{p}\frac{b_{j}}{\sigma_{i}\sigma_{j}},\quad
i=1,2,\dots,p.$ (3.8)

###### Proof.

Express the noise covariance matrix as

$\bm{\Sigma}_{\varepsilon}=\bm{A}\left[\rho\bm{1}_{p}\bm{1}_{p}+(1-\rho)\bm{I}\right]\bm{A}$
(3.9)

where $A_{ii}=\sigma_{i}$ and $A_{ij}=0$ for $i\neq j$. Thus,

$\bm{\Sigma}_{\varepsilon}^{-1}=\frac{1}{1-\rho}\left[\bm{A}^{-2}-\frac{\rho}{1+\rho(p-1)}\bm{A}^{-1}\bm{1}\bm{1}^{\prime}\bm{A}^{-1}\right].$
(3.10)

Now it can be shown that $\bm{w}_{MAF}=\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}$
(see (LABEL:*)eq:24 in Appendix B) and the lemma follows by substitution. ∎

We now consider the special case where all input time series have the same
noise variance. Substitution into (LABEL:*)eq:6 gives the MAF coefficient
vector

$w_{i}\propto\;b_{i}-\frac{\rho}{1+\rho(p-1)}\sum_{j=1}^{p}b_{j},\quad
i=1,2,\dots,p.$ (3.11)

An alternative way to derive the MAF coefficients in this special case of
common noise variance for all time series is to find the eigenvectors of
(LABEL:*)eq:15 analytically. Using this method also illustrates that only in
this special case of common noise variance can the leading MAF coefficient
vector in (3.11) be constructed from a linear combination of the first two PC
coefficient vectors. However, when the noise variances are not all equal, this
will not be the case. In fact, no linear combination of the PCs can be used to
obtain the MAF time series. More details are given in Appendix A.

If furthermore, $\rho=0$, i.e. no cross-correlation between noise time series,
then the MAF and PCA coefficient vectors are the same and are proportional to
signal strength vector $\bm{b}$.

##  4 The MAF methodology on sampled time series

In the following section, we specify how the MAF methodology is implemented on
a given collection of sampled time series. First, an algorithm which
calculates the MAF factors is specified. Second, a set of methods for the
selection of number of MAFs is discussed. Lastly, we investigate how to
quantify the uncertainty in the estimated MAF factors in the S+N model.

###  4.1 Calculation of Maximum Autocorrelation Factors

Consider an input set of $p$ time series, each of which is recorded at $n$
time points. Let the column matrix $\bm{Z}\in\mathbb{R}^{n\times p}$ with
$n>p$ denote this collection of time series. First recall from Equation 1.1
that the MAF factors are defined as the eigenvectors of
$\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}$. The following operations are
implemented on the data matrix $\bm{Z}$ to obtain these eigenvectors.

First, we transform $\bm{Z}$ such that its covariance matrix, $\bm{S}$, is the
identity through spectral decomposition. Working with this transformed matrix
of new time series we compute the first differences in time and compute the
corresponding covariance matrix, $\bm{S}_{\Delta}$. Then, we obtain the
eigenvectors of the differenced covariance matrix via spectral decomposition.
These eigenvectors are in turn transformed back to the original coordinate
system to become the columns of the MAF coefficient matrix, defined as
$\bm{W}_{MAF}$. Finally, $\bm{W}_{MAF}$ is pre-multiplied by $\bm{Z}$ to yield
the MAF factors, and are the orthogonal time series with maximum
autocorrelation. Algorithm 1 formally specifies the calculation of these MAF
factors. For the purpose of this algorithm, the covariance operator is defined
as
$Cov(\bm{Z})=\sum_{i=1}^{n}[\bm{Z}_{i\cdot}-\sum_{k=1}^{n}\bm{Z}_{k\cdot}][\bm{Z}_{i\cdot}-\sum_{k=1}^{n}\bm{Z}_{k\cdot}]^{\prime}$,
where $\bm{Z}_{i\cdot}$ is the $i^{th}$ row of $\bm{Z}$.

Algorithm 1 Calculate MAF factors, $\bm{Y}\in\mathbb{R}^{n\times p}$, and MAF
coefficients $\bm{W}_{MAF}(\bm{Z}).$

0: $\bm{Z}\in\mathbb{R}^{n\times p}$

1: Calculate $\bm{S}_{\bm{Z}}=Cov(\bm{Z})$

2: Decompose $\bm{S}_{Z}$ such that $\bm{S}_{Z}=\bm{U}\bm{D}\bm{U}^{\prime}$
where $\bm{U}\in O(p)$ and $D$ is diagonal with the eigenvalues of
$\bm{\Sigma}_{Z}$.

3: Compute $\bm{X}=\bm{Z}\bm{U}\bm{D}^{-0.5}\bm{U}^{\prime}$.

4: Compute $\Delta\bm{X}_{i\cdot}=\bm{X}_{i\cdot}-\bm{X}_{(i+1)\cdot}$

5: Compute $\bm{S}_{\Delta}=Cov(\Delta\bm{X})$.

6: Decompose $\bm{S}_{\delta}=\bm{V}\bm{K}\bm{V}^{\prime}$ where $\bm{V}\in
O(p)$ and $\bm{K}$ is diagonal with eigenvalues in increasing order,
$K_{11}\leq K_{22}\leq\dots\leq K_{pp}$.

7: Let $\bm{W}_{MAF}(\bm{Z})=\bm{U}\bm{D}^{-0.5}\bm{U}^{\prime}\bm{V}$

8: Compute $\bm{Y}=\bm{Z}\bm{W}_{MAF}(\bm{Z})$.

9: Let $\bm{Y}_{\cdot j}$ be the $j^{th}$ column of $\bm{Y}$. For each $j$,
compute $c_{j}=\text{sign}\left[\sum_{i=1}^{n}i\bm{Y}_{\cdot j}\right]$.

9: $c_{j}\bm{Y}_{\cdot j}$ for each $j$ and $\bm{W}_{MAF}(\bm{Z}).$

###  4.2 Uncertainty quantification

Often there is a need to understand the sampling variability of the estimated
$\bm{w}_{MAF}(\bm{Z})$ and the associated MAF factors in the S+N model. One
natural way to undertake this is to use resampling of the data.

In resampling the time series we seek to preserve the underlying signal while
resampling the noise. As such, an underlying smooth signal estimate is
obtained by smoothing the original time series. Denote the smooth estimate as
$\bm{\tilde{}}{Z_{i}}(t)$, for $i=1,\dots,p$. Possible smoothing techniques
include local regression (Loess) (Cleveland, 1979) and spline smoothers
(Hastie et al., 2009).

The residuals between the original and the smooth time series,
$\hat{\bm{\varepsilon}}(t)=\bm{Z}(t)-\tilde{\bm{Z}}(t)$, is then resampled and
added back to the smooth original time series, $\tilde{\bm{Z}}(t)$. Resampling
can be done in blocks as there is temporal structure in the residuals. Denote
the resampled time series as $\bm{Z}^{\ast}(t)$. With the new set of time
series we recompute the MAF coefficients and MAF factors,
$\bm{w}_{MAF}(\bm{Z}^{\ast})$ and $\bm{Y}^{\ast}(t)$.

The above procedure can be repeated to obtain B instances of $Y_{i}^{\ast}(t)$
that can be aggregated to obtain a pointwise confidence interval around
$Y_{i}(t)$. In the event that the resampled MAFs are not well centered around
the original MAFs $\bm{Y}(t)$, information might have been lost as the noise
vector was resampled, compromising the shape of the MAF in question. This
would suggest that there is too much noise present for any isolation of a
signal. Thus, a significance test could be employed to determine the MAF’s
relevance. The procedure to resample the MAF coefficients and MAF factors is
formally given in Algorithm 2.

Algorithm 2 Resample the MAF factors, $\bm{Y}\in\mathbb{R}^{n\times p}$

0: $\bm{Z}\in\mathbb{R}^{n\times p}$

1: For the set of time series $Z_{i}(t)$ for $i=1,...,p$, calculate
$\bm{w}_{MAF}(\bm{Z})$ and $\bm{Y}(t)$.

2: Create a smooth time series from each original time series $Z_{i}(t)$ and
calculate the residual $\hat{\varepsilon}_{i}(t)=Z_{i}(t)-\tilde{Z}_{i}(t)$.

3: From the set of integers $[1,n]$, draw a sequence $\\{s_{i}\\}$ of $n$
integers with replacement, yielding the sequence
${s_{1}},{s_{2}},\dots,{s_{n}}$. The new set of residuals becomes
$\hat{\varepsilon}_{i}({s_{1}}),\hat{\varepsilon}_{i}({s_{2}}),\dots,\hat{\varepsilon}_{i}({s_{n}})$
which we shall call $\hat{\varepsilon}_{i}^{\ast}(t)$. This step could be
modified to resample blocks of residuals.

4: Let $Z_{i}^{\ast}(t)=\tilde{Z}_{i}(t)+\hat{\varepsilon}_{i}^{\ast}(t)$.

5: Calculate new MAFs from the resampled data,
$\bm{Y}^{\ast}(t)=\bm{Z}^{\ast}\bm{w}_{MAF}(\bm{Z}^{\ast})$.

6: Repeat steps 3-5 B times.

6: B realizations of $\bm{Y}^{\ast}(t)$ and $\bm{w}_{MAF}(\bm{Z}^{\ast})$.

###  4.3 Selection of number of MAFs

In real applications, the number of underlying signals is often unknown.
Determining the number of underlying signals can be done in various ways.
First, one can find that the eigenvalues of
$\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}$ and plot them in an
“autocorrelation scree plot” similar to what is done in PCA. Using this plot,
one can look for the presence of a shoulder to define the number of MAFs one
should retain. Alternatively, one could define a cutoff after some fraction
$\alpha$ (such as $95\%$), of the total autocorrelation that is contained in
preceding MAFs.

A second method employs cross validation to find the number of underlying
signals. By defining a hold-out block one can regress each original time
series in $\bm{Z}$ on the $k$ first MAFs for $k=1,\dots,p$. Then, select $k$
such that the RMSE on the hold-out block is minimized.

A third method involves using the framework of hypothesis testing, the
description of which is deferred to the section on statistical inference in
Subsection 5.3.

##  5 Statistical properties

Having looked at the model properties of MAF and PCA, we now turn to their
sampling properties. This section is divided into three parts. In the first
part, we show that the sample covariance and lagged covariance yield
consistent estimates of MAF coefficients and MAF factors under the S+N model
as the number of time steps grows. PCA estimates are treated similarly. In the
second part, a simulation study is conducted to compare MAF and PCA as signal
recovery techniques. We use the signal cross-correlation with MAF and PC time
series as the metric of comparison, and find that MAF is both more resilient
to increased noise and more suitable when the noise has cross-correlation. In
the third subsection, we introduce a hypothesis testing framework to test if
an underlying time trend extracted by MAF is statistically significant.

###  5.1 Consistency

The following theorem shows that as the number of time steps grows for a
$p$-variate time series, the MAF and PCA coefficients will converge to their
model values.

######  Theorem 2.

Consider a set of time series $Z_{n}(t)\in\mathbb{R}^{p}$, such that

$\displaystyle\bm{Z}_{n}(t)=$ $\displaystyle
f_{n}(t)\bm{b}+\bm{\varepsilon}_{n}(t)\quad t=1,\dots,n$
$\displaystyle\bm{\Delta}\bm{Z}_{n}(t)=$
$\displaystyle\bm{Z}_{n}(t)-\bm{Z}_{n}(t+1)=\Delta
f_{n}(t)\bm{b}+\Delta\bm{\varepsilon}_{n}(t)$ (5.1)

with $f_{n}(t)\in\mathbb{R}\>\>\forall
t=1,2,...,n,\bm{b},\bm{\varepsilon}_{n}(t)\in\mathbb{R}^{p}\>\>\forall
t=1,2,...,n$,
$\Delta\bm{\varepsilon}_{n}=\bm{\varepsilon}_{n}(t)-\bm{\varepsilon}_{n}(t+1)$,
and $\Delta f_{n}(t)=f_{n}(t)-f_{n}(t+1)$. Residual time series
$\bm{\varepsilon}_{n}$ is a weakly stationary $p$-variate time series and the
associated autocovariance is absolutely summable. The signal time series is
such that

$\frac{1}{n}\sum_{t=1}^{n}f_{n}(t)=0,\quad\frac{1}{n}\sum_{t=1}^{n}f_{n}^{2}(t)=1,\quad\frac{1}{n-1}\sum_{t=1}^{n-1}[\Delta
f_{n}(t)-\overline{\Delta f_{n}(t)}]^{2}=a,\quad\forall n,$ (5.2)

where $\Delta f_{n}(t)=\frac{1}{n}(f_{n}(1)-f_{n}(n))$.

Then,

$\displaystyle\bm{S}_{n}\overset{p}{\to}$ $\displaystyle\bm{\Sigma}\text{ as
}n\to\infty,$
$\displaystyle\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}\overset{p}{\to}$
$\displaystyle\bm{\Sigma}^{-1/2}\bm{\Sigma}_{\Delta}\bm{\Sigma}^{-1/2}\text{
as }n\to\infty$

###### Proof.

See Appendix B. ∎

###  5.2 Simulation study

We now undertake a simulation study to compare the MAF and PCA proceedures as
signal recovery techniques. First, we generate $100$ simulations of $p=3$
parallel time series of length $n=150$, using the S+N model of Section 2, viz.

$\bm{Z}(t)=f(t)\bm{b}+\bm{\varepsilon}(t)$ (5.4)

where $f(t)$ is a specified underlying signal time series shown in Figure 4.
This time series is a rescaled and interpolated version of the mean annual
surface time series for the northern hemisphere for the years 1850-2007, taken
from Mann et al. (2008). The vector $\bm{b}$ is the $p$-vector of signal
strengths. $\bm{\varepsilon}(t)$ is an $iid$ zero-mean Gaussian noise
$p$-vector for $t=1,\dots,n$. The $3\times 3$ cross-covariance matrix for
$\bm{\varepsilon}$ has a unit diagonal and common value $\rho$ in the off-
diagonal entries.

Figure 4: Signal used in the signal recovery process, shown with a mean zero
and unit variance.

For each of the 100 realizations of the three parallel time series we compute
the combined MAF$(t)$ time series and the combined PCA$(t)$ time series. The
cross-correlations of the MAF time series and the PCA time series with the
true signal time series as given in Figure 4 are used as a metric for
comparison.

Two specific simulations of the data are shown in Figure 5 with their
associated smoothed MAF and PCA time series on the right. We use a LOESS
filter to smooth the time series222We use local regression to smooth with 60
years in the span and tricubic weighting. The equivalent span as a fraction of
the total time series is $60/150=2/5$ and the tricubic weight go as
$(1-(d/60)^{3})^{3}$ with $d$ the distance from the point of interest.. The
first row of Figure 5 shows the three parallel time series with
$\bm{b}=(0.8,0.4,0.2)$ and cross-correlation $\rho=0.25$. The second row shows
a parallel time series with a weaker signal strength with
$\bm{b}=(0.4,0.2,0.1)$ and $\rho=0.25$. The results of the analysis are
compelling. Namely, the cross-correlation of the MAF time series with the
underlying signal for the first row of Figure 5 is $0.67$ and the PCA
equivalent is $0.52$, while for the second row MAF and PCA cross-correlations
are $0.35$ and $0.08$, respectively.

Figure 5: Top: One realization of the data with the signal shown in bold on
top of each time series. The Signal-to-Noise Ratio (SNR), corresponding to
$\bm{b}=(0.8,0.4,0.2)$, is annotated above each figure and each time series
has been scaled to have unit variance and zero mean. The smoothed MAF1 and PC1
are shown on the right. Bottom: Same as top with different SNR,
$\bm{b}=(0.4,0.2,0.1)$. The noise cross-correlation, $\rho=0.25$.

We proceed to undertake further analysis of this model in order to fully
understand how MAF and PCA perform when the cross-correlation, $\rho$ and the
signal strength vector $\bm{b}$ changes. The full set of scenarios that we
consider in this example are:

  * •

A fixed signal strength vector, $\bm{b}=(0.8,0.4,0.2)$, with changing noise
cross-correlation $\rho$.

  * •

A fixed noise cross-correlation $\rho=0.25$ with changing
$\bm{b}=(0.8c,0.4c,0.2c)$ for $c\in[0.5,2.5]$.

Figure 6 contains plots of signal cross-correlations with MAF and PCA time
series for each of the parameter combination scenarios. Each plotted point
represents an average over 100 simulations. MAF yields higher correlation with
the signal uniformly. It is clear that MAF takes advantage of cross-
correlation in the noise and uses it to amplify the signal, while PCA fails to
exploit this property in the noise and thus under-performs compared to MAF.

The signal information contained jointly across PC1 and PC2 is also less than
the information contained in MAF1 only. This can be seen by regressing the
signal of both PC1 and PC2 and extracting the root of the $R^{2}$ value. This
is the multivariate equivalent to correlation between the a signal and a
signal estimate. This result is also shown in Figure 6 under the legend
PC$1+2$.

Figure 6: Left: Correlation of signal estimate, using MAF or PCA, with true
signal while changing cross-correlation of noise. Right: Correlation of signal
estimate, using MAF or PCA, with true signal while multiplying the signal
strength vector by a factor as shown on the x-axis. The error bars show twice
the standard error from the mean using 100 repetitions.

###  5.3 Hypothesis Testing

Consider the same $p$-variate time series of length $n$ as described in
Equation 5.4. One might want to test whether a time signal is indeed present
in the data or not. We consider the following hypotheses:

$\displaystyle H_{0}:$ $\displaystyle\bm{Z}_{n}(t)=\bm{\varepsilon}(t)$
$\displaystyle H_{A}:$
$\displaystyle\bm{Z}_{n}(t)=f(t)\bm{b}+\bm{\varepsilon}(t),\quad
f(t)\neq\text{constant},$ (5.5)

where $\bm{\varepsilon}(t)$ is an $iid$ zero-mean Gaussian noise p-vector time
series with cross-correlation $\rho$ and unit variance.

To test for the presence of a signal in a MAF we introduce empirical signal-
to-noise ratio, $SNR_{empir}(x(t))$, as a function of an arbitrary time series
$x(t)$ for $t=1,...,n$, using a smooth version of $x(t)$, called
$\tilde{x}(t)$,

$SNR_{empir}(x(t))=\frac{SD(\tilde{x})}{SD(x-\tilde{x})},$ (5.6)

surpressing the $t$ argument for brevity and with
$SD=\sqrt{\sum_{t=1}^{n}\left[x(t)-\sum_{k=1}^{n}x(k)\right]^{2}}$.

Under the null model, there is no signal. So, subtracting out a smooth trend
should not affect the corresponding null distribution. The only difference
would be the slightly reduced degrees of freedom of the $\chi^{2}$ associated
with the residuals after regressing on the smooth. This is accounted for by
inflating the residuals by a factor of $\frac{n}{n-\nu_{i}}$ where $\nu_{i}$
is the degrees of freedom associated with each smoothed time series. We then
resample the resulting residuals recompute the test statistic. Resampling can
be done in blocks if there is temporal structure in the residuals.

From the newly created time series, we can obtain new MAF factors. Through the
test statistic, we see whether the first MAF factor remain after the noise
estimate has been shuffled, an operation that should preserve the MAF factor
if it indeed represents a signal.

Algorithm 3 provides the details to the hypothesis testing procedure.

Algorithm 3 Calculate number of MAF factors, $\bm{Y}\in\mathbb{R}^{n\times p}$

0: $\bm{Z}\in\mathbb{R}^{n\times p}$

1: For each time series $Z_{i}(t)$ for $i=1,...,p$, calculate
$SNR_{empir}(Z_{i}(t))$.

2: Create a smooth time series from each original time series $Z_{i}(t)$ and
calculate the residual $\hat{\varepsilon}_{i}(t)=Z_{i}(t)-\tilde{Z}_{i}(t)$.

3: From the set of integers $[1,n]$, draw a sequence $\\{s_{i}\\}$ of $n$
integers with replacement, yielding the sequence
${s_{1}},{s_{2}},\dots,{s_{n}}$. The new set of residuals becomes
$\hat{\varepsilon}_{i}({s_{1}}),\hat{\varepsilon}_{i}({s_{2}}),\dots,\hat{\varepsilon}_{i}({s_{n}})$
which we shall call $\hat{\varepsilon}_{i}^{\ast}(t)$.

4: Let
$Z_{i}^{\ast}(t)=\sqrt{\frac{n}{n-\nu_{i}}}\hat{\varepsilon}_{i}^{\ast}(t)$.

5: Calculate new MAFs from the resampled data,
$\bm{Y}^{\ast}(t)=\bm{Z}^{\ast}\bm{w}_{MAF}(\bm{Z}^{\ast})$, and their
associated $SNR_{empir}(\bm{Y}_{i}^{\ast}(t))$ for $i=1,...,p$.

6: Repeat steps 3-5 B times.

7: Calculate MAF1’s associated p-value,  $p=\frac{\\{\\#\text{ of
}SNR_{empir}(\bm{Y}_{1}^{\ast}(t))>SNR_{empir}(\bm{Y}_{1}(t))\\}}{B}$ (5.7)

The empirical SNR is used as test statistic since it’s model counterpart
maximizes the expected likelihood under the mode described in Equation 4 and
proved in Lemma 3. Sample autocorrelation was also explored as a test
statistic but was found to be less powerful than empirical SNR.

In resampling the residuals one can sample with or without replacement, where
the former is referred to as the bootstrap. Permuting the residuals, i.e.
resampling without replacement allows for exact type 1 error control because
we sample from the population as opposed to an estimate of the population
which is the case for the bootstrap. Furthermore, the validity of the
bootstrap depends on the empirical distribution’s asymptotic convergence to
the population distribution, but the permutation test does not have this
requirement.

However, if there is autocorrelation present in the residuals, one would
normally the data in blocks to account for the temporal structure. In this
case, permuting the data is less suitable due to the smaller number of
permutations possible. Sampling with replacement does not have a reduction in
the number of possible combinations and might thus be there method of choice.

This hypothesis testing procedure can be extended to multiple signals. Because
MAF solves an eigenvalue/eigenvector problem the MAF factors are orthogonal.
As a corollary, the second MAF maximizes autocorrelation on a dataset that
lies in the space perpendicular to the first MAF. Similarly, the third lies in
the space perpendicular to the first two MAFs. In this vein, each MAF will
produce a signal estimate orthogonal to the other MAFs. In the hypothesis
testing framework, one would test whether each signal estimate is significant
or not. Our method for creating the null distribution outlined in the single-
signal case would still be valid. The only difference would be that multiple
signals are subtracted out of the dataset followed by a permutation of the
residuals. To reflect the multiple signal extension in Algorithm 3, one would
only replace Step 7 by the calculation of the empirical SNR and p-values of
MAFs 1 through $k$ for $k=1,...,p$.

We continue using the two examples presented in Figure 5. The top panels show
three time series where the signal strength vector $\bm{b}_{A}=(0.4,0.2,0.1)$,
while the lower panels contain a weaker signal strength,
$\bm{b}_{B}=(0.4,0.2,0.1)/2$. Furthermore, the smooth line in each panel is
the underlying signal before the noise is added, while the gray lines show the
raw observations used to calculate the MAF transformation.

The MAF SNR distributions are shown in Figure 7, where the solid vertical line
is the SNR of the original observations while the histograms represent the SNR
of 1000 sets of permuted observations. The $p$-value represents the
probability of the observed MAF SNR under the null hypothesis, i.e. the
absence of a signal. In the strong-signal case the $p$-value is 0,
corresponding to the left panel, while the weak-signal case has a $p$-value of
0.893, corresponding to the right panel of Figure 7.

Figure 7: The distributions of the maximized SNR under the null hypothesis,
where the original data has been resampled by permuting the time steps. The
black vertical line shows the original SNR, with the strong-signal example on
the left and the weak-signal example on the right.

We proceed by calculating the power of the test under various signal
strengths. To calculate the power we do the following,

  1. 1.

Simulate B instances from the null model with no signal and calculate the
associated test statistics.

  2. 2.

Find the $(1-\alpha)^{th}$ quantile of the null distribution and call it
$T_{1-\alpha}$.

  3. 3.

Simulate B instances of the alternative and calculate the associated test
statistics.

  4. 4.

Find the area under the curve of the alternative distribution for which the
test statistics are greater than $T_{1-\alpha}$. This area is the power.

  5. 5.

Repeat steps 3-4 for different signal strengths.

A plot of the power as a function of the power for a number of signal
strengths is shown in Figure 8, with $B=5000$, and $\rho=0.5$. The x-axis
represents the coefficient by which the base signal vector,
$\bm{b}=(0.8,0.4,0.2)$, is multiplied. Both SNR and autocorrelation are used
as test statistics and shown in separate panels. Note that SNR has higher
power than autocorrelation.

Figure 8: The power of the hypothesis test as a function of the signal-
multiplication factor. The base signal strength vector $\bm{b}=(0.8,0.4,0.2)$
is multiplied by $c\in[0,1]$.

##  6 Real Data: Application to Tree rings

To illustrate the efficacy of the MAF methodology, an application using the
tree ring data from Western United States is presented in this section. First,
we extract the MAFs to obtain an estimate of the underlying signal(s) present
in the data. Uncertainty of the MAF factors is then estimated. Lastly, we test
for the significance of the underlying signals through the hypothesis testing
framework introduced in previous section.

Overview of the data: The data is obtained from the Mann et al. (2008) and
quantifies the annual growth of tree rings. It has been pre-processed as
described in Mann _et al_ ’s Supplemental Section. We selected 21 concurrent
tree ring time series for the period 1850-1999, of which 4 were already shown
in Figure 1. Figure 9 shows all 21 time series, scaled, centered, and
annotated by their names333The raw data was download from the Supplemental
section from Mann at
http://www.meteo.psu.edu/holocene/public_html/supplements/MultiproxyMeans07/.
Some time series show more temporal coherence than others. The goal here is to
extract a common underlying temporal signal.

Figure 9: 21 tree ring time series from the W77estern United States, centered
and scaled.

MAF estimation: The 3 first MAFs and PCs are shown in Figure 10 where each
time series is annotated by its sample autocorrelation. A smooth version of
each time series is shown in bold with 30 years per knot starting at the last
year. Note that MAF produces time series that are more autocorrelated than PCA
and are sorted in decreasing autocorrelation.

Figure 10: Centered and scaled MAFs and PCs of the tree rings with smoothed
equivalents shown in bold.

Uncertainty quantification: Quantifying uncertainty of the estimated MAF
factors can be obtained through the methodology outlined in subsection 4.2. A
plot of this is shown in Figure 11, where 1000 resampled datasets are created
by doing a block bootstrap with a block size of 5 years.

Each new MAF is created by using normalized MAF coefficients. The smoother
applied is the same LOESS smoother as in Section 5.2. We see a clear signal
present in the first two MAFs with the confidence bands containing the
original smooth MAFs. However, the third MAF time series’ (MAF3) original
estimate can be seen almost outside the confidence interval. This suggests
that MAF3 is mainly composed of noise such that when the tree ring data is
resampled and the MAF is recalculated the trend associated with MAF3
disappears.

Signal concentration: The signal information also seems to be more focused in
fewer MAFs compared with PC. This can be illustrated by obtaining the
canonical correlation between the data set time series and the smooth MAFs and
PCs. The first canonical correlation gives the linear combination of the data
time series most correlated with a linear combination of MAFs/PCs. This can be
interpreted as the correlation with a potential underlying signal.

Figure 12 shows a sampling distribution of the canonical correlations for each
data set previously obtained through resampling and the associated smooth MAFs
and PCs. We see that MAF has a consistently higher canonical correlation until
we include three components, at which point the two methods equalize. The
result is even more pronounced using unsmoothed MAFs/PCs. Notice also that the
MAF distributions are narrower than those of PCA. This is due to the smaller
uncertainty about the MAF factor estimates shown in Figure 11.

Figure 11: A set of 1000 resampled MAF factors using the block bootstrap. The
thick green line shows the original MAF with a smoother applied, while the
grey lines are the un-smoothed resampled MAFs. Dashed lines show the $95^{th}$
confidence bands.  Figure 12: The first canonical correlation between the
resampled datasets and the associated MAF and PCA components 1 to $k$ for
$k=1,2,3$.

##  7 Discussion

We demonstrated advantages of the MAF optimization criterion in comparison
with PCA for the purpose of extracting a common time trend component from
multiple concurrent time series. In particular, under a model where each time
series is a combination of the underlying time trend with additive noise, we
showed that the MAF-optimized linear combination of time series, i.e.,
maximizing autocorrelation, also maximizes the signal-to-noise ratio among all
possible linear combinations. The sub-optimality of PCA can become worse as
the number of available time series grows, as the cross-correlation between
time series increases, and as the noise levels increase. We also investigated
some sampling properties of the MAF analysis and showed through simulations
that the MAF-optimized combined time series can be statistically more stable
than the corresponding PCA-optimized time series obtained from the same set of
concurrent time series data.

We generalize the signal-plus-noise model to include $q$ multiple underlying
signal time series embedded in $p\geq q$ time series. Considering a response
variable comprised of linear combinations of these predictive signals, we
showed that the first $q$ MAFs span the same space as the first $q$ CCFs. And
since CCA by definition maximizes correlation with the signals, the
corresponding $q$ dimensional subspace spanned by the first $q$ CCFs is
optimal when regressing the signal onto these. So, by transitivity, the first
$q$ MAFs then also contain the optimal subspace for replicating the response.
The advantage of MAF is that knowledge of the shape of the signal is not
necessary. So, MAF compresses a $p$-dimensional time series into a $q$
dimensional data set without losing any information.

Lastly, we illustrated some initial applications of MAF applied to combining
21 concurrent annual tree ring time series for a region in western North
America, covering the period 1850-1999, with the goal of extracting common
time trend information. Regional tree ring time series data are believed to be
imperfect proxies for regional weather time series, such as average annual
temperature, and in a subsequent paper we are investigating the calibration
between regional temperature time series and regional tree ring proxy time
series for these and other regions of the globe. An important step in the
calibration is the extraction of common time trend information from the proxy
data.

## References

  * Arbenz and Golub [1988] Peter Arbenz and Gene H Golub.  On the spectral decomposition of hermitian matrices modified by low rank perturbations with applications.  _SIAM Journal on Matrix Analysis and Applications_ , 9(1):40–58, 1988.
  * Briffa et al. [2008] Keith R Briffa, Vladimir V Shishov, Thomas M Melvin, Eugene A Vaganov, Håken Grudd, Rashit M Hantemirov, Matti Eronen, and Muktar M Naurzbaev.  Trends in recent temperature and radial tree growth spanning 2000 years across northwest eurasia.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ , 363(1501):2269–2282, 2008.  doi: 10.1098/rstb.2007.2199.
  * Bunch et al. [1978] R. Bunch, James, P. Nielsen, Christopher, and C. Sorensen, Danny.  Rank-one modification of the symmetric eigenproblem.  _Numerische Mathematik_ , 31(1):31–48–, 1978\.  ISSN 0029-599X.  URL http://dx.doi.org/10.1007/BF01396012.
  * Cleveland [1979] W.S. Cleveland.  Robust locally weighted regression and smoothing scatterplots.  _J. Amer. Statist. Assoc._ , 74(368):829–836, 1979.
  * Doob [1953] Joseph L Doob.  _Stochastic processes_ , volume 101.  New York Wiley, 1953.
  * Durrett [2010] Rick Durrett.  _Probability: theory and examples_.  Cambridge university press, 2010.
  * Gallagher et al. [2014] Neal B. Gallagher, Jeremy M. Shaver, Randall Bishop, Robert T. Roginski, and Barry M. Wise.  Decompositions using maximum signal factors.  _J. Chemometrics_ , 28(8):663–671, August 2014\.  ISSN 1099-128X.  URL http://dx.doi.org/10.1002/cem.2634.
  * Hastie et al. [2009] T. Hastie, R. Tibshirani, and J. Friedman.  _Elements of Statistical Learning_.  Springer, 2nd edition, 2009.
  * Jansen and Rajaratnam [2014] L. Jansen and B. Rajaratnam.  Robust reconstructions with temporal dependencies.  _Journal of the American Statistical Association (in press)_ , 2014\.
  * Li et al. [2007] B. Li, D. W. Nychka, and C.M. Ammann.  The ’hockey stick’ and the 1990s: a statistical perspective on reconstructing hemispheric teperatures.  _Tellus A_ , 59(5):591–598, 2007.
  * Mann et al. [2008] Michael E. Mann, Zhihua Zhang, Malcolm K. Hughes, Raymond S. Bradley, Sonya K. Miller, Scott Rutherford, and Fenbiao Ni.  Proxy-based reconstructions of hemispheric and global surface temperature variations over the past two millennia.  _Proceedings of the National Academy of Sciences_ , 2008.  doi: 10.1073/pnas.0805721105.
  * McShane and Wyner [2011] B.B. McShane and J. Wyner.  A statistical analysis of multiple temperature proxies: Are reconstructions of surface temperatures over the last 1000 years reliable?  _The Annals of Applied Statistics_ , 2011.
  * Muirhead [2005] R.J. Muirhead.  _Aspects of Multivariate Statistical Theory_.  Wiley, 2005.
  * Shapiro and Switzer [1989] D.E. Shapiro and P. Switzer.  Extracting time trends from multiple monitoring sites.  Technical report, Stanford University, 1989.
  * Switzer and Green [1984] P. Switzer and A. A. Green.  Min/max autocorrelation factors for multivariate spatial imagery.  Technical report, Stanford University, 1984.

## Supplemental section

##  Appendix A Getting general MAF coefficients

We present an alternative method for deriving the MAF coefficients under the
general model given in Equation 3.7. To do this, we first develop the case
where all input time series have the same noise level we get the following
covariance for $\bm{Z}(t)$,

$\Sigma_{Z}=\bm{b}\bm{b}^{\prime}+\rho\bm{1}_{p}\bm{1}_{p}^{\prime}+(1-\rho)\bm{I}$
(A.1)

where $\rho$ is the common cross-correlation across all the time series. The
lagged covariance structure is given in 2.11. The PCs are given by the
eigenvectors of $\bm{\Sigma}_{Z}$, whereas the MAFs are the eigenvectors of
$\bm{\Sigma}_{Z}^{-1/2}\bm{\Sigma}_{\Delta Z}\bm{\Sigma}_{Z}^{-1/2}$.

First consider the special case where $\bar{b}=\sum_{i=1}^{p}b_{i}=0$. This
implies that $\bm{b}^{\prime}\bm{1}_{p}=0$ and the eigenvectors are $\bm{b}$,
$\bm{1}_{p}$ and all the vectors perpendicular to these two. The corresponding
eigenvalues are $(|\bm{b}|^{2},\rho p,0)+1-\rho$ where the last eigenvalue is
repeated $p-2$ times. If $|\bm{b}|^{2}>\rho p$, PC1 will be $\bm{b}$. PC2 will
then be $\bm{1}_{p}$, unless $\rho<0$ in which case PC2 will be in the
aforementioned nullspace. Lastly, if $|\bm{b}|^{2}<\rho p$, PC1 will be
$\bm{1}_{p}$.

Another special case if where $\rho=0$. Here the highest eigenvalue and
corresponding eigenvector will be proportional to $\bm{b}$ while all the
others will be perpendicular to $\bm{b}$ for both MAF and PCA.

In the general case where $\bar{b}\neq 0$ and $\rho\neq 0$, notice that
$\bm{b}\bm{b}^{\prime}+\rho\bm{1}_{p}\bm{1}_{p}^{\prime}$ is rank 2. So the
dimensionality of that nullspace is $p-2$. All vectors in this nullspace will
have an eigenvalue of $1-\rho$. The remaining two eigenvectors are found by
assuming a general structure of the eigenvectors,
$\bm{v}=a_{1}\bm{1}_{p}+a_{2}\bm{b}$. It then follows that the
eigenvectors/eigenvalues of $\Sigma_{Z}$ are given by

$\displaystyle v_{1}=$ $\displaystyle\left(\frac{\rho
p-|\bm{b}|^{2}+\Delta}{2\bar{\bm{b}}p}\right)\bm{1}_{p}+\bm{b},\quad$
$\displaystyle\lambda_{1}=\frac{\rho p+|\bm{b}|^{2}+\Delta}{2}+1-\rho$
$\displaystyle v_{2}=$ $\displaystyle\left(\frac{\rho
p-|\bm{b}|^{2}-\Delta}{2\bar{\bm{b}}p}\right)\bm{1}_{p}+\bm{b},\quad$
$\displaystyle\lambda_{2}=\frac{\rho p+|\bm{b}|^{2}-\Delta}{2}+1-\rho$ (A.2)

where

$\displaystyle\Delta^{2}=$ $\displaystyle(|\bm{b}|^{2}-\rho
p)^{2}+4(\bar{\bm{b}}p)^{2}\rho,$ (A.3)

where $|\bm{b}|^{2}$ is the squared sum of the SNRs of the input time series.

The vectors in the nullspace have mean equal to zero. This can be seen by
considering any eigenvector
$\bm{v}\in\mathcal{N}(\bm{b}\bm{b}^{\prime}+\rho\bm{1}_{p}\bm{1}_{p}^{\prime})$,

$\displaystyle(\bm{b}\bm{b}^{\prime}+\rho\bm{1}_{p}\bm{1}_{p}^{\prime})\bm{v}=0$
$\displaystyle(\bm{b}^{\prime}\bm{v})\bm{b}+\rho\bar{v}p\bm{1}_{p}=0.$ (A.4)

Because $\bm{b}$ is in general not equal to $\bm{1}_{p}$, we have
$\bm{b}^{\prime}\bm{v}=\bar{v}=0$

To get the eigenvectors corresponding to the MAFs, consider also the lagged
covariance matrix,

$\displaystyle\Sigma_{\Delta Z}=$ $\displaystyle
k_{f}\bm{b}\bm{b}^{\prime}+k_{\varepsilon}(\rho\bm{1}_{p}\bm{1}_{p}^{\prime}+(1-\rho)\bm{I})$
(A.5)

Letting
$\bm{\Sigma}_{Z}^{-1}=\bm{\Gamma}\bm{D}^{-1}\bm{\Gamma}^{\prime}=\bm{H}\bm{H}^{\prime}$,
where $\bm{H}=\bm{\Gamma}\bm{D}^{-1/2}$. Furthermore, because the optimal SNR
in Equation 2.2 does not depend on $k_{f}$ as long as $k_{f}<k_{\varepsilon}$,
we can set $k_{f}=0$, w.l.o.g. This is because the optimal SNR coefficients
for each time series is equivalent to the MAF1 coefficients, which is the
eigenvector corresponding to the smallest eigenvalue of

$\tilde{\bm{\Sigma}}_{\Delta}=\bm{H}^{\prime}\left(\rho\bm{1}_{p}\bm{1}_{p}^{\prime}+(1-\rho)\bm{I}\right)\bm{H},$
(A.6)

where the coordinate system has been rotated such that
$\bm{\Sigma}_{Z}=\bm{I}$. The MAF coefficients will change but the resulting
MAF factors will not under this rotation, as shown in Lemma 2.

Now, let $\bar{\bm{u}}_{i}$ be the mean of the normalized version of the
vectors in Equation A, $\lambda_{i}$ be the corresponding eigenvectors. By
letting $\Lambda$ be the diagonal matrix with
$\frac{(1-\rho)}{\sqrt{\lambda_{i}}}$ along the diagonal and $c_{i}=\frac{\rho
p\bar{\bm{u}}_{i}}{\sqrt{\lambda_{i}}}$, we can recast this in a more familiar
form,

$(\tilde{\Sigma}_{\Delta})_{ij}=\bm{\Lambda}+\bm{c}\bm{c}^{\prime}.$ (A.7)

A closer look at this matrix will reveal that $c_{i}=0,\,\forall i>2$.
Furthermore, $\bm{\Lambda}_{ii}=1,\,\forall i>2$. This means that we can
decompose the matrix as follows,

$\begin{bmatrix}\bm{A}&\bm{0}_{2\times(p-2)}\\\ \bm{0}_{(p-2)\times
2}&\bm{1}_{(p-2)\times(p-2)}.\end{bmatrix}$ (A.8)

And because $\bm{A}$ is symmetric and $2\times 2$, its
eigenvalues/eigenvectors can be found in closed form. The remaining
eigenvectors can be made the standard basis vectors
$\bm{e}_{i}=(0_{1},...,0_{i-1},1,0_{i+1},...,0_{p}),\,\forall i>2$. In
particular, by solving

$\bm{A}\begin{bmatrix}x\\\
y\end{bmatrix}=\begin{bmatrix}\lambda^{-1}_{1}\left(1-\rho+\rho
p^{2}\bar{\bm{u}}^{2}_{1}\right)&\frac{\rho
p^{2}}{\sqrt{\lambda_{1}\lambda_{2}}}\bar{\bm{u}}_{1}\bar{\bm{u}}_{2}\\\
\frac{\rho
p^{2}}{\sqrt{\lambda_{1}\lambda_{2}}}\bar{\bm{u}}_{1}\bar{\bm{u}}_{2}&\lambda^{-1}_{2}\left(1-\rho+\rho
p^{2}\bar{\bm{u}}^{2}_{2}\right)\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}=\begin{bmatrix}a&b\\\ b&d\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}=\mu\begin{bmatrix}x\\\ y\end{bmatrix},$ (A.9)

we find that the eigenvalues/eigenvectors are

$\begin{bmatrix}x\\\
y\end{bmatrix}=\frac{1}{\sqrt{b^{2}+(\mu-d)^{2}}}\begin{bmatrix}\mu-d\\\
b\end{bmatrix},\quad\mu_{1},\mu_{2}=\frac{a+d\pm\sqrt{(a-d)^{2}+4b^{2}}}{2},$
(A.10)

where we are interested in the smallest eigenvalue, i.e. where we subtract the
term involving the discriminant.

Now, let $\tilde{\bm{w}}$ be the full vector in $\mathbb{R}^{p}$ with zeroes
everywhere except in the first two entries which take the values $x$ and $y$.
To get the values of each coefficient in the basis of the original time
series, we do an inverse transformation,

$\bm{w}=(\bm{H}^{\prime})^{-1}\tilde{\bm{w}}=\bm{H}\tilde{\bm{w}}=\begin{bmatrix}\bm{h}_{1}&\bm{h}_{2}&\ldots&\bm{h}_{p}\end{bmatrix}\begin{bmatrix}x\\\
y\\\ 0\\\ \vdots\\\ 0\end{bmatrix}=x\bm{h}_{1}+y\bm{h}_{2}.$ (A.11)

Note that this expression is a linear combination of PC1 and PC2. Furthermore,
the largest eigenvalue of the two in Equation A.10 is equal to 1, just like
the other degenerate eigenvalues. This leaves only one non-degenerate
eigenvalue, which can be interpreted as there being only one signal present in
different strenths.

The result in Equation A.11 can be used to obtain MAF1 in a more general
setting, where the noise is of unequal variance,

$\bm{z}_{t}=\bm{b}f_{t}+\bm{\varepsilon}_{t},\>\>\text{with}\>\>(\varepsilon_{t})_{i}\sim(0,\sigma_{i}^{2})\quad\forall
t,$ (A.12)

by taking advantage of the fact that MAF is preserved under linear
transformations. We can write the modified covariance matrix as

$\bm{A}^{\prime}\bm{\Sigma}_{Z}\bm{A}=\bm{A}^{\prime}\left(\tilde{\bm{b}}\tilde{\bm{b}}^{\prime}+\rho\bm{1}_{p}\bm{1}_{p}^{\prime}+(1-\rho)\bm{I}\right)\bm{A},$
(A.13)

where $\tilde{\bm{b}}=\bm{b}/{\bm{\sigma}}$ and $\bm{\sigma}^{2}$ is the
vector is noise variances. Similarly for the lagged covariance matrix.

It then follows that the eigenvalue equation to be solved is

$\bm{\Sigma}_{Z}^{-1/2}\bm{\Sigma}_{\Delta}\bm{\Sigma}_{Z}^{-1/2}\tilde{\bm{w}}=\lambda\tilde{\bm{w}},$
(A.14)

where $\tilde{\bm{w}}=\bm{A}\bm{w}$. But $\tilde{\bm{w}}$ is already given in
(A.11), and thus $\bm{w}=\bm{A}^{-1}\tilde{\bm{w}}$, which are the MAF1
coefficients in the original coordinate system.

This this general case with unequal noise variance, the MAF will not be a
linear combination of PC1 and PC2. The reason for this is that the vectors in
the nullspace of the new covariance matrix’ rank-2 update will not in general
be an eigenvector of the full covariance matrix, with the unequal variance
terms in the diagonal. This means that the eigenvalue problem cannot be
rewritten in a form similar to the one given in (A.8). In fact, the PCA
eigenvectors do not even exist in closed form, but must be obtained by solving
a determinant equation for the eigenvalues. This problem is explored in Arbenz
and Golub [1988], Bunch et al. [1978] and the references therein.

##  Appendix B Proofs

###### Proof of Lemma 4.

Under $H_{0}$ and the Gaussian assumption, the log-likelihood

$\ln(L_{0})=-\frac{1}{2}\ln(|\bm{\Sigma}_{\varepsilon}|)-\frac{1}{2}\sum_{t=1}^{n}\bm{z}(t)^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{z}(t)-\frac{p}{2}\ln(2\pi),$
(B.1)

where we assume that the mean of $\bm{z}_{t}$ is zero without loss of
generality.

Now, the likelihood under the alternative hypothesis

$\ln(L_{A})=-\frac{1}{2}\ln(|\bm{\Sigma}_{\varepsilon}|)-\frac{1}{2}\sum_{t=1}^{n}(\bm{z}(t)-\bm{b}f(t))^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}(\bm{z}(t)-\bm{b}f(t))-\frac{p}{2}\ln(2\pi).$
(B.2)

The likelihood ratio is

$\displaystyle\ln(L_{A})-\ln(L_{0})=$
$\displaystyle-\frac{1}{2}\sum_{t=1}^{n}f^{2}(t)\bm{b}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}+\sum_{t=1}^{n}f(t)\bm{z}(t)^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\bm{b}^{\prime}\Sigma_{\varepsilon}^{-1}\bm{b}+\sum_{t=1}^{n}f(t)\bm{z}(t)^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}$
(B.3)

where the last equality holds because of the unit squared sum of $f(t)$. Now,
if we substitute $\bm{z}(t)$ for the alternative model and taking expectations
under either model, we get

$\displaystyle E[\ln(L_{A})-\ln(L_{0})]=$
$\displaystyle-\frac{1}{2}\bm{b}^{\prime}\Sigma_{\varepsilon}^{-1}\bm{b}+E\left[\sum_{t=1}^{n}f(t)(\bm{b}f(t)-\bm{\varepsilon}(t))^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\bm{b}^{\prime}\Sigma_{\varepsilon}^{-1}\bm{b},$
(B.4)

where the term involving $\bm{\varepsilon}$ is zero after taking the
expectation.

Now, we already proved that MAF1 maximizes the model signal-to-noise ratio,

$SNR=\frac{(\bm{w}^{\prime}\bm{b})^{2}}{\bm{w}^{\prime}\bm{\Sigma}_{\varepsilon}\bm{w}}.$
(B.5)

Making the change of coordinates $\bm{u}=\Sigma_{\varepsilon}^{1/2}\bm{w}$,
using the familiar spectral decomposition for the square root, gives the SNR
representation

$\frac{(\bm{u}^{\prime}\bm{\Sigma}_{\varepsilon}^{-1/2}\bm{b})^{2}}{\bm{u}^{\prime}\bm{u}}.$
(B.6)

Using Cauchy-Schwartz theorem, the normalized vector, $\bm{u}$, which
maximizes SNR is parallel to $\bm{\Sigma}_{\varepsilon}^{-1/2}\bm{b}$. And so
in the original coordinate system,

$\bm{w}_{MAF}=\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}$ (B.7)

.

Substituting the expression for these MAF1 coefficients in our definition for
SNR gives

$\displaystyle SNR_{optimal}=$
$\displaystyle\frac{(\bm{b}^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b})^{2}}{\bm{b}^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{\Sigma}_{\varepsilon}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}}$
$\displaystyle=$
$\displaystyle\bm{b}^{\prime}\bm{\Sigma}_{\varepsilon}^{-1}\bm{b},$ (B.8)

which is proportional to the expected likelihood ratio test statistic. ∎

######  Proposition 2.

The MAF-transformation matrix,
$\bm{W}_{MAF}=[\bm{w}_{1},\bm{w}_{2},...,\bm{w}_{p}]$ , contains the
eigenvectors of $\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}$.

###### Proof of Proposition 2.

By definition, the MAF1 vector of coefficients, $\bm{w}_{1}$, minimizes

$\frac{\bm{w}_{1}^{\prime}\bm{S}_{\Delta}\bm{w}_{1}}{\bm{w}_{1}^{\prime}\bm{S}\bm{w}_{1}}.$
(B.9)

By letting $\bm{S}^{1/2}\bm{w}_{1}=\bm{u}_{1}$ we get

$\frac{\bm{u}_{1}^{\prime}\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}\bm{u}_{1}}{\bm{u}_{1}^{\prime}\bm{u}_{1}},$
(B.10)

which is equivalent to minimizing

$\displaystyle\bm{u}_{1}^{\prime}\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}\bm{u}_{1}$
subject to $\displaystyle\bm{u}_{1}^{\prime}\bm{u}_{1}=1.$ (B.11)

Following Muirhead [2005], the minimizing vector is the eigenvector with the
lowest eigenvalue. Furthermore, the eigenvector with the $k^{th}$ smallest
eigenvalue minimizes

$\displaystyle\bm{u}_{k}^{\prime}\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}\bm{u}_{k}$
subject to $\displaystyle\bm{u}_{k}^{\prime}\bm{u}_{k}=1$ and
$\displaystyle\bm{u}_{k}^{\prime}\bm{u}_{i}=0\quad\forall i<k.$ (B.12)

$\bm{u}_{k}$ corresponds to MAF$k$, after MAF$i,\>\forall i<k$, has been
projected out of the data. The linear transformation
$\bm{w}_{k}=\bm{S}^{-1/2}\bm{u}_{k}$ gives the eigenvectors in the original
coordinate system. ∎

######  Theorem 2 (Simplified).

Consider a set of time series $Z_{n}(t)\in\mathbb{R}^{p}$, such that

$\displaystyle\bm{Z}_{n}(t)=$ $\displaystyle
f_{n}(t)\bm{b}+\bm{\varepsilon}_{n}(t)\quad t=1,\dots,n$
$\displaystyle\bm{\Delta}\bm{Z}_{n}(t)=$
$\displaystyle\bm{Z}_{n}(t)-\bm{Z}_{n}(t+1)=\Delta
f_{n}(t)\bm{b}+\Delta\bm{\varepsilon}_{n}(t)$ (B.13)

with $f_{n}(t)\in\mathbb{R}\>\>\forall
t=1,2,...,n,\bm{b},\bm{\varepsilon}_{n}(t)\in\mathbb{R}^{p}\>\>\forall
t=1,2,...,n$,
$\Delta\bm{\varepsilon}_{n}=\bm{\varepsilon}_{n}(t)-\bm{\varepsilon}_{n}(t+1)$,
and $\Delta f_{n}(t)=f_{n}(t)-f_{n}(t+1)$. Residual time series
$\bm{\varepsilon}_{n}$ is a weakly stationary $p$-variate time series and the
associated autocovariance is absolutely summable. Then,

$\displaystyle\bm{S}_{n}\overset{p}{\to}$ $\displaystyle\bm{\Sigma}\text{ as
}n\to\infty,$
$\displaystyle\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}\overset{p}{\to}$
$\displaystyle\bm{\Sigma}^{-1/2}\bm{\Sigma}_{\Delta}\bm{\Sigma}^{-1/2}\text{
as }n\to\infty$

###### Proof of Theorem 2.

Three parts make up this proof: 1. Stationarity of differenced time series, 2.
Convergence in probability of $\bm{S}_{n}$ and $\bm{S}_{\Delta n}$, the sample
cross-correlation and lagged cross-correlation to their model counterparts.
3\. Consistency of MAF and PCA coefficients to their model counterparts.

We begin by recognizing that since $\bm{\varepsilon}_{n}$ is a weakly
stationary $p$-variate time series, we have

$\displaystyle E[\bm{\varepsilon}_{n}(t)]=$ $\displaystyle\bm{0}$
$\displaystyle Cov[\bm{\varepsilon}_{n}(t)]=$
$\displaystyle\bm{\Sigma}_{\varepsilon}$ $\displaystyle
Cov[\Delta\bm{\varepsilon}_{n}(t)]=$
$\displaystyle\bm{\Sigma}_{\Delta\varepsilon},$ (B.15)

with the assumption of lagged summability,

$\sum_{\tau=0}^{\infty}|\bm{\gamma}_{\tau,i}|<\infty.$ (B.16)

where the lagged autocovariance for noise component $i$ is

$\bm{\gamma}_{\tau,i}=Cov[\varepsilon_{n,i}(t),\varepsilon_{n,i}(t+\tau)].$
(B.17)

Furthermore, let

$\displaystyle\bm{S}_{n}=$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}[\bm{Z}_{n}(t)-\overline{\bm{Z}_{n}(t)}][\bm{Z}_{n}(t)-\overline{\bm{Z}_{n}(t)}]^{\prime}$
$\displaystyle\bm{S}_{\Delta n}=$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\bm{[}\Delta\bm{Z}_{n}(t)-\overline{\Delta\bm{Z}_{n}(t)}][\Delta\bm{Z}_{n}(t)-\overline{\Delta\bm{Z}_{n}(t)}]^{\prime},$
(B.18)

where $\overline{\bm{Z}_{n}(t)}=\sum_{i=1}^{n}\bm{Z}_{n}(t)$, similarly for
$\Delta\bm{Z}_{n}(t)$.

Part I: Stationarity of differenced time series:
We now show that
$\Delta\bm{\varepsilon}_{n}(t)-\overline{\Delta\bm{\varepsilon}_{n}}$ is a
zero mean weakly stationary time series, using the weak stationarity of
$\bm{\varepsilon}_{n}(t)$. Let

$Cov[\bm{\varepsilon}(t),\bm{\varepsilon}(t+h)]=\bm{\gamma}(h),$ (B.19)

which is by definition not a function of $t$. Then

$\displaystyle Cov[\Delta\bm{\varepsilon}(t),\Delta\bm{\varepsilon}(t+h)]=$
$\displaystyle
Cov[\bm{\varepsilon}(t),\bm{\varepsilon}(t+h)]-Cov[\bm{\varepsilon}(t+1),\bm{\varepsilon}(t+h)]+$
$\displaystyle\quad
Cov[\bm{\varepsilon}(t+1),\bm{\varepsilon}(t+h+1)]-Cov[\bm{\varepsilon}(t),\bm{\varepsilon}(t+h+1)]$
$\displaystyle=$ $\displaystyle
2\bm{\gamma}(h)-\bm{\gamma}(h+1)-\bm{\gamma}(h-1)\quad\text{for }h>0,$ (B.20)

is not a function of $t$ and is thus also weakly stationary. It is trivial to
show that the differenced time series has zero mean.

Part II: Consistency of $\bm{S}_{n}$ and $\bm{S}_{\Delta n}$
Substituting Equation 2 in Equation 2.15, we get

$\displaystyle\bm{S}_{n}=$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}[(f_{n}(t)-\overline{f_{n}})\bm{b}+\bm{\varepsilon}_{n}(t)-\overline{\bm{\varepsilon}_{n}}][(f_{n}(t)-\overline{f_{n}})\bm{b}+\bm{\varepsilon}_{n}(t)-\overline{\bm{\varepsilon}_{n}}]^{\prime}$
$\displaystyle=$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}[\bm{b}\bm{b}^{\prime}f_{n}^{2}(t)+2f_{n}(t)\bm{b}(\bm{\varepsilon}_{n}(t)-\overline{\bm{\varepsilon}_{n}})^{\prime}+(\bm{\varepsilon}_{n}(t)-\overline{\bm{\varepsilon}_{n}})(\bm{\varepsilon}_{n}(t)-\overline{\bm{\varepsilon}_{n}})^{\prime}].$
(B.21)

As $n\to\infty$, The first term equals $\bm{b}\bm{b}^{\prime}$ by definition,
the second term goes to zero in probability because the $p$-vector of the
time-averaged residuals goes to the zero vector in probability, and the third
term goes to $\bm{\Sigma}_{\varepsilon}$ because $\bm{\varepsilon}_{n}(t)$ is
a weakly stationary time series with zero mean [Doob, 1953, Durrett, 2010].

Thus,

$\bm{S}_{n}=\frac{1}{n}\sum_{t=1}^{n}[\bm{Z}_{n}(t)-\overline{\bm{Z}_{n}(t)}][\bm{Z}_{n}(t)-\overline{\bm{Z}_{n}(t)}]^{\prime}\overset{p}{\to}\bm{b}\bm{b}^{\prime}+\bm{\Sigma}_{\varepsilon}=E[\bm{S}_{n}]$
(B.22)

Now, consider

$\displaystyle\bm{S}_{\Delta n}=$
$\displaystyle\frac{1}{n}\sum_{t}[\bm{b}\bm{b}^{\prime}(\Delta
f_{n}-\overline{\Delta f_{n}})^{2}(t)+2(\Delta f_{n}(t)-\overline{\Delta
f_{n}})\bm{b}(\Delta\bm{\varepsilon}_{n}(t)-\overline{\Delta\bm{\varepsilon}_{n}})^{\prime}(t)+(\Delta\bm{\varepsilon}_{n}(t)-\overline{\Delta\bm{\varepsilon}_{n}})(\Delta\bm{\varepsilon}_{n}(t)-\overline{\Delta\bm{\varepsilon}_{n}})^{\prime}].$
(B.23)

Applying the same arguments and using the weak stationarity of the differenced
time series, we get

$\bm{S}_{\Delta
n}=\frac{1}{n}\sum_{t}[\Delta\bm{Z}_{n}(t)-\overline{\Delta\bm{Z}_{n}(t)}][\Delta\bm{Z}_{n}(t)-\overline{\Delta\bm{Z}_{n}(t)}]^{\prime}\overset{p}{\to}a\bm{b}\bm{b}^{\prime}+\bm{\Sigma}_{\Delta\varepsilon}=E[\bm{S}_{\Delta
n}].$ (B.24)

To wit,

$\displaystyle\bm{S}_{n}\overset{p}{\to}$ $\displaystyle
E[\bm{S}_{n}]=\bm{\Sigma}=\bm{b}\bm{b}^{\prime}+\bm{\Sigma}_{\epsilon}\text{
as }n\to\infty$ $\displaystyle\bm{S}_{\Delta n}\overset{p}{\to}$
$\displaystyle E[\bm{S}_{\Delta
n}]=\bm{\Sigma}_{\Delta}=a\bm{b}\bm{b}^{\prime}+\bm{\Sigma}_{\Delta\epsilon}\text{
as }n\to\infty,$ (B.25)

Part III: Consistency of MAF and PCA coefficients
PCA and MAF coefficients are the eigenvectors of $\bm{S}$ and
$\bm{S}^{-1/2}\bm{S}_{\Delta}\bm{S}^{-1/2}$ respectively, using a spectral
decomposition $\bm{S}^{-1/2}=\bm{H}\bm{L}^{-1/2}\bm{H}^{\prime}$. The
continuous mapping theorem ensures that consistent estimates of the covariance
and lagged covariance matrices implies consistent estimates of MAF and PCA
coefficients, i.e. the coefficients will also converge to their model values,
the eigenvectors of the model covariance matrix, $\bm{\Sigma}$ for PCA and
$\bm{\Sigma}^{-1/2}\Sigma_{\Delta}\bm{\Sigma}^{-1/2}$, with
$\bm{\Sigma}^{-1/2}=\bm{\Gamma}\bm{D}^{-1/2}\bm{\Gamma}^{\prime}$. ∎

###### Proof of Theorem 1.

First we establish that MAF coefficient vectors $1$ to $q$ are linear
combinations of the signal strength vectors $\bm{b}_{i}$. Then we show that
CCA coefficients also has this property. Thus, these two methods have
coefficients that span the same $q$-subspace in $\mathbb{R}^{p}$.

Part I: MAF coefficient vectors
The definition of the first $q$ MAF factors solve the following sequence of
problems,

maximize
$\displaystyle\quad\rho(\bm{a}_{i})=\frac{\bm{a}_{i}^{\prime}\bm{\Sigma}_{\delta
Z}\bm{a}_{i}}{\bm{a}_{i}^{\prime}\bm{\Sigma}_{Z}\bm{a}_{i}}=k_{\varepsilon}\frac{\bm{B}\bm{\Lambda}\bm{B}^{\prime}+\bm{I}}{\bm{B}\bm{B}^{\prime}+\bm{I}}$
subject to $\displaystyle\quad\bm{a}_{i}^{\prime}\bm{a}_{i}=1$
$\displaystyle\quad\bm{a}_{i}^{\prime}\bm{a}_{j}=0,\quad\forall j<i\leq q,$
(B.26)

where we have changed coordinate system such that
$\bm{\Sigma}_{\varepsilon}=\bm{I}$ and where $\Lambda$ is the matrix with
entries $(k_{i})/k_{\varepsilon}=\lambda_{i}$ for $i=1,\dots,q$ in the
diagonal and zero in the off-diagonals. For now assume that the columns of
$\bm{B}$ are linearly independent.

We want to show that the first $q$ maximizing vectors are in the range of
$\bm{B}$.

First, from the spectral theorem

$\bm{B}\bm{B}^{\prime}=\bm{U}\bm{D}\bm{U}^{\prime},$ (B.27)

where $U\in\mathbb{R}^{p\times p}$ whose last $p-q$ columns are perpendicular
to $\bm{b}_{i}$ $\forall i=1,\dots,q$.

So,

$\bm{B}\bm{B}^{\prime}+\bm{I}=\bm{U}(\bm{D}+\bm{I})\bm{U}^{\prime},$ (B.28)

so the eigenvectors are preserved by adding the identity matrix.

Second, note that the first $q$ columns of $\bm{U}$, $\bm{u}_{i}\in R(\bm{B})$
since we can write

$\bm{B}[\bm{B}^{\prime}\bm{U}\bm{D}^{-1}]=\bm{U}.$ (B.29)

Third, note that for any $\bm{x}\in R(\bm{B})$

$0\leq\sum_{i=1}^{q}(\bm{b}_{i}^{\prime}\bm{x})^{2}=\bm{x}\bm{B}\bm{B}^{\prime}\bm{x}\leq\sum_{i=1}^{q}\lambda_{i}(\bm{b}_{i}^{\prime}\bm{x})^{2}=\bm{x}\bm{B}\bm{\Lambda}\bm{B}^{\prime}\bm{x},$
(B.30)

and thus

$\bm{x}^{\prime}(\bm{B}\bm{B}^{\prime}+\bm{I})\bm{x}\leq\bm{x}^{\prime}(\bm{B}\bm{\Lambda}\bm{B}^{\prime}+\bm{I})\bm{x},$
(B.31)

with strict inequality iff $\bm{x}\in R(\bm{B})$.

Now,

$\bm{B}\bm{\Lambda}\bm{B}^{\prime}+\bm{I}=\bm{V}(\bm{E}+\bm{I})\bm{V}^{\prime},$
(B.32)

and since Equation B.31 holds for any vector in the range of $\bm{B}$, any
eigenvector in $\bm{U}$ or $\bm{V}$ with eigenvalue greater than 1 must be in
the range of $\bm{B}$. So,

$\rho(\bm{x})=k_{\varepsilon}\frac{\bm{x}\bm{V}(\bm{E}+\bm{I})\bm{V}^{\prime}\bm{x}}{\bm{x}\bm{U}(\bm{D}+\bm{I})\bm{U}^{\prime}\bm{x}}\geq
k_{\varepsilon}.$ (B.33)

We can thus find $q$ linearly independent vectors in $R(\bm{B})$, call them
$\bm{x}_{B}$ such that $\rho(\bm{x}_{B})>k_{\varepsilon}$. Any vectors in the
null space of $B$ will have $\rho(\bm{x}_{B})=k_{\varepsilon}$, and so these
will appear after the set of vectors in $R(\bm{B})$. If the rank of
$\bm{B}\bm{B}^{\prime}<q$ the proof follows the same arguments with a lower
dimension substituted for $q$. In the original co-ordinate system where
$\bm{\Sigma}_{\varepsilon}\neq\bm{I}$, the vectors $\bm{a}_{i}$ will be in the
range of $\bm{\Sigma}_{\varepsilon}^{-1}\bm{B}$ as seen by a change of
coordinate transform $\bm{a}=\tilde{\bm{a}}\bm{\Sigma}_{\varepsilon}^{-1/2}$.

Part II: CCA coefficient vectors

Now if we consider Canonical Correlation Analysis (CCA), we look for the
linear combination of the columns of $\bm{Z}$ which maximizes a linear
combination of the signals $F=(f_{1}(t),f_{2}(t),...,f_{k}(t))$. Without loss
of generality,

$\Sigma_{F}=\bm{I},$ (B.34)

and thus

$\Sigma_{ZF}=\bm{B}$ (B.35)

By definition the canonical variables of Z have linear weights given by the
eigenvectors of

$\displaystyle\Sigma_{Z}^{-1}\Sigma_{ZF}\Sigma_{F}^{-1}\Sigma_{FZ}$

which, using Equation B.34 and Equation B.35 and $\Sigma_{F}=\bm{I}$ become

$\displaystyle(\bm{B}\bm{B}^{\prime}+\Sigma_{\varepsilon})^{-1}\bm{B}\bm{B}^{\prime}$

which is equivalent to the sequence of problems

maximize
$\displaystyle\quad\phi(\bm{a}_{i})=\frac{1}{1+\frac{1}{\sigma(\bm{a}_{i})}}$
subject to $\displaystyle\quad\bm{a}_{i}^{\prime}\bm{a}_{i}=1$
$\displaystyle\quad\bm{a}_{i}^{\prime}\bm{a}_{j}=0,\quad\forall j<i<q.$ (B.38)

But this is equivalent to

maximize
$\displaystyle\quad\sigma(\bm{a}_{i})=\frac{\bm{a}_{i}^{\prime}\bm{B}\bm{B}\bm{a}_{i}}{\bm{a}_{i}^{\prime}\bm{\Sigma}_{\varepsilon}\bm{a}_{i}}$
subject to $\displaystyle\quad\bm{a}_{i}^{\prime}\bm{a}_{i}=1$
$\displaystyle\quad\bm{a}_{i}^{\prime}\bm{a}_{j}=0,\quad\forall j<i<q.$ (B.39)

Now, making the change of variable
$\tilde{\bm{a}_{i}}=\bm{\Sigma}_{\varepsilon}^{1/2}a_{i}$
$\tilde{\bm{B}}=\bm{\Sigma}_{\varepsilon}^{-1/2}\bm{B}$ our CCA problem
reduces to finding the eigenvalues of $\tilde{\bm{B}}\tilde{\bm{B}}^{\prime}$,
sorted in the diagonal matrix
$\tilde{\bm{D}}=diag(\tilde{d}_{1},\tilde{d}_{2},\dots,\tilde{d}_{q})$. By the
spectral theorem

$\displaystyle\tilde{\bm{B}}\tilde{\bm{B}}^{\prime}=$
$\displaystyle\tilde{\bm{A}}\tilde{\bm{D}}\tilde{\bm{A}}$
$\displaystyle\bm{A}=$
$\displaystyle\bm{\Sigma}_{\varepsilon}^{-1}\bm{B}\bm{B}^{\prime}\bm{A}\tilde{\bm{D}}^{-1}$
$\displaystyle\bm{A}=$
$\displaystyle\bm{\Sigma}_{\varepsilon}^{-1}\bm{B}\bm{U},$ (B.40)

where $\bm{U}=\bm{B}^{\prime}\bm{A}\tilde{\bm{D}}^{-1}$. Now we see that the
eigenvectors are linear combinations of
$\bm{\Sigma}_{\varepsilon}^{-1}\bm{b}_{i},\quad\forall i\leq q$. ∎

###### Proof of Property 1.

A linear combination, $\bm{a}^{\prime}\bm{Z}(t)$, of the $p$-vector time
series $\bm{Z}(t)$ can be expressed as
$\bm{a}^{\prime}\bm{Z}(t)=\bm{b}^{\prime}\bm{f}(t)$, where
$\bm{b}^{\prime}=\bm{a}^{\prime}\bm{B}$, $\bm{Z}(t)=\bm{B}\bm{f}(t)$, $f(t)$
is the unknown underlying $q$-vector factor time series, and A is the unknown
$p\times q$ loading matrix. Let $r(\bm{a})$ denote the unit-lag
autocorrelation of the scalar time series $\bm{a}^{\prime}\bm{Z}(t)$. Then the
orthogonality of the factors $\bm{f}(t)$ yields
$r(\bm{a})=\bm{b}^{\prime}diag(\bm{k})\bm{b}/\bm{b}^{\prime}\bm{b}$ where $r$
is the $q$-vector of the underlying factor autocorrelations in decreasing
order. Since $r(\bm{a})$ is therefore a convex combination of the factor
autocorrelations, $R(\bm{a})$ cannot be greater than the largest of the
underlying factor autocorrelations, $r_{1}$. Therefore MAF-1, which maximizes
$r(\bm{a})$, yields factor loadings $b_{i}=0$ for $i>1$, and a MAF-1 time
series $\bm{a}^{\prime}\bm{Z}(t)$ is proportional to the underlying factor-1
time series $f_{1}(t)$. For $q<p$, the MAF-1 optimizing coefficient
$\bm{a}^{\prime}$ for the linear combination of $\bm{a}^{\prime}\bm{Z}(t)$ is
not unique because the $p\times p$ cross-covariance of the vector time series
$\bm{Z}(t)$ has rank $q$. But every linear combination that maximizes
autocorrelation will yield a time series that is proportional to the
underlying factor time series $f_{1}(t)$. Similarly, successive orthogonal MAF
time series, up to MAF-$q$, will evaluate to the corresponding successive
underlying factor time series $f_{2}(t),...,f_{q}(t)$. ∎