Abstract

Dynamical components of Earth’s ice–ocean–atmosphere system evolve along characteristic trajectories, which make these components partly predictable. This paper reviews several methods for extracting these predictable components from space–time fields. These methods are optimal persistence analysis (OPA), slow feature analysis (SFA), principal trend analysis (PTA), average predictability time decomposition (APTD), and forecastable components analysis (ForeCA). These methods generally find a set of components that are ordered by their predictability, but each method uses a different measure of predictability. Also, a new bootstrap test for investigating the type of predictability exhibited by these components is introduced. This new test is based on an “integrated red noise” hypothesis. The five methods and new test are applied to a dataset of Australian daily near-surface minimum air temperature, spanning 1910–2013. For all five methods, the two leading predictable components are a long-term trend and a low-frequency pattern that decreased in the first half of the twentieth century and increased after that. The third predictable component differs between the methods based on persistence (e.g., OPA) and those based on more general measures of predictability (APTD and ForeCA). In addition, the use of spectral entropy for analyzing time-dependent predictability is investigated. Further research is needed into the application of predictable component methods to specific problems, such as to fields that require regularization (i.e., using ridge regression), to fields with missing values, and to fields with propagating predictable components.

1. Introduction

The analysis of predictable structures in Earth’s ice–ocean–atmosphere system is an important component of climate research (DelSole and Tippett 2009a,b). A system is predictable if , the state of a system at time t, is not independent of , the state(s) at one or more past times. Understanding predictability on intra-annual and interannual time scales is important for resource management, and understanding predictability on multidecadal time scales is important for assessing natural versus forced variation in the climate system. This paper is concerned with the predictability of space–time fields. Several methods for investigating the predictability of space–time fields have been proposed, including optimal persistence analysis (OPA; DelSole 2006), average predictability time decomposition (APTD; DelSole and Tippett 2009b), forecastable components analysis (ForeCA; Goerg 2013), slow feature analysis (SFA; Wiskott and Sejnowski 2002), and principal trend analysis (PTA; Hannachi 2007). The main aim of this paper is to compare these methods by applying them to a daily minimum air temperature field from Australia (1910–2013). The five methods are briefly reviewed in sections 2 and 3, and a new test to examine the type of predictability exhibited by the components is introduced. The Australian daily temperature field is described in section 4, and the five methods and new inferential test are applied to this field in section 5. The results are also compared with the conventional empirical orthogonal function/principal component analysis (EOF/PCA) method (Monahan et al. 2009).

All the methods mentioned above rely on the same general model; that is, each time series in a climate field is the sum of a small number of common signals (or components) (this general model will be detailed in section 2). For different time series in the field, each common signal can have a different amplitude, and its phase is fixed at either 0° or 180° (the latter means the signal is multiplied by −1). Other space–time decomposition methods, such as empirical mode decomposition (EMD; Huang and Wu 2008) and multichannel singular spectral analysis (MSSA; Ghil et al. 2002), use a different model, where each time series in a field is the sum of its own individual components; these individual components do not have to be the same from one time series to the next. For example, in section 1 of the supplementary material the differences between MSSA and the methods used here are described in more detail; owing to these differences, MSSA (and EMD) are not considered in this study. Another method, predictive power analysis (PPA; Schneider and Griffies 1999) is equivalent to performing APTD with a single lag or lead time (see section 3b), so the PPA components will be dependent upon the chosen single lead or lag time. The advantage of the four methods used in this study (OPA, APTD, ForeCA, and PTA) is that these methods find components that are not dependent on a single lead/lag time. The fifth method, SFA, is dependent on a single lead time, but the estimation of its components is derivative-based (which is unique among the methods here); hence, it is used here to provide a useful comparison.

The following section discusses three items that are relevant to all five predictable component methods investigated in this study: the general model, the criterion of interest, and testing the predictability of the components.

2. Review

a. The general model

The general model is to decompose the space–time field into a set of component time series and their associated spatial patterns :

 
formula

The matrix is an field consisting of p spatial points (stations, columns) and n time points (rows); is a matrix of weights; is a matrix of spatial patterns, where ; and is an matrix of residuals (see Table 1 for notation used hereinafter). Equation (1) can be rewritten as follows:

 
formula

where denotes the jth column vector of matrix , and means to stack the column vectors of the respective matrix [e.g., the left-hand side of Eq. (2) could be rewritten as ]. In Eq. (2), each time series in is thought of as being a linear combination of the same set of components . The components each have unit variance, so the column vectors of the matrix contain the amplitudes (linear coefficients) of the components for each time series in . For example, each time series in the near-surface air temperature field may be thought of as being made up of a number of possible signals, including a trend (e.g., global warming), interdecadal variation (e.g., the Atlantic meridional oscillation), interannual variation (e.g., ENSO), intraseasonal variation (e.g., the Madden–Julian oscillation), and so on; but the amplitude (or contribution) of these signals to different time series in can be different. The way in which such signals are defined is often arbitrary (e.g., a trend may be defined as a linear function of time). Alternatively, it should be possible to extract such signals from the data field itself by finding projections that cancel the noise and elucidate the signals. To do this, the signals need to be defined according to some characteristic that makes them “interesting”: for example, a signal’s variance and its predictability are two different properties that may be of interest. Thus, estimating the general model from data requires three steps:

  1. Define a criterion of interest .

  2. Calculate by finding one-dimensional projections of the data that optimize the criterion of interest such that . Analytical solutions for can typically be found by constraining the different signals to be uncorrelated (i.e., ). How to estimate for the five different methods is described in section 3. Section 2 of the supplementary material shows how the different solutions can be derived by optimizing the respective criterion of interest.

  3. Calculate using least squares principles: that is, , since , where is a covariance matrix (of ), and is an identity matrix of rank d.

