An improved projection of climate observations for detection and attribution

An important goal of climate research is to determine the causal contribution of human activity to observed changes in the climate system. Methodologically speaking, most climatic causal studies to date have been formulating attribution as a linear regression inference problem. Under this formulation, the inference is often obtained by using the generalized least squares (GLS) estimator after projecting the data on the r leading eigenvectors of the covariance associated with internal variability, which are evaluated from numerical climate models. In this paper, we revisit the problem of obtaining a GLS estimator adapted to this particular situation, in which only the leading eigenvectors of the noise's covariance are assumed to be known. After noting that the eigenvectors associated with the lowest eigenvalues are in general more valuable for inference purposes, we introduce an alternative estimator. Our proposed estimator is shown to outperform the conventional estimator, when using a simulation test bed that represents the 20th century temperature evolution.

1 Introduction

1.1 Context

An important goal of climate research is to determine the causes of past global warming in general and the responsibility of human activity in particular (Hegerl et al., 2007); this question has thus emerged as a research topic known as detection and attribution (D&A). From a methodological standpoint, D&A studies are usually based on linear regression methods, often referred to in this particular context as optimal fingerprinting, whereby an observed climate change is regarded as a linear combination of several externally forced signals added to internal climate variability (Hasselmann, 1979; Bell, 1986; North et al., 1982; Allen and Tett, 1999; Hegerl and Zwiers, 2011). On the one hand, the latter signals consist of spatial, temporal or space–time patterns of response to external forcings as anticipated by one or several climate models. Internal climate variability, on the other hand, is usually represented by a centred multivariate Gaussian noise. Its covariance Σ, which is also estimated from model simulations, thus describes the detailed space–time features of internal variability. In particular, its eigenvectors can often be interpreted as the modes of variability, whose relative magnitudes are reflected by the associated eigenvalues.

Denoting y the n-dimensional vector of observed climate change, $\mathbf{x}=(x_{\mathrm{1}},\mathrm{\dots },x_p)$ the n×p matrix concatenating the p externally forced signals and ν the internal climate variability noise, the regression equation is as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{l}\mathit{y}=\mathbf{x}\mathit{\beta }+\mathit{\nu }\end{array}.\end{array}$$

The results of the inference on the vector of regression coefficients β, and the magnitude of its confidence intervals, determine whether the external signals “are present in the observations” (Hegerl and Zwiers, 2011) and whether or not the observed change is attributable to each forcing.

Under the assumption that Σ is known, the above inference problem can be solved using the standard generalized least squares (GLS) setting (e.g. Amemiya, 1985). Under this framework, the expression of the best linear unbiased estimator (BLUE) and its variance are given by the following:

$$\begin{array}{}\text{(2)}& \begin{array}{ll}\widehat{\mathit{\beta }}=\left({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}\mathbf{x}{)}^{-\mathrm{1}}\right({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}\mathit{y}),& \text{Var}\left(\widehat{\mathit{\beta }}\right)=({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}\mathbf{x}{)}^{-\mathrm{1}}.\end{array}\end{array}$$

The GLS estimator of Eq. (2) can be obtained as the ordinary least squares (OLS) estimator for the linearly transformed regression equation:

$$\begin{array}{}\text{(3)}& \mathbf{T}\mathit{y}=\mathbf{Tx}\mathit{\beta }+\mathbf{T}\mathit{\nu },\end{array}$$

where the n×n transformation matrix T is such that the components of the transformed noise Tν are independent and identically distributed (IID). This condition yields $\mathbf{T}\mathbf{\Sigma }{\mathbf{T}}^{\prime }=\mathbf{I}$ and an admissible choice for T is thus as follows:

$$\begin{array}{}\text{(4)}& \mathbf{T}=\mathbf{V}{\mathbf{\Delta }}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}^{\prime },\end{array}$$

