Comparing climate time series – Part 3: Discriminant analysis

In parts I and II of this paper series, rigorous tests for equality of stochastic processes were proposed. These tests provide objective criteria for deciding whether two processes differ, but they provide no information about the nature of those differences. This paper develops a systematic and optimal approach to diagnosing differences between multivariate stochastic processes. Like the tests, the diagnostics are framed in terms of vector autoregressive (VAR) models, which can be viewed as a dynamical system forced by random noise. The tests depend on two statistics, one that measures dissimilarity in dynamical operators and another that measures dissimilarity in noise covariances. Under suitable assumptions, these statistics are independent and can be tested separately for significance. If a term is significant, then the linear combination of variables that maximizes that term is obtained. The resulting indices contain all relevant information about differences between data sets. These techniques are applied to diagnose how the variability of annual-mean North Atlantic sea surface temperature differs between climate models and observations. For most models, differences in both noise processes and dynamics are important. Over 40 % of the differences in noise statistics can be explained by one or two discriminant components, though these components can be model dependent. Maximizing dissimilarity in dynamical operators identifies situations in which some climate models predict large-scale anomalies with the wrong sign.

1 Introduction

Comparing time series is a common task in climate studies, yet there are few rigorous criteria for deciding whether two time series come from the same stationary process. Such comparisons are vital in climate studies, particularly for deciding whether climate models produce realistic simulations. Part I of this study (DelSole and Tippett, 2020) derived a test for the hypothesis that two univariate time series come from the same stationary autoregressive process. Part II generalized this test to multivariate time series (DelSole and Tippett, 2021a). However, the decision that two processes differ provides no information about the nature of those differences. The purpose of this paper is to develop optimal diagnostics of differences between multivariate stochastic processes. Such diagnostics may help to identify errors in climate models and suggest ways to correct them.

The starting point of our approach is the objective criterion derived in Part II for deciding whether two multivariate time series come from the same process. This criterion assumes that each multivariate time series comes from a vector autoregressive (VAR) model. Part II derived the likelihood ratio test for equality of VAR parameters. The associated test statistic, called a deviance statistic, depends on all parameters in the VAR models, as it should, but in specific applications it may be dominated by particular aspects of the VAR models, such as differences in noise processes or differences in dynamics. In this paper, we first show that the deviance can be partitioned into two terms such that one term measures differences in noise processes and the other term measures differences in dynamics. Under suitable assumptions, these terms are independent. Then, the individual terms are further decomposed using discriminant analysis techniques. Discriminant analysis finds linear combinations of variables that maximize some quantity.

When applied to diagnose dissimilarity in noise processes, discriminant analysis yields a procedure that is similar to principal component analysis and other related techniques for diagnosing variability (Schneider and Griffies, 1999; DelSole et al., 2011; Wills et al., 2018). When applied to diagnose dissimilarity in dynamics, discriminant analysis yields a procedure that is equivalent to generalized singular value decomposition. Although both procedures have appeared in previous studies, no rigorous significance test accompanied these applications. In this work, these procedures emerge naturally from a rigorous statistical framework, and, as a result, they each have a well-defined significance test. A significance test is possible because each dissimilarity measure is based on prewhitened variables rather than on the original, serially correlated variables.

These new analysis techniques are applied to diagnose how variability of annual-mean North Atlantic sea surface temperature (SST) differs between observations and model simulations. It is found that differences in noise processes and differences in dynamics both play a significant role in explaining the overall differences. The components that dominate these differences are illustrated for a few cases. This paper closes with a summary and discussion.

2 Likelihood ratio test for differences in vector autoregressive models

This section briefly reviews the procedure derived in Part II for deciding whether two time series come from the same stochastic process. This procedure considers only stationary, Gaussian processes. Such processes are described completely by a mean vector and time-lagged covariance matrices. Standard tests for differences in these quantities (see chapter 10 of Anderson, 1984) assume independent data and therefore are not suitable for comparing serially correlated multivariate time series. Our approach is to fit a VAR model to data and then test equality of model parameters. Accordingly, consider a J-dimensional random vector z_t, where t is a time index. A VAR model for z_t is

$$\begin{array}{}\text{(1)}& {\mathit{z}}_t={\mathbf{A}}_{\mathrm{1}}{\mathit{z}}_{t-\mathrm{1}}+\mathrm{\dots }+{\mathbf{A}}_p{\mathit{z}}_{t-p}+\mathit{\mu }+{\mathit{ϵ}}_t,\end{array}$$

where ${\mathbf{A}}_{\mathrm{1}},\mathrm{\dots },{\mathbf{A}}_p$ are constant J×J matrices, μ is a constant vector that controls the mean of z_t, and ϵ_t is a Gaussian white noise process with covariance Γ. We refer to ϵ_t as the noise and ${\mathbf{A}}_{\mathrm{1}},\mathrm{\dots },{\mathbf{A}}_p$ as the AR parameters. The term model parameters refers to the set $\mathit{\{}\mathit{\mu },\mathbf{\Gamma },{\mathbf{A}}_{\mathrm{1}},\mathrm{\dots },{\mathbf{A}}_p\mathit{\}}$, which involves a total of $J+J(J+\mathrm{1})/\mathrm{2}+J^{\mathrm{2}}p$ independent parameters.

Typically, AR parameters characterize dynamics. For instance, some stochastic models are derived by linearizing the equations of motion about a time mean flow and parameterizing the eddy–eddy nonlinear interactions by a stochastic forcing and linear damping (Farrell and Ioannou, 1995; DelSole, 2004). Similar stochastic models can be derived by linear inverse methods (Penland, 1989; Winkler et al., 2001). Majda et al. (2006) showed that such stochastic models can be derived under suitable assumptions of a timescale separation (see the review by Gottwald et al., 2017). Realizations from these stochastic models, sampled at discrete time intervals, are indistinguishable from realizations of VAR models. In each of these models, the AR parameters characterize physical processes such as advection, dissipation, and wave propagation. Accordingly, we say that AR parameters characterize the dynamics of the system and that differences in AR parameters indicate differences in dynamics.

The model parameters can be estimated by the method of least squares, and the distributions of the estimates are asymptotically normal (Lütkepohl, 2005, Sect. 3.2.2). Accordingly, it proves convenient to express the VAR model as a regression model. Suppose Eq. (1) is used to generate N+p realizations ${\mathit{z}}_{\mathrm{1}},\mathrm{\dots },{\mathit{z}}_{N+p}$. Then, these realizations can be modeled by

$${\mathbf{Y}}^{*}={\mathbf{X}}^{*}\mathbf{B}+\mathbf{1}{\mathit{\mu }}^T+\mathbf{E},$$

where 1 is an N-dimensional vector of ones, superscript T denotes the transpose operation, the rows of E are independent realizations of ϵ_t, and

$$\begin{array}{c}{\displaystyle }{\mathbf{Y}}^{*}=\left(\begin{array}{c}{\mathit{z}}_{p+\mathrm{1}}^T\\ \\ {\mathit{z}}_{p+\mathrm{2}}^T\\ \\ \mathrm{⋮}\\ \\ {\mathit{z}}_{N+p}^T\end{array}\right)\in {\mathcal{R}}^{N\times J},\\ {\displaystyle }{\mathbf{X}}^{*}=\left(\begin{array}{cccc}{\mathit{z}}_p^T& {\mathit{z}}_{p-\mathrm{1}}^T& \mathrm{\dots }& {\mathit{z}}_{\mathrm{1}}^T\\ \\ {\mathit{z}}_{p+\mathrm{1}}^T& {\mathit{z}}_p^T& \mathrm{\dots }& {\mathit{z}}_{\mathrm{2}}^T\\ \\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \\ {\mathit{z}}_{N+p-\mathrm{1}}^T& {\mathit{z}}_{N+p-\mathrm{2}}^T& \mathrm{\dots }& {\mathit{z}}_N^T\end{array}\right)\in {\mathcal{R}}^{N\times Jp},\\ {\displaystyle }\mathbf{B}=\left(\begin{array}{c}{\mathbf{A}}_{\mathrm{1}}^T\\ \\ {\mathbf{A}}_{\mathrm{2}}^T\\ \\ \mathrm{⋮}\\ \\ {\mathbf{A}}_p^T\end{array}\right)\in {\mathcal{R}}^{Jp\times J}.\end{array}$$

Inferences based on the above regression model are identical to inferences based on the model $\mathbf{Y}=\mathbf{XB}+\mathbf{E}$, where in the latter model μ is omitted and X and Y contain centered variables in which the mean of each column is subtracted from that column. To simplify the analysis, we hereafter consider regression models of the form $\mathbf{Y}=\mathbf{XB}+\mathbf{E}$ comprising centered variables. To account for centering, 1 degree of freedom per dimension should be removed at the appropriate steps below, but for brevity they will not be discussed further. With this understanding, we hereafter express our two VAR models in the form

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathbf{Y}}_{\mathrm{1}}={\mathbf{X}}_{\mathrm{1}}{\mathbf{B}}_{\mathrm{1}}+{\mathbf{E}}_{\mathrm{1}},\text{(3)}& {\displaystyle }{\mathbf{Y}}_{\mathrm{2}}={\mathbf{X}}_{\mathrm{2}}{\mathbf{B}}_{\mathrm{2}}+{\mathbf{E}}_{\mathrm{2}},\end{array}$$

where

$$\begin{array}{c}{\displaystyle }{\mathbf{Y}}_{\mathrm{1}}\in {\mathcal{R}}^{N_{\mathrm{1}}\times J},{\mathbf{X}}_{\mathrm{1}}\in {\mathcal{R}}^{N_{\mathrm{1}}\times M},{\mathbf{B}}_{\mathrm{1}}\in {\mathcal{R}}^{M\times J},\\ {\displaystyle }{\mathbf{Y}}_{\mathrm{2}}\in {\mathcal{R}}^{N_{\mathrm{2}}\times J},{\mathbf{X}}_{\mathrm{2}}\in {\mathcal{R}}^{N_{\mathrm{2}}\times M},{\mathbf{B}}_{\mathrm{2}}\in {\mathcal{R}}^{M\times J}.\end{array}$$