Table 1.

Notation used.

Notation used.
Notation used.

b. Criterion of interest

There are many possible criteria of interest. For example, for PCA, , where and . However, a signal that explains a lot of variance in a field is not necessarily predictable in time. The interest here is in finding signals that are predictable. Since a signal is predictable if is not independent of , then a suitable may be based on the autocorrelation function . The autocorrelation function describes the correlation between and :

 
formula

where , and τ is the lead or lag time (it can take positive or negative values). Now, needs to be a scalar value, so the simplest that describes a predictable signal would be . Other possibilities are or . Further, since the Fourier transform of the autocorrelation function is the spectrum, other suitable may be derived from the Fourier spectrum. Since all these potential are based on the autocorrelation function (or its Fourier transform), then the leading signals found by optimizing should have nonzero autocorrelation. This is what is meant here by the term “predictable components.” In contrast, choosing a criterion that is based only on or a criterion that is unrelated to the autocorrelation function will not necessarily find components that have nonzero autocorrelation.

c. Predictability of components

It is important to test whether the components are more or less predictable than the components that could be produced from some null hypothesis process. The simplest null hypothesis is a multivariate white noise process. A second null hypothesis is a multivariate red noise process, where the rate of change of a field is forced by a white noise process. An example of this type of process is sea surface temperature driven by atmospheric noise. A third null hypothesis is that the rate of change of a field is forced by another red noise process (Newman et al. 2003; Shakun and Shaman 2009; Di Lorenzo and Ohman 2013). This can occur when the slow components of a system are forced by a process that is already red. An example of this third type of process is sea ice driven by sea surface temperature (which is already red, as in the previous example). A time series produced by this third type of process will have a spectrum that is more red than the spectrum of time series produced by the first two null hypotheses. Employing the third null hypothesis is particularly important for investigating the stochastic nature of a long-term trend or the stochastic nature of a multidecadal oscillation. For example, if a slow component of the climate system is identified using a predictable component method, then it is interesting to ask whether that slow component is a red noise process, or is it due to a process that integrates both atmospheric noise as well as red noise signals from other parts of the climate system that are “slower” than the atmosphere. Thus, the third null hypothesis is a further test about the type of predictability that a predictable component displays.

The third null hypothesis can be expressed as follows:

 
formula

where is the jth column vector of , except in the case of , which means the field lagged by one step in time. The variable is a lag-1 autoregression coefficient matrix; () is a univariate or multivariate time series with a nonwhite spectrum; is a matrix of regression weights; and is a matrix of residuals (not necessarily with a white spectrum). Equation (4) states that each time series in is a linear combination of the previous state of , plus some other forcing and some noise . For example, it is possible that the South Pacific sea surface temperature field () integrates both South Pacific atmospheric anomalies () and tropical Pacific sea surface temperature anomalies () (Shakun and Shaman 2009). In statistical parlance, Eq. (4) is a vector autroregressive model with exogenous variables (VARX), of order 1.

To test the third null hypothesis, it is first necessary to estimate the VARX model from the data, which requires an estimate of . One method of estimating is to assume that relates to a known process (e.g., ENSO) and use an appropriate sea surface temperature index. The matrix could also be estimated directly from the data, for example by extracting using bandpass filtering, assuming that the process operates over particular frequencies. Here, was estimated using the latter method: was filtered using a bandpass Butterworth filter to extract the 2–7-yr signal, and was estimated as the leading three principal components of the filtered field. Three components were chosen by examining the scree plot of eigenvalues, and the leading three principal components explained 66% of the variance in the filtered field. After estimating , the other components of Eq. (4) can be estimated using the standard least squares solution for a VARX model (Tsay 2013).

In this study, the three null hypotheses mentioned above were tested in the following manner.

  1. The first null hypothesis was modified to include both and [see Eq. (4)], so the first null hypothesis is now a test of whether the predictable components have the same predictability as a field with interannual variance plus white noise. This was done for comparative purposes (which will be shown in section 5b). To implement this hypothesis, the model was fitted to the data field , and then the bootstrapped field was calculated as . The realizations were sampled from the estimated multivariate normal distribution of the regression coefficients [i.e., ]. The realizations were generated by randomly selecting (with replacement) the rows of .

  2. The second null hypothesis was implemented by fitting a lag-1 vector autoregressive (VAR) model (without exogenous variables) to the field : that is, . The residual matrix of the VAR model was then randomized by block bootstrapping using a block length of 73 pentads. The was generated using , where and are row vectors of , , is a row vector of , and is the bootstrapped residual matrix from the VAR model. The block bootstrap was used in order to preserve any short-term autocorrelation (pink noise) in the residuals of the VAR model.

  3. The parameters of the third null hypothesis were estimated by fitting the VARX model [Eq. (4)] using , because was already estimated for the first null hypothesis. The bootstrapped field was generated using , where , , and is the block bootstrapped residual matrix from the VARX model. A block length of 73 pentads was also used for this block bootstrap.