where $\mathbf{\Delta }=\text{diag}({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_n)$ with ${\mathit{\lambda }}_{\mathrm{1}}\ge \mathrm{\dots }\ge {\mathit{\lambda }}_n>\mathrm{0}$ is the matrix of eigenvalues of Σ and $\mathbf{V}=({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\mathit{v}}_n)$ is its matrix of column eigenvectors, i.e. $\mathbf{\Sigma }=\mathbf{V}\mathbf{\Delta }{\mathbf{V}}^{\prime }$. In Eq. (3), the data are thus projected on the eigenvectors of the covariance matrix scaled by the inverse square root of their eigenvalues.

In most applications, the covariance Σ is not known and Eq. (2) cannot be implemented. This has given rise to many theoretical developments (e.g. Poloni and Sbrana, 2014; Naderahmadian et al., 2015) which often rely on further assumptions regarding the covariance Σ, in order to allow for its estimation simultaneously with β. In the applicative context of D&A considered here, the situation may be viewed as an in-between between known and unknown Σ because we have access to an ensemble of unforced climate model simulations. This ensemble can be used to extract valuable information about Σ. However, the resulting knowledge of Σ is imperfect for two reasons. First, the size of the ensemble is often very small, typically a few dozens, because of the high cost of numerical experiments requiring intensive computational resources. For the range of dimensions n usually considered, typically a few hundreds or thousands, the empirical covariance S may thus be non-invertible. Second, climate model simulations are known to inadequately represent the modes of internal variability at the smaller scale. Those smaller scale modes have lower magnitude and are expected to be associated with eigenvectors with low eigenvalues. Furthermore, the latter eigenvectors with low eigenvalues are also poorly sampled by construction and therefore are difficult to estimate reliably (North et al., 1982). For these two reasons, it is often assumed that only a small number r of leading eigenvectors of Σ, which are denoted $({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$ and r is typically on the order of 10, can be estimated from the control runs with enough reliability. Accordingly, it is assumed that the remaining n−r eigenvectors and eigenvalues are unknown.

Under this assumption, the inference on β must therefore be obtained based on these r leading eigenvectors only. The conventional approach presently used in D&A studies to tackle this problem follows from a simple and straightforward idea: adapting Eq. (4) to this situation by restricting the image of the transformation matrix T to the subspace generated by the r leading eigenvectors, thus leading to a truncated transformation having the following expression:

$$\begin{array}{}\text{(5)}& {\mathbf{T}}_r={\mathbf{V}}_r{\mathbf{\Delta }}_r^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}_r^{\prime },\end{array}$$

where ${\mathbf{\Delta }}_r=\text{diag}({\mathit{\lambda }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\lambda }}_r)$ and ${\mathbf{V}}_r=({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$. A new estimator thus follows by similarly applying the OLS estimator to the transformed data (T_rx,T_ry):

$$\begin{array}{}\text{(6)}& \begin{array}{l}{\widehat{\mathit{\beta }}}_r=\left({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}_r^{-\mathrm{1}}\mathbf{x}{)}^{-\mathrm{1}}\right({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}_r^{-\mathrm{1}}\mathit{y}),\end{array}\end{array}$$

where ${\mathbf{\Sigma }}_r^{-\mathrm{1}}={\mathbf{V}}_r{\mathbf{\Delta }}_r^{-\mathrm{1}}{\mathbf{V}}_r^{\prime }$ is the pseudo-inverse of ${\mathbf{V}}_r{\mathbf{\Delta }}_r{\mathbf{V}}_r^{\prime }$. Finally, from a mere terminology standpoint, applying the truncated transformation T_r to the data is often referred to in D&A as “projecting” the data onto the subspace spanned by the leading eigenvectors. However, T_r is in fact not a projection from a mathematical standpoint insofar as ${\mathbf{T}}_r^{\mathrm{2}}\ne {\mathbf{T}}_r$ hence we prefer to use the word “transformation”.

1.2 Objectives

The solution recalled above in Eq. (6) is easy to use and may seem to be a natural choice. It emerged as a popular approach in D&A arguably for these reasons. Nevertheless, this approach has two important limitations. Firstly, as pointed out by Hegerl et al. (1996), the leading eigenvectors sometimes fail to capture the signals x. Indeed, the leading eigenvectors capture variability well, but not necessarily the signal, merely because the corresponding physical processes may differ and be associated with different patterns. Secondly, by construction, the subspace spanned by the leading eigenvectors maximizes the amplitude of the noise, thereby minimizing the ratio of signal to noise and affecting the performance of the inference on β.

The first objective of this article is to highlight in detail these two limitations: both will be formalized and discussed immediately below in Sect. 2. Its second objective is to circumvent these limitations by building a new estimator of β that performs better than the conventional estimator ${\widehat{\mathit{\beta }}}_r$ of Eq. (6) while still making the same assumptions and using the same data. Most importantly, the central assumption here is that only the r leading eigenvectors of Σ can be estimated with enough accuracy from the control runs. Furthermore, we assume that r is given and is small compared to n (e.g. r is on the order of 10), and that N control runs are available to estimate the r leading eigenvectors with N≥r.

Section 2 highlights the limitations of the conventional estimator ${\widehat{\mathit{\beta }}}_r$. Section 3 introduces our proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$. Section 4 evaluates the benefit of our approach based on a realistic test bed, and illustrates it on actual observations of surface temperature. Section 5 concludes the paper.

2 Limitations of the conventional estimator

2.1 General considerations

It is useful to start this discussion by returning for a moment to the case where Σ is assumed to be fully known, i.e. its n eigenvectors are all known perfectly, as well as their eigenvalues. Under this assumption, let us evaluate the benefit of a given eigenvector v_k, $k=\mathrm{1},\mathrm{\dots },n$, for the inference of β. For this purpose, we propose to compare the variance $\text{Var}\left(\widehat{\mathit{\beta }}\right)=({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}\mathbf{x}{)}^{-\mathrm{1}}$ of the BLUE estimator of Eq. (2) obtained by using the transformation matrix $\mathbf{T}=\mathbf{V}{\mathbf{\Delta }}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}^{\prime }$, to the variance of the estimator ${\widehat{\mathit{\beta }}}_{-k}$ obtained by using the transformation matrix ${\mathbf{T}}_{-k}={\mathbf{V}}_{-k}{\mathbf{\Delta }}_{-k}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}_{-k}^{\prime }$ where V_−k and Δ_−k exclude the kth eigenvector and the kth eigenvalue respectively. From Eq. (2), these variances can be rewritten after some algebra:

$$\begin{array}{ll}{\displaystyle }\text{Var}\left(\widehat{\mathit{\beta }}\right)& {\displaystyle }=(\sum _{i=\mathrm{1}}^n{\mathbf{x}}^{\prime }{\mathit{v}}_i{\mathit{v}}_i^{\prime }\mathbf{x}/{\mathit{\lambda }}_i{)}^{-\mathrm{1}}\text{Var}({\widehat{\mathit{\beta }}}_{-k})\\ \text{(7)}& {\displaystyle }& {\displaystyle }=(\sum _{i=\mathrm{1},\mathrm{\dots },ni\ne k}{\mathbf{x}}^{\prime }{\mathit{v}}_i{\mathit{v}}_i^{\prime }\mathbf{x}/{\mathit{\lambda }}_i{)}^{-\mathrm{1}}.\end{array}$$

The Sherman–Morrison–Woodbury formula is a well known formula in linear algebra, which shows that the inverse of a rank-one correction of some matrix, can be conveniently computed by doing a rank one correction to the inverse of the original matrix. By applying the Sherman–Morrison–Woodbury formula to Eq. (7), we obtain the following:

$$\begin{array}{}\text{(8)}& \text{Var}\left({\widehat{\mathit{\beta }}}_{-k}\right)=\text{Var}\left(\widehat{\mathit{\beta }}\right)+{\displaystyle \frac{\text{Var}\left({\widehat{\mathit{\beta }}}_{-k}\right){\mathbf{x}}^{\prime }{\mathit{v}}_k{\mathit{v}}_k^{\prime }\mathbf{x}\text{Var}\left({\widehat{\mathit{\beta }}}_{-k}\right)}{{\mathit{\lambda }}_k+{\mathit{v}}_k^{\prime }\mathbf{x}\text{Var}\left({\widehat{\mathit{\beta }}}_{-k}\right){\mathbf{x}}^{\prime }{\mathit{v}}_k}}.\end{array}$$

Unsurprisingly, an immediate implication of Eq. (8) is that removing any eigenvector v_k results in an increase of the variance of the estimator. Indeed, the matrix in the second term of the right-hand side of Eq. (8) is positive definite; it can thus be interpreted as the cost in variance increase of excluding v_k or conversely as the benefit in variance reduction of including v_k. Moreover, the most valuable eigenvectors for inferring β in the above-defined sense are those combining a small eigenvalue (low λ_k) and a good alignment with the signals x (high x^′v_k). Accordingly, for any eigenvector v_k, its benefit is a decreasing function of λ_k, everything else being constant.

A similar conclusion can be reached by using a different approach. Consider this time the transformation ${\mathbf{T}}^{*}={\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}$. Since T^* is a p×n matrix, the data dimension is thus drastically reduced from n to p in the transformed Eq. (3). Furthermore, noting that the noise components T^*ν have covariance ${\mathbf{T}}^{*}\mathbf{\Sigma }{{\mathbf{T}}^{*}}^{\prime }$, and applying the GLS estimator to the transformed data $({\mathbf{T}}^{*}\mathit{y},{\mathbf{T}}^{*}\mathbf{x})$, we obtain an estimator ${\widehat{\mathit{\beta }}}^{*}$ which has the following expression:

$$\begin{array}{}\text{(9)}& {\widehat{\mathit{\beta }}}^{*}=\left({\mathbf{x}}^{\prime }{\mathbf{T}}^{*}\right({\mathbf{T}}^{*}\mathbf{\Sigma }{\mathbf{T}}^{*}{}^{\prime }{)}^{-\mathrm{1}}{\mathbf{T}}^{*}{}^{\prime }\mathbf{x}{)}^{-\mathrm{1}}\left({\mathbf{x}}^{\prime }{\mathbf{T}}^{*}\right({\mathbf{T}}^{*}\mathbf{\Sigma }{\mathbf{T}}^{*}{}^{\prime }{)}^{-\mathrm{1}}{\mathbf{T}}^{*}{}^{\prime }\mathit{y}).\end{array}$$

After some algebra, the above expression of ${\widehat{\mathit{\beta }}}^{*}$ simplifies to the following:

$$\begin{array}{}\text{(10)}& \begin{array}{ll}{\widehat{\mathit{\beta }}}^{*}& =\left({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}\mathbf{x}{)}^{-\mathrm{1}}\right({\mathbf{x}}^{\prime }{\mathbf{\Sigma }}^{-\mathrm{1}}\mathit{y}),\end{array}\end{array}$$

which is exactly the expression of the estimator of Eq. (2). Therefore, remarkably, the BLUE estimator obtained from the full data (y,x) is entirely recovered from the transformed data $({\mathbf{T}}^{*}\mathit{y},{\mathbf{T}}^{*}\mathbf{x})$ when using the subspace spanned by T^*. The latter subspace can thus be considered as optimal, in the sense that it reduces the data dimension to its bare minimum p while still preserving all the available information for inferring β. The p directions of the subspace spanned by T^* are denoted $(u_{\mathrm{1}}^{*},\mathrm{\dots },u_p^{*})$. They have the following expression:

$$\begin{array}{}\text{(11)}& u_k^{*}=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{{\mathit{v}}_i^{\prime }{x}_k}{{\mathit{\lambda }}_i}}{\mathit{v}}_i.\end{array}$$

