Comparing climate time series – Part 1: Univariate test

This paper proposes a new approach to detecting and describing differences in stationary processes. The approach is equivalent to comparing auto-covariance functions or power spectra. The basic idea is to fit an autoregressive model to each time series and then test whether the model parameters are equal. The likelihood ratio test for this hypothesis has appeared in the statistics literature, but the resulting test depends on maximum likelihood estimates, which are biased, neglect differences in noise parameters, and utilize sampling distributions that are valid only for large sample sizes. This paper derives a likelihood ratio test that corrects for bias, detects differences in noise parameters, and can be applied to small samples. Furthermore, if a significant difference is detected, we propose new methods to diagnose and visualize those differences. Specifically, the test statistic can be used to define a “distance” between two autoregressive processes, which in turn can be used for clustering analysis in multi-model comparisons. A multidimensional scaling technique is used to visualize the similarities and differences between time series. We also propose diagnosing differences in stationary processes by identifying initial conditions that optimally separate predictable responses. The procedure is illustrated by comparing simulations of an Atlantic Meridional Overturning Circulation (AMOC) index from 10 climate models in Phase 5 of the Coupled Model Intercomparison Project (CMIP5). Significant differences between most AMOC time series are detected. The main exceptions are time series from CMIP models from the same institution. Differences in stationary processes are explained primarily by differences in the mean square error of 1-year predictions and by differences in the predictability (i.e., R-square) of the associated autoregressive models.

1 Introduction

Climate scientists often confront questions of the following types.

1. Does a climate model realistically simulate a climate index?

2. Do two climate models generate similar temporal variability?

3. Did a climate index change its variability?

4. Are two power spectra consistent with each other?

Each of the above questions requires deciding whether two time series come from the same stochastic process. Although numerous papers in the weather and climate literature address questions of the above types, the conclusions often are based on visual comparison of estimated auto-covariance functions or power spectra without a rigorous significance test. Lund et al. (2009) provide a lucid review of some objective tests for deciding whether two time series come from the same stationary process. An additional test that was not considered by Lund et al. (2009) is to fit autoregressive models to time series and then to test differences in parameters (Maharaj, 2000; Grant and Quinn, 2017). Grant and Quinn (2017) showed that this test has good power and performs well even when the underlying time series are not from an autoregressive process. The latter test has been applied to such problems as speech recognition but not to climate time series. The purpose of this paper is to further develop this test for climate applications.

The particular tests proposed by Maharaj (2000) and Grant and Quinn (2017) test equality of autocorrelation without regard to differences in variances. However, in climate applications, differences in variance often are of considerable importance. In addition, these tests employ a sampling distribution derived from asymptotic theory and therefore may be problematic for small sample sizes. In this paper, the likelihood ratio test is derived for the more restrictive case of equality of noise variances, which leads to considerable simplifications, including a test that is applicable for small sample sizes. Further comments about how our proposed test compares with previous tests will be discussed below.

If the time series are deemed to come from different processes, then it is desirable to characterize those differences in meaningful terms and to group time series into clusters. Piccolo (1990) and Maharaj (2000) propose such classification procedures. Following along these lines, our hypothesis test suggests a natural measure for measuring the distance between two stationary processes that can be used for clustering analysis. For multi-model studies, this distance measure can be used to give a graphical summary of the similarities and differences between time series generated by different models. In addition, we extend the interpretation of such summaries considerably by proposing a new approach to diagnosing differences in stationary processes based on finding the initial condition that optimally separates one-step predictions.

2 Previous methods for comparing time series

In this section we review previous methods for comparing two time series. These methods are based on the theory of stochastic processes and assume that the joint distribution of the values of the process at any collection of times is multivariate normal and stationary. Although non-stationarity is a prominent feature in many climate time series, a stationary framework is a natural starting point for non-stationary generalizations. Stationarity implies that the expectation at any time is constant, and the second-order moments depend only on the difference times (more precisely, this is called weak-sense stationarity). Accordingly, if X_t is a stationary process, then the mean is independent of time,

$$\begin{array}{}\text{(1)}& \mathbb{E}\left[X_t\right]={\mathit{\mu }}_X\text{(a constant)},\end{array}$$

and the time-lagged auto-covariance function depends only on the difference in times,

$$\begin{array}{}\text{(2)}& \mathbb{E}\left[\right(X_{t+\mathit{\tau }}-{\mathit{\mu }}_X\left)\right(X_t-{\mathit{\mu }}_X\left)\right]=c_X\left(\mathit{\tau }\right).\end{array}$$

Because a multivariate normal distribution is fully characterized by its first and second moments, the stochastic process is completely specified by μ_X and the auto-covariance function c_X(τ). Stationarity further implies that the auto-covariance function is an even function of the lag τ. Following standard practice, we consider discrete time series where values are available at N_X equally spaced time steps $X_{\mathrm{1}},X_{\mathrm{2}},\mathrm{\dots },X_{N_X}$.

Now consider another stationary process Y_t, with mean μ_Y and auto-covariance function c_Y(τ). The problem we consider is this: given sample time series $X_{\mathrm{1}},\mathrm{\dots },X_{N_X}$ and $Y_{\mathrm{1}},\mathrm{\dots },Y_{N_Y}$, decide whether the two time series come from the same stationary process. For stationary processes, we often are not concerned with differences in means (e.g., in climate studies, these often are eliminated through “bias corrections”), hence we allow μ_X≠μ_Y. Therefore, our problem is equivalent to deciding whether the two time series have the same auto-covariance function, i.e., deciding whether

$$\begin{array}{}\text{(3)}& c_X\left(\mathit{\tau }\right)=c_Y\left(\mathit{\tau }\right)\text{for all }\mathit{\tau }=\mathrm{0},\mathrm{1},\mathrm{\dots }.\end{array}$$

This problem can be framed equivalently as deciding whether two stationary processes have the same power spectrum. Recall that the power spectrum is the Fourier transform of the auto-covariance function. We define the power spectrum of X_t as

$$\begin{array}{}\text{(4)}& p_X\left(\mathit{\omega }\right)=\sum _{\mathit{\tau }=-\mathrm{\infty }}^{\mathrm{\infty }}{c}_X\left(\mathit{\tau }\right)e^{i\mathit{\omega }\mathit{\tau }}.\end{array}$$

The spectrum of Y_t is defined similarly and denoted as p_Y(ω). Because the Fourier transform is invertible, equality of auto-covariances is equivalent to equality of spectra. Thus, our problem can be framed equivalently as deciding whether

$$\begin{array}{}\text{(5)}& p_X\left(\mathit{\omega }\right)=p_Y\left(\mathit{\omega }\right)\text{for all }\mathit{\omega }\in [\mathrm{0},\mathit{\pi }).\end{array}$$

Estimates of the power spectrum are based on the periodogram (Box et al., 2008).

The above hypothesis differs from hypotheses about a single process that are commonly tested with auto-covariance functions or power spectra. For instance, the hypothesis of vanishing auto-correlation often is assessed by comparing the sample auto-correlation function to $\pm \mathrm{1.96}/\sqrt{N}$, where N is the length of the time series (e.g., Brockwell and Davis, 2002, chap. 1). In the spectral domain, the hypothesis of white noise often is tested based on the Kolmogorov–Smirnov test, in which the standardized cumulative periodogram is compared to an appropriate set of lines (e.g., Jenkins and Watts, 1968, Sect. 6.3.2). These tests consider hypotheses about one process. In contrast, our hypothesis involves a comparison of two processes.

Coates and Diggle (1986) derived spectral-domain tests for equality of stationary processes. The underlying idea of these tests is that equality of power spectra implies that their ratio is independent of frequency and therefore indistinguishable from white noise. This fact suggests that standard tests for white noise can be adapted to periodogram ratios. Coates and Diggle (1986) derive a second parametric test that assumes that the log ratio of power spectra is a low-order polynomial of frequency.

Lund et al. (2009) explored the above tests in the context of station data for temperature. They found that spectral methods have relatively low statistical power – that is, the methods are unlikely to detect a difference in stationary processes when such a difference exists. Our own analysis using the data discussed in Sect. 5 is consistent with Lund et al. (2009) (not shown). The fact that spectral-domain tests have less statistical power than time-domain tests is not surprising. After all, spectral tests are based on comparing periodograms in which the number of unknown parameters grows with sample size. In particular, for a time series of length N, the periodogram combines coefficients for sines and cosines of a Fourier transform into N∕2 coefficients for the amplitude. Thus, a time series of length N=64 yields a periodogram with 32 amplitudes; a time series of length N=512 yields a periodogram with 256 amplitudes; and so on. Because the number of unknowns grows with sample size, the sampling error of the individual periodogram estimates does not decrease with sample size. Typically, sampling errors are reduced by smoothing periodogram estimates over frequencies, but the smoothing makes implicit assumptions about the shape of the underlying power spectrum. In the absence of hypotheses to constrain the power spectra, the large number of estimated parameters results in low statistical power.

Lund et al. (2009) also develop and discuss a time-domain test. This test is based on the sample estimate of the auto-covariance function

$$\begin{array}{}\text{(6)}& {\widehat{c}}_X\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{N_X}}\sum _{t=\mathrm{1}}^{N_X-\left|\mathit{\tau }\right|}\left(X_{t+\mathit{\tau }}-{\widehat{\mathit{\mu }}}_X\right)\left(X_t-{\widehat{\mathit{\mu }}}_X\right),\end{array}$$

where ${\widehat{\mathit{\mu }}}_X$ is the sample mean of X_t based on sample size N_X. The analogous estimate for the auto-covariance function of Y_t is denoted ${\widehat{c}}_Y\left(\mathit{\tau }\right)$. The test is based on differences in auto-covariance functions up to lag τ₀. That is, the test is based on

$$\begin{array}{}\text{(7)}& \mathbf{\Delta }=\left(\begin{array}{c}{\widehat{c}}_X\left(\mathrm{0}\right)-{\widehat{c}}_Y\left(\mathrm{0}\right)\\ {\widehat{c}}_X\left(\mathrm{1}\right)-{\widehat{c}}_Y\left(\mathrm{1}\right)\\ \mathrm{⋮}\\ {\widehat{c}}_X\left({\mathit{\tau }}_{\mathrm{0}}\right)-{\widehat{c}}_Y\left({\mathit{\tau }}_{\mathrm{0}}\right)\end{array}\right).\end{array}$$

In the case of equal sample sizes $N_X=N_Y=N$, Lund et al. (2009) propose the statistic

$$\begin{array}{}\text{(8)}& {\mathit{\chi }}^{\mathrm{2}}={\displaystyle \frac{N}{\mathrm{2}}}{\mathbf{\Delta }}^T{\widehat{\mathbf{W}}}^{-\mathrm{1}}\mathbf{\Delta },\end{array}$$