For the three tests, SFA was applied to each bootstrapped field , and the statistic [Eq. (12)] was calculated for the first 30 SFA components. This procedure was repeated 3000 times, in order to generate 99% confidence limits for the bootstrapped statistic. The reason that SFA was used here to generate the components is that SFA is a fast method (see section 3d).

3. Methods

a. Optimal persistence analysis

Optimal persistence analysis (DelSole 2006) can be expressed as the reduced rank prediction of a weighted average of :

 
formula
 
formula

where is an element of the weight vector ; and are row vectors of ; and is a row vector of the residual matrix . Section 2 of the supplementary material shows how OPA can be cast in the above form. To calculate , is first obtained from the eigenvalue decomposition: , where and are matrices of eigenvectors and eigenvalues. Then . Note that, in the eigenvalue decomposition step, denotes the matrix produced by filtering with weights [Eq. (5a)]. The weights can be obtained from a suitable taper, such as the Parzen taper:

 
formula

where . The criterion of interest that OPA optimizes is

 
formula

where is the autocorrelation function [Eq. (3)]. OPA decomposes into a set of component time series where . Time series that have a high include red noise processes and time series with low-frequency components (e.g., multidecadal trends). Time series that have a low include white noise, blue noise, and any series with an oscillating autocorrelation function. Note that, for any general index of predictability, a blue noise process should be more predictable than a white noise process. So can distinguish some types of predictability from white noise (such as persistent processes like red noise), but it will not identify all types of predictability (such as blue noise).

b. Average predictability time decomposition

APTD (DelSole and Tippett 2009b) can be expressed as the eigenvalue decomposition of the weighted covariance matrix of the predicted fields over multiple lead times:

 
formula
 
formula

where is a matrix of eigenvectors. For APTD, . Choices for include or a suitable taper, such as the Parzen taper [Eq. (6)]. The latter is preferable, because cannot be reliably estimated for long lead times. The criterion of interest that APTD optimizes is

 
formula

(see section 2 of the supplementary material). Red noise, blue noise, and any series with an oscillating autocorrelation function will have a high . Only time series of white noise (i.e., with a uniform spectrum) will have a low .

c. Forecastable components analysis

Forecastable component analysis (Goerg 2013) makes direct use of the multivariate spectrum :

 
formula

where . The spectrum has dimensions for each frequency ω. ForeCA is based on the eigenvalue decomposition of the weighted multivariate spectrum; the weights essentially filter out particular frequencies. A solution can be found using a fixed-point algorithm:

 
formula
 
formula

where is the trailing column vector of (i.e., the column vector associated with the smallest eigenvalue). The two equations are iterated until converges. Outside the fixed-point algorithm, the trailing column vectors are concatenated into a new matrix . Additional components are found by applying the above algorithm to the orthogonal subspace of : . For each new component, the weights are adjusted as . The weights for the original field are . The criterion of interest that ForeCA optimizes is

 
formula

(see section 2 of the supplementary material), where is the scaled spectral density [i.e., ]. The absolute value of the sum term in Eq. (12) is the discrete spectral entropy. For a white noise process, , whereas, for a time series that is perfectly predictable (e.g., a noiseless sinusoid), .

Also, and are both measures of spectral flatness; this is discussed in section 3 of the supplementary material.

d. Slow feature analysis

Slow feature analysis (Wiskott and Sejnowski 2002) is based on the singular value decomposition of the whitened first difference matrix of :

 
formula

where is the first difference matrix of (i.e., ), and and are the left and right singular vectors. In the above decomposition, the SFA components are the trailing components. The is estimated as . SFA is actually a type of OPA, with the weights [Eq. (5)] equal to

 
formula

Note that, in this case [i.e., Eq. (5) with Eq. (14)], the SFA components are the leading components of the eigenvalue decomposition described in section 3a. In both cases [Eq. (13) or Eq. (5) with Eq. (14)], the criterion of interest that is being optimized is the lag-1 autocorrelation:

 
formula

this is shown in section 2 of the supplementary material. Both blue noise and white noise time series have a low , so, like , is not a general measure of predictability.

SFA is a relatively fast method because it based on a simple subtraction, ; hence, it is used in this study as part of the bootstrap tests of the type of predictability.

e. Principal trend analysis

Principal trend analysis (Hannachi 2007; Fischer and Paterson 2014) is designed to extract low-frequency components from a field. It is recast here as a reduced rank regression model (see section 2 of the supplementary material):

 
formula
 
formula