Consistently with Eq. (8), Eq. (11) shows that for each $k=\mathrm{1},\mathrm{\dots },p$, the optimal direction $u_k^{*}$ emphasizes the eigenvectors that best trade off the smallest possible eigenvalue (i.e. small noise) with the best possible alignment between the eigenvector and the signal x_k (i.e. strong signal). It is interesting to notice that the foremost importance of small eigenvalues, which is underlined here, is also quite obvious in the formulation of optimal detection proposed by North and Stevens (1998). However, the importance of small eigenvalues was lost later on, as truncation on leading eigenvectors became a popular practice.

Figure 1 (a) Proportion of the total variance of the projected anthropogenic signal ${\mathbf{V}}_r^{\prime }{x}_{\mathrm{1}}$ (thick blue line) and of the projected noise ${\mathbf{V}}_r^{\prime }\mathit{\nu }$ (thick black line) as functions of r; uniform increase r∕n (light black line). (b) Signal-to-noise ratio of the projected data when using the projection subspace V_r consisting of the r leading eigenvectors (thick green line) and when using the orthogonal subspace ${\mathbf{V}}_r^{⟂}$ (thick red line) as functions of r. Both values are normalized to the signal-to-noise ratio of the full data.

In light of these considerations, the two main limitations of the conventional estimator ${\widehat{\mathit{\beta }}}_r$ of Eq. (6), mentioned in Sect. 1, appear more clearly. Firstly, in agreement with Hegerl et al. (1996), there is no guarantee that the leading eigenvectors will align well with the signals at stake. Secondly, the leading eigenvectors are by definition those with the largest eigenvalues. Therefore, a straightforward truncation onto the leading eigenvectors, such as the one used in ${\widehat{\mathit{\beta }}}_r$, maximizes the noise instead of minimizing it. While the first limitation may be an important pitfall in some cases, it will not be an issue whenever the signals align well with the leading eigenvectors, which nothing prevents in general. In contrast, the second limitation always manifests whenever ${\widehat{\mathit{\beta }}}_r$ is used and is thus arguably a more serious problem of the conventional approach. In the present paper, we thus choose to focus restrictively on this second limitation, and to not address the first one.

2.2 Illustration

