# Lagged teleconnections of climate variables
identified via complex rotated Maximum Covariance Analysis

Niclas Rieger
Centre de Recerca Matemàtica (CRM)
Departament de Física
Universitat Autònoma de Barcelona
Bellaterra, Spain
<EMAIL_ADDRESS>
& Álvaro Corral
Centre de Recerca Matemàtica (CRM)
Bellaterra, Spain
Complexity Science Hub Vienna
Vienna, Austria
<EMAIL_ADDRESS>Estrella Olmedo
Institute of Marine Sciences (ICM-CSIC)
Barcelona Expert Center (BEC)
Barcelona, Spain
<EMAIL_ADDRESS>& Antonio Turiel
Institute of Marine Sciences (ICM-CSIC)
Barcelona Expert Center (BEC)
Barcelona, Spain
<EMAIL_ADDRESS>
###### Abstract

A proper description of ocean-atmosphere interactions is key for a correct
understanding of climate evolution. The interplay among the different
variables acting over the climate is complex, often leading to correlations
across long spatial distances (teleconnections). In some occasions, those
teleconnections occur with quite significant temporal shifts that are
fundamental for the understanding of the underlying phenomena but which are
poorly captured by standard methods. Applying orthogonal decomposition such as
Maximum Covariance Analysis (MCA) to geophysical data sets allows to extract
common dominant patterns between two different variables, but generally
suffers from (i) the non-physical orthogonal constraint as well as (ii) the
consideration of simple correlations, whereby temporally offset signals are
not detected. Here we propose an extension, complex rotated MCA, to address
both limitations. We transform our signals using the Hilbert transform and
perform the orthogonal decomposition in complex space, allowing us to
correctly correlate out-of-phase signals. Subsequent Varimax rotation removes
the orthogonal constraints, leading to more physically meaningful modes of
geophysical variability. As an example of application, we have employed this
method on sea surface temperature and continental precipitation; our method
successfully captures the temporal and spatial interactions between these two
variables, namely for (i) the seasonal cycle, (ii) canonical ENSO, (iii) the
global warming trend, (iv) the Pacific Decadal Oscillation, (v) ENSO Modoki
and finally (vi) the Atlantic Meridional Mode. The complex rotated modes of
MCA provide information on the regional amplitude, and under certain
conditions, the regional time lag between changes on ocean temperature and
land precipitation.

_Keywords_ MCA $\cdot$ rotation $\cdot$ complex $\cdot$ Hilbert transform
$\cdot$ lagged teleconnections $\cdot$ xmca

## Significance Statement

Correlations between time series of different climate variables are often
time-lagged and can appear over long spatial distances. Our goal was to
develop a method that allows the simultaneous identification of the dominant
spatial patterns and their time lags between two different climate variables.
Using sea surface temperatures and continental rainfalls as example, our
method extracts the different dynamics of the seasonal cycle as well as 5
other well-known climate phenomena. Especially for cyclic time series like the
seasonal cycle, the relative time lags at different locations can be
determined precisely, whereas for acyclic time series only qualitative
statements about time lags can be made. In future studies, we expect new
insights into the dynamical structure of the Madden-Julian oscillation thanks
to this method, which is readily available as a Python package.

##  1 Introduction

The Earth’s climate system is extremely complicated and deciphering the web of
interdependencies and influences of different climate subsystems is an
involved challenge. As the quantity and quality of Earth observations increase
thanks to the advances in remote sensing, so too does the amount of data that
needs to be processed. Data-driven dimensionality reduction methods are
therefore crucial for climate studies, as they allow high-dimensional spatio-
temporally resolved signals to be disaggregated into the dominant patterns,
while still capturing the subtle details of higher resolution data. As such,
Principal Component Analysis (PCA), or Empirical Orthogonal Functions (EOF)
analysis as it is often referred to in climate science, allows to identify the
dominant internal structure of the variability as expressed by the variance,
with a variety of different available versions of PCA proving the popularity
of such methods in climate science [24, 27, e.g.].

Climate phenomena with different expression in oceanic and atmospheric
variables, such as the El Niño-Southern Oscillation (ENSO), however, require
the simultaneous analysis of several variables for a more comprehensive
description. In principle, multivariate PCA [40] makes it possible to extract
the patterns of co-variability of more than one variable. However,
multivariate PCA accumulates the variance and the covariance of variables with
very different variability in the same quantities. In consequence this may
mask co-varying patterns as low-variability patterns of one variable can be
erroneously accumulated in very dominant structures of one of the other,
large-variability variables [8].

Maximum Covariance Analysis (MCA)111Sometimes referred to as Singular Value
Decomposition (SVD) analysis. This name is unfortunate and should not be
confused with the actual factorisation technique of a real/complex matrix.
avoids this masking by taking into account only the covariance between two
sets of variables. As such, it bears similarity to Canonical Correlation
Analysis [31, CCA,] which aims at maximising the temporal correlation between
both variables. When the number of grid points (i.e. number of time series) is
higher than the number of observations (i.e. number of time steps) and the
data exhibits multicollinearity, as it is often the case for climate data, CCA
fails as it requires the individual variance matrices to be non-singular
unless regularised [65, 14, 26]. In case the two fields of variables are
identical, MCA reduces to PCA, the former thus being a natural generalisation
of PCA.

Yet the methods discussed above maximise instantaneous correlation and do not
consider time-delayed signals. To gain a deeper understanding of the dynamics
of climate phenomena, however, it is necessary to systematically investigate
time lags. A typical approach to tackle with this problem is to consider one
variable set with a time lag defined a priori followed by a MCA [44, e.g.].
However, this requires knowledge of the time lag which may vary from one
location to another [4].

In this paper, we propose complex rotated MCA to systematically investigate
the phase shift of two variables. We generate complex time series known as the
analytical signal, where the real and imaginary parts are related to each
other by the Hilbert transform, and decompose the covariance matrix in complex
space, in analogy to complex PCA [30, 6]. We also effectively reduce spectral
leakage inherent in the Hilbert transform of non-cyclic signals by using an
extrapolation method of the signal beyond its boundaries. Finally, to relax
the orthogonality constraint of the obtained solutions, we apply Varimax
rotation to the spatial patterns, which leads to more localised solutions and
thus facilitates their physical interpretation [55, 9].

To make the method readily accessible as a tool, we provide it as a Python
package, called xmca [56]. Due to the power and popularity of NumPy [64] and
xarray [32], both packages form the basis of xmca, so that their typical data
format can be used directly as input for analysis. The package is modularised
in a way that provides the user free choice whether standard, complex, rotated
or complex rotated MCA is to be performed. The user can also choose between
Varimax orthogonal rotation as well as Promax oblique rotation. Further, if
desired, standardisation of the input data is computed on the fly. The
different flavours work in the same way for PCA, if one instead of two fields
is provided as input.

The remainder of the article is structured as follows. Section 2 introduces
the methodology, where we briefly discuss MCA (Sec. 2.1), complex MCA (Sec.
2.2) and rotated MCA (Sec. 2.3). Section 3 describes the data used to test the
method using first synthetic data (Sec. 3.1) and then climatic variables (Sec.
3.2). Section 4 presents the results of both the synthetic (Sec. 4.1) and
real-world analysis (Sec. 4.2). We conclude our study and provide directions
for future research in Section 5.

##  2 Methods

###  2.1 Maximum Covariance Analysis

Let us consider two spatio-temporal data fields
$\boldsymbol{X}_{A}\in\mathbb{R}^{m\times n_{A}}$ and
$\boldsymbol{X}_{B}\in\mathbb{R}^{m\times n_{B}}$ representing two different
geophysical fields $s\in\left\\{A,B\right\\}$, both having temporal dimension
$m$ and spatial dimensions $n_{A}$ and $n_{B}$, respectively. Throughout the
text, the index $s$ is used implicitly without further definition to represent
one of the two fields. In the following, we will refer to the temporal
dimensions as the number of observations while we denote the spatial
dimensions by the number of grid points. Assuming each time series to have
zero mean, MCA then aims at maximising

$\displaystyle\boldsymbol{v}_{A}^{T}\boldsymbol{C}\boldsymbol{v}_{B},\quad
s.t.\quad\boldsymbol{v}_{A}^{T}\boldsymbol{v}_{A}=\boldsymbol{v}_{B}^{T}\boldsymbol{v}_{B}=1$
(1)

where $\boldsymbol{C}$ denotes the temporal covariance matrix and
$\boldsymbol{v}_{A},\boldsymbol{v}_{B}$ the spatial patterns, of both fields,
respectively. Mathematically, this can be achieved by applying the singular
value decomposition (SVD) to the covariance matrix,

$\displaystyle\boldsymbol{C}=\frac{1}{m-1}\boldsymbol{X}_{A}^{T}\boldsymbol{X}_{B}=\boldsymbol{V}_{A}\boldsymbol{\Sigma}\boldsymbol{V}_{B}^{T},$
(2)

with the columns of the obtained singular vector matrices
$(\boldsymbol{V}_{s})$ representing pairs of spatial patterns describing the
maximum amount of temporal covariance between both variables. The entries of
$\boldsymbol{\Sigma}\in\mathbb{R}^{n_{A}\times n_{B}}$ along the main
diagonal, the singular values $\sigma_{k}$, represent the covariance of each
spatial pattern pair $k$, providing a mean of estimating the relative
importance of each pair via the covariance fraction222Typically the squared
covariance fraction defined as
$\gamma_{k}^{*}=\sigma_{k}^{2}/\left(\sum_{j=1}^{\min(n_{A},n_{B})}\sigma_{j}^{2}\right)$
is considered for the relative importance of each mode for MCA. However, we
opt for the non-squared covariance fraction since the total explained
covariance is conserved under rotation, i.e. for $r$ rotated modes
$\sum_{i}^{r}\sigma_{i}=\sum_{i}^{r}\sigma_{i}^{*}$ where $\sigma_{i}^{*}$
refers to the covariance associated to mode $i$ after rotation. Furthermore,
this measure is comparable to the solutions obtained by PCA, and in fact it is
equivalent when $\boldsymbol{X}_{A}=\boldsymbol{X}_{B}$, for which MCA reduces
to PCA and the singular values equal the eigenvalues in PCA. $\gamma_{k}$

$\displaystyle\gamma_{k}=\sigma_{k}\left(\sum_{j=1}^{\min(n_{A},n_{B})}\sigma_{j}\right)^{-1}.$
(3)