M=Jp, and E₁ and E₂ are independent and distributed as

$$\begin{array}{c}{\displaystyle }\text{rows of }{\mathbf{E}}_{\mathrm{1}}\stackrel{iid}{\sim }{\mathcal{N}}_J\left(\mathrm{0},{\mathbf{\Gamma }}_{\mathrm{1}}\right),\\ {\displaystyle }\text{rows of }{\mathbf{E}}_{\mathrm{2}}\stackrel{iid}{\sim }{\mathcal{N}}_J\left(\mathrm{0},{\mathbf{\Gamma }}_{\mathrm{2}}\right).\end{array}$$

The null hypothesis is

$$H_{\mathrm{0}}:{\mathbf{B}}_{\mathrm{1}}={\mathbf{B}}_{\mathrm{2}}\text{and}{\mathbf{\Gamma }}_{\mathrm{1}}={\mathbf{\Gamma }}_{\mathrm{2}}.$$

Let B₀ and Γ₀ denote the common regression parameters and noise covariance matrix under H₀, respectively. The total number of parameters estimated under H₀ is $\mathrm{2}J+J(J+\mathrm{1})/\mathrm{2}+MJ$, derived by summing 2J for μ₁ and μ₂, MJ for B₀, and $J(J+\mathrm{1})/\mathrm{2}$ for Γ₀. On the other hand, the case of unconstrained VAR processes is denoted as

$$H_{\mathrm{2}}:\text{no restriction on }{\mathbf{B}}_{\mathrm{1}},{\mathbf{B}}_{\mathrm{2}},{\mathbf{\Gamma }}_{\mathrm{1}},{\mathbf{\Gamma }}_{\mathrm{2}}.$$

The total number of parameters estimated under H₂ is $\mathrm{2}J+J(J+\mathrm{1})+\mathrm{2}MJ$.

As discussed in Part II of this paper series, the bias-corrected likelihood ratio for testing H₀ versus H₂ is

$${\mathrm{\Lambda }}_{\mathrm{0}:\mathrm{2}}={\left({\displaystyle \frac{\left|{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{|}^{{\mathit{\nu }}_{\mathrm{1}}}\right|{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{|}^{{\mathit{\nu }}_{\mathrm{2}}}}{|{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}{|}^{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}}\right)}^{\mathrm{1}/\mathrm{2}},$$

where $|\cdot |$ denotes the determinant and

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\nu }}_{\mathrm{1}}=N_{\mathrm{1}}-M-\mathrm{1},\text{(5)}& {\displaystyle }{\mathit{\nu }}_{\mathrm{2}}=N_{\mathrm{2}}-M-\mathrm{1},\text{(6)}& {\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{1}}={\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}\right)}^{-\mathrm{1}}{\mathbf{X}}_{\mathrm{1}}^T{\mathbf{Y}}_{\mathrm{1}},\text{(7)}& {\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{2}}={\left({\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right)}^{-\mathrm{1}}{\mathbf{X}}_{\mathrm{2}}^T{\mathbf{Y}}_{\mathrm{2}},\text{(8)}& {\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{0}}={\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}+{\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{Y}}_{\mathrm{1}}+{\mathbf{X}}_{\mathrm{2}}^T{\mathbf{Y}}_{\mathrm{2}}\right),\text{(9)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}=\lfloor \lfloor {\mathbf{Y}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{1}}\rfloor \rfloor /{\mathit{\nu }}_{\mathrm{1}},\text{(10)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}=\lfloor \lfloor {\mathbf{Y}}_{\mathrm{2}}-{\mathbf{X}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{2}}\rfloor \rfloor /{\mathit{\nu }}_{\mathrm{2}},\text{(11)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}={\displaystyle \frac{{\displaystyle }\lfloor \lfloor {\mathbf{Y}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{0}}\rfloor \rfloor +\lfloor \lfloor {\mathbf{Y}}_{\mathrm{2}}-{\mathbf{X}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{0}}\rfloor \rfloor }{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}},\end{array}$$

and we have used the following operator for any matrix Z:

$$\lfloor \lfloor \mathbf{Z}\rfloor \rfloor ={\mathbf{Z}}^T\mathbf{Z}.$$

In our applications, ${\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}$ and ${\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}$ are positive definite, and this will be assumed hereafter.

MLEs and unbiased estimates are asymptotically equivalent. Hence, asymptotic theory (Hogg et al., 2019) gives the distribution of the deviance statistic under H₀ as

$$\begin{array}{}\text{(12)}& D_{\mathrm{0}:\mathrm{2}}=-\mathrm{2}\log {\mathrm{\Lambda }}_{\mathrm{0}:\mathrm{2}}\sim {\mathit{\chi }}_{MJ+J(J+\mathrm{1})/\mathrm{2}}^{\mathrm{2}}\text{when }{H}_{\mathrm{0}}\text{ is true}.\end{array}$$

If no difference in VAR models is detected, then the data are assumed to come from the same stochastic process.

Note that covariance matrices Γ₁ and Γ₂ characterize the residuals of the VAR models. If ${\mathbf{A}}_{\mathrm{1}}=\mathrm{\dots }={\mathbf{A}}_p=\mathrm{0}$, then the process is white noise, and the above test reduces to the standard likelihood ratio test for differences in covariance matrices (e.g., see Anderson, 1984). If the process is VAR(p), then Y_i−X_iB_i transforms the correlated process Y_t into an uncorrelated white noise process. Such a transformation is called a prewhitening transformation (Box et al., 2008). In practice, B_i is not known, and hence the sample estimate ${\widehat{\mathbf{B}}}_i$ is substituted, yielding an approximately white noise process. Accordingly, ${\widehat{\mathbf{\Gamma }}}_i$ is the covariance matrix of the prewhitened variables ${\mathbf{Y}}_i-{\mathbf{X}}_i{\widehat{\mathbf{B}}}_i$ and characterizes the noise statistics.

3 Partitioning the deviance statistic

If H₀ is rejected, then a difference in VAR models is said to be detected. However, this result does not tell us the nature of those differences. A natural question is whether the difference can be attributed to a difference in noise statistics or a difference in AR parameters. To answer these questions, we consider a hypothesis that is intermediate between H₀ and H₂, namely,

$$H_{\mathrm{1}}:{\mathbf{\Gamma }}_{\mathrm{1}}={\mathbf{\Gamma }}_{\mathrm{2}}\text{and no restriction on }{\mathbf{B}}_{\mathrm{1}},{\mathbf{B}}_{\mathrm{2}}.$$

As a shorthand, H₁ will be called equality of noise. Under H₁, the maximum likelihood estimates of B₁ and B₂ are Eqs. (6) and (7), and the (bias-corrected) maximum likelihood estimate of the covariance matrix Γ is

$$\begin{array}{}\text{(13)}& {\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}={\displaystyle \frac{{\mathit{\nu }}_{\mathrm{1}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}.\end{array}$$

It follows that the deviance statistic for testing H₀ versus H₁ is

$$\begin{array}{}\text{(14)}& D_{\mathrm{0}:\mathrm{1}}=({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}})\log \left|{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}\right|-({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}})\log \left|{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}\right|,\end{array}$$

and the deviance statistic for testing H₁ versus H₂ is

$$\begin{array}{}\text{(15)}& D_{\mathrm{1}:\mathrm{2}}=({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}})\log \left|{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}\right|-{\mathit{\nu }}_{\mathrm{1}}\log \left|{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}\right|-{\mathit{\nu }}_{\mathrm{2}}\log \left|{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}\right|.\end{array}$$

It is readily verified that the sum of these deviances gives the original deviance:

$$\begin{array}{}\text{(16)}& D_{\mathrm{0}:\mathrm{2}}=D_{\mathrm{0}:\mathrm{1}}+D_{\mathrm{1}:\mathrm{2}}.\end{array}$$

In the next section, we will see that D_1:2 measures differences in noise covariances, and D_0:1 measures differences in AR parameters. Under this interpretation, Eq. (16) decomposes total deviance into a part due to differences in noise statistics and a part due to differences in AR parameters.

Table 1 Table summarizing the hypotheses considered in the hierarchical test procedure.

Testing significance requires knowing the sampling distributions of D_0:1 and D_1:2. Under H₁, standard regression theory indicates that ${\mathit{\nu }}_{\mathrm{1}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}$ and ${\mathit{\nu }}_{\mathrm{2}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}$ each have independent Wishart distributions (Fujikoshi et al., 2010, Theorem 7.2.1). Importantly, these distributions are independent of the model parameters. Since D_1:2 depends only on ${\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}$ and ${\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}$, it too has a distribution that is independent of the model parameters (under H₁). For large sample sizes, this distribution approaches a chi-squared distribution with $J(J+\mathrm{1})/\mathrm{2}$ degrees of freedom (see Hogg et al., 2019, Theorem 6.5.1), which we denote as

$$\begin{array}{}\text{(17)}& D_{\mathrm{1}:\mathrm{2}}\sim {\mathit{\chi }}_{J(J+\mathrm{1})/\mathrm{2}}^{\mathrm{2}}\text{when }{H}_{\mathrm{1}}\text{ is true}.\end{array}$$

In contrast, the distribution of D_0:1 is more complicated. In the special case in which H₁ is true, the hypothesis of equal AR parameters can be tested as a standard problem in general linear hypothesis, which yields a statistic that depends only on ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}$ and ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$. Again, the sampling distribution is independent of the model parameters. Since D_0:1 depends only on ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}$ and ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$, we may infer that the distribution of D_0:1 is also independent of the model parameters. In this case, asymptotic theory indicates that, for large sample sizes, the distribution is

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& D_{\mathrm{0}:\mathrm{1}}\sim {\mathit{\chi }}_{MJ}^{\mathrm{2}}\text{when }{H}_{\mathrm{0}}\text{ is true, and }\\ & H_{\mathrm{1}}\text{ is the alternative hypothesis}.\end{array}\end{array}$$