In order to illustrate the second issue more concretely, we use surface temperature data from Ribes and Terray (2012), described in more detail in Section 4. The signals x considered here consist of the responses to anthropogenic and natural forcings, with p=2 and n=275. Both x and the covariance Σ are obtained from climate model simulations. The r leading eigenvectors ${\mathbf{V}}_r=({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$ are obtained from Σ and are assumed to be known. The proportion of the noise's total variance retained by projecting the data on V_r is denoted ${\mathit{\sigma }}_{\text{n}}^{\mathrm{2}}\left(r\right)$ while the proportion of the signal's total variance associated with the same projection is denoted ${\mathit{\sigma }}_{\text{s}}^{\mathrm{2}}\left(r\right)$. Since the projection matrix associated with projecting on the subspace V_r is equal to V_rV_r^′, we have the following:

$$\begin{array}{ll}{\displaystyle }{\mathit{\sigma }}_{\text{n}}^{\mathrm{2}}\left(r\right)& {\displaystyle }=\mathbb{E}\left(\widehat{\text{Var}}\left({\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathit{\nu }\right)\right)/\mathbb{E}\left(\widehat{\text{Var}}\left(\mathit{\nu }\right)\right)\\ \text{(12)}& {\displaystyle }& {\displaystyle }=\mathbb{E}\left({\mathit{\nu }}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathit{\nu }\right)/\mathbb{E}\left({\mathit{\nu }}^{\prime }\mathit{\nu }\right)=\sum _{i=\mathrm{1}}^r{\mathit{\lambda }}_i/\sum _{i=\mathrm{1}}^n{\mathit{\lambda }}_i,{\displaystyle }{\mathit{\sigma }}_{\text{s}}^{\mathrm{2}}\left(r\right)& {\displaystyle }=\mathbb{E}\left(\widehat{\text{Var}}\left({\mathbf{V}}_r{\mathbf{V}}_r^{\prime }{x}_i\right)\right)/\mathbb{E}\left(\widehat{\text{Var}}\left(x_i\right)\right)\\ \text{(13)}& {\displaystyle }& {\displaystyle }=x_i^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }{x}_i/x_i^{\prime }{x}_i\text{for}i=\mathrm{1},\mathrm{\dots },p,\end{array}$$

where for any vector z, $\widehat{\text{Var}}\left(\mathit{z}\right)$ denotes the empirical variance of z. Figure 1a plots ${\mathit{\sigma }}_{\text{n}}^{\mathrm{2}}\left(r\right)$ and ${\mathit{\sigma }}_{\text{s}}^{\mathrm{2}}\left(r\right)$ computed from Eq. (12) applied to the above-mentioned data. For simplicity, ${\mathit{\sigma }}_{\text{s}}^{\mathrm{2}}\left(r\right)$ is shown only for the anthropogenic signal x₁. As expected from the expression of ${\mathit{\sigma }}_{\text{n}}^{\mathrm{2}}\left(r\right)$ given in Eq. (12), it can be seen from this plot that ${\mathit{\sigma }}_{\text{n}}^{\mathrm{2}}\left(r\right)$ increases with r much faster than r∕n, since by definition the leading eigenvectors concentrate most of the noise's variance. In contrast, ${\mathit{\sigma }}_{\text{s}}^{\mathrm{2}}\left(r\right)$ increases with r roughly like r∕n: this shows that, in this example, the signal projects nearly uniformly on each eigenvector. For instance, when r=20, V_r involves 7 % of the n eigenvectors, and accordingly, projecting on V_r retains 8 % of the signal's total variance. In contrast, as much as 55 % of the noise's total variance is retained by projecting on V_r. Consequently, the signal-to-noise ratio of the projected data (see Fig. 1b) is reduced by a factor of 7 compared to the signal-to-noise ratio of the full data. Hence, it is fair to say that, in this case, the projection on the leading eigenvectors is greatly suboptimal, insofar as it drastically reduces the signal-to-noise ratio.

Let us now consider the quantities ${\mathit{\sigma }}_{\text{n}}^{\mathrm{2}*}\left(r\right)$ and ${\mathit{\sigma }}_{\text{n}}^{\mathrm{2}*}\left(r\right)$ defined as above, but projecting this time on the subspace orthogonal to V_r. By construction,

$$\begin{array}{}\text{(14)}& \begin{array}{l}{\mathit{\sigma }}_{\text{n}}^{\mathrm{2}*}\left(r\right)=\mathrm{1}-{\mathit{\sigma }}_{\text{n}}^{\mathrm{2}}\left(r\right),\\ {\mathit{\sigma }}_{\text{s}}^{\mathrm{2}*}\left(r\right)=\mathrm{1}-{\mathit{\sigma }}_{\text{s}}^{\mathrm{2}}\left(r\right).\end{array}\end{array}$$

As a direct consequence, for r=20, the projection on the subspace orthogonal to V_r captures 92 % = 100 %–8 % of the total variance of the signal, while the variance of the projected noise is only 20 % = 100 %–80 % of the total. Therefore, the signal-to-noise ratio of the projected data (see Fig. 1b) is this time magnified by a factor of 2 compared to the signal-to-noise ratio of the full data, and by a factor of 14 when comparing to the signal-to-noise ratio of the data projected on V_r. It can thus be seen that projecting the data on the orthogonal to V_r greatly magnifies the signal-to-noise ratio by filtering out a large fraction of the noise, while retaining most of the signal. These characteristics arguably make it a much more relevant projection subspace for the purpose of estimating β than the subspace of the leading eigenvectors. We further elaborate on this idea in Sect. 3.

3 An orthogonal estimator

3.1 Description

With these preliminary considerations in hand, we return to the situation of direct interest here: inferring β when only the r leading eigenvectors are available. In an attempt to improve the conventional estimator ${\widehat{\mathit{\beta }}}_r$, we propose the following new estimator:

$$\begin{array}{}\text{(15)}& {\widehat{\mathit{\beta }}}_r^{*}=({\mathbf{x}}^{\prime }\mathbf{x}-{\mathbf{x}}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathbf{x}{)}^{-\mathrm{1}}({\mathbf{x}}^{\prime }\mathit{y}-{\mathbf{x}}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathit{y}).\end{array}$$

Like the conventional estimator ${\widehat{\mathit{\beta }}}_r$, the proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$ only uses the data (y,x) and the r leading eigenvectors ${\mathbf{V}}_r=({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$. Its formulation stems from the above considerations regarding the suboptimality of the subspace T_r of Eq. (5) spanned by the leading eigenvectors, as we now highlight.

Let us denote ${\mathbf{V}}_r^{⟂}$ a $n\times (n-r)$ matrix consisting of n−r orthonormal column vectors spanning the subspace orthogonal to $({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$. In other words, we have the following by construction:

$$\begin{array}{}\text{(16)}& \begin{array}{l}{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }+{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }=\mathbf{I},\\ {{\mathbf{V}}_r^{⟂}}^{\prime }{\mathbf{V}}_r^{⟂}=\mathbf{I}.\end{array}\end{array}$$

Note that infinitely many matrices ${\mathbf{V}}_r^{⟂}$ satisfy Eq. (16), and that it is straightforward to construct one such matrix from the r leading eigenvectors $({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$, for instance by using the Gram–Schmidt algorithm. Furthermore, it is also important to note that, since ${\mathbf{V}}_r^{⟂}$ is orthogonal to V_r, it spans the same subspace as the subspace spanned by the n−r eigenvectors $({\text{v}}_{r+\mathrm{1}},\mathrm{\dots },{\mathit{v}}_n)$ associated with the smallest eigenvalues, because the latter are also orthogonal to $({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$ by definition. Therefore, it is not necessary to have knowledge of these n−r eigenvectors $({\text{v}}_{r+\mathrm{1}},\mathrm{\dots },{\mathit{v}}_n)$ in order to be able to project the data onto the latter subspace: only the knowledge of $({\mathit{v}}_{\mathrm{1}},\mathrm{\dots },{\text{v}}_r)$ is required for that purpose. Thus, with the knowledge of the r leading eigenvectors only, we are able to project the data onto the orthogonal subspace spanned by ${\mathbf{V}}_r^{⟂}$ using the transformation matrix ${\mathbf{T}}_r^{⟂}={\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }=\mathbf{I}-{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }$, thus obtaining $({\mathbf{T}}_r^{⟂}\mathit{y},{\mathbf{T}}_r^{⟂}\mathbf{x})$. The estimator proposed in Eq. (15) consists of using the OLS estimator corresponding to these transformed data:

$$\begin{array}{ll}{\displaystyle }{\widehat{\mathit{\beta }}}_r^{*}& {\displaystyle }=\left({\mathbf{x}}^{\prime }{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }\mathbf{x}{)}^{-\mathrm{1}}\right({\mathbf{x}}^{\prime }{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }\mathit{y})\\ \text{(17)}& {\displaystyle }& {\displaystyle }=({\mathbf{x}}^{\prime }\mathbf{x}-{\mathbf{x}}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathbf{x}{)}^{-\mathrm{1}}({\mathbf{x}}^{\prime }\mathit{y}-{\mathbf{x}}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathit{y}).\end{array}$$

It is interesting to note that the orthogonal projection subspace ${\mathbf{V}}_r^{⟂}$ does not actually need to be obtained to derive ${\widehat{\mathit{\beta }}}_r^{*}$ because only ${\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }$ matters for this purpose, and this matrix happens to be directly known from V_r through the orthogonality Eq. (16). This yields the right-most term of Eq. (17), which corresponds to the expression of ${\widehat{\mathit{\beta }}}_r^{*}$ given above in Eq. (15).

However, the choice of applying the OLS estimator to the projected data as implied by Eq. (17), rather than applying the GLS estimator, may seem surprising at first. The rationale for this choice is the following. Since we only assume knowledge of the r leading eigenvectors, the projection subspace ${\mathbf{V}}_r^{⟂}$ is thus by construction unknown from the standpoint of the covariance of the projected noise. In other words, while the projected noise ${\mathbf{T}}_r^{⟂}\mathit{\nu }$ in the projected regression equation remains Gaussian centred, we do not know anything about its covariance. Thus, we use by default the assumption that ${\mathbf{T}}_r^{⟂}\mathit{\nu }$ is IID with unknown variance σ², which is a commonplace assumption used in regression analysis when nothing is known about the noise's covariance. An immediate implication of this assumption is that the associated log likelihood is the following:

$$\begin{array}{}\text{(18)}& -\mathrm{2}\log \mathrm{\ell }(\mathit{\beta },{\mathit{\sigma }}^{\mathrm{2}})=n\log {\mathit{\sigma }}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}}}(\mathit{y}-\mathbf{x}\mathit{\beta }{)}^{\prime }{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }(\mathit{y}-\mathbf{x}\mathit{\beta }).\end{array}$$

The complete likelihood of Eq. (18) can be concentrated in β by maximizing out σ² to obtain the concentrated likelihood ℓ_c(β):

$$\begin{array}{ll}{\displaystyle }-\mathrm{2}\log {\mathrm{\ell }}_c\left(\mathit{\beta }\right)& {\displaystyle }=n\log \left\{(\mathit{y}-\mathbf{x}\mathit{\beta }{)}^{\prime }{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }(\mathit{y}-\mathbf{x}\mathit{\beta })\right\}\\ \text{(19)}& {\displaystyle }& {\displaystyle }=n\log \left\{(\mathit{y}-\mathbf{x}\mathit{\beta }{)}^{\prime }(\mathbf{I}-{\mathbf{V}}_r{{\mathbf{V}}_r}^{\prime }\left)\right(\mathit{y}-\mathbf{x}\mathit{\beta })\right\}.\end{array}$$

Our proposed OLS estimator ${\widehat{\mathit{\beta }}}_r^{*}$ results from the maximization of the above likelihood. The above considerations also have some implications with respect to the uncertainty analysis on the proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$. Indeed, in the present context, it is straightforward to obtain confidence intervals around ${\widehat{\mathit{\beta }}}_r^{*}$ based on the ratio of likelihoods ${\mathrm{\ell }}_c\left(\mathit{\beta }\right)/{\mathrm{\ell }}_c\left({\widehat{\mathit{\beta }}}_r^{*}\right)$, following a classic approach in statistics that is often used in D&A (e.g. Hannart et al., 2014).

It is important to underline that the latter approach to uncertainty quantification relies on the assumption that the projected noise ${\mathbf{T}}_r^{⟂}\mathit{\nu }$ is IID. Since the unknown true covariance of ${\mathbf{T}}_r^{⟂}\mathit{\nu }$ is in fact equal to ${\mathbf{\Sigma }}_r^{⟂}={\mathbf{V}}_r^{⟂}{\mathbf{\Delta }}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }$, this assumption is incorrect and hence the resulting confidence intervals are prone to be incorrect as well. An initial assessment under test bed conditions described in Section 4 suggested that the proposed approach is conservative in the sense that it leads to an overestimation of uncertainty (not shown). However, a more detailed assessment of the performance of this uncertainty quantification procedure is needed, but is beyond the scope of the present paper.

3.2 Extension to total least squares

The signals x are estimated as the empirical averages of several ensembles of finite size m, which introduces sampling uncertainty. In order to account for the latter, it is commonplace in D&A to reformulate the linear regression model of Eq. (1) as an error-in-variables (EIV) model (Allen and Stott, 2003):

$$\begin{array}{}\text{(20)}& \begin{array}{l}\mathit{y}={\mathbf{x}}^{*}\mathit{\beta }+\mathit{\nu },\\ \mathbf{x}={\mathbf{x}}^{*}+\mathit{\epsilon },\end{array}\end{array}$$

where x^* is now an unobserved latent variable, x is the ensemble average estimator of the true responses x^* and ε is a Gaussian centred noise term with covariance Σ∕m.

It is straightforward to generalize the approach exposed in Sect. 2.3 to this error-in-variables regression model. Likewise, the data (y,x), the latent variable x^*, and the noise terms ν and ε are projected on the orthogonal subspace using ${\mathbf{V}}_r^{⟂}$. The projected residuals are then also assumed IID with unknown variance σ². This yields the following complete likelihood:

$$\begin{array}{ll}{\displaystyle }-\mathrm{2}\log \mathrm{\ell }& {\displaystyle }(\mathit{\beta },{\mathbf{x}}^{*},{\mathit{\sigma }}^{\mathrm{2}})=n(p+\mathrm{1})\log {\mathit{\sigma }}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{2}}}}(\mathit{y}-{\mathbf{x}}^{*}\mathit{\beta }{)}^{\prime }{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }\\ \text{(21)}& {\displaystyle }& {\displaystyle }(\mathit{y}-{\mathbf{x}}^{*}\mathit{\beta })+{\displaystyle \frac{m}{{\mathit{\sigma }}^{\mathrm{2}}}}\text{Tr}\left\{(\mathbf{x}-{\mathbf{x}}^{*}{)}^{\prime }{\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }(\mathbf{x}-{\mathbf{x}}^{*})\right\},\end{array}$$

which can be concentrated in β by maximizing out σ² and x^* to obtain the concentrated likelihood ℓ_c(β):

$$\begin{array}{ll}{\displaystyle }-\mathrm{2}& {\displaystyle }\log {\mathrm{\ell }}_c\left(\mathit{\beta }\right)=n(p+\mathrm{1})\log \\ \text{(22)}& {\displaystyle }& {\displaystyle }\left\{{\displaystyle \frac{(\mathit{y}-\mathbf{x}\mathit{\beta }{)}^{\prime }(\mathbf{I}-{\mathbf{V}}_r{{\mathbf{V}}_r}^{\prime }\left)\right(\mathit{y}-\mathbf{x}\mathit{\beta })}{\mathrm{1}+\frac{\mathrm{1}}{m}{\mathit{\beta }}^{\prime }\mathit{\beta }}}\right\}.\end{array}$$

The so-called total least squares (TLS) estimator results from the maximization of the above concentrated likelihood in β, and is available under a closed form expression:

$$\begin{array}{}\text{(23)}& {\widehat{\mathit{\beta }}}_r^{*}=-{\displaystyle \frac{\sqrt{m}}{z_{p+\mathrm{1}}}}\mathit{z},\end{array}$$

where $\mathit{z}=(z_{\mathrm{1}},\mathrm{\dots },z_p{)}^{\prime }$ is a p-dimensional vector and z_p+1 is a scalar, such that $(z_{\mathrm{1}},\mathrm{\dots },z_{p+\mathrm{1}}{)}^{\prime }$ is the eigenvector associated with the smallest eigenvalue of the positive definite matrix $\mathbf{S}=\left[\sqrt{m}\mathbf{x}\mathit{y}{]}^{\prime }\right(\mathbf{I}-{\mathbf{V}}_r{{\mathbf{V}}_r}^{\prime }\left)\right[\sqrt{m}\mathbf{x}\mathit{y}]$ of size p+1. The estimator of Eq. (23) is a generalization of the OLS estimator of Eq. (17), the latter is indeed a special case of the former for $m=+\mathrm{\infty }$. Furthermore, in the particular case p=1 which is often referred to as the Deming regression, the estimator of Eq. (23) takes the following expression:

$$\begin{array}{}\text{(24)}& {\widehat{\mathit{\beta }}}_r^{*}={\displaystyle \frac{\sqrt{(m{S}_{xx}-S_{yy}{)}^{\mathrm{2}}+\mathrm{4}m{S}_{xy}^{\mathrm{2}}}-(m{S}_{xx}-S_{yy})}{\mathrm{2}{S}_{xy}}},\end{array}$$

where $S_{xx}=x^{\prime }(\mathbf{I}-{\mathbf{V}}_r{{\mathbf{V}}_r}^{\prime })x$, $S_{yy}={\mathit{y}}^{\prime }(\mathbf{I}-{\mathbf{V}}_r{{\mathbf{V}}_r}^{\prime })\mathit{y}$ and $S_{xy}=x^{\prime }(\mathbf{I}-{\mathbf{V}}_r{{\mathbf{V}}_r}^{\prime })\mathit{y}$. Next, following the same approach as in Sect. 2.3, it would be tempting to derive confidence intervals around ${\widehat{\mathit{\beta }}}_r^{*}$ based on the ratio of likelihoods ${\mathrm{\ell }}_c\left(\mathit{\beta }\right)/{\mathrm{\ell }}_c\left({\widehat{\mathit{\beta }}}_r^{*}\right)$. While this approach performs reasonably well in the context of EIV regression with known variance (Hannart et al., 2014), it is known to be problematic in the present context of EIV regression with unknown variance (Casella and Berger, 2008; Sect. 12.2). In this case, more suitable alternatives are available, e.g. an estimator of the variance of ${\widehat{\mathit{\beta }}}_r^{*}$ for p=1 is given by Fuller (1987) Sect. 1.3; it has the following expression:

$$\begin{array}{}\text{(25)}& \widehat{V}\left({\widehat{\mathit{\beta }}}_r^{*}\right)={\displaystyle \frac{\mathrm{1}}{n-r}}\left({\widehat{\mathit{\beta }}}_r^{*}{\displaystyle \frac{S_{xx}}{S_{xy}}}-\mathrm{1}\right)\left(m{\widehat{\mathit{\beta }}}_r^{*}{\displaystyle \frac{S_{xx}}{S_{xy}}}+{\widehat{\mathit{\beta }}}_r^{*\mathrm{2}}\right),\end{array}$$

from which confidence intervals can be obtained. When p > 1 the above variance estimator is applied successively to each component i of ${\widehat{\mathit{\beta }}}_r^{*}$ ($i=\mathrm{1},\mathrm{\dots },p$) by replacing y with $\mathit{y}-{\sum }_{k\ne i}{\widehat{\mathit{\beta }}}_{r,k}^{*}{x}_k$ in Eq. (25). It is worthwhile noting that various approaches are available in the literature to tackle the general problem of deriving confidence intervals for β under the present EIV regression model with unknown variance. While it would be worthwhile to investigate and to compare their performance in practice in a D&A context, such investigation is beyond the scope of this paper.

Summarizing, the proposed solution here to deal with the case of error in variables thus consists of using Eq. (23) to derive an estimate of β and Eq. (25) to derive a confidence interval.

3.3 Continuity considerations

The motivation of this subsection is merely to provide a more in-depth mathematical grounding to the proposed estimator, beyond the general considerations of Sect. 2.1. This subsection can safely be skipped without affecting the understanding of the remainder of this article, in so far as it does not provide any additional results or properties regarding the estimator itself.

The idea here is that, in addition to the general considerations of Sect. 2.1, the proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$ can also be justified by following a continuity extension argument. For any given x and y, let us consider the BLUE estimator of Eq. (2) as a function of the known covariance Σ. For this purpose, we define the function f on the set of positive definite matrices of dimension n, denoted ${\mathcal{S}}_n^{+*}\left(\mathbb{R}\right)$, as follows:

$$\begin{array}{}\text{(26)}& \begin{array}{l}f:{\mathcal{S}}_n^{+*}\left(\mathbb{R}\right)\to {\mathbb{R}}^p,\\ f\left(\mathbf{A}\right)=\left({\mathbf{x}}^{\prime }{\mathbf{A}}^{-\mathrm{1}}\mathbf{x}{)}^{-\mathrm{1}}\right({\mathbf{x}}^{\prime }{\mathbf{A}}^{-\mathrm{1}}\mathit{y}),\end{array}\end{array}$$

so that we have $\widehat{\mathit{\beta }}=f\left(\mathbf{\Sigma }\right)$.

Figure 2 (a) Scatter plot of one realization of the raw data (y,x) using the simulation test bed of Section 4: the observation vector y is plotted on the vertical axis and the first column vector x₁ of the matrix x is plotted on the horizontal axis. (b) Same as (a) for the transformed data Ty and Tx using the optimal transformation matrix $\mathbf{T}=\mathbf{V}{\mathbf{\Delta }}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}^{\prime }$. (c) Same as (b) when using the truncated transformation matrix ${\mathbf{T}}_r={\mathbf{V}}_r{\mathbf{\Delta }}_r^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}_r^{\prime }$ on the leading eigenvectors, with r=20. (d) Same as (b) when using the proposed transformation matrix ${\mathbf{T}}_r^{⟂}={\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }$ with r=20.

In the present context, we assumed that Σ is unknown but that ${\mathbf{\Sigma }}_r={\mathbf{V}}_r{\mathbf{\Delta }}_r{\mathbf{V}}_r^{\prime }$ is known. While it is tempting to apply f in Σ_r, unfortunately Σ_r is not positive definite so that f is not defined in Σ_r. Hence, an estimate of β cannot be obtained by straightforwardly applying the function f in Σ_r. Nevertheless, Σ_r belongs to the set of positive semi-definite matrices, denoted ${\mathcal{S}}_n^{+}\left(\mathbb{R}\right)$, which has the following interesting topological property:

$$\begin{array}{}\text{(27)}& {\mathcal{S}}_n^{+}\left(\mathbb{R}\right)=\stackrel{\mathrm{‾}}{{\mathcal{S}}_n^{+*}\left(\mathbb{R}\right)},\end{array}$$

where for any subset E, the notation $\stackrel{\mathrm{‾}}{E}$ used above in Eq. (27) denotes the adherence of E. Equation (27) thus implies that ${\mathcal{S}}_n^{+*}\left(\mathbb{R}\right)$ is dense in ${\mathcal{S}}_n^{+}\left(\mathbb{R}\right)$, which means that for any positive semi-definite matrix A there exists a sequence of positive definite matrices converging to A. Equation (27) therefore suggests the possibility to define a function $\stackrel{\mathrm{‾}}{f}$ which continuously extends f to the set of positive semi-definite matrices. For this purpose, let us consider the matrix A+εI with ε>0. It is immediate to show that for any $\mathbf{A}\in {\mathcal{S}}_n^{+}\left(\mathbb{R}\right)$ and for any arbitrarily small but positive ε, we have $\mathbf{A}+\mathit{\epsilon }\mathbf{I}\in {\mathcal{S}}_n^{+*}\left(\mathbb{R}\right)$. Taking advantage of this remark, we define $\stackrel{\mathrm{‾}}{f}$ by $\stackrel{\mathrm{‾}}{f}\left(\mathbf{A}\right)={lim}_{\mathit{\epsilon }\to {\mathrm{0}}^{+}}f(\mathbf{A}+\mathit{\epsilon }\mathbf{I})$. With this definition, $\stackrel{\mathrm{‾}}{f}$ can be obtained under a closed form expression:

$$\begin{array}{}\text{(28)}& \begin{array}{l}\stackrel{\mathrm{‾}}{f}:{\mathcal{S}}_n^{+}\left(\mathbb{R}\right)\to {\mathbb{R}}^p,\\ \stackrel{\mathrm{‾}}{f}\left(\mathbf{A}\right)=({\mathbf{x}}^{\prime }\mathbf{x}-{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathbf{x}{)}^{-\mathrm{1}}({\mathbf{x}}^{\prime }\mathit{y}-{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathit{y}).\end{array}\end{array}$$

The proof of Eq. (28) is detailed in Appendix A. The proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$ of Eq. (15) is therefore obtained by applying the extended function $\stackrel{\mathrm{‾}}{f}$ defined above, to the known, positive semi-definite matrix Σ_r:

$$\begin{array}{}\text{(29)}& {\widehat{\mathit{\beta }}}_r^{*}=\stackrel{\mathrm{‾}}{f}\left({\mathbf{\Sigma }}_r\right)=({\mathbf{x}}^{\prime }\mathbf{x}-{\mathbf{x}}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathbf{x}{)}^{-\mathrm{1}}({\mathbf{x}}^{\prime }\mathit{y}-{\mathbf{x}}^{\prime }{\mathbf{V}}_r{\mathbf{V}}_r^{\prime }\mathit{y}).\end{array}$$

4 Simulations and illustration

4.1 Data

We illustrate the method by applying it to surface temperature data and climate model simulations over the 1901–2000 period and over the entire surface of the Earth. For this purpose, we use the data described and studied in Ribes and Terray (2012). The temperature observations are based on the HadCRUT3 merged land–sea temperature data set (Brohan et al., 2006). The simulated temperature data used to estimate the response patterns are provided by results from the CMIP5 archive arising from simulations performed with the CNRM-CM5 model of Météo France under NAT and ANT experiments. For both experiments, the size of the ensembles is m=6. However, estimates of internal variability are based on intra-ensemble variability from 10 CMIP5 models (CNRM-CM5, CanESM2, HadGEM2-ES, GISS-E2-R, GISS-E2-H, CSIRO-Mk3-6-0, IPSL-CM5A-LR, BCC-CSM1-1, NorESM1-M, FGOALS-g2) under NAT and ANT experiments as well as pre-industrial simulations from both the CMIP3 and CMIP5 archives, leading to a large ensemble of 374 members. Standard pre-processing is applied to all observed and simulated data. Data are first interpolated on the 5^∘ × 5^∘ observational grid and then the spatio-temporal observational mask is applied. The dimension of the data set is then further reduced by computing decadal means and projecting the resulting spatial patterns onto spherical harmonics with resolution T4 to obtain vectors of size n=275 (see Ribes et al., 2012 for more details about ensembles and pre-processing). The ANT and NAT individual responses (p=2) are estimated by averaging the corresponding two ensembles to obtain the matrix of regressors x. The covariance matrix Σ is estimated as the empirical covariance of the large ensemble of control runs.

Figure 3 Performance results on simulated data: MSE of several estimators of β (vertical axis) as functions of the number of eigenvectors r (horizontal axis) assumed to be reliably estimable from control runs: conventional estimator ${\widehat{\mathit{\beta }}}_r$ computed based on estimated eigenvectors (thick green line) with N=30, and based on the true eigenvectors (light green line); proposed orthogonal estimator ${\widehat{\mathit{\beta }}}_r^{*}$ computed based on estimated eigenvectors (thick red line) with N=30, and based on the true eigenvectors (light red line); GLS estimator computed based on the true value of the covariance Σ (thick blue line).

4.2 Performance on simulated data

This section evaluates the performance of the estimators described above, by applying it to simulated values of y obtained from the linear regression equation $\mathit{y}=\mathbf{x}\mathit{\beta }+\mathit{\nu }$ assumed by the model of Eq. (1), where the noise ν is simulated from a multivariate Gaussian distribution with covariance Σ and where $\mathit{\beta }=(\mathrm{1},\mathrm{1}{)}^{\prime }$ is used. This setting thus assumes that x is known exactly, therefore it does not include the EIV situation prevailing in Eq. (20) and in Sect. 3.3 where x is measured with an error.

The use of simulated rather than real data aims at verifying that our proposed estimator performs correctly, and at comparing its performance with the conventional procedure, a goal which requires the actual values of β to be known. A sample of N=30 control runs was also simulated from a multivariate Gaussian distribution with covariance Σ, and the r leading eigenvectors and eigenvalues were estimated from these control runs.

Figure 2 shows scatter plots of one realization of the data (x,y) under this test bed, for the raw data (x,y) as well as for the transformed data (Tx,Ty), when using the three main transformations considered above: the optimal transformation $\mathbf{T}=\mathbf{V}{\mathbf{\Delta }}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}^{\prime }$ which assumes a fully known covariance Σ; the truncated transformation matrix ${\mathbf{T}}_r={\mathbf{V}}_r{\mathbf{\Delta }}_r^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}_r^{\prime }$ on the r leading eigenvectors, with r=20, and the proposed projection matrix ${\mathbf{T}}_r^{⟂}={\mathbf{V}}_r^{⟂}{{\mathbf{V}}_r^{⟂}}^{\prime }$ on the orthogonal to the r leading eigenvectors, again with r=20. Figure 2 allows the behavior of these transformations to be visualized. Firstly, the benefit of the optimal transformation $\mathbf{T}=\mathbf{V}{\mathbf{\Delta }}^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{V}}^{\prime }$ appears clearly, as it greatly enhances the signal. Secondly, it is also apparent that the standard truncated transformation using the r leading eigenvectors fails to enhance the signal, actually leading to a cloud of points which is more noisy than the original raw data. In contrast, the proposed transformation on the orthogonal to the r leading eigenvectors accurately captures part of the signal enhancement produced by the optimal transformation.