where $\widehat{\mathbf{W}}$ is an estimate of the covariance matrix between auto-covariance estimates, namely

$$\begin{array}{}\text{(9)}& {\widehat{W}}_{i+\mathrm{1},j+\mathrm{1}}=\sum _{k=-K}^K\left(\widehat{c}\left(k\right)\widehat{c}(k-i+j)+\widehat{c}(k+j)\widehat{c}(k-i)\right),\end{array}$$

and $\widehat{c}\left(\mathit{\tau }\right)$ is the pooled auto-covariance estimate

$$\begin{array}{}\text{(10)}& \widehat{c}\left(\mathit{\tau }\right)={\displaystyle \frac{{\widehat{c}}_X\left(\mathit{\tau }\right)+{\widehat{c}}_Y\left(\mathit{\tau }\right)}{\mathrm{2}}}.\end{array}$$

The reasoning behind this statistic is lucidly discussed in Lund et al. (2009) and follows from standard results in time series analysis (see proposition 7.3.2 in Brockwell and Davis, 1991). Under the null hypothesis of equal auto-covariance functions, and for large N, the statistic χ² has an approximate chi-square distribution with τ₀+1 degrees of freedom. A key parameter in this statistic is the cutoff parameter K in the sum for $\widehat{\mathbf{W}}$. Lund et al. (2009) propose using $K=N^{\mathrm{1}/\mathrm{3}}$ but acknowledge that this rule needs further study.

We have applied Lund et al.'s test to the numerical examples discussed in Sect. 5 and find that $\widehat{\mathbf{W}}$ is not always positive definite. In such cases, χ² is not guaranteed to be positive and therefore does not have a chi-square distribution. The lack of positive definiteness sometimes can be avoided by choosing a slightly different K, but in many of these cases the resulting χ² depends sensitively on K. This sensitivity arises from the fact that $\widehat{\mathbf{W}}$ is close to singular, so changing K by one unit can change χ² by more than an order of magnitude. Conceivably, some modification of $\widehat{\mathbf{W}}$ or some rule for choosing K can remove this sensitivity to K, but without such modification this test is considered unreliable. Further comments about this are given at the end of Sect. 5.

3 Comparing autoregressive models

Several authors have proposed approaches to comparing stationary processes based on assuming that the time series come from autoregressive models of order p, denoted AR(p). Accordingly, we consider the two AR(p) models

$$\begin{array}{}\text{(11)}& {\displaystyle }{X}_t& {\displaystyle }={\mathit{\varphi }}_{\mathrm{1}}^X{X}_{t-\mathrm{1}}+{\mathit{\varphi }}_{\mathrm{2}}^X{X}_{t-\mathrm{2}}+\mathrm{\dots }+{\mathit{\varphi }}_p^X{X}_{t-p}+{\mathit{\gamma }}_X+{\mathit{ϵ}}_t^X,\text{(12)}& {\displaystyle }{Y}_t& {\displaystyle }={\mathit{\varphi }}_{\mathrm{1}}^Y{Y}_{t-\mathrm{1}}+{\mathit{\varphi }}_{\mathrm{2}}^Y{Y}_{t-\mathrm{2}}+\mathrm{\dots }+{\mathit{\varphi }}_p^Y{Y}_{t-p}+{\mathit{\gamma }}_Y+{\mathit{ϵ}}_t^Y,\end{array}$$

where the ϕs are autoregressive parameters, the γs are constants that control the mean, and the ϵ_ts are independent Gaussian white noise processes with zero mean and variances

$$\begin{array}{}\text{(13)}& {\mathit{\sigma }}_X^{\mathrm{2}}=\mathrm{var}\left[{\mathit{ϵ}}_t^X\right]\text{and}{\mathit{\sigma }}_Y^{\mathrm{2}}=\mathrm{var}\left[{\mathit{ϵ}}_t^Y\right].\end{array}$$

By construction, X_t and Y_s are independent for all t and s. Our method can be generalized to handle correlations between X_t and Y_s, specifically by using a vector autoregressive model with coupling between the two processes, but this generalization will not be considered here. For the above models, there exists a one-to-one relation between the first p+1 auto-covariances $c\left(\mathrm{0}\right),\mathrm{\dots },c\left(p\right)$ and the parameters ${\mathit{\varphi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\varphi }}_p,{\mathit{\sigma }}^{\mathrm{2}}$. This relation is expressed through the Yule–Walker equations and a balance equation for variance (Box et al., 2008). Therefore, equality of the first p+1 auto-covariances is equivalent to equality of the AR(p) parameters. Furthermore, estimates of the parameters depend only on the first p+1 auto-covariances. The remaining auto-covariances of an AR(p) process $c(p+\mathrm{1}),c(p+\mathrm{2}),\mathrm{\dots }$ are obtained through recursion formulas that depend only on the first p+1 auto-covariances. Importantly, the auto-covariance function of an AR(p) process does not depend on γ; rather, γ controls the mean of the process. In climate studies, the mean often is removed from time series before any analysis, hence we do not require the means to be equal; i.e., we allow γ_X≠γ_Y.

In previous tests, different authors make different assumptions about the noise variance. Maharaj (2000) allows the noise in the two models to differ in their variances and to be correlated at zero lag. Grant and Quinn (2017) allow the noise variances to differ but assume the noises are independent at zero lag. In the climate applications we have in mind, differences in variance are important, hence we are interested in detecting differences in noise variances. Accordingly, our null hypothesis of equivalent stationary processes is the following:

$$\begin{array}{}\text{(14)}& H_{\mathrm{0}}:\left\{{\mathit{\varphi }}_{\mathrm{1}}^X={\mathit{\varphi }}_{\mathrm{1}}^Y,\mathrm{\dots },{\mathit{\varphi }}_p^X={\mathit{\varphi }}_p^Y\right\}\text{and }{\mathit{\sigma }}_X={\mathit{\sigma }}_Y.\end{array}$$

The alternative hypothesis is that at least one of the above parameters differs between the two processes. Constraining processes to be AR(p) means that each process is characterized by p+1 parameters, in contrast to the unconstrained case in which the auto-covariance function or power spectrum is characterized by an infinite number of parameters.

Table 1 Definition of the variables in the test statistic D_ϕ,σ (Eq. 17).

The likelihood ratio test for hypothesis H₀ is derived in the Appendix. The derivation differs from that of Grant and Quinn (2017) primarily by including the hypothesis σ_X=σ_Y in the null hypothesis, which leads to considerable simplifications in estimation and interpretation. Furthermore, we propose a modification to correct for biases and a Monte Carlo technique for determining significance thresholds. We describe the overall procedure here and direct the reader to the Appendix for details. The test relies on sample estimates of the above parameters, so we define these first. The autoregressive parameters in Eq. (11) are estimated using the method of least squares, yielding the least squares estimates ${\widehat{\mathit{\varphi }}}_{\mathrm{1}}^X,\mathrm{\dots },{\widehat{\mathit{\varphi }}}_p^X,{\widehat{\mathit{\gamma }}}^X$. (The least squares method arises when maximizing the conditional likelihood, which is conditional on the specific realization of $X_{\mathrm{1}},\mathrm{\dots },X_p$ in the data. For large sample sizes, the conditional least squares estimates will be close to the familiar maximum likelihood estimates.) The noise variance for X_t is estimated using the unbiased estimator

$$\begin{array}{}\text{(15)}& {\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}^X{X}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p^X{X}_{t-p}-{\widehat{\mathit{\gamma }}}^X\right)}^{\mathrm{2}}}{{\mathit{\nu }}_X}},\end{array}$$

where ${\mathit{\nu }}_X=N_X-\mathrm{2}p-\mathrm{1}$ is degrees of freedom, computed as the difference between the sample size N_X−p and the number of estimated parameters p+1. This noise variance estimate is merely the residual variance of the AR model and is readily available from standard time series analysis software. Similarly, the method of least squares is used to estimate the AR parameters for Y_t in Eq. (12), yielding the least squares estimates ${\widehat{\mathit{\varphi }}}_{\mathrm{1}}^Y,\mathrm{\dots },{\widehat{\mathit{\varphi }}}_p^Y,{\widehat{\mathit{\gamma }}}^Y$ and the noise variance estimate

$$\begin{array}{}\text{(16)}& {\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_Y}{\left(Y_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}^Y{Y}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p^Y{Y}_{t-p}-{\widehat{\mathit{\gamma }}}^Y\right)}^{\mathrm{2}}}{{\mathit{\nu }}_Y}},\end{array}$$

where ${\mathit{\nu }}_Y=N_Y-\mathrm{2}p-\mathrm{1}$. Then, the test of H₀ is based on the statistic

$$\begin{array}{}\text{(17)}& D_{\mathit{\varphi },\mathit{\sigma }}=D_{\mathit{\sigma }}+D_{\mathit{\varphi }\mid \mathit{\sigma }},\end{array}$$

where