By projecting the data fields on their respective singular vectors, we obtain
the corresponding temporal evolution for each spatial pattern given by the
columns of $\boldsymbol{P}_{s}=\boldsymbol{X}_{s}\boldsymbol{V}_{s}$. Since
the singular vectors are orthonormal, i.e.
$\boldsymbol{V}_{s}^{T}\boldsymbol{V}_{s}=\boldsymbol{1}_{n_{s}}$ with
$\boldsymbol{1}_{n_{s}}$ being the identity matrix of rank $n_{s}$, the
projections of the left and right field are uncorrelated, i.e.
$\boldsymbol{P}_{A}^{T}\boldsymbol{P}_{B}=\boldsymbol{\Sigma}/(m-1)$, while
the projections of the same field are correlated in general, i.e.
$\boldsymbol{P}_{s}^{T}\boldsymbol{P}_{s}$ is not a diagonal matrix. In this
paper, we will refer to the spatial patterns and their corresponding time
projections as empirical orthogonal functions (EOFs) and principal components
(PCs), respectively, according to the usual convention in climate science. The
EOFs and the PCs associated with a specific singular value $\sigma_{k}$ are
denoted as mode $k$.

###  2.2 Complex MCA

Propagating features or lagged signals could be detected by using a complex
representation of the input fields. In analogy to complex PCA [67, 53, 30, 6],
we complexify the real input fields via the Hilbert transform to construct the
analytical signal $\hat{\boldsymbol{X}}_{s}$ defined as

$\displaystyle\hat{\boldsymbol{X}}_{s}=\boldsymbol{X}_{s}+i\altmathcal{H}(\boldsymbol{X}_{s})$
(4)

where $\altmathcal{H}(\cdot)$ denotes the column-wise applied Hilbert
transform. The analytical signal constructed in that way is a unique complex
representation of the real signal, but whether it also represents a physical
reality depends on the frequency spectrum of the analysed signal. By
construction, the frequency components of the Hilbert transform are phase
shifted by $-\nicefrac{{\pi}}{{2}}$ with respect to those of the original
signal. Therefore, for narrow-bandwidth signals, the Hilbert transform has a
simple physical interpretation i.e., it represents a signal which arrives with
a lag of one fourth of the typical period. If the signal consists of multiple
dominating frequencies, however, the interpretation of the phase is more
elusive, as it cannot be simply associated with a single frequency. Thus, the
more narrow the signal bandwidth, the more directly we can relate the phase to
specific timings of its Hilbert transform [7].

Figure 1: Example illustrating the Hilbert transform using the Theta
extension. (a) Input signal $s(t)=sin(t)+ct+\epsilon(t)$ using an arbitrary
constant $c$ and Gaussian white noise $\epsilon$ with zero mean and unit
variance (blue) as well as the extended time series via
forecasting/backcasting using the optimised Theta model [20]. (b) Hilbert
transform $\altmathcal{H(\cdot)}$ of the original and extended signal,
respectively. (c) Comparison of both Hilbert transforms over the domain of the
original series.

A fundamental issue in the computation of the Hilbert transform arises when
non-stationary or drifting signals are processed. Such signals are non-cyclic,
and therefore, when the Fourier coefficients are calculated, strong boundary
effects can occur due to spectral leakage [7] (Fig. 1). This problem can be
circumvented by detrending the time series and considering only integer cycles
as well as by applying window functions to the time series (e.g. Hanning,
Hamming). However, this comes at the cost of information loss. Additionally,
such mitigation techniques are particularly ill suited to deal with non-
stationarities associated to climate change (that include not only noticeable
trends on the mean level, but also increases in the amplitude of some periodic
phenomena). Therefore, it seems important to introduce techniques capable to
dealing with non-stationary signals.

To mitigate spectral leakage across the boundaries of the time series, we
extrapolate the time series at both boundaries, to the past and to the future,
using the optimised Theta model [3, 20], a special case of an autoregressive
integrated moving average model with drift, ARIMA(0,1,1) [35]. The Theta model
is a relatively simple yet well performing extrapolation method. While the
forecast in itself actually works with non-cyclic signals, seasonal features
are considered via multiplicative classical decomposition, thus allowing
cyclic and non-cyclic signals to be extrapolated. This approach of handling
the seasonal structure of a time series requires the user to specify the
dominant period of the signal, $T_{s}$, beforehand. For discrete time series
$T_{s}$ represents the number of time steps needed to complete one cycle, that
is, e.g. $365$ for daily data considering an annual cycle, or $24$ for hourly
data with a daily cycle (for more details we refer the reader to [20]). We
then apply the Hilbert transform to the extended series, so the spectral
leakage is only important on the backward and forward extensions of it.
Finally, we extract the central part (removing the parts corresponding to the
extension), that correctly corresponds to the Hilbert transform of the
original series. Using this approach, we effectively reduce the edge effects
of the Hilbert transform compared to a non-processed time series (Fig. 1).
Notice that it is not necessary that the extrapolation faithfully reproduces
the characteristics of the original series; it just suffices for our purpose
that the extrapolated time series approximately continues the cyclic structure
at the original time series boundaries in order to reduce spectral leakage.
Apart from tracking trends, the exact extrapolation beyond the boundaries is
not essential since its effects on the central part of the Hilbert transform
are very marginal at most.

After the described complexification of the original time series, we follow
the steps of standard MCA, with the difference that the transpose
$(\cdot)^{T}$ incorporates the complex conjugate $\bar{(\cdot)}$ and the
obtained EOFs and projections PCs are complex and unitary. This allows us to
calculate the spatial amplitude $\boldsymbol{\mathcal{A}}_{s}$ and phase
function $\boldsymbol{\Theta}_{s}$ for both fields,

$\displaystyle\boldsymbol{\mathcal{A}}_{s}=$
$\displaystyle\left(\boldsymbol{V}_{s}\odot\boldsymbol{\bar{V}}_{s}\right)^{\odot\frac{1}{2}},$
(5) $\displaystyle\boldsymbol{\Theta}_{s}=$
$\displaystyle\arctan\\!2\left(\text{Im}(\boldsymbol{V}_{s}),\text{Re}~{}(\boldsymbol{V}_{s})\right),$
(6)

where $\odot$ denotes the element-wise multiplication/exponentiation, and
$\arctan\\!2$ refers to the two-argument arctangent which is calculated
element-wise (Appendix 5.1). Although this matrix notation seems somewhat
cumbersome compared to the more direct expression through scalar fields, it
allows us to be coherent with the rest of the paper. The phase function can be
interpreted directly as a time lag if the corresponding (real) PC has a
narrow-band spectrum with just one dominant frequency. If the spectrum is
rather broad-band or has several dominant frequencies, an interpretation of
the phase function is usually not straightforward. We note that the complex
EOFs derived from the SVD are only defined up to a phase shift of
$\exp{i\theta}$ with $\theta\in[0,2\pi]$. However, if complex EOFs are
obtained with a different phase shift $\theta$ (e.g. due to another SVD
algorithm), the change will also be reflected in the projected PCs, so that
taking into account both, PC and phase function, the results are unambiguous.

###  2.3 Rotated MCA

While orthogonality is often a mathematically desirable property, it does not
make a lot of sense from a purely geophysical standpoint. Therefore, standard
EOFs are difficult to interpret in the case of geophysical data. The major
drawbacks of EOFs due to orthogonality are twofold: First, EOFs are sensitive
to the selected spatial domain, that is including or removing some regions may
change large parts of the EOFs. Secondly, EOFs tend to split certain
geophysically meaningful patterns across several consecutive modes [55]. To
relax the orthogonality constraint to better accommodate the geophysical
reality, the EOFs can be rotated, which implies a linear transformation of the
first $r$ loaded333Loaded EOFs are weighted by the square root of the
corresponding singular value. EOFs $\boldsymbol{L}_{s,r}$. This concept, which
was originally developed in the context of PCA, can also be applied to MCA
[9]. For this, we apply the rotation matrix
$\boldsymbol{R}\in\mathbb{C}^{r\times r}$ to the loading matrix
$\boldsymbol{L}_{r}\in\mathbb{C}^{n\times r}$, $n=n_{A}+n_{B}$,

$\displaystyle\boldsymbol{L}_{r}=$
$\displaystyle\begin{pmatrix}\boldsymbol{L}_{A,r}\\\
\boldsymbol{L}_{B,r}\end{pmatrix}=\begin{pmatrix}\boldsymbol{V}_{A,r}\\\
\boldsymbol{V}_{B,r}\end{pmatrix}\boldsymbol{\Sigma}_{r}^{1/2},$ (7)

where $r$ reflects the respective submatrices containing only the first $r$
columns. In addition, $\boldsymbol{\Sigma}_{r}$ is the diagonal submatrix
containing only the first $r$ columns and rows.

There are a number of different criteria for defining the rotation matrix
$\boldsymbol{R}$ [55, e.g.], including the Varimax orthogonal rotation [37]
and the Promax oblique rotation [28], whose general aim is to regroup the
obtained patterns by approximating simple structures [63]. Mathematically,
Varimax rotation seeks to maximise the summed variances of squared loadings
which is achieved by (i) restricting rotated EOFs to be composed by only a few
numbers of grid points with high loadings while the remaining grid points
exhibit near-zero loadings and by (ii) limiting each grid point to contribute
to only one rotated EOF while having near-zero loadings for the other EOFs.
Since non-rotated EOFs are typically dense, that is consisting of mostly non-
zero values, Varimax rotation produces more sparse EOFs containing mostly zero
or close-to-zero values, leading to spatially compact structures which allow a
clearer interpretation. Promax oblique rotation builds upon the Varimax
solution by raising the rotated, normalised EOFs to the power $p\geq 1$ while
retaining the original sign, thus further reducing low loading compared to
high loading of the EOFs. Promax can be understood as an oblique
generalisation, with $p=1$ yielding a Varimax orthogonal solution. [28]
provides a value for $p$, which the authors consider appropriate for most
applications ($p=4$). In the extensive review [55] points out, however, that
the Promax rotation using $p=2$ consistently performs better, which is what we
will use in this paper. In order to keep the paper self-contained, we provide
a brief summary of both rotation criteria in Appendix 5.2.

The main difference between both rotation types is that Promax allows rotated
PCs to be correlated with each other, with higher values of $p$ typically
leading to stronger correlations. In contrast, Varimax solutions yield always
uncorrelated PCs. For both, Varimax orthogonal and Promax oblique rotation,
the obtained EOFs are no longer orthogonal. The question which rotation method
is the most suitable for a given analysis remains unsettled in the literature.
In reality, we do not expect geophysical signals to be perfectly uncorrelated,
which generally argues for applying an oblique rotation. Nevertheless, [19]
showed that Varimax orthogonal and Promax oblique solutions perform similarly,
in particular when the PCs obtained by the oblique solution exhibit low linear
Pearson correlation coefficients. In the presence of simple structures,
however, Promax oblique rotation performs better by effectively reducing the
number of grid points that contribute to each mode, hence further simplifying
the EOFs and increasing correlations among the PCs [19]. Therefore, the
decision on how many EOFs to rotate and which rotation type to perform is a
choice to be taken case-by-case and which we will explore in Sec. 4.2.