where and is obtained from the singular value decomposition: ; is a matrix that contains the time positions of the sorted columns of . For example, if , with time positions , then the x-sorted time positions are . The criterion of interest that PTA optimizes is the rank correlation between a time series x and time itself:

 
formula

(see section 2 of the supplementary material). Thus, is a measure of the monotonocity of a time series; the leading component of PTA typically exhibits a strong monotonic trend, while the trailing components are trendless.

4. Data

Daily land surface temperature

Daily minimum temperature (Tmin) and daily maximum temperature (Tmax) are available for 60 sites for the period 1910–2013 from the Australian Climate Observations Reference Network–Surface Air Temperature dataset (ACORN-SAT; Trewin 2013; ABOM 2014). The focus here is on the daily minimum temperature data. In this dataset, there are eight sites that each have over 2000 missing values, so those eight sites were omitted from the analyses in this paper (Burketown, Kerang, Miles, Moree, Nhill, Normanton, Palmerville, and Tibooburra). For the other 52 sites, the number of missing values ranged from 24 days (Melbourne regional office) to 1987 days (Bathurst). Prior to analysis, the data were processed in three steps (in the following order): deseasonalization, calculating pentad (5-day mean) temperature, and imputation of the missing values. First, the data were deseasonalized by subtracting the time mean and the first two harmonics. Pentads were calculated for the deseasonalized data, and a 6-day average was used for the 12th pentad during leap years. Third, missing values in the pentad data were imputed using the data interpolating empirical orthogonal functions (DINEOF) method (Taylor et al. 2013). Last, the data series were scaled to unit variance.

The methods in section 3 were applied to the pentad data matrix of 7592 time points × 52 stations. Note, the data matrix was not first projected into a principal component subspace, because here (also, the results would be the same if the field was projected into its principal component space with all the components retained). The initial analyses (using the methods in section 3) showed that, for two methods (ForeCA and APTD), the third predictable component was heavily weighted by one station (Longreach), suggesting that the temperature series for this station contains some artificial predictability not found in other stations. Thus, the Longreach series was removed prior to the imputation step above, and the DINEOF method was rerun on the matrix of 51 stations.

5. Results

a. Comparison of methods

The five methods in section 3 and PCA (for comparison) were applied to the imputed pentad data matrix of 7592 time points × 51 stations. For APTD, τ was truncated at pentads (i.e., 10 years). For ForeCA, was estimated using Welch’s overlapped segment averaging (WOSA) method, with block size =n/4, and 50% block overlap. In this subsection, the leading principal components and predictable components are compared, followed by a more detailed analysis of the similarities and differences between the leading predictable components.

To compare the leading components from each method, three different measures of persistence or predictability (, , and ) were calculated for the first 10 leading components of each method, by applying Eqs. (7), (12), and (15) to the individual component time series. The results are shown in Fig. 1. The leading two principal components (Fig. 2, top row) explain 45% of the total variance in the minimum temperature field but have decorrelation times of <500 days (PCA2 has days) (Fig. 1, left column). In contrast, for the predictable component methods, the leading two predictable components explain <16% of the field’s total variance but have decorrelation times of >1000 days (Fig. 1). An Australia-wide average time series was calculated (by averaging all 51 time series in the field), and this time series was correlated with the leading component of each method. The Australia-wide average time series has the best correlation with PCA1 () and the worst correlation with APTD1 and ForeCA1 (both ) (note, APTD1, ForeCA1, etc. denote the first component for the respective method). The second and third principal components are not necessarily physically interpretable (their spatial patterns are shown in Fig. 3). Using gridded 3-month mean minimum and maximum temperature anomalies for Australia (1950–94), Jones (1999) showed that the leading five rotated principal components could be better interpreted (than their unrotated counterparts) as the average time series for five different regions in Australia. The spectrum of the Australian average time series (like PCA1) has power at many frequencies (the spectrum of PCA1 is shown in Fig. 4). In contrast, the leading component for all predictable component methods is dominated by low-frequency variance (Fig. 4, left column). The principal components (PCA) are not discussed further in this section.

Fig. 1.

Three different measures of predictability for the first 10 component time series found using five predictable component methods and PCA. Each row corresponds to a different method (the method names are on the right). The three measures of predictability are (left) the decorrelation time, (center) 1-spectral entropy, and (right) the lag-1 autocorrelation.

Fig. 1.

Three different measures of predictability for the first 10 component time series found using five predictable component methods and PCA. Each row corresponds to a different method (the method names are on the right). The three measures of predictability are (left) the decorrelation time, (center) 1-spectral entropy, and (right) the lag-1 autocorrelation.

Fig. 2.

The three leading component time series () found using five predictable component methods and PCA. Each row corresponds to a different method (the method names are on the right).

Fig. 2.

The three leading component time series () found using five predictable component methods and PCA. Each row corresponds to a different method (the method names are on the right).

Fig. 3.

The spatial patterns of the leading three components for five predictable component methods and PCA. Each row corresponds to a different method (the method names are on the right). The filled circles show the correlation between each respective component time series (Fig. 2) and the field . The circle areas are proportional to the size of the correlation, and the circle color shows the sign of the correlation. Red shading corresponds to positive correlation, and blue shading corresponds to negative correlation.