$$\begin{array}{}\text{(18)}& {\displaystyle }\begin{array}{rl}{D}_{\mathit{\sigma }}& ={\mathit{\nu }}_X\log \left({\mathit{\nu }}_X+{\mathit{\nu }}_Y/F_{\mathit{\sigma }}\right)+{\mathit{\nu }}_Y\log \left({\mathit{\nu }}_X{F}_{\mathit{\sigma }}+{\mathit{\nu }}_Y\right)\\ & -({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\log ({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\end{array}\text{(19)}& {\displaystyle }{D}_{\mathit{\varphi }\mid \mathit{\sigma }}=({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\log \left(\mathrm{1}+{\displaystyle \frac{p{F}_{\mathit{\varphi }\mid \mathit{\sigma }}}{{\mathit{\nu }}_X+{\mathit{\nu }}_Y}}\right),\end{array}$$

and the remaining variables are defined in Table 1. Under the null hypothesis H₀, it is shown in the Appendix that F_σ and F_ϕ∣σ have the following approximate distributions:

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }{F}_{\mathit{\sigma }}\sim F_{{\mathit{\nu }}_X,{\mathit{\nu }}_Y},\text{(21)}& {\displaystyle }& {\displaystyle }{F}_{\mathit{\varphi }\mid \mathit{\sigma }}\sim F_{p,{\mathit{\nu }}_X+{\mathit{\nu }}_Y}.\end{array}$$

Furthermore, the two statistics are independent. Critical values for D_σ and D_ϕ∣σ are obtained individually from the critical values of the F-distribution, taking care to use the correct one-tailed or two-tailed test (see Appendix, particularly Eq. A31). In principle, the exact sampling distribution of D_ϕ,σ can be derived analytically because D_ϕ,σ is the sum of two random variables with known distributions. However, this analytic calculation is cumbersome, whereas the quantiles of $D_{\mathit{\sigma }}+D_{\mathit{\varphi }\mid \mathit{\sigma }}$ can be estimated accurately and quickly by Monte Carlo techniques. Essentially, one draws random samples from $F_{{\mathit{\nu }}_X,{\mathit{\nu }}_Y}$ and $F_{p,{\mathit{\nu }}_X+{\mathit{\nu }}_Y}$, substitutes these into Eqs. (18)–(19), evaluates D_ϕ,σ in Eq. (17), and then repeats this many times (e.g., 10 000 times). Note that the required repetitions do not grow with sample size N_X and N_Y since the F-distributions can be sampled directly. The 5 % significance threshold for D_ϕ,σ, denoted D_crit, is then obtained from the 95th percentile of the Monte Carlo samples. Although a Monte Carlo technique is used to obtain critical values, the test still constitutes a test for small sample sizes.

The above test assumes that time series come from an AR(p) process of known order and are excited by Gaussian noise. Grant and Quinn (2017) use Monte Carlo simulations to show that the test is robust to non-Gaussianity. Since our proposed method assumes normal distributions through Eqs. (20) and (21), its robustness to departures to Gaussian remains an open question and is a topic for future study. If the order of the process is unknown and has to be selected or the time series does not come from an autoregressive model (e.g., the process is a moving average process), then the type I error rate does not match its nominal value. In such cases, Grant and Quinn (2017) propose using some sufficiently large value of p, such as

$$\begin{array}{}\text{(22)}& p^{*}=\lfloor {\left(\log min[N_X,N_Y]\right)}^{\mathit{\nu }}\rfloor ,\end{array}$$

where ⌊k⌋ denotes the largest integer smaller than or equal to k. The resulting AR(p^*) model can capture the first $p^{*}+\mathrm{1}$ auto-covariances regardless of their true origin. Of course, there will be some loss of power when a test based on p^* is applied to a time series that comes from an AR(p) model with $p<p^{*}$. This criterion seems to work well in our examples, as discussed in more detail in Sect. 5.

4 Interpretation of differences in AR processes

If a difference in the AR process is detected, then we would like to interpret those differences. After all, such comparisons often are motivated to validate climate models, so if a discrepancy is found we would like to describe those differences in meaningful terms. The fact that the statistic D_ϕ,σ can be expressed as the sum of two independent terms suggests that it is natural to consider the two terms separately. The first term D_σ depends only on the difference in noise variance estimates ${\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}$ and ${\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}$. The noise variance ${\mathit{\sigma }}_X^{\mathrm{2}}$ is the one-step prediction error of the AR model. After all, if the AR parameters in Eq. (11) were known exactly, then the conditional mean would be

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{S}_X=& \mathbb{E}\left[X_t\mid X_{t-\mathrm{1}},X_{t-\mathrm{2}},\mathrm{\dots },X_{t-p}\right]={\mathit{\varphi }}_{\mathrm{1}}^X{X}_{t-\mathrm{1}}\\ & +{\mathit{\varphi }}_{\mathrm{2}}^X{X}_{t-\mathrm{2}}+\mathrm{\dots }-{\mathit{\varphi }}_p^X{X}_{t-p}+{\mathit{\gamma }}_X,\end{array}\end{array}$$

where “S” stands for “signal”, and the one-step prediction error would be

$$\begin{array}{}\text{(24)}& X_t-\mathbb{E}\left[X_t\mid X_{t-\mathrm{1}},\mathrm{\dots },X_{t-p}\right]={\mathit{ϵ}}_t^X.\end{array}$$

Comparing the variance of one-step prediction errors is more statistically straightforward than comparing variances of the original time series because prediction errors are approximately uncorrelated, since they are the residuals of the AR model, which acts as a pre-whitening transformation. In contrast, comparing variances of time series is not straightforward because of serial correlation. In practice, the residuals are only approximately uncorrelated because the parameter values are only approximate.

The second term D_ϕ∣σ vanishes when ${\widehat{\mathit{\varphi }}}_X={\widehat{\mathit{\varphi }}}_Y$ and is positive otherwise, hence it clearly measures the difference in AR parameters. To further interpret this term, it is helpful to partition the AR model for X_t into two parts, an unpredictable part associated with the noise ${\mathit{ϵ}}_t^X$ and a predictable part S_X. The predictable part is the response to the “initial condition”

$$\begin{array}{}\text{(25)}& \mathit{u}=\left(\begin{array}{c}{X}_{t-\mathrm{1}}\\ X_{t-\mathrm{2}}\\ \mathrm{\dots }\\ X_{t-p}\end{array}\right).\end{array}$$

With this notation, ${\widehat{S}}_X={\mathit{u}}^T{\widehat{\mathit{\varphi }}}_X$ is the estimated predictable response of X_t. The estimated predictable response of Y_t can be defined analogously and denoted ${\widehat{S}}_Y$. The initial conditions u may be drawn from the stationary distribution of X_t, the stationary distribution Y_t, or some mixture of the two. In fact, the covariance matrix Σ_HM defined in Eq. (A40) is proportional to the “harmonic mean” of the sample time-lagged covariance matrices for X_t and Y_t. Assuming that the initial condition u is drawn from a distribution with covariance matrix Σ_HM independently of u, then

$$\begin{array}{}\text{(26)}& \mathrm{var}\left[{\widehat{S}}_X\right]=\mathrm{var}\left[{\mathit{u}}^T{\widehat{\mathit{\varphi }}}_X\right]={\mathit{\varphi }}_X^T{\mathbf{\Sigma }}_{\mathrm{HM}}{\mathit{\varphi }}_X.\end{array}$$

If the initial condition for X_t and Y_t is contrived to be the same u, then

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\mathrm{var}\left[{\widehat{S}}_X-{\widehat{S}}_Y\right]& =\mathrm{var}\left[{\mathit{u}}^T\left({\mathit{\varphi }}_X-{\mathit{\varphi }}_Y\right)\right]\\ & ={\left({\mathit{\varphi }}_X-{\mathit{\varphi }}_Y\right)}^T{\mathbf{\Sigma }}_{\mathrm{HM}}\left({\mathit{\varphi }}_X-{\mathit{\varphi }}_Y\right).\end{array}\end{array}$$

Comparing this expression to F_ϕ∣σ in Table 1 suggests that F_ϕ∣σ is a kind of signal-to-noise ratio, where the “signal” is a difference in predictable responses for the same initial condition.

This difference in predictable responses, and hence F_ϕ∣σ, can be related to the difference in R-squares of the individual time series. To see this, expand Eq. (27) as

$$\begin{array}{}\text{(28)}& \mathrm{var}\left[S_X-S_Y\right]=\mathrm{var}\left[S_X\right]+\mathrm{var}\left[S_Y\right]-\mathrm{2}{\mathit{\varphi }}_X^T{\mathbf{\Sigma }}_{\mathrm{HM}}{\mathit{\varphi }}_Y.\end{array}$$

The Cauchy–Schwarz inequality implies that

$$\begin{array}{}\text{(29)}& {\mathit{\varphi }}_X^T{\mathbf{\Sigma }}_{\mathrm{HM}}{\mathit{\varphi }}_Y\le \sqrt{\mathrm{var}\left[S_X\right]\mathrm{var}\left[S_Y\right]}.\end{array}$$

Hence, the above two expressions imply that

$$\begin{array}{}\text{(30)}& \mathrm{var}\left[S_X-S_Y\right]\ge {\left(\sqrt{\mathrm{var}\left[S_X\right]}-\sqrt{\mathrm{var}\left[S_Y\right]}\right)}^{\mathrm{2}}.\end{array}$$

In the special case of equal noise variances, the above inequality becomes

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{var}\left[S_X-S_Y\right]}{{\mathit{\sigma }}^{\mathrm{2}}}}\ge {\left(\sqrt{{\text{SNR}}_X}-\sqrt{{\text{SNR}}_Y}\right)}^{\mathrm{2}},\end{array}$$

where SNR_X and SNR_Y are the signal-to-noise ratios of the two AR models:

$$\begin{array}{}\text{(32)}& {\text{SNR}}_X={\displaystyle \frac{\mathrm{var}\left[S_X\right]}{{\mathit{\sigma }}_X^{\mathrm{2}}}}\text{and}{\text{SNR}}_X={\displaystyle \frac{\mathrm{var}\left[S_Y\right]}{{\mathit{\sigma }}_Y^{\mathrm{2}}}}.\end{array}$$

Note that pF_ϕ∣σ is an estimate of the left-hand side of Eq. (31). Recall that R-square is defined as one minus the ratio of the noise variance to the total variance:

$$\begin{array}{}\text{(33)}& R_X^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{\mathrm{var}\left[{\mathit{ϵ}}_t^X\right]}{\mathrm{var}\left[X_t\right]}}={\displaystyle \frac{{\text{SNR}}_X}{{\text{SNR}}_X+\mathrm{1}}}.\end{array}$$

Because R-square and signal-to-noise ratio are one-to-one, inequality Eq. (31) implies the following: if the noise variances are identical, then a large difference in predictabilities (i.e., R-squares) necessarily implies a tendency for F_ϕ∣σ to be large. This suggests that it may be informative to compare R-squares. We use the sample estimate

$$\begin{array}{}\text{(34)}& \begin{array}{rl}& {\widehat{R}}_X^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\mathit{\nu }}_X{\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}}{{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-\stackrel{\mathrm{‾}}{X}\right)}^{\mathrm{2}}}}\\ & \text{where}\stackrel{\mathrm{‾}}{X}={\displaystyle \frac{\mathrm{1}}{N_X-p}}\sum _{t=p+\mathrm{1}}^{N_X}{X}_t,\end{array}\end{array}$$

which is always between 0 and 1.

It should be recognized that a difference in R-square is a sufficient, but not necessary, condition for a difference in predictable responses. This fact can be seen from the fact that R-square for an AR(p) model is

$$\begin{array}{}\text{(35)}& R_X^{\mathrm{2}}={\mathit{\varphi }}_{\mathrm{1}}^X{\mathit{\rho }}_X\left(\mathrm{1}\right)+\mathrm{\dots }+{\mathit{\varphi }}_p^X{\mathit{\rho }}_X\left(p\right).\end{array}$$

In particular, two processes may have the same R-square because the combination of ϕs yields the same R-square, but the specific values of the ϕs may differ. Although a difference in R-squares is not a perfect indicator of F_ϕ∣σ, it is at least a sufficient condition and therefore worth examination.

On the other hand, if the R-squares are the same, a difference in process still could be detected and should be explained. Simply identifying differences in ϕs would be unsatisfying because those differences have a complicated relation to the statistical characteristics of the process when p>1. Accordingly, we propose the following approach, which despite its limitations may still be insightful. Our basic idea is to choose initial conditions that maximize the mean square difference in predictable responses. To be most useful, we want this choice of initial condition to account for a multi-model comparison over M models. Accordingly, we define the mean square difference between predictable responses as