##  3 Data

To test our method, we apply it to artificial and real climate data. For the
artificial data sets, we consider complex MCA without rotation, as studies
already exist that demonstrate the better interpretability and lower
sensitivity to sampling errors of the Varimax-rotated solutions compared to
the unrotated EOFs. [45, 55, 10]. By means of two synthetic experiments we
seek to illustrate the advantages and caveats of complex MCA. In a first
experiment (Experiment I), we test the performance of complex MCA compared to
standard MCA considering time-lagged signals. In a second experiment
(Experiment II), we investigate how the Theta extension can improve the result
of complex MCA to non-stationary processes. Finally, we apply complex MCA with
rotation to climatic data that we expect to have intrinsic geophysical cycles
but are also affected by the non-stationarity of climate change.

###  3.1 Synthetic Data

We create two 2D spatio-temporal data fields $X_{A},X_{B}$ with coordinates
representing longitude $\lambda\in[0,359]$ and time $t\in[0,364]$ days. The
data generation model follows

$\displaystyle X_{s}(t,\lambda)=$ $\displaystyle
c_{s}(\lambda)\zeta_{s}(t,\lambda)+\epsilon(t,\lambda)$ (8)

where $\zeta_{s}(t,\lambda)$ depends on the individual experiment design,
$c_{s}(\lambda)=\cos^{2}{a_{s}\lambda}$ represents a "zonal" modulation factor
with scale factor $a_{s}$ and $\epsilon(t,\lambda)$ denotes Gaussian white
noise with zero mean and variance of $1$.
The idea of Experiment I is to highlight the advantage of using complexified
fields compared to standard MCA in the presence of moving patterns and phase-
shifted, stationary fields. Therefore, we define both signals as travelling
waves $\zeta_{s}(t,\lambda)=\cos{(\omega t+k_{s}\lambda+\phi_{s})}$. Our
parameter choices are motivated by the Madden-Julian Oscillation (MJO), which
is an eastwards propagating mode of deep convection and associated zonal wind
circulation in the tropical atmosphere [46, 71]. As such, the MJO is one of
the dominant drivers of intraseasonal variability in the tropics characterized
by a zonal propagation period of about 30-90 days. In general, MJO events tend
to dominate over the Indian Ocean and the western Pacific before they decay
towards the eastern part of the Pacific. Furthermore, when observing the MJO
through different variables, the spatial scale of MJO events can vary. While
the zonal wind circulation typically exhibits a wave number of about $1$, the
convective precipitation patterns may have zonal wave numbers of $1-3$. With
this in mind, we fix the model parameter to $\omega=\nicefrac{{2\pi}}{{56}}$
$\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}^{-1}$ (representing a period of
$56\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$),
$k_{A}=\nicefrac{{-2\pi}}{{360}}$, $k_{B}=\nicefrac{{-6\pi}}{{360}}$
(representing wave numbers 1 and 3, respectively) and phase shifts
$\phi_{A}=0$ and $\phi_{B}=\nicefrac{{\pi}}{{2}}$. Finally, the modulation
factor $c_{s}(\lambda)$ can be thought of describing regions of enhanced and
suppressed MJO activity in the tropics. We therefore choose the scale factor
$a_{A}=\nicefrac{{\pi}}{{360}}$ and $a_{B}=\nicefrac{{2\pi}}{{360}}$
representing 1 and 2 regions of enhanced activity, respectively. (Fig.2a).
In Experiment II, we investigate the response of complex MCA to non-cyclical,
non-stationary signals. Since MCA seeks to maximise covariance through a new
set of linear combinations, we restrict ourselves here to linear trends (Tab.
1). To keep the experiment as simple as possible, we also consider static
fields only. Then, the signal may be simply defined as
$\zeta_{s}(t,\lambda)=\zeta_{s}(t)=h_{s}t$, where
$h_{A}=\nicefrac{{3}}{{365}}$ $\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}^{-1}$
and $h_{B}=\nicefrac{{5}}{{365}}$
$\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}^{-1}$ describe linear trends. In
this case, the modulation factor $c_{s}(\lambda)$ may describe, for instance,
differential heating of the continents and the oceans due to global warming.
The actual values of the scale factors remain unchanged compared to Experiment
I (Fig. 3a). For a better overview, the parameters of both experiments are
summarised in Tab. 1.

Table 1: Models for synthetic data generation of $X_{s}$ for Experiment I and II (Section 33.1). For our final model, we consider the coordinate ranges $t=0,1,\dots,364$ days and $\lambda=0,1\dots,359$ (wihtout unit). | Experiment I | Experiment II
---|---|---
Field | $a_{s}$ |  $\omega$ ($\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}^{-1}$) | $k_{s}$ | $\phi_{s}$ |  $h_{s}$ ($\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}^{-1}$)
$A$ | $\nicefrac{{\pi}}{{360}}$ | $\nicefrac{{2\pi}}{{56}}$ | $-\nicefrac{{2\pi}}{{360}}$ | $0$ | $\nicefrac{{3}}{{365}}$
$B$ | $\nicefrac{{2\pi}}{{360}}$ | $\nicefrac{{2\pi}}{{56}}$ | $-\nicefrac{{6\pi}}{{360}}$ | $\nicefrac{{\pi}}{{2}}$ | $\nicefrac{{5}}{{365}}$

###  3.2 Climate Data

We analyse monthly means of global sea surface temperature (SST) and
continental precipitation using the extended ERA5 data set from 1950-2019 [29,
5] provided by the European Centre for Medium-Range Weather Forecasts (ECMWF)
as a state-of-the-art replacement of the ERA-Interim reanalysis [16]. In
general, trends and low frequency variability of surface temperature and
humidity are represented well for the period 1979-2019 [29, 60]. However, [60]
noticed strong biases in temperature and humidity records for Central and East
Africa in the 1950s. Via inspection, we found that East African precipitation
records in particular are strongly biased during that period towards much
higher values. In order to avoid providing our Theta model extension with an
incorrect starting point of the time series, we remove the first 10 years,
providing us data from 1960 to 2019. Furthermore, the SVD of the covariance
matrix is a rather memory intensive numerical operation, which is why we limit
the domain of interest from $None$ S to $None$ N with a $None$x$None$ spatial
resolution, guaranteeing most of the continents to be included into the
analysis. In total, the fields of SST and continental precipitation cover
$n_{A}=$21\,816$$ and $n_{B}=$14\,184$$ grid points, respectively. To align
the different temporal scales of highly variable precipitation and slow-
varying SST and to filter out the high frequency signals, we smooth both data
sets with a 6-month moving average, for each month taking into account the 3
preceding and the 2 following months. We choose this particular time window to
effectively filter out the semi-annual cycles present in the geophysical
observations. We then normalise both data fields by dividing each grid point
by its temporal standard deviation to ensure equal spatial importance. In a
non-standardised MCA analysis, regions of high rainfall, like the tropics,
would dominate the covariance patterns, since annual rainfall in the mid-
latitudes and subtropics is typically much lower. After normalisation, we
weight the data points located on the regular $None$x$None$ grid according to
their associated area on a sphere [51] by multiplying each grid point with
$\sqrt{\cos(\varphi_{j})}$, $\varphi_{j}$ being the latitude at grid point
$j$.