Fig. 3.

The spatial patterns of the leading three components for five predictable component methods and PCA. Each row corresponds to a different method (the method names are on the right). The filled circles show the correlation between each respective component time series (Fig. 2) and the field . The circle areas are proportional to the size of the correlation, and the circle color shows the sign of the correlation. Red shading corresponds to positive correlation, and blue shading corresponds to negative correlation.

Fig. 4.

The power spectra (blue) of the leading three component time series (the component time series were first prewhitened for this diagram). The solid red line is the theoretical mean spectrum for a prewhitened Gaussian red noise process, and the dashed red line is the 95% confidence interval of the mean spectrum.

Fig. 4.

The power spectra (blue) of the leading three component time series (the component time series were first prewhitened for this diagram). The solid red line is the theoretical mean spectrum for a prewhitened Gaussian red noise process, and the dashed red line is the 95% confidence interval of the mean spectrum.

For the five predictable component methods, the first component is the long-term trend in temperature (Fig. 2), which has a positive pattern for all stations (Fig. 3). For predictable component methods, the first component will reveal if a field has a long-term trend, because a long-term trend optimizes both measures of persistence and measures of predictability (Fig. 1). The rank correlation between the first predictable component and time () is highest for PTA () and lowest for SFA () (note, PCA1 had ), showing that, for PTA, the first component more cleanly separates the long-term trend component from other low-frequency variance.

Components 2 and 3 also showed similarities and differences between the five predictable component methods. To compare these components, components from four methods (OPA, ForeCA, SFA, and PTA) were projected into the subspace formed by components 2 and 3 of the APTD method. This projection was done by regressing the former components against the APTD components. For example, in Fig. 5 the projection of OPA2 (the second OPA component) into the subspace formed by APTD2 and APTD3 gives the following equation:

 
formula

Note that all components have unit variance. The regression coefficients represent the projection (hence, the OPA2 vector falls in the positive quadrant in Fig. 5a), and the is a measure of how close a particular component lies to the subspace into which it is projected (vectors that have a low are pointing out of the page in Fig. 5). Also, the correlation coefficient between any two components shown on the figure is given by the cosine of the angle between the components and/or axes. Components with small angles or opposite angles are positively and negatively correlated (e.g., ForeCA3 and SFA7 are negatively correlated), while components that are perpendicular are independent (e.g., ForeCA2, ForeCA3, and ForeCA4).

Fig. 5.

Biplots showing the relationship between predictable components 2–4 of five different methods. The vectors (arrows) were calculated by projecting each component into the subspace formed by two APTD components. (a) The subspace formed by components APTD2 and APTD3. (b) The subspace formed by components APTD3 and APTD4. In (a) and (b), the percentage for each component in parentheses refers to the variance explained by regressing each component on the respective APTD components [e.g., Eq. (18)]. Only components with have been plotted for clarity. In (a) the vectors for SFA2 and PTA2 overlap; hence, the SFA2 label has been moved to the left.

Fig. 5.

Biplots showing the relationship between predictable components 2–4 of five different methods. The vectors (arrows) were calculated by projecting each component into the subspace formed by two APTD components. (a) The subspace formed by components APTD2 and APTD3. (b) The subspace formed by components APTD3 and APTD4. In (a) and (b), the percentage for each component in parentheses refers to the variance explained by regressing each component on the respective APTD components [e.g., Eq. (18)]. Only components with have been plotted for clarity. In (a) the vectors for SFA2 and PTA2 overlap; hence, the SFA2 label has been moved to the left.

Components OPA2, APTD2, ForeCA2, SFA2, and PTA2 are related, because their vectors generally lie along the APTD2 axis (Fig. 5a). For all five predictable component methods, the second component has a higher index of persistence () and predictability () than all the other components (not including the first component) (Fig. 1). The time series of the second component has a quadratic pattern (Fig. 2), which is partly a consequence of the orthogonality imposed by the methods: the component time series will be like members of an orthogonal polynomial set. The associated spatial patterns are typically positive everywhere (Fig. 3) and have large weights in the southeast region, reflecting that the quadratic time trend is stronger in the southeast than in the north. For the second component time series, their power spectra exhibit some power in the 2–4-yr range, especially SFA2 and PTA2 (Fig. 4). SFA2 is the only second component to show large power in the 0.1–0.3-yr range.