$$\begin{array}{}\text{(36)}& \begin{array}{rl}\mathrm{\Gamma }\left[\mathit{u}\right]& =\sum _{m=\mathrm{1}}^M\sum _{m^{\prime }=\mathrm{1}}^M{∥S_m\left(\mathit{u}\right)-S_{m^{\prime }}\left(\mathit{u}\right)∥}^{\mathrm{2}}\\ & =\sum _{m=\mathrm{1}}^M\sum _{m^{\prime }=\mathrm{1}}^M{∥{\mathit{u}}^T\left({\mathit{\varphi }}_m-{\mathit{\varphi }}_{m^{\prime }}\right)∥}^{\mathrm{2}}={\mathit{u}}^T\mathbf{A}\mathit{u},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(37)}& \mathbf{A}=\sum _{m=\mathrm{1}}^M\sum _{m^{\prime }=\mathrm{1}}^M\left({\mathit{\varphi }}_m-{\mathit{\varphi }}_{m^{\prime }}\right){\left({\mathit{\varphi }}_m-{\mathit{\varphi }}_{m^{\prime }}\right)}^T.\end{array}$$

Our goal is to choose the initial condition u that maximizes Γ[u]. The initial condition u must be constrained in some way, otherwise Γ[u] is unbounded and there is no maximum. If the initial conditions are drawn from a distribution with covariance matrix Σ_M, then an appropriate constraint is to fix the associated Mahalanobis distance:

$$\begin{array}{}\text{(38)}& {\mathit{u}}^T{\mathbf{\Sigma }}_M^{-\mathrm{1}}\mathit{u}=\mathrm{1}.\end{array}$$

The problem is now to maximize u^TAu subject to the constraint ${\mathit{u}}^T{\mathbf{\Sigma }}_M^{-\mathrm{1}}\mathit{u}=\mathrm{1}$. This is a standard optimization problem that is merely a generalization of principal component analysis (also called empirical orthogonal function analysis). The solution is obtained by solving the generalized eigenvalue problem

$$\begin{array}{}\text{(39)}& \mathbf{A}\mathit{u}=\mathit{\lambda }{\mathbf{\Sigma }}_M^{-\mathrm{1}}\mathit{u}.\end{array}$$

This eigenvalue problem yields p eigenvalues and p eigenvectors. The eigenvalues can be ordered from largest to smallest, ${\mathit{\lambda }}_{\mathrm{1}}\ge {\mathit{\lambda }}_{\mathrm{2}}\ge \mathrm{\dots }\ge {\mathit{\lambda }}_p$, and the corresponding eigenvalues can be denoted ${\mathit{u}}_{\mathrm{1}},{\mathit{u}}_{\mathrm{2}},\mathrm{\dots },{\mathit{u}}_p$. The first eigenvector u₁ gives the initial condition that maximizes the sum square difference in predictable responses; the second eigenvector u₂ gives the initial condition that maximizes Γ subject to the condition that ${\mathit{u}}_{\mathrm{2}}^T{\mathbf{\Sigma }}_M^{-\mathrm{1}}{\mathit{u}}_{\mathrm{1}}=\mathrm{0}$; and so on. The eigenvalues ${\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_p$ give the corresponding values of Γ. The eigenvectors can be collected into the p×p matrix

$$\begin{array}{}\text{(40)}& \mathbf{U}=\left[\begin{array}{cccc}{\mathit{u}}_{\mathrm{1}}& {\mathit{u}}_{\mathrm{2}}& \mathrm{\dots }& {\mathit{u}}_p\end{array}\right].\end{array}$$

Because the matrices A and Σ_M are symmetric, the eigenvectors can be chosen to satisfy the orthogonality property

$$\begin{array}{}\text{(41)}& {\mathbf{U}}^T{\mathbf{\Sigma }}_M^{-\mathrm{1}}\mathbf{U}=\mathbf{I}\Rightarrow {\mathbf{\Sigma }}_M={\mathbf{UU}}^T.\end{array}$$

Summing Γ over all eigenvectors gives

$$\begin{array}{}\text{(42)}& \begin{array}{rl}\sum _{k=\mathrm{1}}^p\mathrm{\Gamma }\left[{\mathit{u}}_k\right]& =\sum _{k=\mathrm{1}}^p{\mathit{u}}_k^T\mathbf{A}{\mathit{u}}_k=\mathrm{tr}\left[\mathbf{A}\sum _{k=\mathrm{1}}^p{\mathit{u}}_k{\mathit{u}}_k^T\right]\\ & =\mathrm{tr}\left[{\mathbf{AUU}}^T\right]=\mathrm{tr}\left[\mathbf{A}{\mathbf{\Sigma }}_M\right]\\ & =\sum _{m=\mathrm{1}}^M\sum _{m^{\prime }=\mathrm{1}}^M{\left({\mathit{\varphi }}_m-{\mathit{\varphi }}_{m^{\prime }}\right)}^T\\ & {\mathbf{\Sigma }}_M\left({\mathit{\varphi }}_m-{\mathit{\varphi }}_{m^{\prime }}\right).\end{array}\end{array}$$

Note the similarity of the final expression to Eq. (27). This similarity implies that if only two models are compared and Σ_M=Σ_HM, then the solution exactly matches the variance of signal differences (27). Unfortunately, Σ_HM and the noise variance ${\widehat{\mathit{\sigma }}}^{\mathrm{2}}$ depend on $(m,m^{\prime })$, so it is difficult to generalize the optimal initial condition to explain all pair-wise values of F_ϕ∣σ. We suggest using a covariance matrix that is the harmonic mean across all M time series:

$$\begin{array}{}\text{(43)}& {\mathbf{\Sigma }}_M=M{\left({\widehat{\mathbf{\Sigma }}}_{\mathrm{1}}^{-\mathrm{1}}+{\widehat{\mathbf{\Sigma }}}_{\mathrm{2}}^{-\mathrm{1}}+\mathrm{\dots }+{\widehat{\mathbf{\Sigma }}}_M^{-\mathrm{1}}\right)}^{-\mathrm{1}},\end{array}$$

where ${\widehat{\mathbf{\Sigma }}}_m$ is the sample covariance matrix of the initial condition Eq. (25) for the mth time series. Finally, we note that

$$\begin{array}{}\text{(44)}& {\mathrm{\Gamma }}_{\text{sum}}=\sum _{k=\mathrm{1}}^p\mathrm{\Gamma }\left[{\mathit{u}}_k\right]={\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{2}}+\mathrm{\dots }+{\mathit{\lambda }}_p.\end{array}$$

Thus, λ_k∕Γ_sum is the fraction of sum total variance of signal differences explained by the kth eigenvector. Conceptually, each eigenvector u_k can be interpreted as an initial condition that has been optimized to separate predictable responses.

5 Example: diagnosing differences in AMOC simulations

To illustrate the proposed method, we apply it to an index of the Atlantic Meridional Overturning Circulation (AMOC). This variable is chosen because it is considered to be a major source of decadal variability and predictability (Buckley and Marshall, 2016; Zhang et al., 2019). However, the AMOC is not observed directly with sufficient frequency and consistency to constrain its variability on decadal timescales. As a result, there has been a heavy reliance on coupled atmosphere–ocean models to characterize AMOC variability and predictability. The question arises as to whether simulations of the AMOC by different models can be distinguished and, if so, how they differ.

Figure 1 AMOC index simulated by 10 CMIP5 models under pre-industrial conditions. Each time series is offset by a constant with no re-scaling. The AMOC index is defined as the maximum annual mean meridional overturning streamfunction at 40^∘ N in the Atlantic.

The data used in this study come from pre-industrial control runs from Phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2010). These control simulations lack year-to-year changes in forcing and thereby permit a focus on internal variability without confounding effects due to anthropogenic climate change. We consider only models that have pre-industrial control simulations spanning at least 500 years and contain meridional overturning circulation as an output variable. Only 10 models meet these criteria. An AMOC index is defined as the annual mean of the maximum meridional overturning streamfunction at 40^∘ N in the Atlantic. A third-order polynomial over the 500 years is removed to eliminate climate drift. The AMOC index from the 10 models is shown in Fig. 1. Based on visual comparisons, one might perceive differences in amplitude (e.g., the MPI models tend to have larger amplitudes than other models) and differences in the degree of persistence (e.g., high-frequency variability is more evident in INMCM4 than in other models), but whether these differences are statistically significant remains to be established.

Figure 2 Akaike's information criterion (AIC) as a function of autoregressive model order for AMOC time series from selected CMIP5 models. AIC is computed separately from the first and second 250 years of the data. A box identifies the minimum AIC up to order 11. The actual AIC values are shown – AIC values have not been offset. Models have been selected primarily to avoid line crossings.

To compare autoregressive processes, it is necessary to select the order of the autoregressive model. As discussed in Sect. 3, we use Eq. (22), which for ν=1.0 gives $p^{*}=\mathrm{5}$. One check on this choice is whether the residuals are white noise. We find that AR(5) is adequate for all models except MRI, in the sense that the residuals reveal no serious departures from white noise according to visual inspection and by the Ljung–Box test (Ljung and Box, 1978). Another check is to compute Akaike's information criterion (AIC; Box et al., 2008) for the two halves of the data (i.e., 250 years). Some representative examples are shown in Fig. 2. As can be seen, AIC is a relatively flat function of model order, hence small sampling variations can lead to large changes in order selection. Indeed, the order selected by AIC is sensitive to which half of the data is used. Nevertheless, because AIC is nearly a flat function of order, virtually any choice of order beyond AR(2) can be defended. The highest order selected by AIC is p=11. While we have performed our test for AR(11), this order would be a misleading illustration of the method because one might assume that 11 lags are necessary to identify model differences. Thus, in the results that follow, we choose AR(5) for all cases and bear in mind that comparisons with MRI may be affected by model inadequacy.

Figure 3 Auto-covariance function of the AMOC time series from each CMIP5 model, as estimated from Eq. (6) using the first 250 years of data.

The choice of AR(5) means that all information relevant to deciding differences in stationary processes is contained in the first six values of the sample auto-covariance function $\widehat{c}\left(\mathrm{0}\right),\mathrm{\dots },\widehat{c}\left(\mathrm{5}\right)$. The sample auto-covariance function for each AMOC index is shown in Fig. 3. Recall that the zero-lag auto-covariance $\widehat{c}\left(\mathrm{0}\right)$ is the sample variance. The figure suggests that the time series have different variances and different decay rates, but it is unclear whether these differences are significant. Note that a standard difference-in-variance test cannot be performed here because the time series are serially correlated. One might try to modify the standard F-test by adjusting the degrees of freedom, as is sometimes advocated in the t-test (Zwiers and von Storch, 1995), but the adjustment depends on the autocorrelation function that we are trying to compare. An alternative approach is to pre-whiten the data based on the AR fit and then test differences in variance, but this is exactly equivalent to our test based on F_σ.