Figure 4 Illustration of temperature data with r=20. (a) Regression coefficients β_ANT and β_NAT estimated using the conventional estimator ${\widehat{\mathit{\beta }}}_r$: value of the estimate (green dots) and 95 % confidence interval (thick black lines). (b) Same as (a) under the proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$.

The conventional estimator ${\widehat{\mathit{\beta }}}_r$ and the proposed estimator ${\widehat{\mathit{\beta }}}_r^{*}$ of Eq. (15) were both derived based on the r leading estimators estimated from the N control runs. Then, these two estimators were derived again, this time using the true values of the r leading eigenvectors, for comparison. Finally, the GLS estimator of Equation (1) was derived based on the true value of the covariance Σ. The performance of each five estimators of β thus obtained was evaluated based on average mean squared error (MSE) ${\sum }_{j=\mathrm{1}}^J({\widehat{\mathit{\beta }}}_j-\mathit{\beta }{)}^{\prime }({\widehat{\mathit{\beta }}}_j-\mathit{\beta })/N$, where $j=\mathrm{1},\mathrm{\dots },J$ denotes the simulation number and J=1e4. These calculations were repeated for several values of r ranging from 5 to 30. Results are plotted in Fig. 3. They show an important gap in MSE between the proposed estimator and the conventional approach with the former exceeding the latter by a factor of 1000 for r=5, decreasing to 80 for r=30. This gap suggests a substantial benefit of using ${\widehat{\mathit{\beta }}}_r^{*}$ and strongly emphasizes the relevance of this approach compared to the conventional one.