Climate indices of monthly means, used in this study for the sake of
comparison, are downloaded from the websites of the National Oceanic and
Atmospheric Administration (NOAA;
https://psl.noaa.gov/data/climateindices/list/), namely the Oceanic Niño Index
(ONI) provided by the NOAA Climate Prediction Center, the Pacific Decadal
Oscillation [[, PDO,]]mantua_pacific_2002 and the Atlantic Meridional Mode [[,
AMM,]]chiang_analogous_2004. We calculate the ENSO Modoki Index [[,
EMI,]]ashok_nino_2007 according to
$\text{EMI}=R_{A}-\nicefrac{{1}}{{2}}R_{B}-\nicefrac{{1}}{{2}}R_{C}$ where $R$
denotes the area-averaged ERA5 SST anomalies over the regions $A$ ($None$E -
$None$W, $None$S - $None$N), $B$ ($None$W - $None$W, $None$S - $None$N) and
$C$ ($None$E - $None$E, $None$S - $None$N), respectively. Furthermore, we
define an oceanic warming index (OWI) as the area average of ERA5 SST
anomalies ($None$ S - $None$ N) using the same latitude area weighting as for
the MCA data preparation. Finally, all indices are smoothed using a 6-month
moving average in alignment with our data pre-processing. We further center
and max-normalise all indices for better comparison with our obtained PCs.

##  4 Results

In the following, we discuss the results from the synthetic experiments before
investigating the results of the analysis of SST and continental
precipitation.

###  4.1 Synthetic Data

Before we start the investigation of the synthetic data, we apply standard and
complex MCA as a baseline study on Gaussian white noise. By comparing the
individual singular values as well as their sum obtained for each method (Tab.
2), we observe that complex MCA consistently yields higher covariance both for
each mode individually and overall, as a result of the additional time-lagged
cross-covariance between field $A$ and $B$. Therefore, when comparing standard
and complex MCA in the following, we will make use of the singular values
directly instead of the covariance fraction $\gamma_{k}$ in order to assess
the explained covariance by each mode.

####  4.1.1 Experiment I: Lagged signals

Applying standard MCA (denoted by subscript $std$), the travelling wave is
split into two modes the first explaining $10\%$ of the shared covariance, and
the second $9\%$ (Fig. 2b). This is not surprising since the singular vectors
are orthogonal and the associated PCs are uncorrelated, forcing the signals
into two distinct modes following a sine and a cosine, which are perfectly
uncorrelated for a full number of cycles and thus orthogonal. While the first
two PCS correctly indicate the temporal evolution of the field following a
sine/cosine, the EOFs fail to provide a clear and easy interpretation of the
associated spatial structures. In particular, without further a priori
knowledge of the expected covarying structures, it seems a trying exercise to
deduce the travelling wave from the obtained PCS and EOFs. Apart from that the
remaining modes do not show any more distinct patterns, basically representing
noise.

Figure 2: Results of Experiment I. (a) Illustration of synthetic data fields
$A$ (left) and $B$ (right) showing (top) the temporal evolution
$X_{s}(t,\lambda)$ for some selected longitudes $\lambda$ smoothed via a
weighted running average using a Gaussian kernel with a standard deviation of
$7$ days and (bottom) the associated modulation patterns $c_{s}(\lambda)$ for
both fields, respectively. (b) Standard MCA solution showing mode 1
(explaining $10\%$ of the total covariance) and mode 2 ($9\%$) with (top) the
associated PCs and (bottom). (c) Complex MCA solution showing mode 1 ($16\%$)
with (top) the associated real (Re) and imaginary (Im) parts of the first PC,
(centre) the spatial amplitude functions and (bottom) the phase functions. For
comparison the denoised signal $X_{s}(t,\lambda_{s,0}=0)$ (grey line,
$\lambda_{0}$), the spatial modulation pattern $c_{s}(\lambda)$ (grey line,
$c_{s}$) and the theoretical phase shift $\Phi_{s}=k_{s}\lambda-\phi_{s}$
(grey line, $\Phi_{s}$) are shown for both fields, respectively. The phase
functions indicates the phase shift that has to be applied to the time series
at a specific longitude $\lambda$ in order to obtain the real part of the non-
shifted PC (shown in top panel). Percentages in parenthesis represent the
covariance fraction $\gamma_{k}$. All EOFs and spatial amplitude functions are
max-normalised for the sake of readability.

In comparison, complex MCA (denoted by subscript $com$) allows to represent
the travelling wave by a single mode (Fig. 2c). Additionally, the varying
strength of the signal at different longitudes is clearly captured by the
spatial amplitude for both fields. Another advantage in the interpretation
compared to standard MCA is the spatial phase, which shows the migration of
the signal to higher longitudes. Identical phase values at different
longitudes indicate correlation while a phase shift of $\pi$ represents anti-
correlation. We note that we have not used the Theta model extension here,
since the signals in this experiment are stationary and spectral leakage due
to non-integer number of periods can be assumed to be marginal. Furthermore,
comparing the singular values $\sigma_{com,1}>\sigma_{std,1}+\sigma_{std,2}$
(Tab. 2) indicates that the complex mode 1 accounts for more covariance and
thus captures more of the travelling wave than the combined mode 1 and 2 of
the standard solution. This makes sense, since the contributions of longitudes
to the individual EOFs 1 and 2 of the standard solution will decrease for
phase shifts which are not equal to $0$ (i.e. sine) or $\nicefrac{{\pi}}{{2}}$
(i.e. cosine). As a pitfall of complex MCA, however, it should be noted that
longitudes whose spatial amplitude function is very low tend to have noisy
phase function values. In our experiment, this is obvious for longitudes of
both fields where $c_{s}(\lambda)\approx 0$ and the associated phase function
does no longer follow the positive linear trend. At these longitudes, the max-
normalised spatial amplitude function falls below $0.2$ for both fields,
serving as a general orientation for the consideration of phase function
values in the following analysis. In general, it is therefore advisable not to
consider regions with low amplitudes.

Figure 3: Results of Experiment II. (a) Illustration of synthetic data fields
$A$ (left) and $B$ (right) showing (top) the temporal evolution
$X_{s}(t,\lambda)$ for some selected longitudes $\lambda$ smoothed via a
weighted running average using a Gaussian kernel with a standard deviation of
$7$ days and (bottom) the associated modulation patterns $c_{s}(\lambda)$ for
both fields, respectively. (b) Complex MCA solution showing mode 1 (explaining
$28\%$ of the total covariance) with (top) the associated real and imaginary
parts of the first PC, (centre) the spatial amplitude functions and (bottom)
the phase functions. (c) As in (b) but showing mode 1 ($17\%$) for the Theta
extended complex MCA. For comparison theoretical real and imaginary part of a
infinite linear trend (grey lines, Re∞ and Im∞), the spatial modulation
pattern $c_{s}(\lambda)$ (grey line, $c_{s}$) and the theoretical phase shift
$\Phi_{s}=0$ (grey line, $\Phi_{s}$) are shown for both fields, respectively.
Percentages in parenthesis represent the covariance fraction $\gamma_{k}$. All
spatial amplitude functions are max-normalised for the sake of readability.

####  4.1.2 Experiment 2: Trends and non-cyclic signals

As discussed in Sec. 2.2, non-periodic behaviour usually leads to undesired
boundary effects caused by spectral leakage in the frequency domain. In our
experiment, this effect is clearly evident for both PCs of the complex MCA
(denoted by subscript $com$) (Fig. 3b). Note that although the boundary
effects seem to occur only in the imaginary parts of the PCs, this may not be
the case in general. The Theta extended complex MCA (denoted by subscript
$thc$), on the other hand, successfully mitigates the boundary effects of the
PCs (Fig. 3c), where we have set the Theta period $T_{s}=1$ as the time series
have no seasonality. Nevertheless, the spatial patterns are similar for both
methods and can hardly be distinguished visually. Considering only spatially
static fields in this experiment, the phase function is constant for both
fields, with exception at longitudes of low amplitude values as already
mentioned in Experiment I. Examining the singular values, we notice that
$\sigma_{com,1}>>\sigma_{thc,1}\approx\sigma_{std,1}$ (Tab. 2), indicating the
increased covariance of the standard Hilbert transform due to boundary effects
compared to the Theta model extended Hilbert transform. This implies that
boundary effects created by the Hilbert transforms can strongly affect
correlations and, depending on their magnitude, lead to a severe "inflation"
of the singular values. As a consequence, the boundary effects can appear as
parts of the first modes, thus completely misleading the interpretation of the
results.

More generally, the frequency spectrum of a linear trend on a given interval
is typically broadband and thus the analytical signal constructed by the
Hilbert transform has not a physical interpretation in terms of characteristic
frequencies. Therefore, the phase function cannot be interpreted in terms of a
physical phase shift. The only exception is for a phase shift of
$\theta=0;\pm\pi$ since these relative phase shifts represent correlating and
anti-correlating signals disregarding the mathematical nature of the phase. A
fundamental consequence of this is that for modes whose PC is broadband (e.g.
a trend), only correlating ($\theta\approx 0$) as well as anti-correlating
($\theta\approx\pm\pi$) patterns should be considered.

Table 2: Singular values of MCA (Standard), complex MCA (Complex) and Theta extended complex MCA (Complex + Theta) obtained for purely white noise random fields (Noise) and the synthetic experiments described in Section 4.1 (Experiment I and II). Covariance fraction $\gamma_{k}$ for mode $k$ (in $\%$ shown within parenthesis) is calculated according to Eq. (3). Some exemplary modes (bold) are depicted in Fig. 2 and 3. Total cumulated covariance is shown in the last row. | Noise | Experiment I | Experiment II
---|---|---|---
Mode | Standard | Complex | Standard | Complex | Standard | Complex | Complex + Theta
$\sigma_{1}$ |  $2.58$ ($1.0$) |  $8.37$ ($1.3$) |  $\mathbf{27.43}$ ($9.8$) |  $\mathbf{102.77}$ ($16.1$) |  $89.1$ ($31.4$) |  $\mathbf{179.82}$ ($28.2$) |  $\mathbf{93.69}$ ($17.0$)
$\sigma_{2}$ |  $2.53$ ($0.9$) |  $8.11$ ($1.3$) |  $\mathbf{24.05}$ ($8.6$) |  $6.97$ ($1.1$) |  $1.95$ ($0.7$) |  $6.19$ ($1.0$) |  $6.18$ ($1.1$)
$\sigma_{3}$ |  $2.50$ ($0.9$) |  $7.84$ ($1.2$) |  $2.21$ ($0.8$) |  $6.82$ ($1.1$) |  $1.92$ ($0.7$) |  $5.96$ ($0.9$) |  $5.97$ ($1.1$)
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
$\sum_{j}\sigma_{j}$ | $269.91$ | $629.95$ | $279.96$ | $640.02$ | $284.1$ | $637.97$ | $551.16$

###  4.2 Climate Data: SST & Continental Precipitation

We apply Theta extended, complex MCA to SST and continental precipitation
using a Theta period $T_{s}=12$ to account for the seasonal cycle. The
dimensionality of the problem can be greatly reduced with the first $72$ modes
explaining more than $99\%$ of the lagged covariance (Fig. 4a). In order to
simplify the obtained patterns and to increase the physical meaning of the
results (Sec. 2.3), we rotate the first $150$ modes explaining
$99.82\text{\,}\mathrm{\char 37\relax}$ of the lagged covariance. Our decision
to rotate $150$ modes is motivated by the idea of retaining as much
information as possible without including the noise which dominates the higher
modes. Since we observe a marked drop in singular values at about the mode
number $100$ followed by an exponential decrease in singular values (Fig. 4b),
we guess that there is a negligible information content in higher modes.
Moreover, the quality of the reconstructed signal is rather independent of the
exact number of rotated modes, with $150\pm 50$ yielding basically identical
results for the first 6 modes for both Varimax orthogonal and Promax $p=2$
oblique rotation. Therefore, we will restrict our discussion in the following
to the first 6 modes. In order to address the question of the rotation method
to be chosen, we note that Promax oblique rotation performs better when simple
structures are present and correlations among PCs are high [19]. We expect our
first mode to be dominated by the shared dynamics of the seasonal cycle of
both SST and continental precipitation, which is indeed what we find for both
Promax (not shown here) and Varimax solution (Fig. 5). This mode, however, is
fairly global and as such does not represent a simple structure. Furthermore,
we observe that the correlations among the first six Promax-rotated PCs range
from -0.15 to 0.13 only, underlining that the Promax oblique solution is close
to orthogonal and differences to the Varimax solutions only marginal, at least
for the first six modes. Since the Promax oblique solution seemingly does not
provide a better results, we opt for the somewhat simpler Varimax orthogonal
rotation. In the following, we will discuss the Varimax-rotated modes by
investigating the real part of the PCs, the spatial amplitude and phase
functions for both fields, SST and continental precipitation, respectively.

In our representations of the modes, we remove non-significant, "noisy" phase
values (see Sec. 4.1) by masking out regions in the spatial amplitude and
phase function exhibiting a max-normalised amplitude of $<0.25$.

Mode 1 describes $62.7\text{\,}\mathrm{\char 37\relax}$ of the covariance
between SST and continental precipitation clearly showing the annual cycle
(Figs. 5 and 12a). As expected, the annual cycle shows itself in both
variables on a global scale, with the exception of the equatorial ocean, where
the seasonal SST variations are only weak. The phase function correctly
identifies the anti-correlation between the northern and southern oceans. It
also suggests that the eastern equatorial Pacific and the equatorial Indian
Ocean nevertheless show a weak seasonal signal, which is, however, positively
phase shifted relative to the rest of the southern ocean. On the continents,
precipitation dominates mainly in monsoon areas and the tropics. The phase
function illustrates the division into June to August (white) and December to
February (black) dominated rainfall systems and identifies corresponding
transition zones, as e.g. over the South American rain forest and central
North America. It also highlights some interesting dynamical regions which
stand out of their respective environment, namely, the West and East coast of
North America, the Mediterranean region, and to a lesser degree the East Asian
monsoon region in China.

Figure 4: Described covariance of the first 300 modes of complex MCA applied
to SST and continental precipitation showing the (a) accumulated and (b)
individual share. Dashed line estimates the boundary of suspected "noisy"
modes. Figure 5: Complex Varimax-rotated MCA of (left) SST and (right)
continental precipitation showing mode 1 and its relative importance indicated
by the covariance fraction. (top) The amplitude functions show the regions
predominantly contributing to the mode. (centre) The phase functions depict
the relative phase shifts with respect to (bottom) the corresponding real part
of the PC, where $0$ (blue) means correlation and $\pm\pi$ (red) shows anti-
correlation. For each grid point, the corresponding PC can be computed from
the given phase value by applying the negative phase shift to the PC of phase
$0$ (bottom). The amplitude functions, PCs and indices shown are all max-
normalised for the sake of comparability. In both amplitude and phase function
regions with a max-normalised amplitude below $0.25$ are masked out. Red
circles mark North Pacific ($None$ E, $None$ N), Indian Ocean ($None$ W,
$None$ N) and South Pacific ($None$ E, $None$ S) for SST and Poyang lake,
China ($None$ W, $None$ N), Phnom Penh, Cambodia ($None$ W, $None$ N) and
Darwin, Australia ($None$ W, $None$ S) for precipitation. The phase shifted
PCs of these locations are examined in Fig. 6.

The obtained mode provides an instructive example to highlight the benefits of
the complexified approach due to the mode’s corresponding narrow-band
frequency spectrum (Fig. 12a). Using the dominant periodicity of mode 1,
$T_{1}=12$ months, the interpretation of a phase shift $\theta$, given by the
spatial phase function, as time lag $\tau$ is straightforward and can be
computed via $\tau=\nicefrac{{T_{1}\theta}}{{2\pi}}$. Doing this for some
exemplary locations (denoted by red circles in Fig. 5), we observe that the
seasonal SST maximum of the North Pacific follows the one in the South Pacific
by $196$ days, that is, being almost perfectly anti-correlated (Fig. 6). The
relative time shifts of the SST and precipitation maximum of $81$ days in the
Indian Ocean and $86$ days at Poyang lake, China, is close to $3$ months and
as such translates to a phase shift of about $-\nicefrac{{\pi}}{{2}}$,
something that cannot be picked up as feature within one mode when using
standard MCA. It should also be noted that although the sampling frequency of
SST and precipitation is monthly, the phase function is continuous and thus
allows to infer time lags on shorter time scales. This is, for instance, the
case for the seasonal precipitation maximum in Darwin, Australia, which
precedes the seasonal SST maximum of the South Pacific by $27$ days.

Figure 6: Comparison between 6-month moving average of ERA5 SST and
continental precipitation (dotted lines) and the reconstructed time series
based on mode 1 only (thick lines) for the exemplary locations defined in Fig.
5 covering the years 1980 to 1985. Days in legend refer to the time shift
$\tau$ of the individual locations with respect to the SST of the South
Pacific derived from the spatial phase functions (Fig. 5). Figure 7: As in
Fig. 5 but showing mode 2 corresponding to the Oceanic Niño Index (ONI)
provided by NOAA Climate Prediction Center as described in Section 3.2. The
Spearman correlation between ONI and the real part of SST PC 2 is a $r=0.941$
with $p\text{-value }<0.001$. Figure 8: As in Fig. 5 but showing mode 3
corresponding to the oceanic warming index (OWI) as described in Section 3.2.
The Spearman correlation between the OWI and the real part of SST PC 3 is
$r=0.781$ with a $p\text{-value }<0.001$. Figure 9: As in Fig. 5 but showing
mode 4 corresponding to the Pacific Decadal Oscillation [[,
PDO,]]mantua_pacific_2002 as described in Section 3.2. The Spearman
correlation between the PDO and the real part of SST PC 4 is $r=0.565$ with a
$p\text{-value }<0.001$.

Mode 2 ($7.0\text{\,}\mathrm{\char 37\relax}$) can be clearly associated to
the ocean-atmosphere phenomenon of El Niño - Southern Oscillation (ENSO) (Fig.
7 and 12b). During El Niño, higher SST in the central and eastern Pacific and
lower SST in the western Pacific positively correlate with heavier rainfall
within a narrow band along the west coast and the southeastern coast of South
America [62], the east and west coasts of North America [57], the East Asian
monsoon region [68] and the Horn of Africa [36]. At the same time,
precipitation decreases in northern South America [62], Oceania [15], South
Africa [21] and in the Indian monsoon region [11]. During La Niña (the
counterphase of El Niño), these correlations are reversed. Recently, similar
teleconnections have been identified via event coincidence analysis [69]. For
a current summary of established ENSO-related rainfall patterns during El Niño
and La Niña see [43].

Our result also shows the dynamical link between ENSO in the Pacific Ocean and
the Indian Ocean [39], the South China Sea [38] and the Tropical North
Atlantic [17, 58, 1, 12] in accordance with previous studies. Interestingly,
it was found that the ENSO related SST teleconnections in the remote oceans
often occur with some delays, with the Indian Ocean typically peaking $\sim 3$
months and the South China Sea and tropical North Atlantic $\sim 5$ months
after the SST peak in the Pacific ENSO region [17, 38, 58]. The mechanism
behind these lagged responses, known as the atmospheric bridge, is based on
the characteristic atmospheric circulation during El Niño which causes changes
in cloud cover and evaporation over the remote oceans, leading to increased
net heat flux and SSTs [42, 38]. However, due to the broadband frequency
spectrum of the SST PC (Fig. 12b), with most of the energy contained at four
different peaks around $2.5$, $3.5$, $5$ and $11$ years, the phase cannot
simply be translated into a time shift. Nevertheless, the tropical North
Atlantic clearly exhibits more positive phase values compared to the Pacific
El Niño region, therefore indicating to be positively time shifted relative to
the Pacific.

Mode 3 ($6.9\text{\,}\mathrm{\char 37\relax}$) represents global warming and
the associated changes in precipitation patterns (Fig. 8). The warming SST
patterns clearly emerge in all major ocean basins, although more pronounced in
the northern hemisphere due to the asymmetric response of the northern and
southern trade winds to global warming [70]. We also note a pronounced warming
of the western part of both, the Pacific and the Atlantic basins, both regions
of enhanced ocean heat transport [70, 25]. Similar to these oceanic trends, we
also observe global trends in the precipitation patterns, with decreasing
rainfall over the Mediterranean, South Africa, Australia, South America and
parts of western North America. There seems to be also a decrease of rainfall
over the west Asian monsoon region. In contrast to that, the results suggest
increased precipitation over the Indian monsoon region as well as some
localised regions in Europe, South and North America. These results are
largely in agreement with studies based on observational data [23] and, more
recently, on CMIP5 climate simulations [22]. Finally, it should be stressed,
that both PCs, SST and precipitation, do only provide a meaningful
interpretation for phases $\theta\approx 0;\pm\pi$ (correlation, anti-
correlation), due to their broadband frequency spectra (Fig. 12c). For phase
shifts different from that, no clear conclusions can be drawn, as it is the
case e.g. for equatorial Africa where the phase shift is approximately
$+\pi/2$.

Mode 4 ($1.6\text{\,}\mathrm{\char 37\relax}$) shows a slow-oscillating
pattern of SST in the northeastern Pacific correlating with localized
precipitation patterns distributed over all continents (Fig. 9). The typical
spatial SST pattern is known as the Pacific Decadal Oscillation (PDO) [48] and
a well-established climate index. A combination of various processes
originating in the tropics and extra-tropics has been proposed as the physical
source of the PDO [50], with ENSO and PDO likely responding to the same
forcing function [52]. Our analysis, however, provides a means to disentangle
ENSO and PDO related precipitation patterns, which are often similar for
western North America [33], though reveal important differences e.g. for
Australia, the Indian subcontinent or the African Sahel region. Yet, care must
be taken when interpreting regions which have a phase shift different from
$0,\pm\pi$. Although the PDO exhibits a strong spectral energy at about
$35\text{\,}\mathrm{y}\mathrm{r}$, the mode contains also important features
at about $1\text{\,}\mathrm{y}\mathrm{r}$ to $4\text{\,}\mathrm{y}\mathrm{r}$
(Fig. 12d) making the phase interpretation physically less clear.

Mode 5 ($1.6\text{\,}\mathrm{\char 37\relax}$) describes an oscillating SST
anomaly mainly limited to the central Pacific (Fig. 10 and 12e) describing El
Niño Modoki [2] and represented by the El Niño Modoki Index (EMI). Higher SST
in the central Pacific and lower SST in the east Pacific correlate with
reduced precipitation in the East Asian monsoon region [18], Australia [61],
parts of South America [62] and South Africa [54] and vice versa. Some of
these teleconnections have recently been revealed by event coincidence
analysis [69] though important links to e.g. the East Asian monsoon region
were missing. Although our result suggests that continental rainfall in
certain regions of Africa, Arabia and the Americas are phase-shifted El Niño
Modoki expressions, future work has to show if these weak amplitudes are
significant.

Figure 10: As in Fig. 5 but showing mode 5 corresponding to the ENSO Modoki
Index [[, EMI,]]ashok_nino_2007 as described in Section 3.2. The Spearman
correlation between the EMI and the real part of SST PC 5 is $r=0.506$ with a
$p\text{-value }<0.001$. Figure 11: As in Fig. 5 but showing mode 6
corresponding to the Atlantic Meridional Mode [[, AMM,]]chiang_analogous_2004
as described in Section 3.2. The Spearman correlation between the AMM and the
real part of SST PC 6 is $r=0.741$ with a $p\text{-value }<0.001$.

Mode 6 ($1.1\text{\,}\mathrm{\char 37\relax}$) is characterised by a SST
pattern concentrated to the tropical and subarctic North Atlantic (Fig. 11 and
12f). The SST pattern, representing the Atlantic Meridional Mode (AMM), is the
dominant coupled ocean-atmospheric phenomenon in the tropical Atlantic [59]
and its impact on precipitation of the African Sahel zone, Central America and
the northern South American continent are well known [41, 49]. Only very
recently, [66] found a link between the AMM and Indian summer monsoon
rainfall, which they used to improve precipitation forecast. In addition, our
result suggests more tele-connections between the AMM and precipitation
variability over the Mediterranean, East Africa, the Congo basin, North
America and the East Asian monsoon region, providing much potential of
advancing local rainfall predictions in those areas. Due to the missing
intrinsic time scale with no clear periodicity (Fig. 12f), only correlating
and anti-correlating patterns should be interpreted. However, most of the
regions identified by this mode satisfy being correlated or anti-correlated.

Figure 12: Magnitude spectrum of the first 6 PCs considering only the real
part and using a Hanning window.

##  5 Conclusion

Understanding the intricate climate system is a challenging task, that
requires advanced statistical methods. Finding correlations among a set of
different climate variables is complicated by the frequently present lagged
responses of different variables to the same forcing. We show that complex
rotated MCA provides a practical tool to single out such modes from a high-
dimensional data space, even when the narrow-band assumption of the input
signals is only partially satisfied.

By taking into account the spatial amplitude and phase function of the
obtained complex modes, we obtain a simple approach to examine otherwise
complicated spatial and temporal structures. Our synthetic experiments
highlight that, in the case of phase-shifted signals, complex MCA can capture
a more comprehensive and complete picture of the correlations present.
However, they also show the sensitivity of the Hilbert transform to spectral
leakage caused by boundary conditions of the given time series, e.g. when the
time series clearly consists of a non-integer number of cycles and/or in the
presence of trends. Since the stationarity assumption does often not hold (due
to climate change), time series should generally be pre-processed when
applying complex MCA.

Extending the time series via the optimised Theta model mitigates the effect
of spectral leakage and produces physically reasonable PCs. This procedure
allows us to resolve trends and non-cyclic signals, although for trends the
obtained spatial phase function has a simple physical interpretation only for
correlating and anti-correlating patterns. Nevertheless, this approach
provides a means of applying complex MCA without the need to detrend the time
series of interest. Moreover, excluding the first mode, the original fields
can be reconstructed without the seasonal cycle, providing an advanced tool to
pre-process time series containing non-stationary and non-linear seasonal
features.

A general caveat in complex MCA is the fact that the phase function loses its
interpretation for PCs with a broadband frequency spectrum. But although the
spatial phase function has not always a simple physical interpretation for
most of the modes due to their broadband frequency spectrum, complex MCA
nevertheless can always be interpreted for the correlating and anti-
correlating patterns.

Applying complex rotated MCA to SST and continental precipitation, we clearly
identify the main shared dynamics in both variables, namely (i) the seasonal
cycle, (ii) the canonical ENSO, (iii) the trends associated to global warming
(iv) the PDO, (v) ENSO Modoki and (vi) the AMM. We also retrieve phase shifted
signals between the two climate variables. While for the seasonal cycle these
phase shifts can directly be translated into a time shift, the remaining modes
generally do not lend itself to such a simple interpretation due to their
broadband frequency spectra. But even without a precise equivalent as time
lag, the phase function provides a means to identify regions of lagged
correlations, for instance, between the SST of the Pacific and the tropical
North Atlantic during El Niño events. In addition, by focusing on narrow
frequency bands only, one may detect frequency ranges over which the phase is
not random, thus potentially uncovering more of dynamic structure of the
individual modes.

Many of the obtained correlation patterns in SST and continental precipitation
had already been evidenced by a multiplicity of different, partly regional
studies. The great advantage of complex rotated MCA is that it allows to
obtain all those patterns by a single analysis of the correlation of two
geophysical variables at global scale in a more compact and easy-to-interpret
way. Besides, our results also point out to new ocean-atmospheric
teleconnections, that, to our knowledge, have not been reported, most notably
for the PDO and the AMM.

Regarding future applications of complex rotated MCA, this method has the
potential for shedding light in the investigation of seasonal and sub-seasonal
phenomena, as well as for spatially propagating patterns. As future work, we
plan to analyse the Madden-Julian oscillation. Additionally, complex rotated
MCA could be used to evidence other connections between less studied
variables, as for instance sea surface salinity, sea surface height, soil
moisture, winds, etc, what has the potential of evidence new phenomena and
novel aspects of existing or new teleconnections.

## Acknowledgments

This work is part of the Climate Advanced Forecasting of sub-seasonal Extremes
(CAFE) project and has been prepared in the framework of the doctorate in
Physics of the Autonomous University of Barcelona. The authors gratefully
acknowledge funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement 813844,
and also from the Spanish government through the ‘Severo Ochoa Centre of
Excellence’ accreditation (CEX2019-000928-S). This work is a contribution to
CSIC Thematic Interdiciplinary Platform TELEDETECT.

The authors acknowledge the Copernicus Climate Change Service and the Physical
Sciences Laboratory of NOAA for providing the data. The authors would also
like to thank J. Ballabrera-Poy for his valuable comments on a previous
version of the paper.

## Data Availability Statement

The ERA5 data used in this study are publicly available from the Copernicus
Climate Data Store at https://doi.org/10.24381/cds.f17050d7 as described in
[29]. ERA5 preliminary extension from 1950-1978 is described in [5]. All
climate indices used in this study are openly accessible. ENSO Modoki Index
(EMI) is provided by the website of JAMSTEC at
http://www.jamstec.go.jp/virtualearth/general/en/. The remaining indices can
be found on the website of the Physical Sciences Laboratory of NOAA at
https://psl.noaa.gov/data/climateindices/list/.

## Appendix

###  5.1 Definition of $\arctan\\!2$

Every non-zero complex number in Cartesian coordinates, $z=x+iy$, can be
transformed into polar complex coordinates, $z=ae^{i\theta}$, where
$a=(x^{2}+y^{2})^{1/2}$ is the amplitude, and $\theta$ the phase of $z$. For
each $z\neq 0$, the phase is only defined up to an integer multiple of $2\pi$,
resulting in an infinite number of possible values. In order to construct a
well-defined function $\theta(x,y)$, one typically limits the phase $\theta$
to $\left(-\pi,\pi\right]$. Then, the two-argument arctangent function
$\arctan\\!2(y,x)$ converts the values of $y$ and $x$ to the polar phase via:

$\theta=\arctan\\!2(y,x)=\begin{cases}\arctan(\frac{y}{x})&\quad\text{if
}x>0\\\ \arctan(\frac{y}{x})+\pi&\quad\text{if }x<0\text{ and }y\geq 0\\\
\arctan(\frac{y}{x})-\pi&\quad\text{if }x<0\text{ and }y<0\\\
+\frac{\pi}{2}&\quad\text{if }x=0\text{ and }y>0\\\
-\frac{\pi}{2}&\quad\text{if }x=0\text{ and }y<0\\\
\text{undefined}&\quad\text{if }x=0\text{ and }y=0\end{cases}$

###  5.2 Finding the rotation matrix $\boldsymbol{R}$

Let us assume a complex loading matrix
$\boldsymbol{L}_{r}\in\mathbb{C}^{n\times r}$ containing only the first $r$
columns (modes) which are to be rotated and $n$ denoting the number of grid
points. The number of grid points may be the sum of the grid points of two
different fields $n=n_{A}+n_{B}$ as it is the case for MCA and described by
Equation (7) or simply the total number of grid points if only one field is
considered as it is for PCA. In the following, we will drop the subscript $r$
in order to keep the notation simple and let $\widetilde{\cdot}$ denote the
rotated solutions. Furthermore, ∗ refers to the conjugate transpose of a
matrix and $|\cdot|$ denotes the absolute value of a complex number.

##### Varimax orthogonal rotation

The goal of Varimax rotation is to approximate simple structures [63] of the
EOFs which is achieved by simplifying the columns of $\boldsymbol{L}$ via an
orthogonal rotation $\boldsymbol{R}$. For this purpose, [37] defines the
simplicity $S_{k}$,

$\displaystyle
S_{k}=\frac{1}{n}\sum_{j=1}^{n}\left(|\widetilde{l}_{jk}|^{2}\right)^{2}-\frac{1}{n^{2}}\left(\sum_{j=1}^{n}|\widetilde{l}_{jk}|^{2}\right)^{2},\quad
k=1,\dots,r$ (9)

which measures the variance of the squared amplitude of the rotated loadings
$\widetilde{l}_{jk}$. With increasing variance, the squared rotated amplitudes
$|\widetilde{l}_{jk}|^{2}$ either become low or large, thus increasing
simplicity. The normalised simplicity $S$ then reads

$\displaystyle
S=\sum_{k=1}^{r}\left[\frac{1}{n}\sum_{j=1}^{n}\left(\frac{|\widetilde{l}_{jk}|^{2}}{h_{j}^{2}}\right)^{2}-\frac{1}{n^{2}}\left(\sum_{j=1}^{n}\frac{|\widetilde{l}_{jk}|^{2}}{h_{j}^{2}}\right)^{2}\right],$
(10)

where $h_{j}=\left(\sum_{k=1}^{r}|l_{jk}|^{2}\right)^{1/2}$ represents the
communality of grid point $j$, which is the amount of variance of the $j$th
gird point accounted for by the $r$ retained modes. Subsequently, the
normalised Varimax-rotated EOFs, $\widetilde{\boldsymbol{L}}_{V,norm}$, are
the solution to the Varimax criterion,

$\displaystyle\widetilde{\boldsymbol{L}}_{V,norm}=\boldsymbol{H}^{-1}\boldsymbol{L}\boldsymbol{R}\quad\text{s.t.}\quad\operatorname*{argmax}_{l}(S)\text{
and }\boldsymbol{R}^{*}\boldsymbol{R}=1,$ (11)

with the communality matrix $\boldsymbol{H}\in\mathbb{R}^{n\times n}$ whose
elements are given by $\text{diag}(h_{1},\dots,h_{n})$. Equation (11) can be
solved by an iterative process, in which the EOFs are rotated in pairs in
order to maximise $S$. Finally, the de-normalised Varimax-rotated EOFs can be
computed via
$\widetilde{\boldsymbol{L}}_{V}=\boldsymbol{H}\widetilde{\boldsymbol{L}}_{V,norm}$.

##### Promax oblique rotation

Achieving simple structures with Promax is done via an oblique Procrustes
transformation [34]. Every target matrix of rotated EOFs $\boldsymbol{T}$ can
always be approximated from a base matrix $\boldsymbol{B}$ via a linear
transformation $\boldsymbol{R}$,

$\displaystyle\boldsymbol{T}=\boldsymbol{B}\boldsymbol{R}+\boldsymbol{E},$
(12)

where $\boldsymbol{E}$ is an error matrix. Minimising
$\text{trace}(\boldsymbol{E}^{*}\boldsymbol{E})$ yields the complex Procrustes
equation,

$\displaystyle\boldsymbol{R}=\left(\boldsymbol{B}^{*}\boldsymbol{B}\right)^{-1}\boldsymbol{B}^{*}\boldsymbol{T}.$
(13)

The basic assumption of Promax is that an Varimax orthogonal rotation is a
reasonable approximation to an optimal oblique solution. Therefore, the base
matrix is chosen to be
$\boldsymbol{B}=\widetilde{\boldsymbol{H}}^{-1}\widetilde{\boldsymbol{L}}_{V}$
whose entries are normalised by the Varimax communalities
$\widetilde{h}_{j}=\sum_{k=1}^{r}|\widetilde{l}_{jk}|^{2}$. Then, the Promax
equation defines the elements of the target matrix $\boldsymbol{T}$,

$\displaystyle t_{jk}=|b_{jk}^{+}|^{p+1}/b_{jk}^{+},$ (14)

where + denotes the max-normalised entries given by
$b_{jk}^{+}=b_{jk}/\max_{j}|b_{jk}|$. The power parameter $p$ thus defines the
strength of the Promax operation, while the sign remains unchanged. Using
Equation (13), the de-normalised Promax-rotated EOFs are given by

$\displaystyle\widetilde{\boldsymbol{L}}_{P}=$
$\displaystyle\widetilde{\boldsymbol{L}}_{V}\boldsymbol{R}\boldsymbol{D}$
$\displaystyle=$
$\displaystyle\widetilde{\boldsymbol{H}}^{-1}\widetilde{\boldsymbol{L}}_{V}\left[\left(\widetilde{\boldsymbol{L}}_{V}^{*}\widetilde{\boldsymbol{H}}^{-2}\widetilde{\boldsymbol{L}}_{V}\right)^{-1}\widetilde{\boldsymbol{L}}_{V}^{*}\widetilde{\boldsymbol{H}}^{-1}\boldsymbol{T}\right]\boldsymbol{D},$
(15)

where the normalisation matrix is given by
$\boldsymbol{D}^{2}=\text{diag}\left(\boldsymbol{R}^{*}\boldsymbol{R}\right)^{-1}$.

## References

  * [1] Michael Alexander and James Scott  “The influence of ENSO on air-sea interaction in the Atlantic” _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2001GL014347  In _Geophysical Research Letters_ 29.14, 2002, pp. 46–1–46–4  DOI: https://doi.org/10.1029/2001GL014347
  * [2] Karumuri Ashok et al.  “El Niño Modoki and its possible teleconnection” _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2006JC003798  In _Journal of Geophysical Research: Oceans_ 112, 2007  DOI: 10.1029/2006JC003798
  * [3] V. Assimakopoulos and K. Nikolopoulos  “The theta model: a decomposition approach to forecasting”  In _International Journal of Forecasting_ 16.4, The M3- Competition, 2000, pp. 521–530  DOI: 10.1016/S0169-2070(00)00066-2
  * [4] J. Ballabrera-Poy, R. Murtugudde and A.. Busalacchi  “On the potential impact of sea surface salinity observations on ENSO predictions” _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2001JC000834  In _Journal of Geophysical Research: Oceans_ 107, 2002, pp. SRF 8–1–SRF 8–11  DOI: 10.1029/2001JC000834
  * [5] B. Bell et al.  “ERA5 monthly averaged data on single levels from 1950 to 1978 (preliminary version)”, 2020  DOI: 10.24381/cds.f17050d7
  * [6] Peter Bloomfield and Jerry M. Davis  “Orthogonal rotation of complex principal components” _eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/joc.3370140706  In _International Journal of Climatology_ 14.7, 1994, pp. 759–775  DOI: 10.1002/joc.3370140706
  * [7] B. Boashash  “Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals” Conference Name: Proceedings of the IEEE  In _Proceedings of the IEEE_ 80.4, 1992, pp. 520–538  DOI: 10.1109/5.135376
  * [8] C Bretherton, C Smith and JM Wallace  “An intercomparison of methods for finding coupled patterns in climate data”  In _Journal of climate_ 5.6, 1992, pp. 541–560
  * [9] Xinhua Cheng and Timothy J. Dunkerton  “Orthogonal Rotation of Spatial Patterns Derived from Singular Value Decomposition Analysis” Publisher: American Meteorological Society  In _Journal of Climate_ 8.11, 1995, pp. 2631–2643  DOI: 10.1175/1520-0442(1995)008<2631:OROSPD>2.0.CO;2
  * [10] Xinhua Cheng, Gregor Nitsche and John M. Wallace  “Robustness of Low-Frequency Circulation Patterns Derived from EOF and Rotated EOF Analyses” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 8.6, 1995, pp. 1709–1713  DOI: 10.1175/1520-0442(1995)008<1709:ROLFCP>2.0.CO;2
  * [11] Annalisa Cherchi and Antonio Navarra  “Influence of ENSO and of the Indian Ocean Dipole on the Indian summer monsoon variability” Publisher: Springer  In _Climate dynamics_ 41.1, 2013, pp. 81–103
  * [12] John C.. Chiang and Adam H. Sobel  “Tropical Tropospheric Temperature Variations Caused by ENSO and Their Influence on the Remote Tropical Climate” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 15.18, 2002, pp. 2616–2631  DOI: 10.1175/1520-0442(2002)015<2616:TTTVCB>2.0.CO;2
  * [13] John C.. Chiang and Daniel J. Vimont  “Analogous Pacific and Atlantic Meridional Modes of Tropical Atmosphere–Ocean Variability” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 17.21, 2004, pp. 4143–4158  DOI: 10.1175/JCLI4953.1
  * [14] Raul Cruz-Cano and Mei-Ling Ting Lee  “Fast regularized canonical correlation analysis”  In _Computational Statistics & Data Analysis_ 70, 2014, pp. 88–100  DOI: 10.1016/j.csda.2013.09.020
  * [15] Aiguo Dai and T… Wigley  “Global patterns of ENSO-induced precipitation” _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/1999GL011140  In _Geophysical Research Letters_ 27.9, 2000, pp. 1283–1286  DOI: https://doi.org/10.1029/1999GL011140
  * [16] DP Dee et al.  “The ERA-Interim reanalysis: configuration and performance of the data assimilation system”  In _Quarterly Journal of the Royal Meteorological Society_ 137.656, 2011, pp. 553–597  DOI: 10.1002/qj.828
  * [17] David B. Enfield and Dennis A. Mayer  “Tropical Atlantic sea surface temperature variability and its relation to El Niño-Southern Oscillation” _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/96JC03296  In _Journal of Geophysical Research: Oceans_ 102, 1997, pp. 929–945  DOI: https://doi.org/10.1029/96JC03296
  * [18] Juan Feng, Wen Chen, C.-Y. Tam and Wen Zhou  “Different impacts of El Niño and El Niño Modoki on China rainfall in the decaying phases”  In _International Journal of Climatology_ 31.14, 2011, pp. 2091–2101  DOI: 10.1002/joc.2217
  * [19] Holmes Finch  “Comparison of the Performance of Varimax and Promax Rotations: Factor Structure Recovery for Dichotomous Items” _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1745-3984.2006.00003.x  In _Journal of Educational Measurement_ 43.1, 2006, pp. 39–52  DOI: https://doi.org/10.1111/j.1745-3984.2006.00003.x
  * [20] Jose A. Fiorucci et al.  “Models for optimising the theta method and their relationship to state space models”  In _International Journal of Forecasting_ 32.4, 2016, pp. 1151–1161  DOI: 10.1016/j.ijforecast.2016.02.005
  * [21] Andrea E. Gaughan et al.  “Inter- and Intra-annual precipitation variability and associated relationships to ENSO and the IOD in southern Africa” _eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/joc.4448  In _International Journal of Climatology_ 36.4, 2016, pp. 1643–1656  DOI: https://doi.org/10.1002/joc.4448
  * [22] Filippo Giorgi, Francesca Raffaele and Erika Coppola  “The response of precipitation characteristics to global warming from climate projections” Publisher: Copernicus GmbH  In _Earth System Dynamics_ 10.1, 2019, pp. 73–89  DOI: https://doi.org/10.5194/esd-10-73-2019
  * [23] Guojun Gu and Robert F. Adler  “Spatial Patterns of Global Precipitation Change and Variability during 1901–2010” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 28.11, 2015, pp. 4431–4453  DOI: 10.1175/JCLI-D-14-00201.1
  * [24] A Hannachi, IT Jolliffe and DB Stephenson  “Empirical orthogonal functions and related techniques in atmospheric science: A review”  In _International Journal of Climatology_ 27.9, 2007, pp. 1119–1152  DOI: 10.1002/joc.1499
  * [25] A. Hannachi and N. Trendafilov  “Archetypal Analysis: Mining Weather and Climate Extremes” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 30.17, 2017, pp. 6927–6944  DOI: 10.1175/JCLI-D-16-0798.1
  * [26] Abdel Hannachi  “Regularised empirical orthogonal functions” Publisher: Taylor & Francis  In _Tellus A: Dynamic Meteorology and Oceanography_ 68.1, 2016, pp. 31723
  * [27] Abdelwaheb Hannachi  “Patterns Identification and Data Mining in Weather and Climate”, Springer Atmospheric Sciences  Springer International Publishing, 2021
  * [28] Alan E. Hendrickson and Paul Owen White  “Promax: A Quick Method for Rotation to Oblique Simple Structure” _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2044-8317.1964.tb00244.x  In _British Journal of Statistical Psychology_ 17.1, 1964, pp. 65–70  DOI: 10.1111/j.2044-8317.1964.tb00244.x
  * [29] H. Hersbach et al.  “ERA5 monthly averaged data on single levels from 1979 to present”, 2019  DOI: 10.24381/cds.f17050d7
  * [30] JD Horel  “Complex Principal Component Analysis: Theory and Examples”  In _Journal of Climate and Applied Meteorology_ 23.12, 1984, pp. 1660–1673  DOI: 10.1175/1520-0450(1984)023<1660:CPCATA>2.0.CO;2
  * [31] Harold Hotelling  “Relations Between Two Sets of Variates” Publisher: [Oxford University Press, Biometrika Trust]  In _Biometrika_ 28.3, 1936, pp. 321–377  DOI: 10.2307/2333955
  * [32] Stephan Hoyer and Joe Hamman  “xarray: N-D labeled Arrays and Datasets in Python” Number: 1 Publisher: Ubiquity Press  In _Journal of Open Research Software_ 5.1, 2017, pp. 10  DOI: 10.5334/jors.148
  * [33] Zeng-Zhen Hu and Bohua Huang  “Interferential Impact of ENSO and PDO on Dry and Wet Conditions in the U.S. Great Plains” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 22.22, 2009, pp. 6047–6065  DOI: 10.1175/2009JCLI2798.1
  * [34] John R. Hurley and Raymond B. Cattell  “The Procrustes program: Producing direct rotation to test a hypothesized factor structure” Publisher: University of Michigan, Mental Health Research Institute  In _Behavioral science_ 7.2, 1962, pp. 258
  * [35] Rob J. Hyndman and Baki Billah  “Unmasking the Theta method”  In _International Journal of Forecasting_ 19.2, 2003, pp. 287–290  DOI: 10.1016/S0169-2070(01)00143-1
  * [36] Matayo Indeje, Fredrick HM Semazzi and Laban J. Ogallo  “ENSO signals in East African rainfall seasons” Publisher: Wiley Online Library  In _International Journal of Climatology: A Journal of the Royal Meteorological Society_ 20.1, 2000, pp. 19–46
  * [37] Henry F. Kaiser  “The varimax criterion for analytic rotation in factor analysis”  In _Psychometrika_ 23.3, 1958, pp. 187–200  DOI: 10.1007/BF02289233
  * [38] Stephen A. Klein, Brian J. Soden and Ngar-Cheung Lau  “Remote Sea Surface Temperature Variations during ENSO: Evidence for a Tropical Atmospheric Bridge” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 12.4, 1999, pp. 917–932  DOI: 10.1175/1520-0442(1999)012<0917:RSSTVD>2.0.CO;2
  * [39] V. Krishnamurthy and Ben P. Kirtman  “Variability of the Indian Ocean: Relation to monsoon and ENSO” _eprint: https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1256/qj.01.166  In _Quarterly Journal of the Royal Meteorological Society_ 129.590, 2003, pp. 1623–1646  DOI: https://doi.org/10.1256/qj.01.166
  * [40] John E. Kutzbach  “Empirical eigenvectors of sea-level pressure, surface temperature and precipitation complexes over North America”  In _Journal of Applied Meteorology and Climatology_ 6.5, 1967, pp. 791–802
  * [41] Peter J. Lamb, Randy A. Peppler and Stefan Hastenrath  “Interannual variability in the tropical Atlantic” Number: 6076 Publisher: Nature Publishing Group  In _Nature_ 322.6076, 1986, pp. 238–240  DOI: 10.1038/322238a0
  * [42] Ngar-Cheung Lau and Mary Jo Nath  “The role of the “atmospheric bridge” in linking tropical Pacific ENSO events to extratropical SST anomalies”  In _Journal of Climate_ 9.9, 1996, pp. 2036–2057
  * [43] Nathan J.. Lenssen, Lisa Goddard and Simon Mason  “Seasonal Forecast Skill of ENSO Teleconnection Maps” Publisher: American Meteorological Society Section: Weather and Forecasting  In _Weather and Forecasting_ 35.6, 2020, pp. 2387–2406  DOI: 10.1175/WAF-D-19-0235.1
  * [44] L Li, R Schmitt, CC Ummenhofer and KB Karnauskas  “North Atlantic salinity as a predictor of Sahel rainfall”  In _Science Advances_ 2.5, 2016, pp. e1501588  DOI: 10.1126/sciadv.1501588
  * [45] Tao Lian and Dake Chen  “An Evaluation of Rotated EOF Analysis and Its Application to Tropical Pacific SST Variability”  In _Journal of Climate_ 25.15, 2012, pp. 5361–5373  DOI: 10.1175/JCLI-D-11-00663.1
  * [46] Roland A. Madden and Paul R. Julian  “Detection of a 40–50 Day Oscillation in the Zonal Wind in the Tropical Pacific” Publisher: American Meteorological Society Section: Journal of the Atmospheric Sciences  In _Journal of the Atmospheric Sciences_ 28.5, 1971, pp. 702–708  DOI: 10.1175/1520-0469(1971)028<0702:DOADOI>2.0.CO;2
  * [47] Nathan J. Mantua and Steven R. Hare  “The Pacific decadal oscillation” Publisher: Springer  In _Journal of oceanography_ 58.1, 2002, pp. 35–44
  * [48] Nathan J. Mantua et al.  “A Pacific Interdecadal Climate Oscillation with Impacts on Salmon Production*” Publisher: American Meteorological Society Section: Bulletin of the American Meteorological Society  In _Bulletin of the American Meteorological Society_ 78.6, 1997, pp. 1069–1080  DOI: 10.1175/1520-0477(1997)078<1069:APICOW>2.0.CO;2
  * [49] Elinor R. Martin, Chris Thorncroft and Ben B.. Booth  “The Multidecadal Atlantic SST—Sahel Rainfall Teleconnection in CMIP5 Simulations” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 27.2, 2014, pp. 784–806  DOI: 10.1175/JCLI-D-13-00242.1
  * [50] Matthew Newman et al.  “The Pacific Decadal Oscillation, Revisited” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 29.12, 2016, pp. 4399–4427  DOI: 10.1175/JCLI-D-15-0508.1
  * [51] G North, T L., R Cahalan and F J.  “Sampling Errors in the Estimation of Empirical Orthogonal Functions”  In _Monthly Weather Review_ 110, 1982  DOI: 10.1175/1520-0493(1982)110<0699:SEITEO>2.0.CO;2
  * [52] David W. Pierce  “The Role of Sea Surface Temperatures in Interactions between ENSO and the North Pacific Oscillation” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 15.11, 2002, pp. 1295–1308  DOI: 10.1175/1520-0442(2002)015<1295:TROSST>2.0.CO;2
  * [53] Eugene M. Rasmusson, Phillip A. Arkin, Wen-Yuan Chen and John B. Jalickee  “Biennial variations in surface temperature over the United States as revealed by singular decomposition”  In _Monthly Weather Review_ 109.3, 1981, pp. 587–598
  * [54] J.. Ratnam, S.. Behera, Y. Masumoto and T. Yamagata  “Remote effects of El Niño and Modoki events on the austral summer precipitation of southern Africa”  In _Journal of Climate_ 27.10, 2014, pp. 3802–3815
  * [55] Michael B. Richman  “Rotation of principal components” Publisher: Wiley Online Library  In _Journal of climatology_ 6.3, 1986, pp. 293–335
  * [56] Niclas Rieger  “xmca v0.3.3: Maximum Covariance Analysis for Climate Science”, 2021  DOI: 10.5281/zenodo.4751735
  * [57] C.. Ropelewski and M.. Halpert  “North American Precipitation and Temperature Patterns Associated with the El Niño/Southern Oscillation (ENSO)” Publisher: American Meteorological Society Section: Monthly Weather Review  In _Monthly Weather Review_ 114.12, 1986, pp. 2352–2362  DOI: 10.1175/1520-0493(1986)114<2352:NAPATP>2.0.CO;2
  * [58] R. Saravanan and Ping Chang  “Interaction between Tropical Atlantic Variability and El Niño–Southern Oscillation” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 13.13, 2000, pp. 2177–2194  DOI: 10.1175/1520-0442(2000)013<2177:IBTAVA>2.0.CO;2
  * [59] Jacques Servain, Ilana Wainer, Julian P. McCreary and Alain Dessier  “Relationship between the equatorial and meridional modes of climatic variability in the tropical Atlantic” _eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/1999GL900014  In _Geophysical Research Letters_ 26.4, 1999, pp. 485–488  DOI: https://doi.org/10.1029/1999GL900014
  * [60] Adrian Simmons et al.  “Low frequency variability and trends in surface air temperature and humidity from ERA5 and other datasets”, 2021  DOI: 10.21957/ly5vbtbfd
  * [61] Andréa S. Taschetto and Matthew H. England  “El Niño Modoki impacts on Australian rainfall”  In _Journal of Climate_ 22.11, 2009, pp. 3167–3174
  * [62] Renata G. Tedeschi, Iracema FA Cavalcanti and Alice M. Grimm  “Influences of two types of ENSO on South American precipitation” Publisher: Wiley Online Library  In _International Journal of Climatology_ 33.6, 2013, pp. 1382–1400
  * [63] L. Thurstone  “The simple structure concept” Publisher: The University of Chicago Press Chicago  In _Multiple Factor Analysis: A Development and Expansion of The Vectors of Mind_ , 1947, pp. 319–346
  * [64] Stefan Van Der Walt, S. Colbert and Gael Varoquaux  “The NumPy array: a structure for efficient numerical computation” Publisher: IEEE  In _Computing in science & engineering_ 13.2, 2011, pp. 22–30
  * [65] Hrishikesh D. Vinod  “Canonical ridge and econometrics of joint production” Publisher: Elsevier  In _Journal of econometrics_ 4.2, 1976, pp. 147–166
  * [66] H. Vittal, Gabriele Villarini and Wei Zhang  “Early prediction of the Indian summer monsoon rainfall by the Atlantic Meridional Mode”  In _Climate Dynamics_ 54.3, 2020, pp. 2337–2346  DOI: 10.1007/s00382-019-05117-0
  * [67] John M. Wallace and Robert E. Dickinson  “Empirical Orthogonal Representation of Time Series in the Frequency Domain. Part I: Theoretical Considerations” Publisher: American Meteorological Society  In _Journal of Applied Meteorology_ 11.6, 1972, pp. 887–892  DOI: 10.1175/1520-0450(1972)011<0887:EOROTS>2.0.CO;2
  * [68] Na Wen, Zhengyu Liu and Laurent Li  “Direct ENSO impact on East Asian summer precipitation in the developing summer”  In _Climate Dynamics_ 52.11, 2019, pp. 6799–6815  DOI: 10.1007/s00382-018-4545-0
  * [69] Marc Wiedermann, Jonatan F. Siegmund, Jonathan F. Donges and Reik V. Donner  “Differential Imprints of Distinct ENSO Flavors in Global Patterns of Very Low and High Seasonal Precipitation” Publisher: Frontiers  In _Frontiers in Climate_ 3, 2021  DOI: 10.3389/fclim.2021.618548
  * [70] Shang-Ping Xie et al.  “Global Warming Pattern Formation: Sea Surface Temperature and Rainfall” Publisher: American Meteorological Society Section: Journal of Climate  In _Journal of Climate_ 23.4, 2010, pp. 966–986  DOI: 10.1175/2009JCLI3329.1
  * [71] Chidong Zhang  “Madden-Julian Oscillation” eprint: https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2004RG000158  In _Reviews of Geophysics_ 43.2, 2005  DOI: 10.1029/2004RG000158