In the third set of components (OPA3, APTD3, ForeCA3, SFA3, and PTA3), there is more variation between the components that capture persistence and those that capture predictability. OPA3 is a persistent feature, with power at interannual and interdecadal frequencies (Fig. 4) and relatively large anomalies on the coastal fringe compared to the continental inland (Fig. 3). SFA3 has less power at interdecadal frequencies, and its spatial pattern is dominated by large weights in the north. SFA3 is uncorrelated with the subspace formed from APTD2, APTD3, and APTD4 (Fig. 5). In contrast, the subspace formed from ForeCA3 and ForeCA4 contains basically the same information as the APTD3–APTD4 subspace. APTD and ForeCA are the only methods for which the third component has an (Fig. 1). The spatial patterns for APTD3 and ForeCA3 are dominated by large weights in only a few stations (Perth in the southwest, and Bourke and Charleville in the central east) (Fig. 3). Their spectra show power at harmonics of the annual cycle (~60, 120, and 365 days) (Fig. 4), which suggests that those sites have atypical annual cycles that include higher-order harmonics [the first two harmonics were removed in the deseasonalization step (section 4)]. While ForeCA and APTD both use general measures of predictability, the small differences between the actual components (as in Fig. 5b) are probably because of the differences in the methods. APTD estimates predictability using finite lags in the time domain, while ForeCA uses the frequency domain. Last, ForeCA and APTD seem prone to finding local (or even site-specific) predictable components (this could be remedied by using truncated principal components before ForeCA or APTD), whereas OPA and SFA are better at finding patterns that relate to large-scale factors. Nevertheless, APTD and ForeCA may be useful for identifying atypical sites in a climate observation network.

b. Interpretation of components

The type of predictability exhibited by the observed predictable components was investigated using the null hypotheses outlined in section 2c. For these three tests, SFA was used to calculate the components of the bootstrapped field (see section 2c). The SFA components of the observed field can be straightforwardly related to the components produced by the other four methods. For example, the first and second components produced by all five methods are similar (Figs. 1 and 2). Also, APTD3, ForeCA3, and SFA7 are related (Fig. 5b).

Figure 6 shows the results of the three bootstrap tests about the type of predictability. Note that the predictability of the SFA components of the three null hypotheses increases from left to right; that is, the VARX model components are more predictable than the VAR model components, which are more predictable than the model consisting of interannual variance and white noise. The observed SFA components are clearly much more predictable than white noise, even when there is interannual variance () embedded in the white noise (Fig. 6a). Also, the observed components SFA1 and SFA7 are more predictable than a red noise process (the lag-1 VAR model; Fig. 6b). Further, SFA1 and SFA7 are more predictable than a first-order Markov process that was itself driven by red noise (Fig. 5c). Instead, the predictability in SFA1 and SFA7 may be due to other factors, such as a higher-order Markov process, external forcing, or a combination of these factors. Unlike a first-order Markov process where the interaction is one-way, a higher-order Markov process requires a two-way interaction between the atmosphere and oceans (e.g., as in the Madden–Julian oscillation or the El Niño–Southern Oscillation). Such higher-order Markov processes may also be partly synchronized to the annual solar cycle (Krishnamurti and Chakraborty 2005; Stein et al. 2011). The predictability in SFA1 probably relates to the trend in greenhouse gases, while the predictability in SFA7 may be due to a higher-order Markov process. In contrast, the predictability of SFA2 lies within the confidence limits for the second and third null hypothesis (Figs. 6b,c). Thus, SFA2 may have been produced by a first-order Markov process driven by pink or red noise.

Fig. 6.

Three bootstrap tests of the type of predictability of the SFA components. The three null hypotheses are (left) interannual variance + white noise, (center) a first-order vector autoregressive model, and (right) a first-order VARX model forced by interannual variance and noise (see section 2c). The measure of predictability used here is [Eq. (12)]. The red bars show the 99% confidence intervals for each hypothesis. The SFA components (first 30 components) are plotted as circles; solid circles lie outside the null hypothesis confidence intervals. Note, the y axis is log scale.

Fig. 6.

Three bootstrap tests of the type of predictability of the SFA components. The three null hypotheses are (left) interannual variance + white noise, (center) a first-order vector autoregressive model, and (right) a first-order VARX model forced by interannual variance and noise (see section 2c). The measure of predictability used here is [Eq. (12)]. The red bars show the 99% confidence intervals for each hypothesis. The SFA components (first 30 components) are plotted as circles; solid circles lie outside the null hypothesis confidence intervals. Note, the y axis is log scale.

To further investigate the predictability in the second and third components, the component time series ForeCA2 and ForeCA3 were prewhitened using a first-order autoregressive [AR(1)] model and analyzed using a Morlet wavelet transform. Also, the local spectra of the wavelet transform were used to calculate the predictability () for each time point: that is, if is the wavelet power spectrum at time t and frequency ω, and , then the local predictability can be calculated by substituting into Eq. (12). This statistic is henceforth referred to as the wavelet entropy, because it is a local measure of entropy (and predictability). The wavelet spectrum for the prewhitened ForeCA2 series shows some significant power at interannual scales (2–10 years) and multidecadal scales (Fig. 7b). Note that the prewhitening did not completely remove the quadratic trend in the ForeCA2 time series (Fig. 2) (hence the large variance at low frequencies in Fig. 7b). The wavelet entropy (Fig. 7a) exhibited no significant seasonality (see below). Thus, the wavelet spectrum in (Fig. 7b) is characteristic of a persistent process driven by red and white noise. Such a process could arise when slow components of the climate system (e.g., soil moisture, snow cover, or the subsurface ocean) are forced by interannual variation in sea surface temperature.