When comparing the performance of ${\widehat{\mathit{\beta }}}_r^{*}$ using estimated vs. actual values of the r leading eigenvectors, a slight benefit of using actual values is found. This benefit increases with r, as eigenvectors with lower eigenvalues tend to be more difficult to estimate. In contrast, when comparing the performance of ${\widehat{\mathit{\beta }}}_r$ using estimated vs. actual values of the r leading eigenvectors, a slight benefit of using the estimated values is found. This superiority of the estimated, and thus incorrect, values may seem surprising at first, but is in fact natural and in line with our main point. Indeed, it can be explained by the fact that the estimated leading eigenvectors ${\widehat{\mathbf{V}}}_r$ match to a large extent with V_r but also span to some extent the true orthogonal subspace ${\mathbf{V}}_r^{⟂}$, which, according to our main point here, is much more suited for estimating β. Unsurprisingly, the estimator based on the estimated leading eigenvectors ${\widehat{\mathbf{V}}}_r$ thus performs better than its counterpart based on actual values.

4.3 Illustration on real data

The method is finally illustrated by applying it to actual observations of surface temperature y, with the same regressors x and covariance Σ as described above. The observed temperature data y is based on the HadCRUT3 merged land–sea temperature data set (Brohan et al., 2006). The method is applied for r=20. Further, we use a number of control runs N that is small compared to n, which is a common situation, even though we have access under the present test bed to 374 control runs. Thus, N=30 control runs are randomly selected from the 374 control runs available and the eigenvectors are estimated from these 30 runs. For the inference of β and confidence intervals, we use the solution described in Sect. 3.2 which consists of using Eq. (23) to derive an estimate of β and Eq. (25) to derive a confidence interval. Results are shown in Fig. 4. For the scaling factor corresponding to anthropogenic forcings, both methods yield an estimate which is close to 1 and significantly positive. However, the level of uncertainty obtained with the proposed estimator is 3 times smaller than when using the conventional method. An even wider uncertainty gap is observed for the coefficient corresponding to natural forcings. In this case, the confidence interval includes zero under the conventional approach, but it does not under the proposed orthogonal approach. However, as already stated in Sect. 3, a caveat in the latter result is that the uncertainty quantification procedure used here has not been evaluated in detail, and thus we do not know for sure that the obtained uncertainty estimate is well calibrated.