Figure 4 Power spectra of AMOC time series from CMIP5 models. Spectra are estimated from the first and second 250 years of each time series using Daniell's estimator. The 95 % confidence interval and bandwidth are indicated by the error bars in the top right corner. The power spectra have been offset by a multiplicative constant to reduce overlap (the y axis is log scale). The longest resolved period is 60 years.

An alternative approach to comparing stationary processes is to compare power spectra. Power spectra of the AMOC time series are shown in Fig. 4. The spectra have been offset by a multiplicative constant to reduce overlap (otherwise the different curves would obscure each other). While many differences can be seen, the question is whether those differences are significant after accounting for sampling uncertainties. The two spectral tests discussed in Sect. 2 find relatively few differences between CMIP5 models, although they do indicate that time series from INMCM4 and MPI differ from those of other CMIP5 models (not shown). Presumably, these differences arise from the fact that the INMCM4 time series is closer to white noise and that the MPI time series have larger total variance than those from other CMIP5 models.

Figure 5 A measure of the “distance” between AMOC time series between the first and second halves of CMIP5 pre-industrial control simulations. The distance is measured by the bias-corrected deviance statistic D_ϕ,σ using an fifth-order AR model. The light and dark gray shadings show values exceeding the 5 % and 1 % significance levels, respectively (the threshold values are 12.7 and 17.0, respectively).

To illustrate our proposed method, we perform the following analysis. First, the AMOC index from each CMIP5 model is split into equal halves, each 250 years long. Then, each time series from the first half is compared to each time series in the second half. Some of these comparisons will involve time series from the same CMIP5 model. Our expectation is that no difference should be detected when time series from two different halves of the same CMIP5 model are compared. To summarize the comparisons, we show in Fig. 5 a matrix of the bias-corrected deviance statistic D_ϕ,σ for every possible model comparison. This statistic is a measure of the “distance” between stationary process. The two shades of gray indicate a significant difference in AR process at the 5 % or 1 % level. Values along the diagonal correspond to comparisons of AMOC time series from the same CMIP5 model. No significant difference is detected when time series come from the same model. In contrast, significant differences are detected in most of the cases when the time series come from different CMIP5 models. Interestingly, models from the same institution tend to be indistinguishable from each other. For instance, the MPI models are mostly indistinguishable from each other, and the CCSM and CESM-BGC models (from the National Center for Atmospheric Research, NCAR) are mostly indistinguishable from each other. We say “mostly” because the difference depends on which half of the simulation is used in the comparison. For instance, a difference is detected when comparing the first half of the CESM1-BGC time series to the second half to the CCSM4 time series, but not vice versa. The models have been ordered to make natural clusters easy to identify in this figure. Aside from a few exceptions, the pattern of shading is the same for AR(2) and AR(10) models (not shown), demonstrating that our conclusion regarding significant differences in AR process is not sensitive to model order. Also, the pattern of shading is similar if the comparison is based on only from the first half of the simulations, or only from the second half of the simulations (not shown), indicating that our results are not sensitive to sampling errors.

Figure 6 A set of points in a two-dimensional Cartesian plane whose pair-wise Euclidean distances most closely approximate the pair-wise distances between autoregressive processes of AMOC time series. The difference between AR processes is measured by the bias-corrected deviance statistic D_ϕ,σ and the points are identified using multidimensional scaling techniques. There are 20 points corresponding to 10 CMIP5 models, each model time series being split in half. Time series from the same model have the same color and are joined by a line segment. Circles around selected points enclose models whose time series are statistically indistinguishable from that of the center model at the 5 % significance level.

To visualize natural clusters more readily, we identify a set of points in a cartesian plane whose pair-wise distances match the above distances as closely as possible. These points can be identified using a procedure called multidimensional scaling (Izenman, 2013). The procedure is to compute the deviance statistic between every possible pair of time series. Because there are 10 CMIP models and time series from each model is split in half, there are 20 time series being compared. Thus, the deviance statistic for every possible pair can be summarized in a 20×20 distance matrix. From this matrix, multidimensional scaling can find the set of points in a 20-dimensional space that has this distance matrix. Moreover, it can find the points in a two-dimensional space whose distance matrix most closely approximates the original distance matrix in some objective sense. In our case, 87 % of original distance matrix can be represented by two-dimensional points. These points are shown in Fig. 6. Although 87 % of the distance matrix is represented in this figure, isolated points may have relatively large discrepancies. The average discrepancy is about 2 units, with the largest discrepancy being about 3.6 units between MRI and MPI-LR. An attractive feature of this representation is that the decision rule for statistical significance, $D_{\mathit{\varphi },\mathit{\sigma }}>D_{\text{crit}}$, is approximately equivalent to drawing a circle of radius $\sqrt{D_{\text{crit}}}$ around a point and identifying all the points that lie outside that circle. This geometric procedure would be exact in a 20-dimensional space, but is only approximate in the 2-dimensional space shown in Fig. 6. The figure suggests that the MPI models and INMCM4 form their own natural clusters. For other models, a particular choice of clusters is shown in the figure, but this is merely an example, and different clusters with different groupings could be justified. Note that all line segments connecting results from the same CMIP model are shorter than the circle radius, indicating no significant difference between time series from the same CMIP model.

It is interesting to relate the above differences to the auto-covariance functions shown in Fig. 3. It is likely that INMCM4 differs from other models because its auto-covariance function decays most quickly to zero. The MPI models are distinguished from the others by their large variances. Interestingly, note that the auto-covariance function for NorESM1-M is intermediate between that of two MPI models, yet the test indicates that NorESM1-M differs from the MPI models. Presumably, the MPI models are indistinguishable because their auto-covariance functions have the same shape, including a kink at lag 1, whereas the NorESM1-M model has no kink at lag 1. This example illustrates that the test does not behave like a mean square difference between auto-covariance functions, which would cluster NorESM1-M with the MPI models.

It is worth clarifying how the above clustering technique differs from those in previous studies. Piccolo (1990) proposed a distance measure based on the Euclidean norm of the difference in AR parameters, namely

$$\begin{array}{}\text{(45)}& d(X,Y)={\left\{\sum _{\mathit{\tau }=\mathrm{1}}^{\mathrm{\infty }}{\left({\mathit{\varphi }}_{\mathit{\tau }}^X-{\mathit{\varphi }}_{\mathit{\tau }}^Y\right)}^{\mathrm{2}}\right\}}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

In contrast, our hypothesis test uses a Mahalanobis norm for measuring differences in AR parameters, where the covariance matrix is based on the sample covariance matrices of the time-lagged data (see Eq. A44). While the Euclidean norm Eq. (45) does have some attractive properties as discussed in Piccolo (1990), it is inconsistent with the corresponding hypothesis test for differences in AR parameters. As a result, the resulting cluster may emphasize differences with large Euclidean norms that are insignificant, or may downplay differences in small Euclidean norms that are significant. In contrast, the distance measure used in our study is consistent with a rigorous hypothesis test. Maharaj (2000) proposes a classification procedure that is consistent with hypothesis testing, but that hypothesis test does not account for differences in noise variances.

Figure 7 Noise variance versus the R-square of AR(5) models estimated from 250-year segments of AMOC time series from 10 CMIP5 models. Estimates from the same CMIP5 model have the same color and are joined by a line segment. The error bar in the bottom left corner shows the critical distance for a significant difference in noise variance at the 5 % level. The x axis is on a log scale so that equal variance ratios correspond to equal linear distances. The error bar in the upper left shows a 95 % confidence interval for R-square using the method of Lee (1971) and computed from the MBESS package in R. The dashed line is the 5 % significance threshold for the R-square.

Having identified differences between stationary processes, it is of interest to relate those differences to standard properties of the time series. Recall that our measure of the difference between stationary processes D_ϕ,σ is the sum of two other measures, namely D_σ and D_ϕ∣σ. Measure D_σ depends only on the noise variances, that is, it depends on the mean square error of a 1-year prediction. In contrast, D_ϕ∣σ measures the difference in predictable responses of the process. As discussed in Sect. 4, we suggest examining differences in R-square. A graph of the noise variance plotted against R-square of each autoregressive model is shown in Fig. 7. The error bars show the critical distance for a significance difference in noise variance (lower left) and for a difference in R-square (upper right). First note that the projections of line segments onto the x axis or y axis are shorter than the respective error bars, indicating that differences in noise and differences in predictability are insignificant when estimated from the same CMIP5 model. Second, note that the relative positions of the dots are similar to those in Fig. 6. Thus, the noise variance and R-square appear to be approximate “principal coordinates” for distinguishing univariate autoregressive processes. Third, using noise variances alone, the MPI models would be grouped together, then INMCM4 would be grouped by itself, while at the bottom end CanESM2 and CCSM4 would be grouped together. These groupings are consistent with the clusters identified above using the full distance measure D_ϕ,σ, suggesting that differences in noise variances explain the major differences between stationary processes. Fourth, the AR models have R-square values mostly between 0.25 and 0.5. Time series from INMCM4 have the smallest R-square values while time series from NorESM1 have the largest R-square values.

Figure 8 Predictions from each estimated AR(5) model using the optimal initial condition derived from Eq. (39). The optimal initial condition is the five black dots joined by lines, and the resulting predictions are the colored lines. Each CMIP model has two predictions corresponding to the two AR(5) models estimated from two non-overlapping 250-year time series from the same CMIP model.

An alternative approach to describing differences in AR parameters is to show differences in response to the same initial condition. We use the optimization method discussed in Sect. 4 to find the initial condition that maximizes the sum square difference in responses. The result is shown in Fig. 8. Essentially, the models separate by predictability – the order of the models from top to bottom closely tracks the order of the models based on R-square seen in Fig. 7. For this initial condition, INMCM4 damps nearly to zero in one time step, whereas NorESM1-M decays the slowest among the models, consistent with expectations from R-square. This initial condition explains 82 % of the differences in response. Because optimal initial conditions form a complete set, an arbitrary initial condition can be represented as a linear combination of optimal initial conditions. If the AR models were used to predict an observational time series with covariance matrix Σ_M, then most differences between model predictions could be explained by the differences in response to a single optimal initial condition.

We end with a brief summary of our exploration of the χ² statistic proposed by Lund et al. (2009). As mentioned earlier, the result is sensitive to the choice of K. However, instead of summing based on the cutoff parameter K, we summed over all possible lags, and set all sample auto-covariances beyond lag p to zero. Using this approach, we find that the resulting χ² statistic gives results very similar to ours for p=5, including clustering the MPI models with each other, and separating them from NorEMS1-M. We have not investigated this procedure sufficiently to propose it as a general rule but mention it to suggest the possibility that some alternative rule for computing the sum Eq. (9) may yield reasonable results.