As expected, for consistency, the sum of the degrees of freedom in Eqs. (17) and (18) equals that in Eq. (12).

However, if H₁ is not true (i.e., the noise covariances differ), then the distribution of D_0:1 depends on how the noise covariances differ. A similar situation arises with the t-test. Recall that the t-test assesses equality of means assuming equal variances. However, if the variances are unequal, the t-test is inappropriate because the distribution of the t-statistic depends on the ratio of variances. This is the central issue in the Behrens–Fisher problem, for which there are multiple tests with various probabilistic interpretations (Scheffé, 1970; Kim and Cohen, 1998). Our procedure would reduce to the t-test if the hypothesis μ₁=μ₂ were tested in the case of J=1 and vanishing AR parameters. Since the intercept μ is merely a special case of a regression matrix B, one can anticipate that testing equality of regression parameters would also lead to a sampling distribution that depends on the difference in noise covariances. Therefore, testing equality of AR parameters without assuming equality of noise covariances would lead to complications analogous to those associated with the Behrens–Fisher problem. To avoid these complications, we test equality of AR parameters under the assumption of equal noise covariances. This restriction is of the same kind as that invoked routinely in the t-test.

Given the above restriction, the hypotheses in our problem are summarized in Table 1 and nested as $H_{\mathrm{0}}\subseteq H_{\mathrm{1}}\subseteq H_{\mathrm{2}}$. A natural procedure is to first test equality of noise covariances and then, if that hypothesis is accepted, test equality of regression parameters under the restriction of equal noise covariances. This stepwise testing procedure is a special case of a general procedure discussed by Hogg (1961). Specifically, Hogg (1961) showed that if the maximum likelihood estimates of a set of nested hypotheses are completely sufficient statistics, then the resulting test statistics are stochastically independent, and the most restrictive hypothesis can be tested by testing a sequence of intermediate hypotheses. Hogg's procedure leads to the following stepwise procedure. First, D_1:2 is tested for significance based on Eq. (17). If it is significantly large, then we decide that the noise covariances differ and stop testing. If D_1:2 is not significant, then D_0:1 is tested for significance based on Eq. (18). If it is significant, then we decide that the noise covariances are equal, but the AR parameters differ. If D_0:1 is not significant, then we conclude that no significant difference in VAR models is detected. This procedure involves two steps performed at the α100 % level. Because the tests are stochastically independent, the family-wise error rate (FWER) for testing H₀ is

$$\text{FWER}=\mathrm{1}-(\mathrm{1}-\mathit{\alpha }{)}^{\mathrm{2}}.$$

Thus, if each step is performed at the α=5 % level, then the corresponding FWER is 9.75 %.

4 Diagnosing differences in noise processes

If D_1:2 is significantly large, then the noise covariances are said to differ. However, this conclusion does not tell us the nature of those differences. The question arises as to how to describe the difference between two covariance matrices in an efficient manner. Flury (1986) suggested using the linear combination of variables with extreme ratios. This approach is a generalization of principal component analysis (PCA), but instead of maximizing variance itself, a ratio of variances is optimized. We call this approach covariance discriminant analysis (CDA). A review of CDA is given in DelSole and Tippett (2022). CDA is a generalization of the discriminant analysis used in classification studies. We refer to both methods as discriminant analysis.

Discriminant analysis begins by considering the linear combination of prewhitened variables,

$$\begin{array}{}\text{(19)}& {\mathit{r}}^{\left(k\right)}=\left({\mathbf{Y}}_k-{\mathbf{X}}_k{\widehat{\mathbf{B}}}_k\right){\mathit{q}}_Y\text{for }k=\mathrm{1},\mathrm{2},\end{array}$$

where q_Y is a vector containing the coefficients of the linear combination. The resulting vector r^(k) is a time series called a variate. The corresponding sample variances are