Fig. 7.

Wavelet analysis of ForeCA2 and ForeCA3. For ForeCA2, the following are shown: (b) the wavelet analysis of the prewhitened component time series and (c) the global wavelet spectrum. (a) A time series of instantaneous predictability, calculated using the local spectral estimates in (b); that is, the time series in (a) is derived from (b) [not (b) from (a)]. (d)–(f) As in (a)–(c), but for ForeCA3. The inset diagram [to the right of (d)] shows there is a seasonal cycle in the instantaneous predictability () of ForeCA3. Note (c) and (f) use the same y axis as for (b) and (e), respectively. The solid red lines in (c) and (f) show the theoretical mean spectrum (white noise) and 95% confidence interval.

Fig. 7.

Wavelet analysis of ForeCA2 and ForeCA3. For ForeCA2, the following are shown: (b) the wavelet analysis of the prewhitened component time series and (c) the global wavelet spectrum. (a) A time series of instantaneous predictability, calculated using the local spectral estimates in (b); that is, the time series in (a) is derived from (b) [not (b) from (a)]. (d)–(f) As in (a)–(c), but for ForeCA3. The inset diagram [to the right of (d)] shows there is a seasonal cycle in the instantaneous predictability () of ForeCA3. Note (c) and (f) use the same y axis as for (b) and (e), respectively. The solid red lines in (c) and (f) show the theoretical mean spectrum (white noise) and 95% confidence interval.

In contrast, the wavelet spectrum for the prewhitened ForeCA3 series shows significant power at interannual and intraseasonal time scales (Fig. 7e). The global wavelet spectrum of the prewhitened ForeCA3 series has no significant power at the annual frequency (Fig. 7f), but the wavelet entropy (Fig. 7d) shows a seasonal pattern whereby predictability is greatest in August–September. This suggests that the prewhitened ForeCA3 series has a cyclostationary spectrum, characteristic of a process linked to the annual cycle, explaining the significance of SFA7 in Fig. 6.

6. Conclusions

To date, several methods for investigating the predictability of space–time fields have been proposed, which have been reviewed here. These methods find optimal projection vectors that capture predictable components in a field. The advantages of the different methods are as follows: PTA more cleanly separates the long-term trend from other low-frequency variation and white noise. APTD and ForeCA use general measures of predictability such that their measure will be high for blue noise and red noise processes and low for white noise processes. The projection vectors for APTD can be found analytically, while the projection vectors for ForeCA are found by iterative decomposition. However, ForeCA can potentially be used for fields with missing values by calculating the multivariate spectrum using the Lomb–Scargle estimate (Pardo-Igúzquiza and Rodríguez-Tovar 2012). The spectral entropy (used by ForeCA) also naturally leads to the idea of instantaneous predictability (and instantaneous entropy), which can be calculated using a local spectral estimate, such as the wavelet transform. The instantaneous predictability allows the investigation of time-variant changes in predictability. OPA and SFA optimize persistence rather than a general measure of predictability, but these methods are potentially useful because they can be expressed in terms of reduced rank regression. This means that OPA, SFA, and PTA can be used with fields where by regularizing the covariance matrix of the field using a ridge parameter () in the reduced rank regression model that underlies these methods. SFA is also fast, which is useful for simulation-based statistical tests. OPA, PTA, and SFA all seem to be very useful for investigating low-frequency variability in the climate system, while APTD and ForeCA seem particularly suitable for investigating components with predictability in the mid- to high-frequency range.

A new test has been developed for investigating the type of predictability exhibited by the components. Previous tests for predictable components have considered only the stationary white noise null hypothesis (e.g., DelSole 2006; DelSole and Tippett 2009b). The new test is based on an expanded null hypothesis that allows different types of stochasticity to be distinguished. Predictable components of the climate system can arise not only from processes that integrate white noise, but also from processes that integrate other red noise processes. The new test is implemented by bootstrapping a first-order VARX model, where the exogenous variable is estimated as the leading principal component (or components) of the interannual bandpass-filtered field. Time series that can beat this null hypothesis must have been produced from some other type of predictability, such as a higher-order Markov process, or a process driven by external forcing with a regular trend (e.g., greenhouse gases) or oscillation (e.g., solar forcing).

The analysis discussed here can be easily expanded in several ways. The ForeCA method can be adapted not only to fields with missing values, but also to examining the predictability between fields (such as between ocean and land temperature fields) by using the Lomb–Scargle cross spectrum. Second, some methods here can be expanded to investigate complex signal fields (i.e., “complex predictable components” analogous to complex principal components) (Schreier 2008). This would allow the investigation of predictable components that propagate within a field. Third, all the methods can be used to investigate the predictable components in a variance field (i.e., time series of absolute values of a detrended anomaly field). This will be useful for answering questions such as the following: Are there predictable patterns in temperature variance, and do they arise from natural variation? Finally, these methods will be applied to an Australian regional climate reanalysis (in development) to better understand the physical mechanisms behind the predictable components.

Acknowledgments

The author wishes to thank Dr. Adrian Paterson and Prof. John Dodson for comments on a draft manuscript.