5 Conclusions

We have introduced a new estimator of the vector β of linear regression coefficients, adapted to the context where only the r leading eigenvectors of the noise's covariance Σ are known, which is an assumption often relevant in climate change attribution studies. General considerations have first shown that when the covariance is known, its most relevant eigenvectors for projection are those trading off the smallest possible eigenvalue with the best possible alignment to the signal. The optimal direction is thus associated in general with the smallest eigenvalues, not with the highest ones. Our proposed estimator builds upon this finding and is based on projecting the data on a subspace which is orthogonal to the r leading eigenvectors.

When applied on a simulation test bed that is realistic with respect to D&A applications, we find that the proposed estimator outperforms the conventional estimator by several orders of magnitude for small values of r. Our proposal is thus prone to significantly affect D&A results by decreasing the estimated level of uncertainty.

Substantial further work is needed to evaluate the performance of the proposed uncertainty quantification, in particular in an EIV context. Such an evaluation was beyond the scope of the present work. Its primary focus was to demonstrate that the choice of the leading eigenvectors as a projection subspace is a vastly suboptimal one for D&A purposes; and that projecting on a subspace which is orthogonal to the r leading eigenvectors yields a significant improvement.

Appendix A: Proof of equation (28)

Let $\mathbf{A}\in {\mathcal{S}}^{+}\left(\mathbb{R}\right)$ be a non-invertible positive matrix. Denote Δ_A the matrix of non-zero eigenvalues of A, V_A the matrix of eigenvectors associated with Δ_A and ${\mathbf{V}}_{\mathrm{A}}^{⟂}$ the matrix of eigenvectors associated with the eigenvalue zero. For ε>0, we have the following:

$$\begin{array}{}\text{(A1)}& \begin{array}{ll}(\mathbf{A}+\mathit{\epsilon }\mathbf{I}{)}^{-\mathrm{1}}& =({\mathbf{V}}_{\mathrm{A}}{\mathbf{\Delta }}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }+\mathit{\epsilon }\mathbf{I}{)}^{-\mathrm{1}}\\ & =({\mathbf{V}}_{\mathrm{A}}{\mathbf{\Delta }}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }+\mathit{\epsilon }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }+\mathit{\epsilon }{\mathbf{V}}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }{)}^{-\mathrm{1}}\\ & ={\left({\mathbf{V}}_{\mathrm{A}}({\mathbf{\Delta }}_{\mathrm{A}}+\mathit{\epsilon }\mathbf{I}){\mathbf{V}}_{\mathrm{A}}^{\prime }+\mathit{\epsilon }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }\right)}^{-\mathrm{1}}\\ & ={\mathbf{V}}_{\mathrm{A}}({\mathbf{\Delta }}_{\mathrm{A}}+\mathit{\epsilon }\mathbf{I}{)}^{-\mathrm{1}}{\mathbf{V}}_{\mathrm{A}}^{\prime }+\frac{\mathrm{1}}{\mathit{\epsilon }}{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }.\end{array}\end{array}$$

Therefore,

$$\begin{array}{ll}{\displaystyle }f(\mathbf{A}+\mathit{\epsilon }\mathbf{I})& {\displaystyle }={\left({\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}({\mathbf{\Delta }}_{\mathrm{A}}+\mathit{\epsilon }\mathbf{I}{)}^{-\mathrm{1}}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathbf{x}+{\displaystyle \frac{\mathrm{1}}{\mathit{\epsilon }}}{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }\mathbf{x}\right)}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }\left({\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}({\mathbf{\Delta }}_{\mathrm{A}}+\mathit{\epsilon }\mathbf{I}{)}^{-\mathrm{1}}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathit{y}+{\displaystyle \frac{\mathrm{1}}{\mathit{\epsilon }}}{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }\mathit{y}\right)\\ {\displaystyle }& {\displaystyle }={\left({\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}{\mathbf{B}}_{\mathit{\epsilon }}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathbf{x}+{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }\mathbf{x}\right)}^{-\mathrm{1}}\\ \text{(A2)}& {\displaystyle }& {\displaystyle }\left({\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}{\mathbf{B}}_{\mathit{\epsilon }}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathit{y}+{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }\mathit{y}\right),\end{array}$$

with ${\mathbf{B}}_{\mathit{\epsilon }}=\mathit{\epsilon }({\mathbf{\Delta }}_{\mathrm{A}}+\mathit{\epsilon }\mathbf{I}{)}^{-\mathrm{1}}$. We have ${lim}_{\mathit{\epsilon }\to {\mathrm{0}}^{+}}{\mathbf{B}}_{\mathit{\epsilon }}=\mathrm{0}$, hence,

$$\begin{array}{}\text{(A3)}& \begin{array}{ll}{lim}_{\mathit{\epsilon }\to {\mathrm{0}}^{+}}& f(\mathbf{A}+\mathit{\epsilon }\mathbf{I})={\left({\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathbf{A}}^{⟂}{{\mathbf{V}}_{\mathbf{A}}^{⟂}}^{\prime }\mathbf{x}\right)}^{-\mathrm{1}}\left({\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}^{⟂}{{\mathbf{V}}_{\mathrm{A}}^{⟂}}^{\prime }\mathit{y}\right)\\ & =({\mathbf{x}}^{\prime }\mathbf{x}-{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathbf{x}{)}^{-\mathrm{1}}({\mathbf{x}}^{\prime }\mathit{y}-{\mathbf{x}}^{\prime }{\mathbf{V}}_{\mathrm{A}}{\mathbf{V}}_{\mathrm{A}}^{\prime }\mathit{y}),\end{array}\end{array}$$

which proves Eq. (28).