$${\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}_k}}{\left({\mathit{r}}^{\left(k\right)}\right)}^T{\mathit{r}}^{\left(k\right)}={\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_k{\mathit{q}}_Y,$$

and the corresponding deviance statistic D_1:2 for the scalar variables r⁽¹⁾ and r⁽²⁾ is

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{d}_{\mathrm{1}:\mathrm{2}}& =\left({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}\right)\log \left({\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_Y\right)\\ & -{\mathit{\nu }}_{\mathrm{1}}\log \left({\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_Y\right)-{\mathit{\nu }}_{\mathrm{2}}\log \left({\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{\mathit{q}}_Y\right),\end{array}\end{array}$$

where Eq. (13) is used to define

$$\begin{array}{}\text{(21)}& {\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_Y=({\mathit{\nu }}_{\mathrm{1}}{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_Y+{\mathit{\nu }}_{\mathrm{2}}{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{\mathit{q}}_Y)/({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}).\end{array}$$

Substituting Eq. (21) into Eq. (20) and using elementary properties of the log function yields

$$d_{\mathrm{1}:\mathrm{2}}\left(\mathit{\theta }\right)={\mathit{\nu }}_{\mathrm{1}}\log \left({\displaystyle \frac{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}\mathit{\theta }}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}\right)+{\mathit{\nu }}_{\mathrm{2}}\log \left({\displaystyle \frac{{\mathit{\nu }}_{\mathrm{1}}/\mathit{\theta }+{\mathit{\nu }}_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}\right),$$

where θ is the ratio of noise variances given by the Rayleigh quotient

$$\begin{array}{}\text{(22)}& \mathit{\theta }={\displaystyle \frac{{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{\mathit{q}}_Y}{{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_Y}}.\end{array}$$

The stationary values of d_1:2 are found by using the chain rule:

$${\displaystyle \frac{\partial d_{\mathrm{1}:\mathrm{2}}}{\partial {\mathit{q}}_Y}}={\displaystyle \frac{\partial d_{\mathrm{1}:\mathrm{2}}}{\partial \mathit{\theta }}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial {\mathit{q}}_Y}}.$$

The function d_1:2(θ) is U-shaped with a minimum at θ=1, but the stationary value at θ=1 is not of interest since it is a minimum. For θ>1 the only stationary values occur at $\partial \mathit{\theta }/\partial {\mathit{q}}_Y=\mathrm{0}$, and for θ<1 the only stationary values occur at $\partial \mathit{\theta }/\partial {\mathit{q}}_Y=\mathrm{0}$. Both cases lead to the same condition $\partial \mathit{\theta }/\partial {\mathit{q}}_Y=\mathrm{0}$, which in turn leads to the generalized eigenvalue problem

$$\begin{array}{}\text{(23)}& {\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{\mathit{q}}_Y=\mathit{\theta }{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_Y\end{array}$$

(see Seber, 2008, result 6.59). This eigenvalue problem has at most J distinct solutions. Because d_1:2(θ) is a U-shaped function, large deviances arise from large and small values of θ. This is sensible because deviance is a measure of the difference in variance, and two variances may differ because one is much greater or much less than the other. We order the eigenvalues by their value of d_1:2(θ_i), which means that the eigenvalues ${\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\theta }}_J$ generally are not ordered from largest to smallest; for instance, the leading component might have the smallest eigenvalue. Let the corresponding eigenvectors be ${\mathit{q}}_{Y,\mathrm{1}},\mathrm{\dots },{\mathit{q}}_{Y,J}$. Then, the maximum possible deviance is obtained from q_Y,1, which yields time series Eq. (19) and has deviance D_1:2(θ₁).

When examining the leading eigenvalue θ₁, it is informative to compare to a significance threshold under H₁. This threshold can be estimated by Monte Carlo methods. Specifically, one generates samples from Eqs. (2) and (3) using ${\mathbf{\Gamma }}_{\mathrm{1}}={\mathbf{\Gamma }}_{\mathrm{2}}=\mathbf{I}$. The sampling distribution is independent of B₁ and B₂, and hence they can be set to zero in the Monte Carlo model. From these samples, one computes the noise covariance matrices, solves the generalized eigenvalue problem Eq. (23), and then repeats this procedure numerous times to build an empirical distribution of θ₁. The resulting test is equivalent to the union-intersection test (Mardia et al., 1979). We emphasize that this test is done only as a reference for judging the significance of a single eigenvalue and is not meant to replace the above test based on D_1:2, which depends on all the eigenvalues.

It can be shown that the variates ${\mathit{r}}_{\mathrm{1}}^{\left(k\right)},\mathrm{\dots },{\mathit{r}}_J^{\left(k\right)}$ corresponding to the eigenvectors are uncorrelated (see Appendix A). Recall that eigenvectors are unique up to a multiplicative constant. It is conventional to normalize the eigenvectors of Eq. (23) so that the denominator in Eq. (22) is one, or equivalently, that the variates ${\mathit{r}}_{\mathrm{1}}^{\left(\mathrm{1}\right)},\mathrm{\dots },{\mathit{r}}_J^{\left(\mathrm{1}\right)}$ each have unit variance. Under this normalization, the prewhitened variables may be decomposed as

$$\begin{array}{}\text{(24)}& {\mathbf{Y}}_k-{\mathbf{X}}_k{\widehat{\mathbf{B}}}_k={\mathit{r}}_{\mathrm{1}}^{\left(k\right)}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{\mathrm{1}}^T+\mathrm{\dots }+{\mathit{r}}_J^{\left(k\right)}{\stackrel{\mathrm{̃}}{\mathit{p}}}_J^T,\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathit{p}}}_j$ is a loading vector defined as

$${\stackrel{\mathrm{̃}}{\mathit{p}}}_j={\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_j.$$

We call the individual term ${\mathit{r}}_i^{\left(k\right)}{\stackrel{\mathrm{̃}}{\mathit{p}}}_i^T$ a noise discriminant. Equation (24) shows that prewhitened variables can be decomposed into a sum of independent noise discriminants, where ${\mathit{r}}_i^{\left(k\right)}$ is a time series describing the temporal structure of the noise discriminant, and ${\stackrel{\mathrm{̃}}{\mathit{p}}}_i$ is a loading vector describing the spatial structure of the noise discriminant. Appendix A also shows that

$$\begin{array}{}\text{(25)}& D_{\mathrm{1}:\mathrm{2}}=\sum _j{d}_{\mathrm{1}:\mathrm{2}}\left({\mathit{\theta }}_j\right).\end{array}$$

That is, D_1:2 is the sum of deviances of the individual noise discriminants. Thus, the decomposition simultaneously decomposes the prewhitened variables and the deviance D_1:2 into a sum of independent components. If a few noise discriminants dominate the total deviance (either because the associated θ_j is much greater than or much less than one), then those discriminants provide a low-dimensional description of differences in noise variances.

5 Diagnosing differences in dynamics

If D_1:2 is sufficiently small, then no difference in noise processes is detected and we turn to D_0:1. If D_0:1 is significantly large, then we conclude there is a difference in dynamics, but the nature of those differences is unknown. Again, we seek linear combinations that maximize D_0:1, consistent with discriminant analysis. Under H₀, the prewhitening transformation depends on ${\widehat{\mathbf{B}}}_{\mathrm{0}}$ rather than ${\widehat{\mathbf{B}}}_{\mathrm{1}}$ and ${\widehat{\mathbf{B}}}_{\mathrm{2}}$, and therefore we introduce the prewhitened variables under H₀,

$${\widehat{\mathbf{E}}}_{\mathrm{0}}=\left(\begin{array}{c}{\mathbf{Y}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{0}}\\ {\mathbf{Y}}_{\mathrm{2}}-{\mathbf{X}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{0}}\end{array}\right).$$

We continue to use $\mathit{q},\mathit{p},\mathit{r}$ to denote a projection vector, loading vector, and variate, respectively, but the interpretation differs from that in Sect. 4 because here we are maximizing D_0:1 rather than D_1:2. A linear combination of ${\widehat{\mathbf{E}}}_{\mathrm{0}}$ is

$${\mathit{r}}^{\left(\mathrm{0}\right)}={\widehat{\mathbf{E}}}_{\mathrm{0}}{\mathit{q}}_Y,$$

which has variance

$${\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}{\left({\mathit{r}}^{\left(\mathrm{0}\right)}\right)}^T{\mathit{r}}^{\left(\mathrm{0}\right)}={\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}{\mathit{q}}_Y.$$

The corresponding deviance statistic from Eq. (14) is then

$$\begin{array}{}\text{(26)}& d_{\mathrm{0}:\mathrm{1}}\left(\mathit{\lambda }\right)=({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}})\log \mathit{\lambda },\end{array}$$

where

$$\begin{array}{}\text{(27)}& \mathit{\lambda }={\displaystyle \frac{{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}{\mathit{q}}_Y}{{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_Y}}.\end{array}$$

Because d_0:1(λ) is a monotone function of λ, maximizing d_0:1(λ) is equivalent to maximizing λ, which leads to the eigenvalue problem

$$\begin{array}{}\text{(28)}& {\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}{\mathit{q}}_Y=\mathit{\lambda }{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_Y.\end{array}$$

Let the eigenvalues be ordered from largest to smallest, i.e., ${\mathit{\lambda }}_{\mathrm{1}}>\mathrm{\dots }>{\mathit{\lambda }}_J$, and let the corresponding eigenvectors be denoted as ${\mathit{q}}_{Y,\mathrm{1}},\mathrm{\dots },{\mathit{q}}_{Y,J}$. The rest of the procedure is similar to that in Sect. 4, and for brevity we simply state the final results. Specifically, the variates ${\mathit{r}}_{\mathrm{1}}^{\left(\mathrm{0}\right)},\mathrm{\dots },{\mathit{r}}_J^{\left(\mathrm{0}\right)}$ are uncorrelated, and the loading vectors are

$${\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,j}={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_{Y,j},$$

and if the eigenvectors are normalized as ${\mathit{q}}_{Y,j}^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_{Y,j}=\mathrm{1}$, then the prewhitened variables are decomposed as

$$\begin{array}{}\text{(29)}& {\widehat{\mathbf{E}}}_{\mathrm{0}}={\mathit{r}}_{\mathrm{1}}^{\left(\mathrm{0}\right)}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,\mathrm{1}}^T+\mathrm{\dots }+{\mathit{r}}_J^{\left(\mathrm{0}\right)}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,J}^T\end{array}$$

and the deviance is decomposed as

$$\begin{array}{}\text{(30)}& D_{\mathrm{0}:\mathrm{1}}=\sum _j{d}_{\mathrm{0}:\mathrm{1}}\left({\mathit{\lambda }}_j\right).\end{array}$$

We call ${\mathit{r}}_i^{\left(\mathrm{0}\right)}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,i}^T$ a difference-in-dynamics component. Identity Eq. (30) shows that D_0:1 can be decomposed as a sum of deviances of the individual components. The significance threshold of λ₁ under H₀ can be estimated by the same Monte Carlo method discussed in Sect.  4, except that Eq. (28) is solved instead of Eq. (23).

Although we have decomposed ${\widehat{\mathbf{E}}}_{\mathrm{0}}$ and D_0:1, it is not obvious how this elucidates differences in dynamics. After all, dynamics are contained in ${\widehat{\mathbf{B}}}_{\mathrm{1}}$ and ${\widehat{\mathbf{B}}}_{\mathrm{2}}$, but these quantities do not appear explicitly in the above decomposition. In Appendix C, we show that the loading vectors ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,\mathrm{1}}^T,\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,J}^T$ appear in the generalized singular value decomposition (GSVD) of B₁−B₂,

$$\begin{array}{}\text{(31)}& {\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}=s_{\mathrm{1}}{\mathit{q}}_{X,\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,\mathrm{1}}^T+\mathrm{\dots }+s_J{\mathit{q}}_{X,J}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,J}^T,\end{array}$$

where $s_{\mathrm{1}},\mathrm{\dots },s_J$ are singular values, ${\mathit{q}}_{X,\mathrm{1}},\mathrm{\dots },{\mathit{q}}_{X,J}$ are left singular vectors, and ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,j},\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,J}$ are right singular vectors. In contrast to the standard SVD, the left and right singular vectors satisfy the generalized orthogonality conditions

$$\begin{array}{}\text{(32)}& {\mathit{q}}_{X,i}^T{\mathbf{\Sigma }}_H{\mathit{q}}_{X,j}={\mathit{\delta }}_{ij}\text{and}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,i}^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{-\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,j}={\mathit{\delta }}_{ij}.\end{array}$$

This decomposition is equivalent to an SVD of a suitably transformed B₁−B₂ (Appendix C; see also Abdi and Williams, 2010). The matrix Σ_H is a covariance matrix defined in Appendix B (see B3). Note that if Σ_H=I and ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}=\mathbf{I}$, then Eq. (32) would be Euclidean norms and Eq. (31) would be the standard SVD. We call s_i a difference-in-dynamics singular value. Appendix C further shows that D_0:1 can be decomposed as

$$\begin{array}{}\text{(33)}& D_{\mathrm{0}:\mathrm{1}}=\sum _i{d}_{\mathrm{0}:\mathrm{1}}(\mathrm{1}+s_i^{\mathrm{2}}).\end{array}$$

s_i can be interpreted as measuring the magnitude of the difference in dynamics, and ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,i}$ describes the manifestation of these differences in the prewhitened variables. Furthermore, Eq. (33) shows that the deviance D_0:1 can be decomposed as a sum of the deviances of the individual components. Although a standard SVD would also decompose ${\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}$, it would not simultaneously decompose D_0:1 into a sum of terms, each of which depends only on a single component. If the sum in Eq. (33) is dominated by the first few components, then differences in AR parameters can be described by those components. We call this decomposition a difference-in-dynamics SVD.

To further clarify the decomposition Eq. (31), note that the form of Eqs. (2) and (3) implies that q_X,j projects X onto the jth difference-in-dynamics component. The corresponding loading vector is

$${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}={\mathbf{\Sigma }}_H{\mathit{q}}_{X,j}.$$

It follows from this and the first orthogonality constraint in Eq. (32) that the loading vectors ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}$ and projection vectors q_X,j form a bi-orthogonal set satisfying ${\mathit{q}}_{X,j}^T{\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,k}={\mathit{\delta }}_{jk}$. Therefore, multiplying the transpose of Eq. (31) by ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}$ yields

$$\begin{array}{}\text{(34)}& ({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}{)}^T{\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}=s_j{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,j}.\end{array}$$

This equation indicates that ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,i}$ defines an initial condition, while the vector ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,i}$ describes the corresponding difference in response, and s_i measures the amplitude of the difference in response. In Appendix C, we show that

$$\begin{array}{}\text{(35)}& s_j^{\mathrm{2}}={\displaystyle \frac{\Vert ({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}{)}^T{\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}{\Vert }_{{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}}^{\mathrm{2}}}{\Vert {\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}{\Vert }_{{\mathrm{\Sigma }}_H}^{\mathrm{2}}}},\end{array}$$

where $\Vert \cdot {\Vert }_{\mathbf{A}}$ is a distance measure defined as

$$\Vert \mathit{a}{\Vert }_{\mathbf{A}}^{\mathrm{2}}={\mathit{a}}^T{\mathbf{A}}^{-\mathrm{1}}\mathit{a},$$

for any vector a and positive definite matrix A. Thus, given the initial condition ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,j}$, the squared singular value $s_j^{\mathrm{2}}$ is the ratio of the dissimilarity in response over the initial amplitude. The dissimilarity in response is measured with respect to a norm based on the average noise covariance matrix ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$. Thus, the numerator in Eq. (35) is a signal-to-noise ratio, which is an attractive basis for measuring differences in responses. Because $s_j^{\mathrm{2}}$ satisfies the same optimization properties as λ_j, it follows that ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,\mathrm{1}}$ maximizes Eq. (35), or equivalently, it maximizes the dissimilarity in response subject to $\Vert {\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,\mathrm{1}}{\Vert }_{{\mathrm{\Sigma }}_H}^{\mathrm{2}}=\mathrm{1}$ (Noble and Daniel, 1988, Theorem 10.22). The latter condition defines a surface of constant probability density under 𝒩(0,Σ_H). Thus, ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{X,\mathrm{1}}$ produces the biggest difference in response among all “equally likely” vectors under 𝒩(0,Σ_H).

