1

Figure 8.4. Wald-Bartlett procedure. The bivariate sample {x(i), y(i)}™=1 is divided into three groups of same size according to the size of the x values; if n is not divisible by 3, then take the closest grouping. Let j index the size-sorted sample. Let the averages of {x(j)}™=1 and {y(j)}™=1, denoted as X1 and y1, define the first group's centre (P1, cross), and let the averages of {x(j)}"=2n/3+1 and {y(j)}"=2n/3+1, denoted as X3 and y3, define the third group's centre (P3). The line P1 P3 (long-dashed) defines the Wald-Bartlett regression estimate of the slope, J = (y3 — yj1)/(X3 — X1). The centre of the complete sample (P) is defined via the averages of {x(j)}n=1 and {y(j)j=1, denoted as x and y. The Wald-Bartlett intercept estimate, f30 = y — /31x, completes the linear fit (solid line).

is made on the noise-influenced observations. Monte Carlo simulations, similar to those in Section 8.3, reveal that the inconsistency leads to an inacceptably poor coverage performance of bootstrap CIs (for /o and /31) (not shown). This limits severely the applicability of the Wald-Bartlett procedure to real-world climatological problems, where the Xtrue(i) are usually unknown.

Wald (1940: p. 298 therein) notes that if prior knowledge exists on the standard deviation ratio, then a consistent estimation could be constructed. This situation is similar to WLSXY estimation (Section 8.1.2.1).

The calculation of classical CIs (Wald 1940; Bartlett 1949) via the Student's t distribution assumed prior knowledge to be available, allowing a consistent estimation, and the errors, Xnoise(i) and Ynoise(i), to be serially independent and of Gaussian shape.

8.2 Bootstrap confidence intervals

Classical CIs are based on the PDF of an estimator (Chapter 3). The PDF can be analytically determined unless the situation (estimation

problem, noise properties) becomes too complex. The construction of classical CIs for the linear errors-in-variables model (Wald 1940; Bartlett 1949; York 1966; Fuller 1987) made a number of assumptions from the following:

1. Gaussian distributional shapes of the noise components, Xnoise(i) and

Ynoise(i);

2. absence of autocorrelation in the noise components;

3. absence of correlation between X(i) and Xnoise(i) as well as between

4. absence of correlation between Xnoise(i) and Ynoise(i).

Some authors treat the correlation effects (points 3 and 4) and non-Gaussian errors (point 1), see the background material (Section 8.7). However, allowance for autocorrelations (point 2) seems to have been made by none.

Here we are interested in linearly relating two climate processes, X(i) and Y(i), and our sample, {t(i),x(i),y(i)}™=1, contains the time. The previous chapters document that non-Gaussian distributions and persistence phenomena are typical of climate processes. We cannot therefore expect the classical method to yield accurate results for climate data. This is, as in previous chapters, the reason to consider the bootstrap method. An additional point is incomplete knowledge about the noise components. Often we have no or only limited information about SY(i), SX(i) or their (squared) ratio, A. Such incomplete knowledge, which may widen the CI, is quantifiable using bootstrap resampling (Booth and Hall 1993).

One resampling algorithm is the pairwise-MBB, which has been found useful in the context of correlation estimation (Algorithm 7.2).

The other algorithm, introduced here for the purpose of enhancing the coverage performance in the context of fitting errors-in-variables models, is called pairwise-moving block bootstrap resampling of residuals or pairwise-MBBres. It is based on the observation that the (linear) errors-in-variables regression (Eq. 8.1) is a model with a deterministic (linear) component. Since pairwise resampling seems to be handicapped in the presence of deterministic components (Chapter 4), the idea of the pairwise-MBBres algorithm is to take the fit and the regression residuals, apply pairwise-MBB to the residuals and add the resampled residuals to the fit. The new approach is that the residuals, eX(i) and eY(i), are in two dimensions (Fig. 8.5). The pairwise-MBBres is given as Algorithm 8.1.

Figure 8.5. Pairwise-MBBres algorithm, definition of residuals. The line La measures the distance from a data point to the fit line (Fig. 8.3). The residuals (dashed lines) are given by eX (i) = [Jo + /J • x(i) — y(i)]/[A//3i + /i and eY(i) = —A• eX (i)//?i.

8.2.1 Simulating incomplete prior knowledge

Assume for the convenience of exposition homoscedastic noise components, SX(i) = SX and SY(i) = SY. For achieving an identifiable problem, OLSBC estimation requires information, not contained in the sample, about SX (Section 8.1.1.2); analogously, WLSXY estimation requires information about both Sx and Sy, or about their ratio, 5 = A1/2 = Sy/Sx (Section 8.1.2.1).

In practical applications such prior knowledge is not always exact. SX or 5 are then described by random variables. Bootstrap resampling can be augmented by a simulation step, where random numbers are drawn from the distribution of Sx or 5. This increases the uncertainty of the OLSBC or WLSXY estimates, leading to wider bootstrap CIs compared to a situation with exact prior knowledge (Booth and Hall 1993).

In the Monte Carlo experiments studying incomplete prior knowledge, we use the model

where Eu[1-a; 1+A](i) is an IID random process with a uniform distribution over the interval [1 - A; 1 + Aj. For example, A = 0.5 specifies a situation where we only know 5 to lie between 0.5 and 1.5 times its true value. Other models are possible.

The construction of bootstrap CIs (Algorithm 8.1) is adapted at Step 8, the calculation of the replications. Instead of applying to the resam-pled data the same estimation procedure that is used for the original data, an adapted estimation is performed (Steps 8a and 8b).

 Step 1 Bivariate time series {i(i),x(i),»(i)>n=1 Step 2 Parameter estimates from OLSBC, WLSXY or Wald-Bartlett procedure ?0,?1 Step 3 Residuals (Fig. 8.5) ex (i), ey (i) Step 4 Fit values xfit (i) = x(i) - ex (i), yfit(i) = y(i) - ey (i) Step 5 Bias-corrected AR(1) parameters, aX = ?y estimated on residuals, block length selection I after Eqs. (7.31) and (7.32) Step 6 Resampled residuals, pairwise-MBB with I { ex (i), eyb (iU n=1 (b, counter) Step 7 Resample x*b(i) = xfit(i)+ eX (i), y*b(i) = yfit(i) + eyb(i),i = 1,...,n Step 8 Bootstrap replications ?o*b,?ib Step 9 Bootstrap prediction y*b(n + 1) = /?*b + ?Obx(n + 1) Step 10 Go to Step 6 until b = B (usually B = 2000) replications exist {?b}f=i, {?*b}f=i, {?b(n + 1)}f=i Step 11 Calculate CIs (Section 3.4)

Algorithm 8.1. Construction of bootstrap confidence intervals for parameters of the linear errors-in-variables regression model, pairwise-MBBres resampling, even spacing. In case of uneven spacing, Step 5 uses ?X = ?y. Step 8 can be adapted as follows for taking incomplete prior knowledge into account. Step 8a: simulate A* , SX ; Step 8b: use OLSBC with S*X instead of Sx (Eq. 8.6), use WLSXY with A* instead of A (Eq. 8.10). Steps 9 and 10: ?ob(n + 1) refers to prediction (Section 8.5).

The first group of experiments (Section 8.3.1) adopts an "easy" setting, where distributional shapes are Gaussian and prior knowledge is exact. The results confirm the success of block resampling methods in preserving serial dependence. The second group of experiments (Section 8.3.2) shows that for realistic settings, with non-Gaussian distributions or incomplete knowledge, the results are less exact. It appears that then WLSXY estimation combined with pairwise-MBBres resampling yields the relatively best results. The third group (Section 8.3.3) quantifies coverage accuracy and RMSE in dependence on the accuracy of the prior knowledge (standard deviation ratio). It demonstrates that even with n ^ to, the RMSE (and the CI length) for /?i (slope) does not go to zero, but rather approaches finite values. On the other hand (intercept), with n ^ to does RMSE^o ^ 0. The last group (Section 8.3.4), finally, explores what happens when we mis-specify the degree of how accurately we know the standard deviation ratio.

8.3.1 Easy setting

The easy setting (Gaussian shapes of Xnoise(i) and Ynoise(i), complete prior knowledge) is further simplified when no autocorrelation resides in the noise components. Table 8.1 exhibits excellent coverage performance of bootstrap CIs with pairwise-MBBres resampling already for sample sizes as small as 20. The excellent performance regards both parameters (intercept and slope) and all estimation procedures (OLSBC, WLSXY and Wald-Bartlett). Similar results for OLSBC and WLSXY were obtained using pairwise-MBB resampling (results not shown).

If there exists autocorrelation, then pairwise-MBBres resampling successfully preserves it, where it does not matter whether the AR(1) parameters are known or, which is more realistic, have to be estimated. However, there is one exception: the Wald-Bartlett estimation of the intercept fails completely, independent of the sample size. Also OLSBC and WLSXY estimations of the intercept are of reduced CI coverage accuracy, but become acceptable for n above, say, 200. It is not clear why intercept estimation is more problematic than slope estimation.

We remark that the inacceptable performance of the Wald-Bartlett procedure occurred even in the presence of knowledge of the size of the true predictor values, which in turn enabled a perfect grouping (independent of the predictor noise). The three points of the (1) inacceptable performance of CIs for the intercept, (2) not better (compared with OLSBC and WLSXY) performance of CIs for the slope and (3) rather strong requirement of knowledge of the size of true predictor values—

Table 8.1. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal shape and complete prior knowledge: CI coverage performance. nsim = 47,500 random samples were generated from Xtrue(i) = En(0, 1)(i) and Y (i) = /0 +/1Xtrue(i)+SY •Ynoise(i), i = 1,...,n, with /0 = 1.0 and /1 = 2.0. Predictor noise was subsequently added, X(i) = Xtrue(i) + SX • Xnoise(i). The Xnoise(i) and Ynoise(i) are mutually independent Gaussian AR(1) processes for even spacing (Eq. 2.1) with parameters aX and aY, respectively. Construction of bootstrap CIs used pairwise-MBBres resampling (Algorithm 8.1), block length selection after Eqs. (7.31) and (7.32), the Student's t interval type (v = n — 2), B = 2000 and confidence level 95%. Prior knowledge of SX = 0.25, SY = 0.5 and the size of the {xtrue(i)}n=1 was exact and utilized in the estimations; AR(1) parameters are in two cases known, in one case unknown and estimated with bias correction.

Estimation method Estimation method

 aX = aY = 0.0 (known ) 10 0.933 0.928 0.933 0.955 0.928 0.938 0.950 20 0.939 0.939 0.942 0.949 0.938 0.942 0.950 50 0.944 0.947 0.948 0.949 0.946 0.947 0.950 100 0.944 0.946 0.947 0.949 0.947 0.946 0.950 200 0.945 0.949 0.950 0.949 0.947 0.945 0.950 500 0.946 0.949 0.949 0.950 0.946 0.946 0.950 1000 0.945 0.948 0.949 0.949 0.948 0.945 0.950 aX = aY = 0.3 (known ) 10 0.839 0.831 0.785 0.949 0.919 0.912 0.950 20 0.863 0.862 0.816 0.947 0.936 0.926 0.950 50 0.895 0.896 0.827 0.948 0.945 0.930 0.950 100 0.909 0.911 0.827 0.947 0.945 0.933 0.950 200 0.921 0.924 0.836 0.947 0.945 0.939 0.950 500 0.931 0.933 0.840 0.948 0.944 0.941 0.950 1000 0.935 0.938 0.844 0.950 0.949 0.945 0.950
 10 0.807 0.797 0.841 0.922 0.892 0.935 0.95 20 0.871 0.869 0.853 0.943 0.932 0.947 0.95 50 0.896 0.897 0.852 0.946 0.943 0.947 0.95 100 0.911 0.913 0.854 0.948 0.947 0.949 0.95 200 0.921 0.923 0.853 0.949 0.948 0.949 0.95 500 0.93 0.932 0.85 0.949 0.945 0.948 0.95 1000 0.934 0.937 0.849 0.949 0.948 0.949 0.95

a Standard errors of Y/j0 and 7/ are nominally 0.001. b Wald—Bartlett procedure.

a Standard errors of Y/j0 and 7/ are nominally 0.001. b Wald—Bartlett procedure.

Table 8.2. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal shape and complete prior knowledge: CI coverage performance (continued). Design is identical to the previous experiment (Table 8.1), with the exception that autocorrelation is ignored at CI construction.

Estimation method Estimation method

OLSBC WLSXY WBb OLSBC WLSXY WBb aX = aY =0.3 (ignored)

OLSBC WLSXY WBb OLSBC WLSXY WBb aX = aY =0.3 (ignored)

 10 0.846 0.838 0.846 0.955 0.928 0.939 0.95 20 0.842 0.84 0.845 0.948 0.937 0.942 0.95 50 0.845 0.847 0.849 0.948 0.944 0.946 0.95 100 0.842 0.846 0.846 0.947 0.946 0.946 0.95 200 0.845 0.848 0.849 0.947 0.946 0.945 0.95 500 0.846 0.849 0.849 0.948 0.945 0.947 0.95 1000 0.846 0.848 0.849 0.946 0.948 0.946 0.95

a Standard errors of and 7g are nominally 0.001. b Wald—Bartlett procedure.

a Standard errors of and 7g are nominally 0.001. b Wald—Bartlett procedure.

provide enough support to exclude the Wald-Bartlett procedure from consideration in the further experiments (which have more realistic settings).

One experiment under the easy setting studied what happens when autocorrelation is ignored (Table 8.2). This was achieved by prescribing positive AR(1) parameters, aX and aY, of both noise components and resetting their estimate values to zero, aX = 0 and aY = 0. The detrimental effect, again on /a but not /?i, was an underestimated bootstrap standard error, which led to too narrow CIs and too low coverages. (Similar results where found when replacing pairwise-MBBres by pairwise-MBB resampling.)

The major findings of the first experiment on CI coverage accuracy (Table 8.1) are reflected in the results on empirical RMSE (Table 8.3). Autocorrelation increases the estimation error of the intercept, but not of the slope. A larger data size means a smaller estimation error of the intercept and the slope. For n ^ to, both RMSE^o and RMSE^ appear to go to zero.

OLSBC and WLSXY estimation methods perform similarly; only for small n and slope estimation does WLSXY seem to give smaller error bars than OLSBC.

 P0 RMSEb pi Estimation method Estimation method OLSBC WLSXY OLSBC WLSXY aX = ay = 0.0 10 0.237 0.245 0.645 0.289 20 0.161 0.164 0.190 0.178 50 0.099 0.101 0.110 0.105 100 0.070 0.071 0.076 0.073 200 0.049 0.050 0.053 0.052 500 0.031 0.032 0.033 0.033 1000 0.022 0.022 0.023 0.023 ax = ay = 0.3 10 0.300 0.310 5.759 0.275 20 0.211 0.216 0.188 0.174 50 0.134 0.137 0.110 0.105 100 0.094 0.096 0.076 0.073 200 0.067 0.068 0.053 0.052 500 0.042 0.043 0.034 0.033 1000 0.030 0.031 0.024 0.023

a Empirical RMSE^, given by b Empirical RMSE^ , given by

a Empirical RMSE^, given by b Empirical RMSE^ , given by

8.3.2 Realistic setting: incomplete prior knowledge

The setting becomes more complex, or realistic, when the prior knowledge about the standard deviations of the measurement noise is not complete. We study (Table 8.4) a situation where the true ratio is 5 = SY/SX = 2.0 but one knows only that 5 is between 1.0 and 3.0. The adapted bootstrap CI construction (WLSXY with A* = (5*)2), for both /0 and /31, yields acceptable accuracies for normal shape and n ^ 200— under the condition that pairwise-MBBres resampling is employed. (For climatological purposes, a 95% CI may be "acceptable" if the true coverage is between, say, 92 and 98%.) The pairwise-MBBres resampling method (Fig. 8.5) is clearly superior to the pairwise-MBB method.

Bootstrap CI construction for OLSBC estimates failed to achieve the accuracies for WLSXY—for both resampling methods. Rather large

Table 8.4. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal/lognormal shape and incomplete prior knowledge: CI coverage performance. Design is identical to the first experiment (Table 8.1), with the following exceptions: (1) autocorrelation parameters are unknown (and estimated with bias correction) and (2) Ynoise(i) has normal or lognormal shape. Estimation and CI construction is identical to the first experiment (Table 8.1), with the following exceptions: (1) the Wald-Bartlett procedure is omitted; (2) prior knowledge of SX = 0.25, SY = 0.5 (5 = 2.0) is incomplete after Eqs. (8.11) and (8.12) with A = 0.5; and (3) CI construction is adapted accordingly (Section 8.2.1).

Nominal

Estimation method

Estimation method

Estimation method

OLSBC

WLSXY

OLSBC

WLSXY

— Y noise

(i) normal shape, pairwise-MBBres

10

0.809

0.802

0.941

0.895

0.950

20

0.871

0.871

0.956

0.931

0.950

50

0.898

0.899

0.955

0.940

0.950

100

0.909

0.913

0.952

0.942

0.950

200

0.921

0.924

0.948

0.943

0.950

500

0.930

0.934

0.939

0.947

0.950

1000

0.937

0.939

0.927

0.952

0.950

2000

0.936

0.939

0.919

0.958

0.950

5000

0.941

0.944

0.909

0.960

 10 0.856 0.875 0.98 0.947 0.95 20 0.864 0.867 0.972 0.944 0.95 50 0.862 0.862 0.958 0.943 0.95 100 0.866 0.866 0.955 0.943 0.95 200 0.865 0.865 0.954 0.948 0.95 500 0.871 0.871 0.943 0.947 0.95 1000 0.873 0.871 0.932 0.952 0.95 2000 0.88 0.88 0.922 0.957 0.95 5000 0.891 0.889 0.915 0.962 0.95
 ax = 0.3, ay = 0.8, Ynoise (i) lognormal shape, pairwise- MBBres 10 0.689 0.673 0.942 0.871 0.950 20 0.738 0.738 0.956 0.902 0.950 50 0.788 0.793 0.961 0.904 0.950 100 0.817 0.824 0.960 0.897 0.950 200 0.849 0.856 0.959 0.897 0.950 500 0.880 0.887 0.949 0.902 0.950 1000 0.894 0.901 0.936 0.917 0.950 2000 0.912 0.919 0.925 0.929 0.950 5000 0.924 0.931 0.917 0.947 0.950

Standard errors of Yß and Yß, are nominally 0.001.

Standard errors of Yß and Yß, are nominally 0.001.

sample sizes (n = 2000 and 5000) reveal the "worrisome" behaviour of Y^ for the OLSBC estimates: they do not saturate and approach the nominal value of 0.95 but seem rather to drift away for large n.

It becomes clear that for realistic settings (autocorrelation, incomplete prior knowledge), WLSXY estimation combined with pairwise-MBBres resampling is the only one-loop option to achieve acceptable levels of CI accuracy. A second loop of resampling (calibration or bootstrap-t) may in principle improve the accuracy, also for errors-in-variables regression (Booth and Hall 1993).

The combination of WLSXY and pairwise-MBBres performed well (Table 8.4) also for a rather difficult setting (stronger, unequal autocorrelations, lognormal shape). It is interesting to note that slope estimation yielded more accurate results than intercept estimation. The data size requirements, however, become rather strong (Table 8.4). Obtaining accurate results for data sizes in the range of 500 and below may require calibration methods.

8.3.3 Dependence on accuracy of prior knowledge

In practical situations, our prior knowledge about the measurement standard errors or their ratio, 5 = SY/SX, may depend to a considerable degree on how good we know the measurement devices (calibration standards, replication analyses, etc.) or the archives "containing" the data (sampling error). The accuracy of that knowledge, parameterized here in form of A (Eqs. 8.11 and 8.12), should influence the estimation RMSE and possibly also the CI accuracy. This is explored by means of a set of simulation experiments (Tables 8.5 and 8.6), where A is varied.

The selection of the other setting parameters follows the previous Monte Carlo experiments in this section: intermediate sizes of autocorrelation, Gaussian shape and a true standard deviation ratio of 5 = 2.0. The data size may take relatively large values (n = 2000 and 5000) because also the limiting behaviour is of interest. We employ the WLSXY estimation and Student's t CI constructed by means of pairwise-MBBres resampling.

The resulting coverages (Table 8.5) approach with increasing data size the nominal value—as they should. In general, the levels are acceptable from n above, say, 200 (slope estimation) or 500 (intercept estimation). In the case of slope estimation, a highly inaccurate prior knowledge (A = 0.9) may require more data points for achieving a coverage level similar to values found for smaller inaccuracies (A < 0.7).

The resulting RMSE values (Table 8.6) for intercept estimation approach zero with increasing data size. The rate of this convergence seems not to depend on the accuracy of the prior knowledge (A). The RMSE

Table 8.5. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal shape: influence of accuracy of prior knowledge on CI coverage performance. Design is identical to the previous experiment (Table 8.4), with the following fixed setting: (1) autocorrelation parameters are ax = ay = 0.3, (2) both noise components have normal shape. Estimation and CI construction is identical to the previous experiment (Table 8.4), with the following exceptions: (1) only WLSXY estimation with pairwise-MBBres resampling is considered; (2) prior knowledge of Sx = 0.25, Sy = 0.5 (S = 2.0) is incomplete after Eq. (8.11) with various A values.

n Ya Nominal

Accuracy of prior knowledge n Ya Nominal

Accuracy of prior knowledge

 A = 0.1 A = 0.3 A = 0.5 A = 0.7 A = 0.9 Intercept estimation 10 0.774 0.772 0.802 0.773 0.782 0.950 20 0.868 0.867 0.871 0.869 0.873 0.950 50 0.895 0.898 0.899 0.897 0.902 0.950 100 0.913 0.912 0.913 0.914 0.915 0.950 200 0.924 0.924 0.924 0.925 0.928 0.950 500 0.931 0.933 0.934 0.934 0.934 0.950 1000 0.936 0.936 0.939 0.939 0.941 0.950 2000 0.940 0.939 0.939 0.942 0.944 0.950 5000 0.945 0.945 0.944 0.944 0.946 0.950 Slope estimation 10 0.873 0.877 0.895 0.877 0.875 0.950 20 0.933 0.931 0.931 0.925 0.916 0.950 50 0.944 0.942 0.940 0.936 0.916 0.950 100 0.948 0.945 0.942 0.936 0.918 0.950 200 0.947 0.948 0.943 0.938 0.920 0.950 500 0.947 0.948 0.947 0.942 0.924 0.950 1000 0.948 0.947 0.952 0.945 0.925 0.950 2000 0.946 0.950 0.958 0.946 0.924 0.950 5000 0.950 0.963 0.960 0.945 0.923 0.950

a Standard errors of y for intercept and slope estimations are nominally 0.001.

a Standard errors of y for intercept and slope estimations are nominally 0.001.

values for slope estimation show an interesting behaviour: they do not vanish with increasing data size but rather approach a finite value. The reason is that the inaccurate prior knowledge about the measurement standard errors (nonzero A) persists to influence the slope estimation— an error source independent of the data size. Similar behaviours were found also for OLSBC estimation of intercept and slope (results not shown). The saturation value of RMSE^ depends on the accuracy of the prior knowledge (A), seemingly in a close-to-linear relation.

Table 8.6. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal shape: influence of accuracy of prior knowledge on RMSE. The experiment is the same as described in Table 8.5.

n RMSE

 A = 0.1 A = 0.3 A = 0.5 A = 0.7 A = 0.9 Intercept estimationa 10 0.310 0.313 0.315 0.329 0.348 20 0.217 0.219 0.220 0.227 0.241 50 0.137 0.137 0.139 0.143 0.151 100 0.096 0.097 0.098 0.100 0.106 200 0.068 0.069 0.069 0.071 0.075 500 0.043 0.043 0.044 0.045 0.048 1000 0.031 0.031 0.031 0.032 0.034 2000 0.022 0.022 0.022 0.023 0.024 5000 0.014 0.014 0.014 0.014 0.015 Slope estimation 10 0.279 0.300 0.428 0.303 0.339 20 0.173 0.177 0.181 0.191 0.223 50 0.105 0.107 0.114 0.126 0.166 100 0.074 0.077 0.085 0.099 0.146 200 0.052 0.056 0.066 0.084 0.135 500 0.033 0.040 0.052 0.073 0.127 1000 0.024 0.032 0.047 0.069 0.125 2000 0.018 0.028 0.043 0.067 0.123 5000 0.012 0.025 0.042 0.066 0.124

a Empirical RMSE^, given by b Empirical RMSE^ , given by a Empirical RMSE^, given by b Empirical RMSE^ , given by

To summarize, measurement error in the predictor requires to modify the OLS method to yield a bias-free slope estimation: OLSBC or WLSXY. These modified estimation methods require prior knowledge about the size of the measurement error. If this knowledge is not exact, which is a typical situation in the climatological practice, then it contributes to the estimation error of the slope (RMSE and CI length). This contribution persists even when the data size goes to infinity.

8.3.4 Mis-specified prior knowledge

What happens if we make a wrong specification of the accuracy of our prior knowledge? We study (Table 8.7) a situation where (1) the

Table 8.7. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal shape: influence of mis-specified prior knowledge on CI coverage performance. Design and estimation (WLSXY) are identical to the previous experiment (Table 8.5). CI construction (via pairwise-MBBres resampling) is identical to that in the previous experiment, with the following exceptions: (1) prior knowledge of Sx = 0.25, Sy = 0.5 (5 = 2.0) is incomplete after Eq. (8.11) with A = 0.5; (2) the adaptive Steps 8a and 8b of Algorithm 8.1 are allowed to mis-specify A.

Table 8.7. Monte Carlo experiment, linear errors-in-variables regression with AR(1) noise of normal shape: influence of mis-specified prior knowledge on CI coverage performance. Design and estimation (WLSXY) are identical to the previous experiment (Table 8.5). CI construction (via pairwise-MBBres resampling) is identical to that in the previous experiment, with the following exceptions: (1) prior knowledge of Sx = 0.25, Sy = 0.5 (5 = 2.0) is incomplete after Eq. (8.11) with A = 0.5; (2) the adaptive Steps 8a and 8b of Algorithm 8.1 are allowed to mis-specify A.

 n Y§0 Nominal True A = 0.5 True A = 0.5 Specified A Specified A 0.3 0.5 0.7 0.3 0.5 0.7 10 0.801 0.802 0.803 0.893 0.895 0.898 0.950 20 0.870 0.871 0.871 0.928 0.931 0.935 0.950 50 0.899 0.899 0.900 0.932 0.940 0.950 0.950 100 0.912 0.913 0.913 0.924 0.942 0.960 0.950 200 0.923 0.924 0.925 0.908 0.943 0.970 0.950 500 0.933 0.934 0.934 0.870 0.947 0.986 0.950 1000 0.939 0.939 0.940 0.827 0.952 0.994 0.950 2000 0.939 0.939 0.940 0.783 0.958 0.998 0.950 5000 0.944 0.944 0.944 0.744 0.960 1.000 0.950

a Standard errors of 7^ and are nominally 0.001.

a Standard errors of 7^ and are nominally 0.001.

true standard deviation ratio is 5 = SY/SX = 2.0, (2) the estimation on the sample is done with an incomplete knowledge of 5, modelled as a uniform distribution over the interval between 1.0 and 3.0 (A = 0.5), and (3) the bootstrap CI construction is allowed to mis-specify the incomplete knowledge by letting 5* be uniformly distributed over the intervals between 1.4 and 2.6 (specified A = 0.3) or between 0.6 and 3.4 (specified A = 0.7). Specifying A = 0.3 (instead of the correct A = 0.5) constitutes a case of overestimation of the accuracy of the prior knowledge, A = 0.7 means an underestimation and A = 0.5 is an unbiased estimation.

The first result is that such a mis-specification has no effect on the accuracy of CIs for the intercept. Table 8.7 displays results (for A = 0.3, 0.5 and 0.7) that are, within the bounds of the "simulation noise," indistinguishable.

The second result is that mis-specified prior knowledge has a clear effect on the accuracy of CIs for the slope. Table 8.7 shows results for A = 0.3 and 0.7 to deviate from those for the correct value of A = 0.5. If we underestimate the accuracy of the prior knowledge about the size of the measurement standard deviations (A = 0.3 instead of 0.5), then the

CIs become too narrow and the coverage is reduced; if we overestimate the accuracy (A = 0.7 instead of 0.5), then the CIs become too wide and the coverage is inflated.

8.4 Example: climate sensitivity

The effective climate sensitivity, denoted here as A-1, is a parameter that relates changes in annual-mean surface temperature to changes in the radiative forcing (greenhouse gases, etc.) of the climate system. Its units are °C (or K) per Wm-2. Climate sensitivity may vary with forcing history and climatic state, reflecting the influence of varying feedback mechanisms (amplifying or attenuating) in the climate system (Mitchell et al. 1987). The lack of an accurate knowledge of A-1 in the recent past (since, say, 1850) is one of the major obstacles for making accurate projections of future temperatures by means of AOGCMs (Forster et al. 2007).

The traditional estimation method for A-1 seems to be via perturbed climate models experiments, where the temperature response of the system is studied for a range of variations of model parameters and forcing scenarios (Forster et al. 2007). Due to the limited performance of climate models, it may be helpful to consider estimations that are based entirely on direct observations. We therefore relate variable Y(i), the observed temperature changes from 1850 to 2001, to variable X(i), the radiative forcing variations. The time series with standard errors are shown in Fig. 8.6. Since the predictor (forcing) has been determined with error, our model is the linear errors-in-variables regression (Eq. 8.1). The estimation objective is the slope, = A-1.

The result (Fig. 8.7) from WLSXY estimation of A-1 is 0.85 K W-1m2. The 95% CI, a Student's t interval obtained from pairwise-MBBres resampling with B = 2000 (Algorithm 8.1), is [0.47 K W-1m2; 1.24 K W-1m2]. The 90% CI, a level often used in the IPCC-WG I Report's chapter on radiative forcing (Forster et al. 2007), is [0.53 K W-1m2; 1.17 KW-1m2]. The climate literature often uses the "equilibrium climate sensitivity," which is defined as the temperature change that would be approached in a (hypothetical) equilibrium following a doubling of the atmospheric "equivalent carbon dioxide concentration" (representing all greenhouse gases). This other sensitivity value is around (4 W-1 m2)A-1, at least in the climate world of the E-R AOGCM of the National Aeronautics and Space Administration Goddard Institute for Space Studies, New York (Foster et al. 2008). Thus, the WLSXY result suggests that a CO2 doubling will lead to a temperature increase of 3.4 K.

What are the effects of autocorrelation? The block bootstrap resampling took into account the relatively strong memory of temperature and

Figure 8.6. Northern hemisphere temperature anomalies and climate forcing, 1850-2001: data. a The temperature time series, y(i), is shown (solid line) as deviation from the 1961-1990 average (n = 152). The annual-mean composite was derived using instrumental data from several thousand stations on land and sea (HadCRUT3 data set). The temperature standard error, sY(i), is shown (shaded band) as ±2sY(i) interval around y(i); it reflects following sources of uncertainty (Brohan et al. 2006): measurements, reporting, inhomogeneity correction, sampling, station coverage and bias correction of sea-surface temperatures. b The radiative forcing time series, x(i), is shown (solid line) with ±2sX (i) uncertainty band (shaded); it comprises following components thought to influence temperature changes (Hegerl et al. 2006; Forster et al. 2007): changes of atmospheric concentrations of greenhouse gases, solar activity variations (Fig. 2.12) and changes of sulfate and other aerosol constituents in the troposphere (lower part of the atmosphere). c Standard deviation ratio, S = sY (i)/sx (i). (Data from (a) Brohan et al. (2006) and (b) Hegerl et al. (2006).)

Year

Figure 8.6. Northern hemisphere temperature anomalies and climate forcing, 1850-2001: data. a The temperature time series, y(i), is shown (solid line) as deviation from the 1961-1990 average (n = 152). The annual-mean composite was derived using instrumental data from several thousand stations on land and sea (HadCRUT3 data set). The temperature standard error, sY(i), is shown (shaded band) as ±2sY(i) interval around y(i); it reflects following sources of uncertainty (Brohan et al. 2006): measurements, reporting, inhomogeneity correction, sampling, station coverage and bias correction of sea-surface temperatures. b The radiative forcing time series, x(i), is shown (solid line) with ±2sX (i) uncertainty band (shaded); it comprises following components thought to influence temperature changes (Hegerl et al. 2006; Forster et al. 2007): changes of atmospheric concentrations of greenhouse gases, solar activity variations (Fig. 2.12) and changes of sulfate and other aerosol constituents in the troposphere (lower part of the atmosphere). c Standard deviation ratio, S = sY (i)/sx (i). (Data from (a) Brohan et al. (2006) and (b) Hegerl et al. (2006).)

Figure 8.7. Northern hemisphere temperature anomalies and climate forcing, 1850-2001: fit. WLSXY estimation yields a straight regression line (solid) with a slope (i.e., effective climate sensitivity) of f31 = 0.85 KW-1m2. Also shown is OLS regression line (dashed).

Figure 8.7. Northern hemisphere temperature anomalies and climate forcing, 1850-2001: fit. WLSXY estimation yields a straight regression line (solid) with a slope (i.e., effective climate sensitivity) of f31 = 0.85 KW-1m2. Also shown is OLS regression line (dashed).

forcing noise components (a'X = a'Y = 0.82) by selecting a block length of l = 18. Ignoring autocorrelation (setting l = 1) would make the CI too narrow; for example, the 90% CI would become [0.56 K W-1m2; 1.14 K W-1m2].

The estimate and, more, the CI for A-1 should be assessed, however, with caution.

■ CI construction (Algorithm 8.1) used pairwise-MBBres resampling with an assumed constant standard deviation ratio of 5 = 0.66 (time-average). This was done because of the absence of Monte Carlo tests of adaptions of pairwise-MBBres resampling with respect to het-eroscedastic errors. Instead we imposed an uncertainty of 5 measured by the "incomplete prior knowledge" parameter A (Eq. 8.11). The employed value of A = 0.5 may have been too small and produced a too narrow CI. Particularly, unrecognized temperature variations not caused by measurement error or forcing changes, that is, "internal temperature variability," may let 5 increase and reduce the sensitivity estimate (Laepple T 2010, personal communication). Note that WLSXY estimation itself recognized heteroscedasticity.

■ The HadCRUT3 temperature data (Brohan et al. 2006) are down-biased between about 1940 and the mid-1960s because of an unrecognized change in 1945 in the sea-surface measurement techniques

(Thompson et al. 2008). Since this interval is short relative to the total observation interval, the influence of the inhomogeneity on the A-1 estimate should be small.

■ The tropospheric aerosol component of the forcing is known only with a "low" to "medium-low" scientific understanding (Forster et al. 2007). The aerosol contribution to X(i) and SX(i) may be large in error. Consequently, the error in Sx (i) and 5 may be large, and the parameter A may be larger than 0.5 (or even another model of the incomplete prior knowledge required). We stress that the large error of the predictor necessitates fitting an errors-in-variables regression model. Ignoring this error (i.e., using OLS estimation) would strongly underestimate the climate sensitivity (Fig. 8.7).

■ Volcanic eruptions, providing large negative forcing components (cooling) have been ignored in the estimation (because of the many unknowns), although the observed temperature time series (Fig. 8.6a) includes this effect. Since the number of large eruptions during the 151-year interval (Hegerl et al. 2006) is assessed as relatively small (about 8 eruptions with < -2.0 Wm-2 in the northern hemisphere), this omission should have a minor influence on the A-1 estimate.

■ Ocean heat uptake has similarly been ignored, although observed temperatures may show this influence. Assuming that it cannot be neglected would (1) increase the A-1 estimate and (2) widen its CI.

■ The analysis focused on the temperature of the northern hemisphere, while the concept of climate sensitivity applies to the globe. The superiority of temperature data quality for the northern part (more stations) suggested this restriction. Obviously, other geographic parts, including the globe, may be analysed in an analogous manner.

8.5 Prediction

A prediction is a statement about an uncertain event. In climate sciences the events lie often in the future (forecast) but frequently also in the past (hindcast), see the introductory examples (p. 3). In the context of the present chapter, we wish to predict an unobserved value, y(n + 1), given a sample, {t(i),x(i),y(i)}™=1, and a new observation, x(n + 1), of the predictor variable made at time t(n + 1).

A typical situation is when a relation between a climate variable, Y(i), and a proxy variable, X(i), is to be established. Suppose we observed {t(i),x(i),y(i)}™=1 over a time interval [t(1);t(n)] but have available a longer proxy time series, {t(i),x(i)}m=1 with m > n (often m » n). If t(i) denotes age and t(i) > t(n) for i > n, then we wish "to hindcast"

y(i) for i > n. An example is as precipitation proxy; y(i) is precipitation, x(i) is ¿18O from a speleothem, [t(1); t(n)] = [0a; 50a] is the interval for which we have instrumental measurements of y(i) (the past 50 years) and [t(1); t(m)] = [0 a; 500 a] is the interval covered by the speleothem samples (the past 500 years). If t(i) denotes time, then we wish to forecast. An example is climate model projections; y(i) is precipitation, x(i) is modelled precipitation (AOGCM), [t(1); t(n)] = [1950; 2010] is the interval for which we have instrumental measurements of y(i) and [t(1); t(m)] = [1950; 2100] is the interval analysed by means of the climate model (a typical value for the upper bound used by IPCC-WG I (Houghton et al. 2001; Solomon et al. 2007) in its reports).

Prediction can be performed by fitting a regression model and utilizing the estimated regression parameters. In the linear case (Fuller 1987: Section 1.6.3 therein):

where /0 and /?1 have been estimated using the sample {t(i), x(i), y(i) }n=1 and the new observation is x(n + 1).

Which method is suitable for estimating /0 and /31? Fuller (1987: pp. 75-76 therein) explains that usage of OLS, ignoring measurement errors of the predictor, is justified when x(n + 1) is drawn from the same distribution that generated {x(i)}™=1. This means effectively that two conditions have to be met:

1. SX(n + 1) ■ Xnoise(n + 1) has the same properties (range, shape, etc.) as Sx(i) ■ Xnoise(i), i = 1,..., n;

2. Xtrue(n + 1) has the same properties as Xtrue(i), i = 1,..., n.

Fuller advises further to take measurement error into account when not both conditions are satisfied. This can be done, for example, by using WLSXY (or OLSBC) estimation. Treating the regression estimates as if they were known parameters, Fuller (1987: p. 76 therein) gives the following expression for the prediction standard error:

We argue that in climatology the above conditions are almost exclusively not satisfied, and we advise to use WLSXY (or OLSBC) as a more conservative approach. In the majority of applications, t(n + 1), the time value related to the new measurement, is outside of [t(1); t(n)], and x(n +1) does not necessarily originate from a random drawing from the process Xtrue(i),i = 1,...,n. The new measurement may rather

0 0