6 Conclusions

This paper examined tests for differences in stationary processes and proposed a new approach to characterizing those differences. The basic idea is to fit each time series to an autoregressive model and then test for differences in parameters. The likelihood ratio test for this comparison was derived in Maharaj (2000) and Grant and Quinn (2017). We have modified the test to correct for certain biases and to include a test for differences in noise variance. The latter test is of major interest in climate applications and leads to considerable simplifications in estimation and interpretation. Furthermore, the proposed test is applicable to small samples. In addition, we propose new approaches to interpreting and visualizing differences in stationary processes. The procedure was illustrated on AMOC time series from pre-industrial control simulations of 10 models in the CMIP5 archive. Based on time series 250 years in length, the procedure was able to distinguish about four clusters of models, where time series from the same CMIP5 model are grouped together in the same cluster. These clusters are identified easily using a multidimensional technique. The clusters obtained by this method were not sensitive to the order of the autoregressive model, although the number of significant differences decreases with AR order due to the larger number of parameters being estimated. Further analysis shows that these clusters can be explained largely by differences in 1-year prediction errors in the AR models and differences in R-square.

The proposed method can be used to compare any stationary time series that are well fit by an autoregressive model, which includes most climate time series. Thus, this method could be used to decide whether two climate models generate the same temporal variability. A natural question is whether this approach can be generalized to compare multivariate time series. This generalization will be developed in Part 2 of this paper. The method also could be used to compare model simulations to observations, provided that the stationarity assumption is satisfied. If non-stationarity is strong, then this method would need to be modified to account for such non-stationarity, such as adding exogenous terms to the AR model. The likelihood-ratio framework can accommodate such extensions and will be developed in Part 3 of this paper.

Appendix A: Derivation of the test

In this Appendix, we derive a test for equality of parameters of autoregressive models based on the likelihood ratio test. The derivation is similar to that in Maharaj (2000) and Grant and Quinn (2017), except modified to test for differences in noise variance and to correct for known biases in maximum likelihood estimates. In addition, we show how the bias-corrected likelihood ratio can be partitioned into two independent ratios and derive the sampling distributions for each.

We consider only the conditional likelihood, which is the likelihood function conditioned on the first p values of the process. The conditional likelihood approach is reasonable for large sample sizes and has the advantage that the estimates can be obtained straightforwardly from the method of least squares. Suppose we have a time series of length N_X for X_t. Then the conditional likelihood function for X_t is

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& L_X={\left(\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_X\right)}^{-(N_X-p)/\mathrm{2}}\\ & \exp \left[-{\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-{\mathit{\varphi }}_{\mathrm{1}}^X{X}_{t-\mathrm{1}}-\mathrm{\dots }-{\mathit{\varphi }}_p^X{X}_{t-p}-{\mathit{\gamma }}^X\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_X^{\mathrm{2}}}}\right]\end{array}\end{array}$$

(see Eq. A7.4.2b of Box et al., 2008). The maximum likelihood estimates of the parameters, denoted by ${\widehat{\mathit{\varphi }}}_{\mathrm{1}}^X,\mathrm{\dots },{\widehat{\mathit{\varphi }}}_p^X,{\widehat{\mathit{\gamma }}}^X$, are obtained from the method of least squares. The maximum likelihood estimate of ${\mathit{\sigma }}_X^{\mathrm{2}}$ is

$$\begin{array}{}\text{(A2)}& {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_X^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}^X{X}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p^X{X}_{t-p}-{\widehat{\mathit{\gamma }}}^X\right)}^{\mathrm{2}}}{N_X-p}}.\end{array}$$

Substituting this into the likelihood function Eq. (A1) and then taking −2 times the logarithm of the likelihood function gives

$$\begin{array}{}\text{(A3)}& -\mathrm{2}\log {\stackrel{\mathrm{‾}}{L}}_X=\left(N_X-p\right)\left(\log {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_X^{\mathrm{2}}+\log \mathrm{2}\mathit{\pi }+\mathrm{1}\right).\end{array}$$

Similarly, suppose we have a time series of length N_Y for Y_t. Then, the corresponding likelihood function is

$$\begin{array}{rl}& L_Y={\left(\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_Y\right)}^{-(N_Y-p)/\mathrm{2}}\\ & \exp \left[-{\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_Y}{\left(Y_t-{\mathit{\varphi }}_{\mathrm{1}}^Y{Y}_{t-\mathrm{1}}-\mathrm{\dots }-{\mathit{\varphi }}_p^Y{Y}_{t-p}-{\mathit{\gamma }}^Y\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_Y^{\mathrm{2}}}}\right],\end{array}$$

the maximum likelihood estimate of ${\mathit{\sigma }}_Y^{\mathrm{2}}$ is

$$\begin{array}{}\text{(A4)}& {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_Y^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_Y}{\left(Y_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}^Y{Y}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p^Y{Y}_{t-p}-{\widehat{\mathit{\gamma }}}^Y\right)}^{\mathrm{2}}}{N_Y-p}},\end{array}$$

and −2 times the logarithm of the likelihood function is

$$\begin{array}{}\text{(A5)}& -\mathrm{2}\log {\stackrel{\mathrm{‾}}{L}}_Y=\left(N_Y-p\right)\left(\log {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_Y^{\mathrm{2}}+\log \mathrm{2}\mathit{\pi }+\mathrm{1}\right).\end{array}$$

Because X_t and Y_t are independent, the likelihood function for all the data is the product L_XL_Y.

Under hypothesis H₀, the likelihood function has the same form as L_XL_Y, except there is only a single set of autoregressive parameters ${\mathit{\varphi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\varphi }}_p$ and a single noise variance σ. The corresponding likelihood function is therefore

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& L_{\mathit{\sigma },\mathit{\varphi }}={\left(\mathrm{2}\mathit{\pi }\mathit{\sigma }\right)}^{-(N_X-p)/\mathrm{2}}{\left(\mathrm{2}\mathit{\pi }\mathit{\sigma }\right)}^{-(N_Y-p)/\mathrm{2}}\\ & \exp \left[{\displaystyle \frac{{\displaystyle \begin{array}{c}-{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-{\mathit{\varphi }}_{\mathrm{1}}{X}_{t-\mathrm{1}}-\mathrm{\dots }-{\mathit{\varphi }}_p{X}_{t-p}-{\mathit{\gamma }}^X\right)}^{\mathrm{2}}\\ -{\sum }_{t=p+\mathrm{1}}^{N_Y}{\left(Y_t-{\mathit{\varphi }}_{\mathrm{1}}{Y}_{t-\mathrm{1}}-\mathrm{\dots }-{\mathit{\varphi }}_p{Y}_{t-p}-{\mathit{\gamma }}^Y\right)}^{\mathrm{2}}\end{array}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right].\end{array}\end{array}$$

The maximum likelihood estimates of ${\mathit{\varphi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\varphi }}_p$, denoted ${\widehat{\mathit{\varphi }}}_{\mathrm{1}},\mathrm{\dots },{\widehat{\mathit{\varphi }}}_p$, are obtained from the least squares estimates of the pooled sample. Again, these estimates are obtained easily by the method of least squares. The maximum likelihood estimate of the common variance σ² is

$$\begin{array}{}\text{(A7)}& {\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}^{\mathrm{2}}={\displaystyle \frac{{\displaystyle \begin{array}{c}{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}{X}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p{X}_{t-p}-{\widehat{\mathit{\gamma }}}^X\right)}^{\mathrm{2}}\\ +{\sum }_{t=p+\mathrm{1}}^{N_Y}{\left(Y_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}{Y}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p{Y}_{t-p}-{\widehat{\mathit{\gamma }}}^Y\right)}^{\mathrm{2}}\end{array}}}{N_X+N_Y-\mathrm{2}p}}.\end{array}$$

The corresponding log-likelihood is

$$\begin{array}{}\text{(A8)}& -\mathrm{2}\log {\stackrel{\mathrm{‾}}{L}}_{\mathit{\sigma },\mathit{\varphi }}=(N_X+N_Y-\mathrm{2}p)\left(\log {\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}^{\mathrm{2}}+\log \mathrm{2}\mathit{\pi }+\mathrm{1}\right).\end{array}$$

Finally, we compute the likelihood ratio, or equivalently, the difference in the log-likelihood functions. This difference (multiplied by −2) is called the deviance statistic and is

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{D}}_{\mathit{\sigma },\mathit{\varphi }}=& \mathrm{2}\log {\stackrel{\mathrm{‾}}{L}}_X+\mathrm{2}\log {\stackrel{\mathrm{‾}}{L}}_Y-\mathrm{2}\log {\stackrel{\mathrm{‾}}{L}}_{\mathit{\sigma },\mathit{\varphi }}\\ & =\log \left({\displaystyle \frac{{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}^{\mathrm{2}(N_X+N_Y-\mathrm{2}p)}}{{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_X^{\mathrm{2}(N_X-P)}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_Y^{\mathrm{2}(N_Y-p)}}}\right).\end{array}\end{array}$$

It is well known that maximum likelihood estimates of variance are biased. Accordingly, we define bias-corrected deviance statistics before proceeding. This can be done by replacing sample sizes by degrees of freedom. Care must be exercised in such substitutions because replacing MLEs with unbiased versions will yield likelihoods that are no longer maximized, and therefore may yield deviance statistics that are negative. Our goal is to define a bias-corrected deviance statistic that is non-negative

To compute the degrees of freedom, note that the sample size of X_t is N_X−p, because the first p time steps are excluded from the conditional likelihood, and p+1 parameters are being estimated, so the degrees of freedom is the difference, $N_X-\mathrm{2}p-\mathrm{1}$. Similarly, the degrees of freedom for Y_t is $N_Y-\mathrm{2}p-\mathrm{1}$. Let the degrees of freedom be denoted

$$\begin{array}{}\text{(A10)}& {\mathit{\nu }}_X=N_X-\mathrm{2}p-\mathrm{1}\text{and}{\mathit{\nu }}_Y=N_Y-\mathrm{2}p-\mathrm{1}.\end{array}$$

Accordingly, an unbiased estimate of ${\mathit{\sigma }}_X^{\mathrm{2}}$ is

$$\begin{array}{}\text{(A11)}& {\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_X}{\left(X_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}^X{X}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p^X{X}_{t-p}-{\widehat{\mathit{\gamma }}}^X\right)}^{\mathrm{2}}}{{\mathit{\nu }}_X}},\end{array}$$

and an unbiased estimate of ${\mathit{\sigma }}_Y^{\mathrm{2}}$ is