SVD has been used in past studies of empirical dynamical models (e.g., Penland and Sardeshmukh, 1995; Winkler et al., 2001; Shin and Newman, 2021), so it is of interest to discuss how the GSVD differs from past applications. For p=1, the VAR model with time step τ is mathematically equivalent to a linear inverse model (LIM) integrated to time τ. In LIM studies, the SVD is used to identify the “optimal initial condition” and the “maximum amplification curve”, which are derived by applying a standard SVD directly to the propagator $\widehat{\mathbf{B}}$. Thus, one obvious difference is that the SVD in this paper diagnoses differences in dynamics rather than dynamics itself. A more fundamental difference is that the standard SVD uses the Euclidean norm to constrain the initial and final vectors, which yields singular values that depend on the coordinate system. In contrast, difference-in-dynamics singular values are invariant to invertible linear transformations of Y and therefore are independent of the coordinate system (this is proven in Appendix D).

6 Application to North Atlantic sea surface temperatures

We now apply the above procedure to diagnose dissimilarities in internal variability of annual mean North Atlantic SST (NASST; $\mathrm{0}-{\mathrm{60}}^{\circ }$N) between observations and climate models. Our choices for climate models, observations, and VAR model were discussed in Part II. We first summarize these relevant results from Part II.

6.1 Summary of NASST analysis in Part II

In Part II, we applied the difference-in-VAR test to pre-industrial control runs from phase 5 of the Coupled Model Intercomparison Project (CMIP5 Taylor et al., 2012). These simulations contain only internal variability and do not contain forced variability on interannual timescales. We analyzed 36 control simulations of length 165 years.

For dimension reduction, NASST fields were projected onto the leading Laplacian eigenvectors in the Atlantic basin. The number of Laplacian eigenvectors J and the order of the VAR model p were chosen using an objective criterion called the Mutual Information Criterion (MIC; DelSole and Tippett, 2021b). This criterion is a generalization of Akaike's information criterion to random predictors and selection of predictands. As discussed in Part II, this criterion selected different values of p and J for different CMIP5 models. To quantify differences using our difference-in-VAR test, a common VAR model must be selected. Which model should be selected? Choosing p too small may result in residuals that contain serial correlations, contradicting our distributional assumptions. Choosing p too large may result in overfitting. Of the two problems, overfitting is less serious because it is taken into account in the statistical test. The main negative impact of overfitting is the loss of statistical power. In this study, loss of power is not a concern because differences grow rapidly with p and J. For this reason, we advocate selecting the maximum p and J among the models selected by MIC. As documented in Part II, these considerations lead to selecting J=7 Laplacian eigenvectors and order p=1.

To assess model adequacy, we performed the multivariate Ljung–Box test (Lütkepohl, 2005, Sect. 4.4.3). After a Bonferroni correction, only 1 of the 36 CMIP5 models had detectable serial correlations, consistent with the expected rate, suggesting that VAR(1) with J=7 was adequate for our data. Our choice of p=1 leads to a model that is indistinguishable from an LIM sampled at discrete time steps. Numerous studies have shown that LIMs provide reasonable models of monthly-mean and annual-mean SSTs (Alexander et al., 2008; Vimont, 2012; Zanna, 2012; Newman, 2013; Huddart et al., 2016; Dias et al., 2018) and can predict tropical SSTs with a skill close to that of state-of-the-art coupled atmosphere–ocean models (Newman and Sardeshmukh, 2017).

For observations, we used version 5 of the extended reconstructed SST data set (ERSSTv5; Huang et al., 2017) and considered the 165-year period 1854–2018. To remove forced variability, we regressed out either a second- or ninth-order polynomial function of time. For the sake of simplicity, we focus only on residuals from a quadratic. The CMIP5 and observation time series were split in half (82 and 83 years), and time series from the earlier half were compared to those in the latter half. The motivation behind this analysis was that if a time series comes from the same source, then the null hypothesis is true and the method should detect a difference at a rate consistent with the significance threshold, which was borne out by the results. Virtually all CMIP5 models were found to be inconsistent with observations. However, the nature of these differences was not diagnosed. The goal of this section is to diagnose these differences.

Figure 1 The deviance statistics D_1:2 (blue) and D_0:1 (red) for comparing ERSSTv5 1854–1935 to 82-year segments from 36 CMIP5 pre-industrial control simulations and to ERSSTv5 1936–2018. The 5 % critical values for D_1:2 and D_0:1 are indicated by the horizontal blue and red lines, respectively. The dots indicate H₁ is not rejected; i.e., equality of noise covariance matrices is not rejected.

6.2 Testing equality of noise covariances

First, we decompose the deviance as $D_{\mathrm{0}:\mathrm{2}}=D_{\mathrm{0}:\mathrm{1}}+D_{\mathrm{1}:\mathrm{2}}$. The individual terms D_0:1 and D_1:2 are shown in Fig. 1. We test each term at the 5 % significance level. The corresponding critical values are indicated as horizontal lines in Fig. 1. In most cases, differences in both noise covariances and AR parameters contribute substantially to the total deviance. We reject H₁ (i.e., equality of noise covariances) when D_1:2 exceeds its critical value (i.e., when blue bars cross the blue line). The figure shows that D_1:2 is significant for all but one model (GFDL-CM3), indicating that for most models the noise covariances differ between observations and CMIP5 models.

Figure 2 Range of noise variance ratios for maximizing the noise deviance D_1:2 between ERSSTv5 1854–1935 and 82-year segments from 36 CMIP5 pre-industrial control simulations. The horizontal dotted lines show the 2.5 % and 97.5 % significance thresholds under hypothesis H₁. The dots identify the noise variance ratio that dominates the deviance (deviance measures the difference from one but does not distinguish between values greater than or less than one).

Figure 3 Percent of deviance D_1:2 explained by the first noise discriminant, for the variance ratios shown in Fig.  2.

Figure 4 Regression patterns between local SST and the discriminant associated with the largest model-to-observation variance ratio. Only regression patterns associated with significant maximum variance ratios are shown. Each regression pattern is scaled to have a maximum absolute coefficient of one and is then plotted using the same color scale (red and blue are opposite signs). The CMIP5 model and percent of deviance D_1:2 explained by the discriminant are indicated in the title of each panel.

Figure 5 Same as Fig. 4 except for discriminants associated with the smallest model-to-observation variance ratio (i.e., CMIP5 models whose stochastic noises are smaller than those in observations).

6.3 Diagnosing differences in noise covariances

To diagnose differences in noise covariances, we perform the discriminant analysis described in Sect. 4 to find the components that maximize D_1:2. Specifically, observations in the first period are compared to models and observations in the second period. The corresponding ratios of noise variances θ are shown in Fig. 2. The ratio is defined as model over observations, so a ratio exceeding one indicates that the discriminant has a larger variance in the model than in observations. The horizontal dotted lines show the 2.5 % and 97.5 % significance thresholds for extreme ratios. The results show that most CMIP5 models have too little variance in their leading discriminant (where “leading” means the maximum value of θ₁ or $\mathrm{1}/{\mathit{\theta }}_{\mathrm{1}}$).

It is worth mentioning that a significant extreme ratio is not always expected when a significant deviance is detected. Discriminant analysis is designed to compress as much deviance as possible into the fewest number of components. The extent of compression depends on the data. If little compression is possible (e.g., the deviances are spread uniformly in phase space), then the most extreme ratio may not be significant, but nevertheless the overall deviance may be significant because individual deviances accumulate. The fact that in some CMIP5 models the most extreme noise variance ratio is insignificant, despite the noise deviance being significant, simply means that the deviance cannot be compressed significantly into a single component.

The percent of deviance D_1:2 explained by the leading discriminant is shown in Fig. 3. The leading discriminant explains more than 40 % of the deviance in about a half of the models. In these cases, the leading discriminant can provide an efficient basis for describing differences in noise covariances. Recall that these discriminants pertain to prewhitened variables. The time-lagged autocorrelation function of discriminant time series ${\mathit{r}}_j^{\left(k\right)}$ reveals little serial correlation (not shown).

To diagnose the spatial structure, we derive regression coefficients in which the observed local SST is the independent variable and ${\mathit{r}}_j^{\left(k\right)}$ in Eq. (19) is the dependent variable. Since the discriminants differ from model to model, differences in amplitude have complex interpretations, and hence for simplicity each regression pattern is scaled to have a maximum absolute coefficient of one. Figure 4 shows regression patterns for the case in which prewhitened variables in CMIP5 models have too much variance compared to ERSSTv5. As can be seen, the spatial structure of these discriminants is similar across CMIP5 models, suggesting a common model deficiency. The regression patterns associated with discriminants with too little (prewhitened) variance are shown in Fig.  5. Here, there is no obvious common structure across different models, although some patterns tend to be similar for CMIP5 models from the same institution.

Figure 6 The leading difference-in-dynamics singular value for each CMIP5 model. The 5 % significance threshold under H₁ is shown as the dotted horizontal line. Significant and insignificant singular values under H₁ are indicated by filled and open circles, respectively.

6.4 Diagnosing differences in dynamics

We now turn to diagnosing differences in AR parameters. Equality of AR parameters is rejected when D_0:1 exceeds its critical value (i.e., when the red bars cross the red line), provided H₁ is not rejected. However, H₁ is rejected for virtually all models, so for these models, strictly speaking, the testing procedure halts before testing equality of AR parameters. Having noted this fact, we nevertheless proceed with discriminant analysis of D_0:1 for the purpose of illustrating the procedure.

Figure 7 The initial condition (a, c, e) that maximizes the difference in 1-year response of VAR models estimated from HadGEM2-ES (a, b) and ERSSTv5 (c, d). The resulting response of the initial condition is shown in (b), (d), and (f). Panel (e) shows panels (a) minus (c); panel (f) shows panels (b) minus (d). The color scale is the same across panels but is not shown because the amplitude is arbitrary for linear initial value problems.

Figure 8 Same as Fig. 7 except for MIROC-ESM.