REFERENCES

REFERENCES
ABOM
,
2014
: Australian Climate Observations Reference Network - Surface Air Temperature (ACORN-SAT). Australian Bureau of Meteorology (ABOM). Accessed 21 Feb 2014. [Available online at http://www.bom.gov.au/climate/change/acorn-sat.]
DelSole
,
T.
,
2006
:
Low-frequency variations of surface temperature in observations and simulations
.
J. Climate
,
19
,
4487
4507
, doi:.
DelSole
,
T.
, and
M. K.
Tippett
,
2009a
:
Average predictability time. Part I: Theory
.
J. Atmos. Sci.
,
66
,
1172
1187
, doi:.
DelSole
,
T.
, and
M. K.
Tippett
,
2009b
:
Average predictability time. Part II: Seamless diagnoses of predictability on multiple time scales
.
J. Atmos. Sci.
,
66
,
1188
1204
, doi:.
Di Lorenzo
,
E.
, and
M. D.
Ohman
,
2013
:
A double-integration hypothesis to explain ocean ecosystem response to climate forcing
.
Proc. Natl. Acad. Sci. USA
,
110
,
2496
2499
, doi:.
Fischer
,
M. J.
, and
A. W.
Paterson
,
2014
:
Detecting trends that are nonlinear and asymmetric on diurnal and seasonal time scales
.
Climate Dyn.
,
43
,
361
374
, doi:.
Ghil
,
M.
, and Coauthors
,
2002
:
Advanced spectral methods for climatic time series
.
Rev. Geophys.
,
40
,
3-1
3-41
, doi:.
Goerg
,
G. M.
,
2013
:
Forecastable components analysis
.
J. Mach. Learn. Res. Workshop Conf. Proc.
,
28
,
64
72
. [Available online at http://jmlr.org/proceedings/papers/v28.]
Hannachi
,
A.
,
2007
:
Pattern hunting in climate: A new method for finding trends in gridded climate data
.
Int. J. Climatol.
,
27
,
1
15
, doi:.
Huang
,
N. E.
, and
Z.
Wu
,
2008
:
A review on Hilbert–Huang Transform: Method and its applications to geophysical studies
.
Rev. Geophys.
,
46
,
RG2006
, doi:.
Jones
,
D. A.
,
1999
:
Characteristics of Australian land surface temperature variability
.
Theor. Appl. Climatol.
,
63
,
11
31
, doi:.
Krishnamurti
,
T. N.
, and
D. R.
Chakraborty
,
2005
:
The dynamics of phase locking
.
J. Atmos. Sci.
,
62
,
2952
2964
, doi:.
Monahan
,
A. H.
,
J. C.
Fyfe
,
M. H. P.
Ambaum
,
D. B.
Stephenson
, and
G. R.
North
,
2009
:
Empirical orthogonal functions: The medium is the message
.
J. Climate
,
22
,
6501
6514
, doi:.
Newman
,
M.
,
G. P.
Compo
, and
M. A.
Alexander
,
2003
:
ENSO-forced variability of the Pacific decadal oscillation
.
J. Climate
,
16
,
3853
3857
, doi:.
Pardo-Igúzquiza
,
E.
, and
F. J.
Rodríguez-Tovar
,
2012
:
Spectral and cross-spectral analysis of uneven time series with the smoothed Lomb–Scargle periodogram and Monte Carlo evaluation of statistical significance
.
Comput. Geosci.
,
49
,
207
216
, doi:.
Schneider
,
T.
, and
S. M.
Griffies
,
1999
:
A conceptual framework for predictability studies
.
J. Climate
,
12
,
3133
3155
, doi:.
Schreier
,
P. J.
,
2008
:
A unifying discussion of correlation analysis for complex random vectors
.
IEEE Trans. Signal Process.
,
56
,
1327
1336
, doi:.
Shakun
,
J. D.
, and
J.
Shaman
,
2009
: Tropical origins of North and South Pacific decadal variability. Geophys. Res. Lett.,36, L19711, doi:.
Stein
,
K.
,
A.
Timmermann
, and
N.
Schneider
,
2011
: Phase synchronization of the El Niño–Southern Oscillation with the annual cycle. Phys. Rev. Lett.,107, 128501, doi:.
Taylor
,
M.
,
M.
Losch
,
M.
Wenzel
, and
J.
Schröter
,
2013
:
On the sensitivity of field reconstruction and prediction using empirical orthogonal functions derived from gappy data
.
J. Climate
,
26
,
9194
9205
, doi:.
Trewin
,
B.
,
2013
:
A daily homogenized temperature data set for Australia
.
Int. J. Climatol.
,
33
,
1510
1529
, doi:.
Tsay
,
R. S.
,
2013
:
Multivariate Time Series Analysis: With R and Financial Applications.
Wiley, 520 pp.
Wiskott
,
L.
, and
T. J.
Sejnowski
,
2002
:
Slow feature analysis: Unsupervised learning of invariances
.
Neural Comput.
,
14
,
715
770
, doi:.

Footnotes

*

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-14-00713.s1.

Supplemental Material