$$\begin{array}{}\text{(A12)}& {\displaystyle }{\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}={\displaystyle \frac{{\sum }_{t=p+\mathrm{1}}^{N_Y}{\left(Y_t-{\widehat{\mathit{\varphi }}}_{\mathrm{1}}^Y{Y}_{t-\mathrm{1}}-\mathrm{\dots }-{\widehat{\mathit{\varphi }}}_p^Y{Y}_{t-p}-{\widehat{\mathit{\gamma }}}^Y\right)}^{\mathrm{2}}}{{\mathit{\nu }}_Y}}.\end{array}$$

As will be shown below, the appropriate bias correction for ${\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}^{\mathrm{2}}$ is

$$\begin{array}{}\text{(A13)}& {\widehat{\widehat{\mathit{\sigma }}}}^{\mathrm{2}}={\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}^{\mathrm{2}}\left({\displaystyle \frac{N_X+N_Y-\mathrm{2}p}{{\mathit{\nu }}_X+{\mathit{\nu }}_Y}}\right).\end{array}$$

Also, it proves helpful to define the unbiased estimate of the common variance σ² as

$$\begin{array}{}\text{(A14)}& {\widehat{\mathit{\sigma }}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\nu }}_X{\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}+{\mathit{\nu }}_Y{\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}}{{\mathit{\nu }}_X+{\mathit{\nu }}_Y}}.\end{array}$$

Then, the bias-corrected deviance statistic Eq. (A9) can be defined as

$$\begin{array}{}\text{(A15)}& D_{\mathit{\varphi },\mathit{\sigma }}=D_{\mathit{\sigma }}+D_{\mathit{\varphi }\mid \mathit{\sigma }},\end{array}$$

where