The leading difference-in-dynamics singular value s₁ for each model is shown in Fig. 6. The model with the largest singular value is HadGEM2-ES. The results of difference-in-dynamics SVD are shown in Fig. 7. The left-hand-side panels show the initial condition p_X,1 that maximizes the difference in responses (which is the same for the CMIP5 model and observations), and the right-hand-side panels show the response p_Y,1 of the respective VAR model to these initial conditions (along with the associated difference). We see that HadGEM2-ES gives a reasonable 1-year evolution over the subtropical–tropical areas but is much too persistent in the northernmost part of the domain compared to observations. To illustrate the diagnostic further, we show one more example. The results of difference-in-dynamics SVD for MIROC-ESM is shown in Fig. 8. As can be seen, the responses are dramatically different. Specifically, the VAR model from ERSSTv5 is predominantly a damping, whereas the VAR model from MIROC-ESM produces anomalies of the opposite sign in the North Atlantic.

7 Conclusions

This paper proposed techniques for optimally diagnosing differences between serially correlated, multivariate time series. The starting point is to assume each time series comes from a vector autoregressive model and then test the hypothesis that the model parameters are equal. The likelihood ratio test for this hypothesis was derived in Part II of this paper series, which leads to a deviance statistic that measures the difference in VAR models. However, merely deciding that two processes differ is not completely satisfying because it does not tell us the nature of those differences. For instance, a VAR can be viewed as a dynamical system driven by random noise, but do differences arise because of different dynamics, different noise statistics, or both? We show that the deviance between VAR models can be decomposed into two terms, one measuring differences in noise statistics and another measuring differences in dynamics (i.e., differences in AR parameters). Under suitable assumptions, the two terms are independent and have asymptotic chi-squared distributions, which leads to straightforward hypothesis testing. Then, the linear combination of variables that maximizes each deviance is obtained. One set of combinations identifies components whose noise covariances differ most strongly between the two processes. The other set of combinations effectively finds the initial conditions that maximize the difference in response from each VAR model.

The proposed method was applied to diagnose differences between annual-mean North Atlantic sea surface temperature variability between climate models and observations. The results show that differences in dynamics and differences in noise covariances both play a significant role in the overall differences. In fact, only one model was found to be consistent with observations. However, this consistency is not robust. In further analyses (not shown), we changed the half of observations used to diagnose differences and found that this model no longer was consistent with observations. In fact, no model was found to be consistent with both halves of the observational period. For half the models, the leading discriminant of noise variance explains 40 % or more of the deviance in noise covariances. For most models, the leading discriminant underpredicts the noise variance, and the associated spatial structures are highly model dependent. In contrast, models that have leading discriminants with too much noise variance tend to have similar spatial structures.

Discriminant analysis was also applied to identify differences in dynamics, or more precisely, differences in AR parameters. These discriminants can be interpreted as finding suitably normalized initial conditions that maximize the difference in response from the two VAR models. Examples from this analysis revealed dramatic errors in the 1-year predictions of some models. For instance, the VAR model trained on HadGEM2-ES gives realistic 1-year predictions at the mid to low latitudes of the North Atlantic but are much too persistent at the high latitudes. The VAR model trained on MIROC-ESM produces 1-year predictions with the wrong sign in the North Atlantic.

Although previous climate studies have diagnosed differences in variability using discriminant techniques (Schneider and Griffies, 1999; DelSole et al., 2011; Wills et al., 2018), the above procedure differs from previous analysis in important ways. Specifically, here discriminant techniques are applied to a measure of dissimilarity that has a simple significance test. This means that the dissimilarity measure can be tested for significance before embarking on a diagnosis of the dissimilarities. The simplicity of the test is a consequence of measuring dissimilarities of prewhitened variables rather than dissimilarities in the original, serially correlated variables. In general, prewhitened variables are closer to white noise than the original data, and hence their asymptotic distributions are simpler than those of the original process.

Diagnosing differences in dynamics leads to generalized SVD, which differs from the SVD used in previous studies to diagnose dynamics of linear inverse models (e.g., Penland and Sardeshmukh, 1995; Winkler et al., 2001; Shin and Newman, 2021). In particular, standard SVD uses the Euclidean norm to constrain the initial and final vectors, which yields singular vectors that depend on the coordinate system. This dependence might be averted by asserting a norm to be some physical quantity like energy or enstrophy, but these assertions usually have no compelling justification and lead to statistics with no simple significance test. In contrast, difference-in-dynamics SVD emerges from a rigorous statistical framework, and, as a result, it leads to simple significance tests, and the norms chosen for maximization are dictated by the statistical framework and lead to results that are independent of the arbitrary coordinate system.

One reason for diagnosing differences in stochastic processes is to gain insight into the causes of the differences. For instance, dissimilarities in noise processes indicate differences in the one-step prediction errors of the associated VAR models. Comparison of prediction errors may give clues as to the locations or processes that dominate the errors, which may be helpful for improving a climate model. Similarly, differences in dynamics indicate differences in the memory or in the evolution of anomalies between observations and climate models. The diagnostics derived here may give clues as to the physical source of the differences, which in turn may give insight into how to improve climate models.

Appendix A: Decomposition of D_1:2

This Appendix proves the key results of decomposing D_1:2, namely, Eqs. (24) and (25). We begin with the eigenvalue problem Eq. (23). Let us collect the eigenvectors ${\mathit{q}}_{Y,\mathrm{1}},\mathrm{\dots },{\mathit{q}}_{Y,J}$ into the matrix

$$\begin{array}{}\text{(A1)}& {\mathit{Q}}_Y=\left(\begin{array}{cccc}{\mathit{q}}_{Y,\mathrm{1}}& {\mathit{q}}_{Y,\mathrm{2}}& \mathrm{\dots }& {\mathit{q}}_{Y,J}\end{array}\right).\end{array}$$

Recall that eigenvectors are unique up to a multiplicative factor, and hence we may normalize them such that ${\mathit{q}}_{Y,i}^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_{Y,i}=\mathrm{1}$. Because ${\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}$ and ${\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}$ are symmetric, the eigenvectors satisfy the following properties (see Seber, 2008, result 6.67):

$$\begin{array}{}\text{(A2)}& {\mathit{Q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{Q}}_Y=\mathit{I}\text{and}{\mathit{Q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{\mathit{Q}}_Y=\mathbf{\Theta },\end{array}$$

where Θ is a diagonal matrix with diagonal elements ${\mathit{\theta }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\theta }}_J$. For each eigenvector q_Y,i, a time series ${\mathit{r}}_i^{\left(k\right)}$ can be computed from Eq. (19). The variates for all the eigenvectors may be collected into a single matrix as

$${\mathit{R}}^{\left(k\right)}=\left(\begin{array}{cccc}{\mathit{r}}_{\mathrm{1}}^{\left(k\right)}& {\mathit{r}}_{\mathrm{2}}^{\left(k\right)}& \mathrm{\dots }& {\mathit{r}}_J^{\left(k\right)}\end{array}\right).$$

Then, Eq. (19) implies

$$\begin{array}{}\text{(A3)}& {\mathit{R}}^{\left(k\right)}=\left({\mathbf{Y}}_k-{\mathbf{X}}_k{\widehat{\mathbf{B}}}_k\right){\mathit{Q}}_Y,\end{array}$$

and this combined with Eqs. (9), (10), and (A2) implies that the associated covariance matrices are

$${\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}_{\mathrm{1}}}}{\left({\mathit{R}}^{\left(\mathrm{1}\right)}\right)}^T{\mathit{R}}^{\left(\mathrm{1}\right)}=\mathit{I}\text{and}{\displaystyle \frac{\mathrm{1}}{{\mathit{\nu }}_{\mathrm{2}}}}{\left({\mathit{R}}^{\left(\mathrm{2}\right)}\right)}^T{\mathit{R}}^{\left(\mathrm{2}\right)}=\mathbf{\Theta }.$$

Diagonal covariance matrices imply that the variates are uncorrelated and therefore independent under normality.

Because the eigenvectors are distinct and form a complete basis set, Q_Y is non-singular, and we may define

$$\begin{array}{}\text{(A4)}& {\stackrel{\mathrm{̃}}{\mathit{P}}}_Y={\mathit{Q}}_Y^{-\mathrm{1}T}={\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{Q}}_Y,\end{array}$$

where the first equality is a definition and the second equality follows from Eq. (A2). Multiplying Eq. (A3) on the right by ${\mathit{Q}}_Y^{-\mathrm{1}}$ and using Eq. (A4) yields

