=Paper=
{{Paper
|id=Vol-1649/102
|storemode=property
|title=Statistical Modelling in Climate Science
|pdfUrl=https://ceur-ws.org/Vol-1649/102.pdf
|volume=Vol-1649
|authors=Nikola Jajcay, Milan Paluš
|dblpUrl=https://dblp.org/rec/conf/itat/JajcayP16
}}
==Statistical Modelling in Climate Science ==
https://ceur-ws.org/Vol-1649/102.pdf
ITAT 2016 Proceedings, CEUR Workshop Proceedings Vol. 1649, pp. 102–109
http://ceur-ws.org/Vol-1649, Series ISSN 1613-0073, c 2016 N. Jajcay, M. Paluš
Statistical modelling in climate science
Nikola Jajcay1,2 and Milan Paluš1
1 Dept. of Nonlinear Dynamics and Complex Systems, Institute of Computer Science, Academy of Sciences of the Czech Republic
2 Dept. of Atmospheric Physics, Faculty of Mathematics and Physics, Charles University in Prague
Abstract: When it comes to modelling in atmospheric and the model errors [3] and these are intrinsic. The problem
climate science, the two main types of models are taken with initial errors is usually tackled by considering an en-
into account – dynamical and statistical models. The for- semble of model forecasts (instead of just one realization -
mer ones have a physical basis: they utilize discretized dif- integration from single initial state), starting with slightly
ferential equations with a set of conditions (boundary con- different initial conditions. The model errors are intrinsi-
ditions + present state as an initial condition) and model cally connected with the exponential error growth emerg-
the system’s state by integrating the equations forward in ing from the chaotic behaviour related to nonlinearities in
time. Models of this type are currently used e.g. as a nu- discretized equations [4]. This limits the predictability of
merical weather prediction models. The statistical models such GCMs to 6-10 days maximum (e.g. [5]).
are considerably different: they are not based on physical
mechanisms underlying the dynamics of the modelled sys-
tem, but rather derived from the analysis of past weather
patterns. An example of such a statistical model based 1.1 Statistical models
on the idea of linear inverse modelling, is examined for
modelling the El Niño – Southern Oscillation phenomenon
with a focus on modelling cross-scale interactions in the The second kind of models used in climate science are sta-
temporal sense. Various noise parameterizations and the tistical models. In their design, they are considerably dif-
possibility of using a multi-variable model is discussed ferent than the dynamical models in the sense that they
among other characteristics of the statistical model. The are not based on physical mechanism underlying the dy-
prospect of using statistical models with low complexity namics of the modelled system, but rather derived from
as a surrogate model for statistical testing of null hypothe- the analysis of past weather patterns. Probably the most
ses is also discussed. used concept is that of inverse stochastic model [6], where
the model is designed, then estimated using past data and,
finally, stochastically integrated forward in time to obtain
1 Modelling in climate science the prediction. The disadvantages connected to this type of
models consist of the selection of variables that capture the
Climate models, which rely on the use of quantitative system we are trying to model. Other possible issue could
methods to simulate interactions in the climate system, are be the non-stationarity of the modelled system - since the
one of the most important tools to predict and asses future statistical model does not involve the underlying physi-
climate projections or to study the climate of the past. In cal mechanisms, just the interaction between subsystems
general, two types of models are mainly used: dynamical (ignoring hidden variables), the model estimated on some
models and statistical models. The base for a dynamical subset of the past data may not correctly capture all possi-
model is a set of discretized differential equations which ble states of the system. In other words, the training period
are integrated forward in time from the present state, pos- of the past data used to estimate the statistical model may
ing as an initial condition. The most prominent example of not capture the full phase space of the modelled system.
the usage of dynamical models is without doubt a general
circulation model (GCM hereafter). It employs a mathe- The motivation for building a statistical model for par-
matical model of circulation of the planetary atmosphere ticular phenomenon, apart from its forecasting, would be
and oceans, therefore it uses the Navier-Stokes equations to scale down the complexity of the problem. When we
on a rotating sphere (describing a motion of viscous fluid) find some e.g. nonlinear interactions in the observed data,
with thermodynamic terms for energy sources and sinks. and we are interested in uncovering the mechanisms, con-
The above described model is used in numerical weather structing a models of different complexity and seeking
prediction, to infer the reanalysis datasets of the past cli- such interactions in them would help to expose the mech-
mate and for future climate projections in climate model anisms and shed some light on the problem.
intercomparison projects CMIP3 [1] and CMIP5 [2]. In the following sections, the inverse stochastic model
The uncertainties of the forecast arisen from the GCM for forecasting the El Niño - Southern Oscillation (ENSO
models are usually classified into two types: the first one is hereafter) phenomenon is built following [7], with the fo-
related to the initial errors (errors in determining the “true” cus on various noise parametrizations and possible use of
present state of the climate), while the second one is due to multiple variables.
Statistical Modelling in Climate Science 103
Figure 1: ENSO phenomenon, its phases and mechanisms: (left) neutral, (center) positive and (right) negative. Figures
taken from [10].
2 Data-based ENSO model (Fig. 1 right) is acting to strengthen the Walker circula-
tion, to move the area of persistent precipitation even more
The ENSO phenomenon exhibits strong interannual cli- westward, the differences in surface pressure is now larger
mate signal and has a great economic and societal impact. and the thermocline is even more tilted. The ENSO tends
It originates from the coupled ocean-atmosphere dynam- to naturally oscillate between these three phases without a
ics of the tropical Pacific [8], but has a strong influence on distinct period (there is no distinct peak in ENSO signal’s
circulation and air-sea interaction also outside the tropical spectrum) and the reasons why are still largely unknown.
belt through teleconnections associated with it [9]. The important aspect of ENSO is that its positive phase
The ENSO phenomenon expresses itself as a sea surface - El Niño is generally characterized by a larger magnitude
temperature (SST hereafter) anomaly and exists in three than its negative phase - La Niña. This statistical skewness
distinct phases - the neutral, positive (El Niño) and nega- is one of the indicators that, at least to some extent, the dy-
tive (La Niña). The basic physical mechanisms for each of namics of ENSO involves nonlinear processes [11]. At the
the phases are depicted in Fig. 1. The normal state of the same time, the most detailed numerical dynamical mod-
equatorial Pacific (Fig. 1 left) is warm SST in the western els seem to severely underestimate this nonlinearity [12],
basin, near Australia and cold SST in the eastern basin, hence the quality of the forecast is not satisfactory.
near the coast of Peru. Above the warm water in the west, From the reviews of statistical models for ENSO fore-
the deep convection takes place, where warm and moist air casting before 2000 [13] it was clear, that majority of mod-
is ascending to the border of troposphere, creating an area els were still linear, but lately the nonlinear models are
of low atmospheric pressure and area of persistent precip- getting more attention (e.g. [14]). In the following, we de-
itation. From the upper part of the troposphere, the air is scribe easy-to-interpret nonlinear model for ENSO fore-
moving eastward and then it descends already as cold and casting.
dry, creating an area of high atmospheric pressure above
the eastern equatorial Pacific. From this basin, the air is 2.1 Inverse models
blowing westward on the surface, in agreement with the The concept of inverse stochastic models are used as the
trade winds, finishing the circulation loop known as the starting point in developing the ENSO model. Let x(t) be
Walker circulation. The easterly surface air flow triggers the state vector of anomalies, so x(t) = X(t) − X, where
the oceanic surface current to flow poleward, effectively X(t) is the climate state vector (could be multi- or univari-
removing water from the surface, thus the water needs to ate climate observations e.g. temperature, pressure etc.
be replaced and this is due to the upwelling, where in the or a PCA time series from eigen-decomposition of some
equatorial area, the water is upwelled from roughly 50 me- climate field) and X is its time-mean. The evolution of
ters depth to the surface. Since the thermocline (a border anomalies could be expressed as
between cold deep ocean and warm surface ocean) is lo-
cated below 50 meters in the west, the upwelled water is ẋ = Lx + N(x) (1)
warm, but in the eastern Pacific the thermocline level is
where L is a linear operator, N represents the nonlinear
above the 50 meters, thus the upwelled water is cold, cre-
terms and dot denotes time derivative.
ating the cold SST in the east and warm SST in the west.
The simplest type of inverse models is linear inverse
The warm phase of ENSO (Fig. 1 center) creates a warm
models (LIM, [6]). By assuming, in eq. (1), that N(x)dx ≈
SST anomaly in the eastern Pacific, acting to weaken the
Tx dt + dr(0) , where T is the matrix describing linear feed-
Walker circulation, to move the area of persistent precip-
backs of unresolved (hidden) processes on x and dr(0) is a
itation eastward, to diminish the difference between east-
white-noise process, eq. (1) could be written as
ern and western Pacific surface pressures and to level off
the thermocline. Reversely, the negative phase of ENSO dx = B(0) xdt + dr(0) , B(0) = L + T. (2)
104 N. Jajcay, M. Paluš
The matrix B(0) and the covariance matrix of the noise The eqs. (4) and (5) describe a wide variety of processes
Q ≡ hr(0) r(0)T i can be directly estimated from the ob- in a fashion that explicitly accounts for the modeled pro-
served statistics of x by multiple linear regression [15]. cess x feeding back on the noise statistics. The linear mul-
The state vector x, or predictor-variable vector, consists of tilevel model is obtained by assuming Ai ≡ 0 and c(0) ≡ 0
amplitudes of corresponding principal components (PCA in eq. (4). Details of the methodology and further discus-
analysis [16] yields spatial patterns - empirical orthogonal sion could be found in [7].
functions and its respective time series - principal com- It is well known, that the extreme ENSO events tend to
ponents), while the vector of response variables contains occur in boreal winter. From several ways to include this
their tendencies ẋ. phase locking to the annual cycle, the alternative approach
used here is to include seasonal dependence in the dynam-
ical part of the first level. Namely, we assume the matrix
2.2 Nonlinear multilevel model
B(0) and vector c(0) to be periodic, with period T = 12
The assumptions of linear, stable dynamics and of additive months:
white-noise used to construct LIMs are only valid to cer- B(0) = B0 + Bs sin(2πt/T ) + Bc cos(2πt/T ),
tain degree of approximation. In particular, the stochastic (0)
forcing dr(0) typically involves serial correlations, and, in c = c0 + cs sin(2πt/T ) + cc cos(2πt/T ) (6)
addition, the matrices B(0) and Q obtained from the data In this case, the whole record is used to estimate four
exhibit substantial dependence on the lag, that was used to seasonal-dependent coefficients. The model is trained
fit them [17]. The two modifications of the basic inverse in the leading EOF (empirical orthogonal function)
model, that address both nonlinearity and serial correla- space [16] of tropical Pacific SST anomalies. The opti-
tions are taken into account, as in [18]. mal number of state-vector components and the degree of
The first modification is obtained by assuming polyno- nonlinearity has to be assessed by cross-validation. The
mial, rather than linear form of N(x) in eq. (1), in par- parameters in this paper were used as in [7].
ticular, a quadratic dependence. The ith component Ni (x)
could be written as
3 Results
(0) (0)
Ni (x) ≈ xT Ai x + ti x + ci dt + dri (3)
In this section, the brief results are presented of how the
The matrices Ai represent the blocks of a third-order ten- statistical model is able to simulate the ENSO signal. The
(0)
sors, while the vectors bi = li + ti are the rows of the skill of the model is determined in the sense of basic linear
matrix B = L + T (as in eq. (2)). These objects, as well
(0) ENSO metrics such as the amplitude of the ENSO signal,
as components of the vector c(0) , are estimated by multiple the seasonality (since the seasonality is important aspect
polynomial regression [19]. of ENSO dynamics) and finally, the power spectrum of
ENSO signal. The model is employed as described in the
The second modification, considering the serial correla-
previous section, the matrices and vectors are estimated
tions in residual forcing, is due to the multilevel structure
from the previous data and then the model is integrated to
of our model. In particular, consider the ith component of
obtain the time series of same length as the training data.
the first, main level of the inverse stochastic model
Since the model is stochastic (forced by a white noise),
(0) (0) (0) we employed an ensemble of 20 members. Each member
dxi = xT Ai x + bi + ci dt + dri , (4)
is integrated with slightly different initial conditions and
these members are referred to as realizations.
where x = {xi } is the state vector and matrices Ai , vec-
(0) (0) The basic ENSO metric is its amplitude, which could
tors bi and the components ci of the vector c(0) as well be characterized by the standard deviation of SST anoma-
(0)
as the components ri of the residual forcing vector r(0) lies averaged over Nino3.4 box (bounded by 5◦ S - 5◦ N
are determined by the least squares. The additional model and 120◦ W - 170◦ W). In Fig. 2 we can see the ENSO
level is added to express the known increments dr(0) as a amplitude as derived from the Nino3.4 index [20] (thick
linear function of an extended state vector [x, r(0) ]. We black line), along with 20 realizations from the data-based
estimate this level’s residual forcing again by the least ENSO model, both linear and quadratic (gold for linear,
squares. More levels are added the same way, until the red for quadratic).
Lth level’s residual, r(L+1) , becomes white in time, and its As can be seen, the linear model slightly overestimates
lag-0 correlation matrix converges to constant, hence the ENSO amplitude, while the quadratic model slightly
underestimate the ENSO amplitude. From the spread of
(0) (1) (1)
dri = bi [x, r(0) ]dt + ri dt, the ensemble members we could infer that the model is
(1) (2) (2) sensitive to initial conditions and the forcing. Still, the
dri = bi [x, r(0) , r(1) ]dt + ri dt,
ensemble averages for both models are within reasonable
... distance from the data borderline, therefore in this aspect
(L) (L+1) (L+1)
dri = bi [x, r(0) , . . . , r(L) ]dt + ri dt (5) the model performs adequately.
Statistical Modelling in Climate Science 105
0.95
ENSO amplitude tive models are more flat in this area of frequencies. In
0.90
the higher frequencies (around annual frequency and less)
the power spectra are in agreement. In general, the spectra
0.85
of modelled time series could be said to copy the actual
STD of SSTA -- Nino3.4
0.80 Nino3.4 time series. The power spectra were computed
0.75 using the Welch method [21].
0.70
10 8
0.65
0.60
10 7
linearan linear
linear model 1
linear model 2
linear model 3
linear model 4
linear model 5
linear model 6
linear model 7
linearar model 98
linear model 10
linear model 11
linear model 12
linear model 13
linear model 14
linear model 15
linear model 16
linear model 17
linear model 18
linear model 19
l 20
quadratic mod .
quadratic model 1
quadratic model 2
quadratic model 3
quadratic model 4
quadratic model 5
quadratic model 6
quadratic model 7
tic mo el 8
quadratic moddel 9
quadratic model 10
quadratic model 11
quadratic model 12
quadratic model 13
quadratic model 14
quadratic model 15
quadratic model 16
quadratic model 17
quadratic model 18
tic mo el 19
del 20
quadra ean quad
line model
mode
me
m
quadra
10 6
power
Figure 2: ENSO amplitude as standard deviation of SST
10 5
anomalies in data (black line) and in 20 realizations of lin-
ear (gold) and quadratic (red) model.
10 4
Other metric connected with ENSO amplitude is its sea-
sonality. As written above, the ENSO phenomenon ex- 10 3
10 -1 10 0
hibits seasonal changes in variance, with elevated variance period [1/yr]
in winter months and lower variance in spring and summer
months. This can be also seen in Fig. 3, where the monthly
Figure 4: ENSO power spectra estimated using the Welch
variance is plotted for the data and for both models. Both
method in data (black curve) and in 20 realizations of lin-
models are capable of modelling higher variance in winter
ear (gold) and quadratic (red) model. Thicker lines repre-
months and drop in variance through spring and summer,
sent the mean over 20 realizations in the respective model.
although the difference in variance is higher in data than in
both models. Still, the ensemble averages are reasonably
close to the data.
4 Noise parametrization in the model
1.2
The statistical model, once estimated, is integrated for-
1.1
ward in time and forced by a noise - usually a realization
1.0 of spatially correlated random process. In the most in-
tuitive and basic case, the last level residuals’ covariance
0.9 matrix is estimated and decomposed using Cholesky fac-
STD of SSTA
0.8
torization yielding a lower triangular matrix R. When the
model is integrated, the random realization of white noise
0.7 is multiplied by the matrix R, yielding spatially corre-
lated white noise which is used as a random forcing in the
0.6
model. The results for quadratic and linear ENSO models
0.5 from the previous section were obtained using this simple
noise parametrization, and the question is whether looking
0.4 deeper into the residuals’ structure could aid the model’s
Jan
b
r
r
y
Jun
Jul
g
p
t
v
c
performance.
Oc
Ma
Ap
De
Ma
No
Fe
Au
Se
month
Figure 3: ENSO seasonality as standard deviation per 4.1 Dependence on the system’s state
month in data (black curve) and in 20 realizations of linear First refinement for the noise parametrization arises from
(gold) and quadratic (red) model. Thicker lines represent the concept of modelling climate processes which exhibit
the mean over 20 realizations in the respective model. low-frequency variability (LFV). In this method, we find
and select noise samples, snippets, from the past noise
The last metric taken into account was the power spec- (residuals) which have forced the system during short time
trum of Nino3.4 time series. The spectrum for the Nino3.4 intervals that resemble the LFV phase just preceding the
data and both linear and quadratic model realizations can currently observed state, and then use these snippets (or in-
be seen in Fig. 4. The main peak in data occurs at roughly formation contained in them) to drive the current state into
5 year period, but still the ensemble averages for respec- the future. For full methodology and discussion, see [22].
106 N. Jajcay, M. Paluš
The found past noise snippets can be used in two differ- 0.90
ENSO amplitude
ent ways. The first one (as used in [22]) seeks various snip-
pets from the past observations and then directly uses them
0.85
to force the model as an ensemble. When e.g. we find 4
STD of SSTA -- Nino3.4
intervals which resemble the LFV phase, we integrate the
model 4 times using all 4 noise snippets directly and than 0.80
average over them. The second version (as used in our
study) is to find, say, 100 samples of the past noise clos-
0.75
est to the current state of the system, cluster them together
and create covariance matrix from them. Afterwards, the
linlienaerar
linear 1
linear 2
linear 3
linear 4
linear 5
linear 6
linear 7
line 8
linearar 9
linear 10
linear 11
linear 12
linear 13
linear 14
linear 15
linear 16
linear 17
linear 18
linear 19
20
conditan cond.
conditioned 1
conditioned 2
conditioned 3
conditioned 4
conditioned 5
conditioned 6
conditioned 7
con ioned 8
conditditioned 9
conditioned 10
conditioned 11
conditioned 12
conditioned 13
conditioned 14
conditioned 15
conditioned 16
conditioned 17
conditioned 18
conditioned 19
20
Cholesky decomposition is used to obtain the matrix R and
ioned
mean
finally, the random white noise realization is multiplied by
me
the matrix R. Using this matrix, the spatial covariance of 1.2
the forcing is dependent on the current state of the sys-
tem. In both noise parametrizations, the current system 1.1
state could be estimated in multiple ways: either using
correlation of the SSA time series, or using the Euclidean 1.0
distance in the subspace spanned by first few EOFs.
As can be seen in Fig. 5, although the amplitude statis-
STD of SSTA
0.9
tics are not substantially shifted, the transient from high-
variance winter period to low-variance spring and summer 0.8
are better captured by the later model, with noise forcing
0.7
conditioned on system’s state. The power spectra for both
models are practically the same (not shown).
0.6
4.2 Seasonal dependence of the forcing 0.5
Jan
b
r
r
y
Jun
Jul
g
p
t
v
c
Oc
Ma
Ap
De
Ma
No
Fe
Au
Se
Although the seasonal dependence of the model is cap- month
tured in model’s dynamics by fitting the seasonally depen-
dent matrices B(0) and c(0) (recall eq. (6)), our analysis Figure 5: ENSO amplitude (upper) and seasonality (bot-
showed, that the last level’s residuals still exhibit season- tom) in data (black curve) and in 20 realizations of linear
ally dependent amplitude. To address this issue, we com- (gold) and linear with conditioned noise on the system’s
puted the standard deviations for each month from the last state (red) model. Thicker lines represent the mean over
level’s residuals, then fitted the 5 harmonics of the annual 20 realizations in the respective model.
cycle to capture the seasonal dependence, removed this
dependence from the residuals, then estimated covariance
matrix and subsequently the matrix R and finally gener- The two latter modifications bring just a slight improve-
ated spatially correlated white noise realization which was ments into ENSO metrics (not shown), but could have
multiplied back by the requisite seasonal amplitude to ac- more substantial advancements in modelling different at-
count for the seasonally dependent amplitude of the forc- mospheric phenomena.
ing. The fitted harmonics of the annual cycle were selected
as
5 Synchronization and causality in the
Pi = cos(2πit/T ) + sin(2πit/T ), i = 1, . . . , 5 (7) observed and modelled data
and then regressed on the seasonally varying standard de-
viation of the last level’s residuals. Better understanding of the complex dynamics of the at-
mosphere and climate is one of the challenges for con-
temporary science. Considering the climate system as
4.3 Using extended covariance matrix
a complex network of interacting subsystems [23] is a
The last modification to the noise is to use the extended co- new paradigm bringing new data analysis methods help-
variance matrix instead of lag-0 covariance matrix. When ing to detect, describe and predict atmospheric phenom-
evaluating system’s state we do not take just the state clos- ena [24]. A crucial step in constructing climate networks
est to the current state of the model, but, say 5 consecu- is inference of network links between climate subsys-
tive months and construct the extended matrix out of this tem [25]. Directed links determine which subsystems in-
snippet. Then the matrix is decomposed using Cholesky fluence other subsystems, i.e. uncover the drivers of at-
factorization and used as a spatial correlation matrix R is mospheric phenomena. Inference of causal relationships
random forcing generation. from climate data is an intensively developing research
Statistical Modelling in Climate Science 107
field, e.g. [26, 27]. Typically, a causal relation is sought
between different variables or modes of atmospheric vari-
ability.
Paluš [28] has open another view at the complexity
of atmospheric dynamics by uncovering causal relations
or information flow between dynamics on different time
scales in the same variable. Recently, phase-phase and
also phase-amplitude interactions between dynamics on
different temporal scales were observed in the ENSO dy-
namics (captured by the Nino3.4 index) using the ap-
proach as in [28]. Shortly, we use the continuous wavelet
transform to the time series for particular time scales to
Figure 6: Phase synchronization (left) and phase-
obtain the instantaneous phase and amplitude of the oscil-
amplitude causality (right) in Nino3.4 time series. Shown
latory mode as
is the significance (over 95th percentile against 500 Fourier
ψ(t) = s(t) + iŝ(t) = A(t)e(iφ (t)) , (8) transform surrogates) of k-nearest neighbours estimate of
mutual information and conditional mutual information.
ŝ(t)
φ (t) = arctan , (9)
s(t)
q to uncover the mechanisms of these interactions and shed
A(t) = s2 (t) + ŝ2 (t). (10) more light onto the dynamics of ENSO in general. We
constructed the ENSO model and repeated the above anal-
Then the time series of phase and / or amplitude are used
ysis to modeled ensemble of the Nino3.4 time series.
to study the interactions. We adopt measures from infor-
mation theory, namely mutual information and conditional
mutual information, where the mutual information could
be expressed as
p(x, y)
I(X;Y ) = ∑ ∑ p(x, y) log , (11)
x∈X y∈Y p(x)p(y)
where p(·) is the probability distribution or joint probabil-
ity distribution and X and Y are our time series of either
phase or amplitude derived from the ENSO SST data. Fi-
nally, the measures we are interested in could be written
as: Figure 7: Phase synchronization (left) and phase-
amplitude causality (right) in modeled Nino3.4 time series
• phase synchronization – I(φ1 (t); φ2 (t)), by the data-based model. Shown is the aggregate of 5 real-
• phase-phase causality – I(φ1 (t); φ2 (t + τ) − izations of k-nearest neighbours estimate of mutual infor-
φ2 (t)|φ2 (t)), mation and conditional mutual information. Significance
against 500 Fourier transform surrogate data.
• phase-amplitude causality – I(φ1 (t); A2 (t +
τ)|A2 (t), A2 (t − η), A2 (t − 2η)), As seen from the analysis of modelled data (Fig. 7),
the main phase synchronization bands (annual cycle with
quasi-biennal cycle and combination frequencies) are also
5.1 Interactions in the data
captured by the modelled data, while the phase - ampli-
As can be seen from Fig. 6, in ENSO dynamics captured tude interactions are not very well captured. This might
by the Nino3.4 index, the synchronization of annual cycle arise from the low complexity of the model, or the absence
with quasi-biennal and combination frequencies (frequen- of some nonlinear interactions in the model design (apart
cies that arise from the interactions between annual and from quadratic).
the most prominent ENSO period) is observed. Also, the
4-6 year cycle of phase in ENSO dynamics influence the
quasi-biennal range of the amplitude time series. 6 Modelling surrogate data with statistical
model
5.2 Interactions in the model
Surrogate data (or analogous data) is a method to generate
Our goal was to simulate the nonlinear cross-scale inter- synthetic data set (time series) that preserve some of the
actions in the model. This is important since it might help statistical properties, while omitting the others. One way
108 N. Jajcay, M. Paluš
of using them, is to test statistical significance by contra-
diction. This involves posing a null hypothesis describing
some kind of a process and then generating an ensemble
of surrogate data according to null hypothesis using Monte
Carlo methods. One of the most used technique for gener-
ating surrogate data is the Fourier transform surrogate [29]
(FT surrogates), which preserve the linear correlations in
the data (periodogram or spectrogram, including autocor-
relation) of the time series, but omits any other interactions
in them.
As an example, consider two intertwined Lorenz sys-
tems, where one of them drives the other. Now, using the
Figure 8: Phase synchronization (left) and phase-
time series in one dimension, say the x dimension from
amplitude causality (right) in modelled Nino3.4 time se-
both Lorenz systems, we can use some method for detect-
ries by the data-based model. Shown is the aggregate of
ing causality, e.g. conditional mutual information between
5 realizations of k-nearest neighbours estimate of mutual
the two time series of two Lorenz systems. We get the
information and conditional mutual information. Signifi-
value of conditional mutual information, but this is still
cance against 500 surrogate time series created with data-
not enough to interpret it in the means of whether there
based model.
is a causal relationship between them or the result arose
by chance. For this purpose, we construct an ensemble of
Fourier transform surrogate data (which qualitatively pre- minimum. This way, we can get better idea of the statisti-
serves properties of the time series, but allows no causal cal significance of the interactions between subsystems, in
relationship between them) and repeat the analysis using particular the nonlinear ones, since we are testing against
the very same method on this ensemble and finally com- the model with just linear interactions.
pare the value for actual data with the histogram of values
obtained from the ensemble of surrogate data. When the 7 Conclusions
value from the data exceeds some percentile (e.g. 95th ) of
the surrogate data distribution, we say that the causal re- Statistical modelling in climate science is continuously
lationship is significant in comparison with e.g. 500 FT getting more attention, since their usage is not limited to
surrogates. forecast some of the phenomena of interest (like ENSO),
When studying nonlinear cross-scale interactions in but could also be used to infer some of the statistical prop-
time series using the above method, the statistical test in- erties and relationships among different subsystems. Since
volves creating an ensemble of surrogate, synthetic time the statistical models live in phase space of particularly re-
series and repeat the analysis for the whole ensemble. duced dimensionality, when we could observe the inter-
Then we computed the percentile, where the observed in- actions of interest, the identification of their sources will
teractions could not arose by random chance. Of course, become more feasible.
one could use Fourier transform method to generate the We showed that the statistical model with the right set-
surrogate time series, effectively posing a null hypothe- tings, which were selected based on careful inspection of
sis of a linear process which has the same spectrum to the modelled system, could generate synthetic time series
that of an observed data. On the other hand, one can cre- of interest, copying the desired properties of the system -
ate a more sophisticated null hypothesis by exploiting the both linear and nonlinear statistics. Since the stochastic-
options of a data-based model: when one consider just a ity is the important aspect of the data-based model, var-
linear model, omit the dynamical seasonal dependence in ious parametrization techniques exist to correctly model
B(0) and c(0) terms (as in eq. (6)) and use the simplest noise the system’s external forcing. Finally, the possibility of
parametrisation (just consider the spatial covariance struc- usage of the low complexity model as surrogate data was
ture), the model will omit the nonlinear interactions and discussed, showing advantages of usage of such technique
could pose as a surrogate data model copying the basic to infer statistical significance.
statistical properties of a modelled time series. This way, The outlook for future work combines various differ-
the analysis would show whether the cross-scale interac- ent paths which appeared. One direction would be fo-
tions are arising from the seasonal dependent dynamics, or cusing on statistical modelling itself, experimenting with
from nonlinear (e.g. quadratic) interactions between sub- various variable model, with input time series and their
systems and so on. preprocessing and so on and so forth. Other direction
When comparing Fig. 6 (testing against 500 Fourier would be connecting the statistical models with dynami-
transform surrogates) and Fig. 8 (testing against 500 data- cal ones, in the sense, that statistical models could be used
based model surrogates), the significant interactions are for parametrization of e.g. sub-grid phenomena (micro-
virtually the same, expect in the latter, the “fluctuations” physics of clouds, local convection etc.) in large coupled
(or they might be false positives as well) are attenuated to atmospheric-oceanic models.
Statistical Modelling in Climate Science 109
References [19] McCullagh, P., and J. A. Nelder: Generalized Linear Mod-
els. Chapman and Hall (1989)
[20] Rayner N. A., D. E. Parker, E. B. Horton, C. K. Folland,
[1] Meehl, G. A., C. Covey, T. Delworth, M. Latif, B. McA- L. V. Alexander, D. P. Rowell, E. C. Kent and A. Kaplan:
vaney, J. F. B. Mitchell, R. J. Stouffer, and K. E. Taylor: Global analyses of sea surface temperature, sea ice, and
The WCRP CMIP3 multi-model dataset: A new era in cli- night marine air temperature since the late nineteenth cen-
mate change research. A Bull. Amer. Meteor. Soc. 88 (2007) tury. J. Geophys. Res. 108 (2003) 4407
1383–1394 [21] Welch, P. D.: The use of Fast Fourier Transform for the
[2] Taylor, K.E., R.J. Stouffer and G.A. Meehl: An Overview estimation of power spectra: A method based on time aver-
of CMIP5 and the experiment design. A Bull. Amer. Meteor. aging over short, modified periodograms. IEEE Trans. Audio
Soc. 93 (2012) 485–498 AU-15 (1967) 70–73
[3] Bjerknes, V.: Dynamic meteorology and hydrology, Part [22] Chekroun, M. D., D. Kondrashov and M. Ghil: Predicting
II. Kinematics. Gibson Bros., Carnegie Institute, New York. stochastic systems by noise sampling, and application to the
(1911) El Niño-Southern Oscillation. P. Natl. Acad. Sci. USA 108
[4] Lorenz, E. N.: Deterministic nonperiodic flow. J. Atmos. Sci. (2011) 11766–11771
20 (1963) 130–141 [23] A. A. Tsonis and P. J. Roebber: The architecture of the
[5] Van den Dool, H. M.: Long-range weather forecasts through climate network. Physica A: Statistical Mechanics and its
numerical and empirical methods. Dyn. Atmos. Oceans 20 Applications, 333 (2004) 497–504
(1994) 247–270 [24] S. Havlin, D. Y. Kenett, E. Ben-Jacob, A. Bunde, R. Cohen,
[6] Penland, C.: Random forcing and forecasting using princi- H. Hermann, J. Kantelhardt, J. Kertész, S. Kirkpatrick, J.
pal oscillation pattern analysis. Mon. Weat. Rev. 117 (1989) Kurths, et al.: Challenges in network science: Applications
2165–2185 to infrastructures, climate, social systems and economics.
[7] Kravtsov, S., D. Kondrashov, and M. Ghil: Multilevel regres- The European Physical Journal Special Topics, 214 (2012)
sion modeling of nonlinear processes: Derivation and appli- 273–293
cations to climate variability. J. Climate 18 (2005) 4404— [25] M. Paluš, D. Hartman, J. Hlinka, and M. Vejmelka: Dis-
4424 cerning connectivity from dynamics in climate networks.
[8] Philander, S. G. H.: El Niño, La Niña, and the Southern Nonlinear Processes in Geophysics 18 (2011) 751–763
Oscillation. Academic Press (1990) [26] Ebert-Uphoff and Y. Deng: Causal discovery for climate re-
[9] Alexander, M. A., I. Bladé, M. Newman, J. R. Lanzante, N.- search using graphical models. Journal of Climate 25 (2012)
C. Lau, and J. D. Scott: The atmospheric bridge: The in- 5648–5665
fluence of ENSO teleconnections on air–sea interaction over [27] Y. Deng and I. Ebert-Uphoff: Weakening of atmospheric
the global oceans. J. Climate 15 (2002) 2205–2231 information flow in a warming climate in the community cli-
[10] WIKIPEDIA . ORG: El Niño–Southern Oscillation mate system model. Geophysical Research Letters 41 (2014)
https://en.wikipedia.org/wiki/El_Ni%C3%B1o% 193–200
E2%80%93Southern_Oscillation, downloaded June 27, [28] Paluš, M.: Multiscale atmospheric dynamics: Cross-
2016. frequency phase-amplitude coupling in the air temperature.
[11] Ghil, M. and A. W. Robertson: Solving problems with Phys. Rev. Lett. 112 (2014) 1–5
GCMs: General circulation models and their role in the cli- [29] Theiler, J., S. Eubank, A. Longtin, B. Galdrikian, and J.
mate modeling hierarchy. Academic Press, (2000) 285–325 Doyne Farmer: Testing for nonlinearity in time series: The
[12] Hannachi, A., D. B. Stephenson, and K. R. Sperber: method of surrogate data. Physica D 58 (1992) 77–94
Probability-based methods for quanfifying nonlinearity in
the ENSO. Clim. Dyn. 20 (2003) 241–256
[13] Ghil, M. and N. Jiang: Recent forecast skill for the El Niño
/ Southern Oscillation. Geophys. Res. Lett. 25 (1998) 171-
–174
[14] Timmermann, A., H. U. Voss, and R. Pasmanter: Empir-
ical dynamical system modeling of ENSO using nonlinear
inverse techniques. J. Phys. Oceanogr. 31 (2001) 1579–1598
[15] Wetherill, G. B.: Regression Analysis with Applications.
Chapman and Hall (1986)
[16] Hannachi, A., I. T. Jolliffe and D. B. Stephenson: Empirical
orthogonal functions and related techniques in atmospheric
science: A review. Int. J. Climatol. 27 (2007) 1119–1152
[17] Penland, C. and M. Ghil: Forecasting Northern Hemi-
sphere 700-mb geopotential height anomalies using empir-
ical normal modes. Mon. Wea. Rev. 121 (1993) 2355—2372
[18] Kondrashov, D., S. Kravtsov, A. W. Robertson and M. Ghil:
A Hierarchy of Data-Based ENSO Models. J. Climate 18
(2005) 4425–4444