$$\begin{array}{}\text{(A16)}& {\displaystyle }{D}_{\mathit{\sigma }}& {\displaystyle }=\log \left({\displaystyle \frac{{\widehat{\mathit{\sigma }}}^{\mathrm{2}{\mathit{\nu }}_X+\mathrm{2}{\mathit{\nu }}_Y}}{{\widehat{\mathit{\sigma }}}_X^{\mathrm{2}{\mathit{\nu }}_X}{\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}{\mathit{\nu }}_Y}}}\right)\text{(A17)}& {\displaystyle }{D}_{\mathit{\varphi }\mid \mathit{\sigma }}& {\displaystyle }=({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\log \left({\displaystyle \frac{{\widehat{\widehat{\mathit{\sigma }}}}^{\mathrm{2}}}{{\widehat{\mathit{\sigma }}}^{\mathrm{2}}}}\right).\end{array}$$

We now prove the D_σ and D_ϕ∣σ are non-negative, independent, and have sampling distributions related to the F-distribution. To show that D_σ is non-negative, note that the ratio in Eq. (A16) is the weighted arithmetic mean over the weighted geometric mean of ${\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}$ and ${\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}$. Consequently, we may invoke the AM-GM inequality to show that D_σ is non-negative and vanishes if and only if ${\widehat{\mathit{\sigma }}}_X={\widehat{\mathit{\sigma }}}_Y$. In this sense, D_σ measures the “distance” between ${\widehat{\mathit{\sigma }}}_X$ and ${\widehat{\mathit{\sigma }}}_Y$. Incidentally, had we used the uncorrected likelihoods, the resulting deviance statistic would vanish at ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_X^{\mathrm{2}}={\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_Y^{\mathrm{2}}$, which would give a biased measure of deviance.

To prove the remaining properties of D_σ and D_ϕ∣σ, we adopt a vector notation that is better suited to the task. Accordingly, let the AR(p) model Eq. (11) be denoted

$$\begin{array}{}\text{(A18)}& {\mathit{w}}_X={\mathbf{Z}}_X{\mathit{\varphi }}_X+\mathit{j}{\mathit{\gamma }}_X+{\mathit{ϵ}}_X,\end{array}$$

where w_X is an (N_X−p)-dimensional vector, Z_X is a $(N_X-p)\times p$ matrix, j is a (N_X−p)-dimensional vector of ones, γ_x is a scalar, and the remaining terms have been defined previously.

$$\begin{array}{}\text{(A19)}& \begin{array}{rl}{\mathit{w}}_X=& \left[\begin{array}{c}{X}_{p+\mathrm{1}}\\ X_{p+\mathrm{2}}\\ \mathrm{⋮}\\ X_{N_X}\end{array}\right],\\ & {\mathbf{Z}}_X=\left(\begin{array}{cccc}{X}_p& X_{p-\mathrm{1}}& \mathrm{\dots }& X_{\mathrm{1}}\\ X_{p+\mathrm{1}}& X_p& \mathrm{\dots }& X_{\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ X_{N_X-\mathrm{1}}& X_{N_X-\mathrm{2}}& \mathrm{\dots }& X_{N_X-p}\end{array}\right),\\ & {\mathit{\varphi }}_X=\left(\begin{array}{c}{\mathit{\varphi }}_{\mathrm{1}}^X\\ {\mathit{\varphi }}_{\mathrm{2}}^X\\ \mathrm{⋮}\\ {\mathit{\varphi }}_p^X\end{array}\right).\end{array}\end{array}$$

Since H₀ does not restrict γ_X, it proves convenient to eliminate this parameter by projecting onto the complement of j. Therefore, we multiply both sides of the equation by the projection matrix

$$\begin{array}{}\text{(A20)}& \mathbf{H}=\mathbf{I}-{\displaystyle \frac{\mathrm{1}}{N_X-p}}\mathit{j}{\mathit{j}}^T.\end{array}$$

This multiplication has the effect of eliminating the γ_X term and centering each column of w_X and Z_X. Henceforth, we assume that each column of w_X and Z_X has been centered. One should remember that the degrees of freedom associated with estimating the noise variance ${\mathit{\sigma }}_X^{\mathrm{2}}$ should be reduced by one to account for this pre-centering of the data. Similarly, the corresponding model for Y_t in Eq. (12) is written as

$$\begin{array}{}\text{(A21)}& {\mathit{w}}_Y={\mathbf{Z}}_Y{\mathit{\varphi }}_Y+\mathit{j}{\mathit{\gamma }}_Y+{\mathit{ϵ}}_Y,\end{array}$$

where w_Y is an (N_Y−p)-dimensional vector, and the remaining terms are analogous to those in Eq. (A18). As for X, we assume that each column of w_Y and Z_Y is centered. The least squares estimates for ϕ_X and ϕ_Y are

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}& {\widehat{\mathit{\varphi }}}_X={\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X\right)}^{-\mathrm{1}}{\mathbf{Z}}_X^T{\mathit{w}}_X\\ & \text{and}{\widehat{\mathit{\varphi }}}_Y={\left({\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}{\mathbf{Z}}_Y^T{\mathit{w}}_Y,\end{array}\end{array}$$

and the corresponding sum square errors are

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}& {\text{SSE}}_X=\Vert {\mathit{w}}_X-{\mathbf{Z}}_X{\widehat{\mathit{\varphi }}}_X{\Vert }^{\mathrm{2}}\\ & \text{and}{\text{SSE}}_Y=\Vert {\mathit{w}}_Y-{\mathbf{Z}}_Y{\widehat{\mathit{\varphi }}}_Y{\Vert }^{\mathrm{2}},\end{array}\end{array}$$

where $\Vert \mathit{a}\Vert ={\mathit{a}}^T\mathit{a}$ denotes the Euclidean norm of vector a. The sum square errors are related to the estimated variances as

$$\begin{array}{}\text{(A24)}& {\displaystyle }{\text{SSE}}_X& {\displaystyle }=(N_X-p){\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_X^{\mathrm{2}}={\mathit{\nu }}_X{\widehat{\mathit{\sigma }}}_X^{\mathrm{2}},\text{(A25)}& {\displaystyle }{\text{SSE}}_Y& {\displaystyle }=(N_Y-p){\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_Y^{\mathrm{2}}={\mathit{\nu }}_Y{\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}.\end{array}$$

Following the standard theory of least squares estimation for the general linear model, Brockwell and Davis (1991; Sect. 8.9) suggest that the sum square errors have the following approximate chi-squared distributions:

$$\begin{array}{}\text{(A26)}& {\displaystyle \frac{{\text{SSE}}_X}{{\mathit{\sigma }}_X^{\mathrm{2}}}}\sim {\mathit{\chi }}_{{\mathit{\nu }}_X}^{\mathrm{2}},\text{and}{\displaystyle \frac{{\text{SSE}}_Y}{{\mathit{\sigma }}_Y^{\mathrm{2}}}}\sim {\mathit{\chi }}_{{\mathit{\nu }}_Y}^{\mathrm{2}}.\end{array}$$

Because X_t and Y_t are independent, the associated sum square errors are independent, hence under H₀, the statistic

$$\begin{array}{}\text{(A27)}& F_{\mathit{\sigma }}={\displaystyle \frac{{\widehat{\mathit{\sigma }}}_X^{\mathrm{2}}}{{\widehat{\mathit{\sigma }}}_Y^{\mathrm{2}}}}\end{array}$$

has an F-distribution with (ν_X,ν_Y) degrees of freedom, as stated in Eq. (20). Furthermore, if H₀ is true, then

$$\begin{array}{}\text{(A28)}& \begin{array}{rl}{\text{SSE}}_{\mathit{\sigma }}& ={\text{SSE}}_X+{\text{SSE}}_Y=(N_X+N_Y-\mathrm{2}p){\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{2}}\\ & =({\mathit{\nu }}_X+{\mathit{\nu }}_Y){\widehat{\mathit{\sigma }}}^{\mathrm{2}}\end{array}\end{array}$$

has a (scaled) chi-squared distribution with ν_X+ν_Y degrees of freedom; i.e.,

$$\begin{array}{}\text{(A29)}& {\displaystyle \frac{{\text{SSE}}_{\mathit{\sigma }}}{{\mathit{\sigma }}^{\mathrm{2}}}}\sim {\mathit{\chi }}_{{\mathit{\nu }}_X+{\mathit{\nu }}_Y}^{\mathrm{2}}.\end{array}$$

Straightforward algebra shows that D_ϕ in Eq. (A16) can be written as a function of F_σ as

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}{D}_{\mathit{\sigma }}\left(F_{\mathit{\sigma }}\right)=& {\mathit{\nu }}_X\log \left({\mathit{\nu }}_X+{\mathit{\nu }}_Y/F_{\mathit{\sigma }}\right)+{\mathit{\nu }}_Y\log \left({\mathit{\nu }}_X{F}_{\mathit{\sigma }}+{\mathit{\nu }}_Y\right)\\ & -({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\log ({\mathit{\nu }}_X+{\mathit{\nu }}_Y).\end{array}\end{array}$$

This is a U-shaped function with a minimum value of zero at F_σ=1 and monotonic on either side of F_σ=1. Note that if X_t and Y_t were swapped, then the X and Y labels would be swapped and $F_{\mathit{\sigma }}\to \mathrm{1}/F_{\mathit{\sigma }}$, which yields the same value of D_σ. Let $F_{\mathit{\alpha },{\mathit{\nu }}_X,{\mathit{\nu }}_Y}$ denote the critical value such that $F>F_{\mathit{\alpha },{\mathit{\nu }}_X,{\mathit{\nu }}_Y}$ has probability α when F has an F-distribution with (ν_X,ν_Y) degrees of freedom. Then, the α100 % significance threshold for rejecting H₀ is

$$\begin{array}{}\text{(A31)}& D_{\mathit{\sigma },\mathit{\alpha }}=D_{\mathit{\sigma }}\left(F_{\mathit{\alpha }/\mathrm{2},{\mathit{\nu }}_X,{\mathit{\nu }}_Y}\right),\end{array}$$

where α on the right-hand side is divided by 2.

The least squares estimate of the common regression parameters ϕ is

$$\begin{array}{}\text{(A32)}& \begin{array}{rl}\widehat{\mathit{\varphi }}& ={\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X+{\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}\left({\mathbf{Z}}_X^T{\mathit{w}}_X+{\mathbf{Z}}_Y^T{\mathit{w}}_Y\right)\\ & ={\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X+{\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}\left(\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X\right){\widehat{\mathit{\varphi }}}_X\right.\\ & \left.+\left({\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right){\widehat{\mathit{\varphi }}}_Y\right),\end{array}\end{array}$$

and the corresponding sum square error is

$$\begin{array}{}\text{(A33)}& {\text{SSE}}_{\mathit{\varphi },\mathit{\sigma }}=\Vert {\mathit{w}}_X-{\mathbf{Z}}_X\widehat{\mathit{\varphi }}{\Vert }^{\mathrm{2}}+\Vert {\mathit{w}}_Y-{\mathbf{Z}}_Y\widehat{\mathit{\varphi }}{\Vert }^{\mathrm{2}}.\end{array}$$

We now express this in terms of parameter estimates of the individual AR models. To do this, we invoke the standard orthogonality relation

$$\begin{array}{}\text{(A34)}& {\displaystyle }\Vert {\mathit{w}}_X-{\mathbf{Z}}_X\widehat{\mathit{\varphi }}{\Vert }^{\mathrm{2}}& {\displaystyle }=\Vert {\mathit{w}}_X-{\mathbf{Z}}_X{\widehat{\mathit{\varphi }}}_X{\Vert }^{\mathrm{2}}+\Vert {\mathbf{Z}}_X({\widehat{\mathit{\varphi }}}_X-\widehat{\mathit{\varphi }}){\Vert }^{\mathrm{2}},\text{(A35)}& {\displaystyle }\Vert {\mathit{w}}_Y-{\mathbf{Z}}_Y\widehat{\mathit{\varphi }}{\Vert }^{\mathrm{2}}& {\displaystyle }=\Vert {\mathit{w}}_Y-{\mathbf{Z}}_Y{\widehat{\mathit{\varphi }}}_Y{\Vert }^{\mathrm{2}}+\Vert {\mathbf{Z}}_Y({\widehat{\mathit{\varphi }}}_Y-\widehat{\mathit{\varphi }}){\Vert }^{\mathrm{2}}.\end{array}$$

It follows that

$$\begin{array}{}\text{(A36)}& \begin{array}{rl}{\text{SSE}}_{\mathit{\varphi },\mathit{\sigma }}=& \left({\text{SSE}}_X+{\text{SSE}}_Y\right)+\Vert {\mathbf{Z}}_X({\widehat{\mathit{\varphi }}}_X-\widehat{\mathit{\varphi }}){\Vert }^{\mathrm{2}}\\ & +\Vert {\mathbf{Z}}_Y({\widehat{\mathit{\varphi }}}_Y-\widehat{\mathit{\varphi }}){\Vert }^{\mathrm{2}}.\end{array}\end{array}$$

The difference between least squares estimates ${\widehat{\mathit{\varphi }}}_X$ and $\widehat{\mathit{\varphi }}$ is

$$\begin{array}{}\text{(A37)}& \begin{array}{rl}{\widehat{\mathit{\varphi }}}_X-\widehat{\mathit{\varphi }}=& {\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X+{\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}\left({\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)\\ & \left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right),\end{array}\end{array}$$

while that between ${\widehat{\mathit{\varphi }}}_Y$ and $\widehat{\mathit{\varphi }}$ is

$$\begin{array}{}\text{(A38)}& \begin{array}{rl}{\widehat{\mathit{\varphi }}}_Y-\widehat{\mathit{\varphi }}=& -{\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X+{\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X\right)\\ & \left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right).\end{array}\end{array}$$

Substituting these expressions into Eq. (A36) and invoking Eq. (A28) yields

$$\begin{array}{}\text{(A39)}& {\text{SSE}}_{\mathit{\varphi },\mathit{\sigma }}={\text{SSE}}_{\mathit{\sigma }}+{\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}^T{\mathbf{\Sigma }}_{\mathrm{HM}}\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right),\end{array}$$

where

$$\begin{array}{}\text{(A40)}& {\mathbf{\Sigma }}_{\mathrm{HM}}={\left({\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X\right)}^{-\mathrm{1}}+{\left({\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}\right)}^{-\mathrm{1}}.\end{array}$$

Σ_HM is proportional to the harmonic mean of two covariance matrices. This matrix arises naturally if we recall the fact that the sample estimates of the parameters have the following distributions:

$$\begin{array}{}\text{(A41)}& \begin{array}{rl}& {\widehat{\mathit{\varphi }}}_X\sim \mathcal{N}\left({\mathit{\varphi }}_X,{\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X\right)}^{-\mathrm{1}}{\mathit{\sigma }}_X^{\mathrm{2}}\right)\\ & \text{and}{\widehat{\mathit{\varphi }}}_Y\sim \mathcal{N}\left({\mathit{\varphi }}_Y,{\left({\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}{\mathit{\sigma }}_Y^{\mathrm{2}}\right).\end{array}\end{array}$$

Therefore, under H₀,

$$\begin{array}{}\text{(A42)}& {\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\sim \mathcal{N}\left(\mathbf{0},\left({\left({\mathbf{Z}}_X^T{\mathbf{Z}}_X\right)}^{-\mathrm{1}}+{\left({\mathbf{Z}}_Y^T{\mathbf{Z}}_Y\right)}^{-\mathrm{1}}\right){\mathit{\sigma }}^{\mathrm{2}}\right),\end{array}$$

hence

$$\begin{array}{}\text{(A43)}& {\displaystyle \frac{{\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}^T{\mathbf{\Sigma }}_{\mathrm{HM}}\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}{{\mathit{\sigma }}^{\mathrm{2}}}}\sim {\mathit{\chi }}_p^{\mathrm{2}}.\end{array}$$

Again, by analogy with the standard theory of least squares estimation for the general linear model, SSE_σ and the quadratic form in Eq. (A39) are approximately independent of each other. Thus, if H₀ is true, then

$$\begin{array}{}\text{(A44)}& F_{\mathit{\varphi }\mid \mathit{\sigma }}={\displaystyle \frac{\mathrm{1}}{p}}{\displaystyle \frac{{\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}^T{\mathbf{\Sigma }}_{\mathrm{HM}}\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}{{\widehat{\mathit{\sigma }}}^{\mathrm{2}}}}\end{array}$$

has an approximate F-distribution with $(p,{\mathit{\nu }}_X+{\mathit{\nu }}_Y)$ degrees of freedom, as indicated in Eq. (21). Using the identity

$$\begin{array}{}\text{(A45)}& {\text{SSE}}_{\mathit{\varphi },\mathit{\sigma }}=(N_X+N_Y-\mathrm{2}p){\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}}^{\mathrm{2}}=({\mathit{\nu }}_X+{\mathit{\nu }}_Y){\widehat{\widehat{\mathit{\sigma }}}}^{\mathrm{2}},\end{array}$$

it follows that

$$\begin{array}{}\text{(A46)}& {\widehat{\widehat{\mathit{\sigma }}}}^{\mathrm{2}}={\widehat{\mathit{\sigma }}}^{\mathrm{2}}+{\displaystyle \frac{{\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}^T{\mathbf{\Sigma }}_{\mathrm{HM}}\left({\widehat{\mathit{\varphi }}}_X-{\widehat{\mathit{\varphi }}}_Y\right)}{{\mathit{\nu }}_X+{\mathit{\nu }}_Y}},\end{array}$$

and therefore

$$\begin{array}{}\text{(A47)}& \begin{array}{rl}{D}_{\mathit{\varphi }\mid \mathit{\sigma }}\left(F_{\mathit{\varphi }\mid \mathit{\sigma }}\right)& =({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\log \left({\displaystyle \frac{{\widehat{\widehat{\mathit{\sigma }}}}^{\mathrm{2}}}{{\widehat{\mathit{\sigma }}}^{\mathrm{2}}}}\right)\\ & =({\mathit{\nu }}_X+{\mathit{\nu }}_Y)\log \left(\mathrm{1}+{\displaystyle \frac{p{F}_{\mathit{\varphi }\mid \mathit{\sigma }}}{{\mathit{\nu }}_X+{\mathit{\nu }}_Y}}\right).\end{array}\end{array}$$

It is clear that F_ϕ∣σ is non-negative, and therefore D_ϕ∣σ is non-negative. Furthermore, D_ϕ∣σ is a monotonic function of F_ϕ|σ. Therefore, the α100 % significance threshold for rejecting H₀ is

$$\begin{array}{}\text{(A48)}& D_{\mathit{\varphi }\mid \mathit{\sigma },\mathit{\alpha }}=D_{\mathit{\varphi }\mid \mathit{\sigma }}\left(F_{\mathit{\alpha },\mathit{\varphi }\mid \mathit{\sigma }}\right).\end{array}$$

It remains to define the critical value for rejecting H₀ based on D_ϕ,σ. D_σ and D_ϕ∣σ are independent because D_σ depends on the ratio of ${\mathit{\chi }}_X^{\mathrm{2}}$ and ${\mathit{\chi }}_Y^{\mathrm{2}}$, while D_ϕ∣σ depends on the sum ${\mathit{\chi }}_X^{\mathrm{2}}+{\mathit{\chi }}_Y^{\mathrm{2}}$, and the ratio and sum are independent as a consequence of Lukacs' proportion-sum independence theorem (Lukacs, 1955). Therefore, we may sample from the F-distributions (20) and (21) and then use Monte Carlo techniques to estimate the upper α100th percentile of $D_{\mathit{\varphi },\mathit{\sigma }}=D_{\mathit{\sigma }}+D_{\mathit{\varphi }\mid \mathit{\sigma }}$ under H₀.