$$\begin{array}{}\text{(A5)}& {\mathit{R}}^{\left(k\right)}{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T=\left({\mathbf{Y}}_k-{\mathbf{X}}_k{\widehat{\mathbf{B}}}_k\right).\end{array}$$

If the column vectors of ${\stackrel{\mathrm{̃}}{\mathit{P}}}_Y$ are denoted as ${\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,\mathrm{1}},\mathrm{\dots },{\stackrel{\mathrm{̃}}{\mathit{p}}}_{Y,J}$, then Eq. (A5) can be written equivalently as Eq. (24).

Importantly, the above decomposition also decomposes the deviance. To see this, note that Eqs. (A2) and (A4) imply

$$\begin{array}{c}{\displaystyle }{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T,{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y\mathbf{\Theta }{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T,\\ {\displaystyle }\text{and}{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y\left({\displaystyle \frac{{\mathit{\nu }}_{\mathrm{1}}\mathit{I}+{\mathit{\nu }}_{\mathrm{2}}\mathbf{\Theta }}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}\right){\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T.\end{array}$$

Substituting these expressions into Eq. (15) and using standard properties of the determinant yields

$$\begin{array}{c}{\displaystyle }{D}_{\mathrm{1}:\mathrm{2}}={\mathit{\nu }}_{\mathrm{1}}\log \left({\displaystyle \frac{\left|{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}\right|}{\left|{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}\right|}}\right)+{\mathit{\nu }}_{\mathrm{2}}\log \left({\displaystyle \frac{\left|{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}\right|}{\left|{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}\right|}}\right)\\ {\displaystyle }={\mathit{\nu }}_{\mathrm{1}}\log {\displaystyle \frac{\left|\frac{{\mathit{\nu }}_{\mathrm{1}}\mathit{I}+{\mathit{\nu }}_{\mathrm{2}}\mathbf{\Theta }}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}\right|}{\left|\mathit{I}\right|}}+{\mathit{\nu }}_{\mathrm{2}}\log {\displaystyle \frac{\left|\frac{{\mathit{\nu }}_{\mathrm{1}}\mathit{I}+{\mathit{\nu }}_{\mathrm{2}}\mathbf{\Theta }}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}\right|}{\left|\mathbf{\Theta }\right|}}=\sum _j{d}_{\mathrm{1}:\mathrm{2}}\left({\mathit{\theta }}_j\right),\end{array}$$

which is Eq. (25). This completes the proofs of Eqs. (24) and (25).

Appendix B: Decomposition of ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}$

In this Appendix, we prove that the matrix ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}$ defined in Eq. (11) can be decomposed as

$$\begin{array}{}\text{(B1)}& {\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}+\mathbf{\Delta },\end{array}$$

where

$$\begin{array}{}\text{(B2)}& \mathbf{\Delta }={\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\mathbf{\Sigma }}_H\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)\end{array}$$

and

$$\begin{array}{}\text{(B3)}& {\mathbf{\Sigma }}_H={\left({\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}\right)}^{-\mathrm{1}}+{\left({\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\right)}^{-\mathrm{1}}/\left({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}\right).\end{array}$$

Identity Eq. (B1) clarifies that D_0:1 defined in Eq. (14) measures differences in AR parameters, because if ${\widehat{\mathbf{B}}}_{\mathrm{2}}={\widehat{\mathbf{B}}}_{\mathrm{1}}$, then D_0:1 vanishes (because then Δ=0 and therefore ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$); otherwise, D_0:1 is positive (because ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$ is positive definite). This identity will be needed in Appendix C.

Starting from the definition Eqs. (6)–(11), we note that the MSEs satisfy the standard orthogonality relations

$${\mathbf{X}}_{\mathrm{1}}^T\left({\mathbf{Y}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)=\mathrm{0}\text{and}{\mathbf{X}}_{\mathrm{2}}^T\left({\mathbf{Y}}_{\mathrm{2}}-{\mathbf{X}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{2}}\right)=\mathrm{0}.$$

It follows from these orthogonality relations that

$$\begin{array}{c}{\displaystyle }\lfloor \lfloor {\mathbf{Y}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{0}}\rfloor \rfloor =\lfloor \lfloor {\mathbf{Y}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{1}}\rfloor \rfloor +\lfloor \lfloor {\mathbf{X}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{1}}-{\widehat{\mathbf{B}}}_{\mathrm{0}}\right)\rfloor \rfloor \\ {\displaystyle }={\mathit{\nu }}_{\mathrm{1}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}+\lfloor \lfloor {\mathbf{X}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{1}}-{\widehat{\mathbf{B}}}_{\mathrm{0}}\right)\rfloor \rfloor ,\\ {\displaystyle }\lfloor \lfloor {\mathbf{Y}}_{\mathrm{2}}-{\mathbf{X}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{0}}\rfloor \rfloor =\lfloor \lfloor {\mathbf{Y}}_{\mathrm{2}}-{\mathbf{X}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{2}}\rfloor \rfloor +\lfloor \lfloor {\mathbf{X}}_{\mathrm{2}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{0}}\right)\rfloor \rfloor \\ {\displaystyle }={\mathit{\nu }}_{\mathrm{2}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}+\lfloor \lfloor {\mathbf{X}}_{\mathrm{2}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{0}}\right)\rfloor \rfloor .\end{array}$$

Therefore, the MLE of Γ₀ in Eq. (11) can be written equivalently as

$$\begin{array}{}\text{(B4)}& {\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}={\displaystyle \frac{{\mathit{\nu }}_{\mathrm{1}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}+\mathbf{\Delta }={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}+\mathbf{\Delta },\end{array}$$

where

$$\mathbf{\Delta }={\displaystyle \frac{\lfloor \lfloor {\mathbf{X}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{1}}-{\widehat{\mathbf{B}}}_{\mathrm{0}}\right)\rfloor \rfloor +\lfloor \lfloor {\mathbf{X}}_{\mathrm{2}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{0}}\right)\rfloor \rfloor }{{\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}}}.$$

To show that this Δ equals Eq. (B2), note that

$$\begin{array}{c}{\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{0}}={\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}+{\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{Y}}_{\mathrm{1}}+{\mathbf{X}}_{\mathrm{2}}^T{\mathbf{Y}}_{\mathrm{2}}\right)\\ {\displaystyle }={\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}+{\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\left(\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}\right){\widehat{\mathbf{B}}}_{\mathrm{1}}+\left({\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right){\widehat{\mathbf{B}}}_{\mathrm{2}}\right).\end{array}$$

To simplify notation, we define

$${\mathbf{A}}_{\mathrm{1}}={\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}\text{and}{\mathbf{A}}_{\mathrm{2}}={\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}.$$

Then,

$$\begin{array}{c}{\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{0}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}={\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathbf{A}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{2}}\right)-{\widehat{\mathbf{B}}}_{\mathrm{1}}\\ {\displaystyle }={\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathbf{A}}_{\mathrm{1}}{\widehat{\mathbf{B}}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}{\widehat{\mathbf{B}}}_{\mathrm{2}}-\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right){\widehat{\mathbf{B}}}_{\mathrm{1}}\right)\\ {\displaystyle }={\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right)}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{2}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right).\end{array}$$

Similarly, we have

$${\widehat{\mathbf{B}}}_{\mathrm{0}}-{\widehat{\mathbf{B}}}_{\mathrm{2}}={\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right)}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{1}}-{\widehat{\mathbf{B}}}_{\mathrm{2}}\right).$$

Thus,

$$\begin{array}{rl}\lfloor \lfloor {\mathbf{X}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{0}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)\rfloor \rfloor & ={\left({\widehat{\mathbf{B}}}_{\mathrm{0}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\mathbf{A}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{0}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)\\ & ={\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\mathbf{A}}_{\mathrm{2}}{\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\\ & {\mathbf{A}}_{\mathrm{1}}{\left({\mathbf{A}}_{\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}\right)}^{-\mathrm{1}}{\mathbf{A}}_{\mathrm{2}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)\\ & ={\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\left({\mathbf{A}}_{\mathrm{1}}^{-\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\\ & {\mathbf{A}}_{\mathrm{1}}^{-\mathrm{1}}{\left({\mathbf{A}}_{\mathrm{1}}^{-\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right).\end{array}$$

Similarly, we have

$$\begin{array}{rl}\lfloor \lfloor {\mathbf{X}}_{\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{0}}-{\widehat{\mathbf{B}}}_{\mathrm{2}}\right)\rfloor \rfloor & ={\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\left({\mathbf{A}}_{\mathrm{1}}^{-\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\\ & {\mathbf{A}}_{\mathrm{2}}^{-\mathrm{1}}{\left({\mathbf{A}}_{\mathrm{1}}^{-\mathrm{1}}+{\mathbf{A}}_{\mathrm{2}}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right).\end{array}$$

Adding these together, dividing by ν₁+ν₂, and returning to the original notation gives

$$\begin{array}{rl}\mathbf{\Delta }& ={\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\left({\left({\mathbf{X}}_{\mathrm{1}}^T{\mathbf{X}}_{\mathrm{1}}\right)}^{-\mathrm{1}}+{\left({\mathbf{X}}_{\mathrm{2}}^T{\mathbf{X}}_{\mathrm{2}}\right)}^{-\mathrm{1}}\right)}^{-\mathrm{1}}\\ & \left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)/({\mathit{\nu }}_{\mathrm{1}}+{\mathit{\nu }}_{\mathrm{2}}).\end{array}$$

which is Eq. (B2), as desired.

Appendix C: Decomposition of D_0:1

This Appendix proves the key results of decomposing D_0:1, namely, Eqs. (31), (33), (34), and (35). Our starting point is the solution to the eigenvalue problem Eq. (28). Let the eigenvalues from Eq. (28) be ${\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_J$ and the associated eigenvectors be ${\mathit{q}}_{Y,\mathrm{1}},\mathrm{\dots },{\mathit{q}}_{Y,J}$. The eigenvectors may be collected into the matrix Q_Y, analogous to Eq. (A1). Furthermore, because the matrices in Eq. (28) are symmetric, the eigenvectors may be normalized to satisfy (see Seber, 2008, result 6.67)

$$\begin{array}{}\text{(C1)}& {\mathit{Q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{Q}}_Y=\mathit{I}\text{and}{\mathit{Q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}{\mathit{Q}}_Y=\mathbf{\Lambda },\end{array}$$

where Λ is a diagonal matrix with diagonal elements ${\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_J$. If the eigenvalues are distinct, then the eigenvectors form a complete basis set and Q_Y is non-singular, and hence we may define

$$\begin{array}{}\text{(C2)}& {\stackrel{\mathrm{̃}}{\mathit{P}}}_Y={\mathit{Q}}_Y^{-\mathrm{1}T}={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{Q}}_Y,\end{array}$$

where the second equality follows from the first identity in Eq. (C1). Using Eq. (C2), the identity Eq. (C1) can be written equivalently as

$$\begin{array}{}\text{(C3)}& {\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T\text{and}{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y\mathbf{\Lambda }{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T.\end{array}$$

Now, consider the Rayleigh quotient Eq. (27). The matrices ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}$ and ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$ are related to each other through Eq. (B1). Substituting Eq. (B1) into Eq. (27) yields

$$\begin{array}{}\text{(C4)}& \mathit{\lambda }=s^{\mathrm{2}}+\mathrm{1},\end{array}$$

where s² is the Rayleigh quotient

$$\begin{array}{}\text{(C5)}& s^{\mathrm{2}}={\displaystyle \frac{{\mathit{q}}_Y^T\mathbf{\Delta }{\mathit{q}}_Y}{{\mathit{q}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}{\mathit{q}}_Y}}.\end{array}$$

Let ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{\mathrm{1}/\mathrm{2}}$ denote the matrix square root of ${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}$. Then, Eq. (C5) can be written equivalently as

$$\begin{array}{}\text{(C6)}& s^{\mathrm{2}}={\displaystyle \frac{{\mathit{v}}^T{\mathbf{\Phi }}^T\mathbf{\Phi }\mathit{v}}{{\mathit{v}}^T\mathit{v}}},\end{array}$$

where

$$\begin{array}{}\text{(C7)}& \mathit{v}={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{\mathrm{1}/\mathrm{2}}{\mathit{q}}_Y\end{array}$$

and

$$\begin{array}{}\text{(C8)}& \mathbf{\Phi }={\mathbf{\Sigma }}_H^{\mathrm{1}/\mathrm{2}}\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right){\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{-\mathrm{1}/\mathrm{2}}.\end{array}$$

A standard result in linear algebra states that extrema of Eq. (C6) are given by the eigenvectors of Φ^TΦ or equivalently the singular vectors of Φ. Let the singular value decomposition (SVD) of Φ be

$$\begin{array}{}\text{(C9)}& \mathbf{\Phi }={\mathbf{USV}}^T,\end{array}$$

where U and V are orthogonal matrices and S is a diagonal matrix with non-negative diagonal elements called singular values. Because $\mathit{\lambda }=s^{\mathrm{2}}+\mathrm{1}$, maximizing s² in Eq. (C6) is equivalent to maximizing λ in Eq. (27). Therefore, the squared singular values $s_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots },s_J^{\mathrm{2}}$ must equal the eigenvalues ${\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_J$, and the eigenvectors ${\mathit{q}}_{Y,\mathrm{1}},\mathrm{\dots },{\mathit{q}}_{Y,J}$ must give the right singular vectors ${\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\mathit{v}}_J$ through Eq. (C7). Let $s_{\mathrm{1}},\mathrm{\dots },s_M$ be ordered from largest to smallest.

Comparing Eqs. (C8) and (C9), we can anticipate that U describes initial conditions, V describes differences in responses, and diagonals of S measure response amplitudes. To be more precise, define

$$\begin{array}{}\text{(C10)}& \begin{array}{cc}{\mathit{Q}}_X={\mathbf{\Sigma }}_H^{-\mathrm{1}/\mathrm{2}}\mathbf{U},& {\mathit{Q}}_Y={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{-\mathrm{1}/\mathrm{2}}\mathbf{V},\\ & \\ {\stackrel{\mathrm{̃}}{\mathit{P}}}_X={\mathbf{\Sigma }}_H^{\mathrm{1}/\mathrm{2}}\mathbf{U},& {\stackrel{\mathrm{̃}}{\mathit{P}}}_Y={\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{\mathrm{1}/\mathrm{2}}\mathbf{V}.\end{array}\end{array}$$

Note that Eqs. (C2) and (C3) hold for the above definitions, as they should since the eigenvalue problem Eq. (28) and SVD Eq. (C9) perform the same decomposition. Furthermore, Eq. (C10) implies the additional identities

$$\begin{array}{}\text{(C11)}& \begin{array}{rl}{\mathit{Q}}_X^T{\mathbf{\Sigma }}_H{\mathit{Q}}_X& =\mathit{I},\\ {\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{-\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y& =\mathit{I},\text{and}\\ {\stackrel{\mathrm{̃}}{\mathit{P}}}_X^T{\mathbf{\Sigma }}_H^{-\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{P}}}_X& =\mathit{I},\end{array}\end{array}$$

the first two of which give Eq. (32), and the last implies $\Vert {\stackrel{\mathrm{̃}}{\mathit{p}}}_X{\Vert }_{{\mathrm{\Sigma }}_H}^{\mathrm{2}}=\mathrm{1}$. However, Eq. (C9) can be rearranged into the form

$$\begin{array}{}\text{(C12)}& {\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}={\mathit{Q}}_X\mathbf{S}{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T,\end{array}$$

which, after expanding Q_X and ${\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T$ into column vectors, can be written equivalently as Eq. (31). Decomposition Eq. (C12) and the first two identities in Eq. (C11) are a weighted form of the generalized singular value decomposition (Abdi and Williams, 2010).

Just as the column vectors of V maximize the Rayleigh quotient Eq. (C6), the column vectors of U maximize

$$\begin{array}{}\text{(C13)}& s^{\mathrm{2}}={\displaystyle \frac{{\mathit{u}}^T\mathbf{\Phi }{\mathbf{\Phi }}^T\mathit{u}}{{\mathit{u}}^T\mathit{u}}}={\displaystyle \frac{\Vert {\widehat{\mathbf{B}}}_{\mathrm{2}}^T{\stackrel{\mathrm{̃}}{\mathit{p}}}_X-{\widehat{\mathbf{B}}}_{\mathrm{1}}^T{\stackrel{\mathrm{̃}}{\mathit{p}}}_X{\Vert }_{{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}}^{\mathrm{2}}}{\Vert {\stackrel{\mathrm{̃}}{\mathit{p}}}_X{\Vert }_{{\mathrm{\Sigma }}_H}^{\mathrm{2}}}},\end{array}$$

which is Eq. (35). In addition, the identity Φ^TU=VS^T implies

$${\left({\widehat{\mathbf{B}}}_{\mathrm{2}}-{\widehat{\mathbf{B}}}_{\mathrm{1}}\right)}^T{\stackrel{\mathrm{̃}}{\mathit{P}}}_X={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y{\mathbf{S}}^T,$$

which gives Eq. (34). An additional property of ${\stackrel{\mathrm{̃}}{\mathit{P}}}_Y$ can be derived from the identity Φ^TΦ=VS²V^T, which implies

$$\mathbf{\Delta }={\stackrel{\mathrm{̃}}{\mathit{P}}}_Y{\mathbf{S}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{P}}}_Y^T,$$

which in turn can be written equivalently as

$$\mathbf{\Delta }=s_{\mathrm{1}}^{\mathrm{2}}{\mathit{p}}_{Y,\mathrm{1}}{\mathit{p}}_{Y,\mathrm{1}}^T+\mathrm{\dots }+s_J^{\mathrm{2}}{\mathit{p}}_{Y,J}{\mathit{p}}_{Y,J}^T.$$

This identity shows that Δ can be decomposed into a sum of terms involving p_Y,i without any cross terms. Finally, Eqs. (30) and (C4) imply

$$\begin{array}{}\text{(C14)}& D_{\mathrm{0}:\mathrm{1}}=\sum _i{d}_{\mathrm{0}:\mathrm{1}}\left(\mathrm{1}+s_i^{\mathrm{2}}\right),\end{array}$$

which shows that the SVD of Φ partitions the deviance into a sum of components ordered by their contribution to D_0:1.

Appendix D: Invariance to linear transformation

In this section, we prove that the eigenvalues from Eqs. (23) and (28) are invariant to invertible linear transformations of the form

$${\mathbf{Y}}^{\prime }=\mathbf{YG},$$

where primes indicate transformed quantities and G is a square, non-singular matrix. Using straightforward matrix algebra, Eqs. (6)–(13) are transformed as

$$\begin{array}{}\text{(D1)}& {\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{1}}^{\prime }={\widehat{\mathbf{B}}}_{\mathrm{1}}\mathbf{G},\text{(D2)}& {\displaystyle }{\widehat{\mathbf{B}}}_{\mathrm{2}}^{\prime }={\widehat{\mathbf{B}}}_{\mathrm{2}}\mathbf{G},\text{(D3)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}^{\prime }={\mathbf{G}}^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}\mathbf{G},\text{(D4)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}^{\prime }={\mathbf{G}}^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}\mathbf{G},\text{(D5)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}^{\prime }={\mathbf{G}}^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}\mathbf{G},\text{(D6)}& {\displaystyle }{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{\prime }={\mathbf{G}}^T{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}\mathbf{G}.\end{array}$$

Substituting transformed quantities into the eigenvalue problem Eq. (23) yields

$$\begin{array}{rl}{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}^{\prime }{\mathit{q}}_Y^{\prime }& =\mathit{\theta }{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}^{\prime }{\mathit{q}}_Y^{\prime },\\ {\mathbf{G}}^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}\mathbf{G}{\mathit{q}}_Y^{\prime }& =\mathit{\theta }{\mathbf{G}}^T{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}\mathbf{G}{\mathit{q}}_Y^{\prime },\\ {\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}\mathbf{G}{\mathit{q}}_Y^{\prime }& =\mathit{\theta }{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}\mathbf{G}{\mathit{q}}_Y^{\prime },\\ {\widehat{\mathbf{\Gamma }}}_{\mathrm{1}}{\mathit{q}}_Y& =\mathit{\theta }{\widehat{\mathbf{\Gamma }}}_{\mathrm{2}}{\mathit{q}}_Y,\end{array}$$

where the second line follows by substituting Eqs. (D3) and (D4), the third line follows by multiplying both sides on the left by the inverse of G^T (which is possible because G is invertible), and the fourth line follows by defining

$$\begin{array}{}\text{(D7)}& {\mathit{q}}_Y=\mathbf{G}{\mathit{q}}_Y^{\prime }.\end{array}$$

This shows that linear transformations do not affect the eigenvalues of Eq. (23), while they transform the eigenvectors as Eq. (D7).

Similarly, substituting transformed quantities into the eigenvalue problem Eq. (28) yields

$${\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{0}}}^{\prime }{\mathit{q}}_Y^{\prime }=\mathit{\lambda }{\widehat{\mathbf{\Gamma }}}_{H_{\mathrm{1}}}^{\prime }{\mathit{q}}_Y^{\prime }.$$

By the same steps as above, this eigenvalue problem can be transformed to Eq. (28). This shows that the eigenvalues of Eq. (28) are invariant to invertible linear transformations. Because of Eq. (C4), it follows that the difference-in-dynamics singular values s_j are also invariant to invertible linear transformations.