Correlations between Surface Temperature Trends and Socioeconomic Activity: Toward a Causal Interpretation Short title: Surface Trends and Socioeconomic Activity: a Causal Interpretation Ross McKitrick Department of Economics University of Guelph Guelph ON Canada N1G 2W1 Tel. 519-824-4120 x52532 Fax. 519-763-8497 [email protected] and Nicolas Nierenberg Nierenberg Foundation REVISED September 6, 2009 to incorporate corrected Schmidt data REVISED January 25, 2010 to insert additional material in Sct. 4.1 and new Table 7 Abstract Evidence of surface data contamination and trend bias due to industrialization has recently been disputed on three main grounds: spatial autocorrelation of the temperature field undermines significance claims from regression models, substitution of modelgenerated data known to be free of contamination yields apparently significant effects, implying such results are spurious, and use of RSS lower tropospheric retrievals rather than UAH retrievals yields insignificant results. However these assertions have not been based on statistical tests. Within a formal testing framework we show that each of the arguments are unfounded. Spatial autocorrelation is not detected in the model residuals when observational data are used. It is present when modeled data are used, however, and correcting for it eliminates the apparent significance of the coefficients. In addition, the results in the analysis using modeled data predict a pattern opposite to that generated with observations, suggesting a causal explanation of the original findings. Use of RSS tropospheric data reduces but does not eliminate significance of the original results. We present evidence that the matrix of socioeconomic variables is essential for a wellspecified model of the surface trend field, further suggesting a causal relationship.

1

Research Supported by Social Sciences and Humanities Research Council of Canada Grant Number 430002 Key words: global warming, data quality, spatial autocorrelation, economic activity

2

Correlations between Surface Temperature Trends and Socioeconomic Activity: Toward a Causal Interpretation

1 1.1

Introduction Background

Empirical climatology relies on the assumption that surface temperature data over land have been adjusted to remove effects due to local non-climatic influences, such as population growth, urbanization, equipment changes, data quality problems in developing countries, variations in local air pollution levels, etc. This assumption is essential to a climatic interpretation of the data (see., e.g., Jun et. al. 2008, p. 935). McKitrick and Michaels (2004a,b, herein MM04, and 2007, herein MM07) tested the assumption by regressing the observed 1979-2002 trends in 440 surface grid cells on a vector of climatological variables (lower tropospheric temperature trends and fixed factors such as latitude, mean air pressure and coastal proximity) augmented with a vector of socioeconomic variables, including education, income and population growth, Gross Domestic Product (GDP) per square km, etc. They rejected, at very high significance levels, independence of the surface temperature field and the socioeconomic variables, thus concluding that the surface climatic data are influenced by the effects of industrialization on local temperature records. They estimated that the non-climatic effects could account for between one-third and one-half of the post-1979 average warming trend over land in the temperature data.

Schmidt (2009, herein S09) critiqued this interpretation on four grounds. First, he noted that warming is observed in numerous data sets. This is not under dispute, instead the accuracy of the land-based data is

3

our focus. Second, he argued that the surface temperature field exhibits spatial autocorrelation (SAC), which reduces the effective degrees of freedom in the sample and biases the test statistics towards overrejection of the null (no correlation) hypothesis. Third, he argued that use of the lower troposphere series from Remote Sensing Systems (Mears et al. 2003, denoted RSS) rather than the University of AlabamaHuntsville series (Spencer and Christy 1990, denoted UAH), reduces the significance of the coefficients, indicating a lack of robustness of the conclusions. Fourth, he argued that the results were spurious on the basis of a comparison with results obtained by swapping the observed surface and tropospheric trends with model-generated data from NASA’s Goddard Institute of Space Studies (GISS) model E, denoted herein as GISS-E. These model-generated data are by construction uncontaminated by industrializationinduced surface changes. Schmidt’s hypothesis was that if the GISS-E data yield the same regression coefficients as the observational data in MM07, it would indicate that the seeming correlations between patterns of warming and patterns of industrialization were a fluke. This is not what the S09 GISS-E runs showed however (as we explain below), but S09 also proposed a more general argument that if any significant correlations appeared, this would imply the results of MM07 were spurious.

Focusing on the latter three points we can show that they do not hold up on close scrutiny. In explaining why not, we present evidence of a structural mismatch between models and observed data that has hitherto escaped notice. Tropospheric trends as generated by GISS-E are not able to condition the observed surface trend field and leave an independent residual, whereas the MM07 regression model does. This implies that there is a structural match between the MM07 regression model and the observed trend field, which is missing from the GISS general circulation model. An additional implication is that SAC of the residuals emerges when the GISS-E data are substituted for observations, something

4

overlooked in the S09 analysis. Upon treating this problem the model-based coefficients reported in S09 become insignificant, as well as failing to match the signs and magnitudes of the observed values. Consequently the argument in S09 concerning the swapping in of GISS-E data not only fails, but points to a conclusion opposite to that expressed in S09. Specifically, a regression model that omits the surface process and inhomogeneity measures fails to explain the observed trend field, and generates a misspecified residual (i.e. one exhibiting significant SAC). But when the surface process and inhomogeneity measures are included in the regression model, the match to the observations becomes quite significant and the residual misspecification disappears. This implies a causal interpretation of the regression results in MM07, since both the geographical and socioeconomic data are necessary to yield a well-specified model of the process that generated the trend field, whereas use of data generated by a climate model run on the assumption that non-climatic biases have been removed from the observations does not.

We do find that use of RSS data rather than UAH data weakens the MM07 coefficients. However S09 failed to present the joint significance tests on which the main conclusions were based, and using RSS data these still uphold the MM07conclusions, albeit at reduced significance. For instance the null hypothesis of no-socio-economic effects rejects at P = 1.81 × 10 −6 using RSS data, rather than P = 9.92 × 10 −13 using UAH data, when the dependent variable is the CRUTEM2v trend field; also some individual coefficients lose size and significance.

In the next section we explain the data sets used throughout this paper. Then in Section 2 we examine the spatial autocorrelation topic, providing some theoretical intuition and then detailed results for the data

5

configurations of interest. In section 3 we explore the mismatch between the regression results from model-generated and observed data. Section 4 discusses the RSS-UAH contrasts, and Section 5 concludes.

1.2

Data sets

All data sets except the one discussed in Section 4.2 are taken from MM07 and S09.1 Readers should consult both these papers for detailed explanations; only a brief summary will be presented herein.

MM07 estimated the regression equation

θ i = β 0 + β 1TROPi + β 2 PRESS i + β 3 DRYi + β 4 DSLPi + β 5WATERi + β 6 ABSLATi + β 7 p i + β 8 mi + β 9 y i + β 10 ci + β 11 ei + β 12 g i + β 13 xi + u i

(1)

where θ i is the 1979-2002 trend in gridded surface climate data, TROPi is the time trend of Microwave Sounding Unit (MSU)-derived temperatures in the lower troposphere in the same grid cell as θ i over the same time interval, PRESS i is the mean sea level air pressure in grid cell i, DRYi is a dummy variable denoting when a grid cell is characterized by predominantly dry conditions (which is indicated by the mean dewpoint being below 0 oC). DSLPi is DRYi × PRESSi , WATERi is a dummy variable indicating the grid cell contains a major coastline, ABSLATi denotes the absolute latitude of the grid cell, pi is local population change from 1979 to 2002, mi is per capita income change from 1979 to 2002, yi is total

6

Gross Domestic Product (GDP) change from 1979 to 2002, ci is coal consumption change from 1979 to 2002, g i is GDP density (national Gross Domestic Product per square kilometer) as of 1979, ei is the average level of educational attainment, and xi is the number of missing months in the observed temperature series and u i is the regression residual. Equation (1) was estimated using the generalized least squares routine in Stata 8.0 with corrections for error clustering and heteroskedasticity.

Summary statistics are in Table 1. The MM07 data set has 440 records, one for every 5x5 degree grid cell over land for which adequate observations were available in the Climatic Research Unit (CRU) data archive to identify a trend over the 1979-2002 interval. Each record contains the linear surface trend expressed as degrees C per decade and the corresponding linear trend from the University of AlabamaHuntsville lower tropospheric record of Spencer and Christy (1990), denoted UAH.

The S09 data set comprises surface and tropospheric gridcell trends like those in MM07, except the surface trends are from later CRU compilations and the tropospheric trends are from Remote Sensing Systems (RSS) (Mears et al. 2003). S09 provides trends derived from the CRUTEM2v and CRUTEM3v data sets (Brohan et al. 2006). For brevity they are denoted herein as CRU2v and CRU3v. As is clear in Table 1 these data sets are very similar to one another. CRU3v is the most recent but has slightly lower geographical coverage compared to CRU2v.

1

Available respectively at http://www.uoguelph.ca/~rmckitri/research/jgr07/jgr07.html and http://www.giss.nasa.gov/staff/gschmidt/supp_data_Schmidt09.zip.

7

S09 also provided synthetic trends from GISS-E. For a description of this model see CCSP (2008 Sct. 2.5.3) and Schmidt et al. (2008). The climate model was run five times, and the mean over the five runs was taken as the ensemble average. The mean trends in the GISS-E surface data are denoted herein as GISSES, and the mean trends in the GISS-E lower troposphere data are denoted GISSET.

The average GISS-E surface trend is 0.14 oC/decade, well below the observed trend of 0.30 oC/decade in the CRU3v compilation. Applying cosine-latitude weighting of the gridcells the difference between GISS-E and CRU3v is 0.14 oC/decade. GISS model runs for a complete surface grid, taken from the coupled model version available at http://aom.giss.nasa.gov/lptime.html (applying both greenhouse gas and sulfate aerosol forcings) exhibits a trend of 0.32 oC/decade over the 1979-2002 interval, matching more closely the observed surface trends. It is not clear from S09 why the modeled surface warming trends are so much lower.

2 2.1

Spatial autocorrelation of the trend field Preliminary theory

S09 and Benestadt (2004) raised the issue of SAC as a weakness of the findings in MM07 and MM04. Both point out, correctly, that the surface temperature field is spatially autocorrelated, and argue that this can, in principle, bias the inferences from regressions on the spatial trend field. Neither author presents test statistics, though S09 presents variograms of the dependent and some independent variables from MM07. In this section we will confirm S09’s assertion that the trend field exhibits SAC. However, S09 was mistaken to conclude that the inferences in MM07 and related analyses are therefore biased. An

8

additional step in the argument is required, namely showing that the regression residuals also exhibit SAC. As we will show, they do not in the case of the observations used in MM07, but they do in the case of the model-generated data regressions in S09. It is important to interpret this contrast correctly.

A standard linear statistical model takes the form:

Random dependent variable = function of non-random explanatory variables + independent errors.

This is represented in a regression framework as

Random dependent variable = linear regression model + residuals.

Inferences concerning the coefficients in a regression model are based on the statistical properties of the residuals, not the dependent variable: the dependent variable does not even appear in the equation of the variance-covariance matrix estimator, since it is the residuals that contribute the stochastic component.

The importance of distinguishing between the dependent variable and the residuals is illustrated in Figure 1. We consider a case of time-series autocorrelation, which is easier to visualize graphically, but the point is the same. Figure 1 uses an artificial time series y consisting of a third-degree time polynomial plus independent random noise. In the top panel of Figure 1, y is modeled using a right hand side array denoted f(X), which in this case is a constant plus a linear term. The associated residuals e(X) are clearly

9

autocorrelated, indicating that f(X) is misspecified. Faced with this situation a researcher could do one of two things: either adopt an error-term model that takes the autocorrelation into account, or correct the model misspecification.

In time series analysis the former approach usually takes the form of an ARMA(p,q) residual model (e.g. Davidson and MacKinnon 2004 ch. 13), though in applied climatology a popular expedient is to apply an “effective degrees of freedom adjustment,” as mentioned in S09. An ARMA(p,q) model breaks the residuals down into two components:

residuals = structural component + independent errors.

The simplest example of a structural component is a lag-one autoregressive model. A properly-specified structural component will leave independent errors. The inferences for the coefficients of the regression model are then based on the properties of those errors, taking into account the degrees of freedom used up in the structural component.

The problem with this approach is that the structural component does not actually explain the dependent variable in the regression model, it simply isolates an independent residual for the purpose of generating asymptotically valid coefficient variances. That is, while the residual model yields asymptotically correct inferences on the coefficients in f(X), it does not correct the model misspecification. It would therefore be preferable to find an improved regression specification that leaves an independently-distributed residual. This would increase confidence that the model properly represents the actual process that

10

generated y. In Figure 1, bottom panel, the regression model is enhanced with higher-order terms, yielding a right hand side denoted f(Z), and an associated set of residuals e(Z) that are independently distributed. Note that y is still autocorrelated, but the residuals now are not, because the regression model is well-specified. Consequently, autocorrelation in y would not, in this case, bias the inferences on the coefficients in f(Z).

There are two points to emphasize. First, autocorrelation (including SAC) in a dependent variable does not bias t-statistics if the regression model explains it and leaves an independent residual. Second, when comparing two regression models with a spatially autocorrelated dependent variable, if one leaves an autocorrelated residual and the other leaves an independent residual, this could be interpreted as evidence that the latter model is a better specification of the underlying process that generated the dependent variable.

With these points in mind we can now turn to formal diagnostics of the models in MM07 and S09. A standard test for residual spatial dependence is implemented as follows. The regression model (1) can be rewritten in matrix notation as

T = Xb + u

(2)

where T is the 440x1 vector of linear trends in the temperature series for each of 440 surface grid cells, X is the 440xk matrix of climatic and socioeconomic covariates, b is the kx1 vector of least-squares slope

11

coefficients and u is the 440x1 residual vector. Spatial autocorrelation in the residual vector can be treated using

u = λ Wu + e

(3)

where λ is the autocorrelation coefficient, W is a symmetric n × n matrix of weights that measure the influence of each location on the other, and e is a vector of homoskedastic Gaussian disturbances, (Pisati 2001). n equals 440 except in some regressions where grid cells are missing, as noted below.

A test of H 0 : λ = 0 measures whether the error term in (1) is spatially dependent. Anselin et al. (1996) point out that if the alternative model allows for possible spatial dependence of T, i.e.

T = φZT + Xb + e

(4),

where Z is a matrix of spatial weights for T and may not be identical to W, then conventional tests of

λ = 0 assuming an alternative model of the form y = Xβ + e is a misspecification and will be severely biased towards over-rejection of the null. They derive a χ 2 (1) Lagrange Multiplier (LM) test of λ = 0

robust to possibly nonzero φ in (4), which has substantially superior performance in Monte Carlo evaluations compared to the non-robust LM test.

Hypothesis tests, and any subsequent parameter estimations, are conditional on the assumed form of the spatial weights matrix W in (3). We examine three forms herein. Denote the great circle distance

12

between the grid cell centers from which observation i and observation j are drawn as g ij . Inverse square root weights are 1 / g ij , inverse linear weights are 1 / g ij and inverse square weights are 1 / g ij2 . These

weights allow the relative influence of one cell on adjacent cells to decline by relatively slower to relatively faster rates, respectively.

In each case the rows of W are standardized to sum to one. We will show that evidence of SAC is absent across a wide range of specifications so we will not double the discussion by including estimations without row standardization. In general we found row standardization improves the speed of optimization but yields increased variance estimates, hence the results reported herein are conservative. In almost every case the conclusions are unaffected, however there is one case in which row standardization is an influential assumption: see Section 4.1 below.

2.2

Results

Table 2 presents the results of SAC hypothesis tests on both the dependent variable and residuals for six different model configurations, reporting for each of them the results based on three different spatial weighting schemes. A result common to all blocks is that the log-likelihood values are lowest for the inverse-root weighting scheme 1 / g ij and highest for the inverse-square weighting scheme 1 / g ij2 . This means that the data give progressively more support to the spatial weighting schemes that decay most quickly with distance.

13

The first block of results refers to the original configuration in MM07: the CRU gridded trends regressed on the UAH tropospheric trends and the rest of the MM07 model (1), and the second block shows the results with the dependent variable updated to CRU2v. Under the three spatial weighting schemes for both configurations we reject the null hypothesis of spatial independence in the dependent variable: indeed it is rejected in all cases shown in Table 2, amply confirming the conjecture in S09. But in all six cases, the corresponding hypothesis is not rejected for the residuals. In other words the regression model explains the SAC pattern and we find no evidence of spatial dependence in the residuals. Hence there is no evidence SAC biases the regression inferences reported in MM07.

The next two blocks report the results for CRU and CRU2v with the RSS data substituted in for UAH. Using the original CRU data the dependent variable exhibits SAC, as do the residuals under two of the three weighting schemes. But when the dependent variable is updated to CRU2v the SAC pattern in the residuals disappears. Thus while RSS appears to have a somewhat weaker ability to explain the SAC pattern in the surface data compared to UAH, on the more updated surface data set we have evidence of a spatially independent residual. In section 4 we will discuss the RSS results in more detail, including presenting SAC tests with the CRU3v data, in which one of the three weighting schemes yields evidence of SAC in the residuals.

Having established the general ability of tropospheric data to yield an uncorrelated residual we then asked whether the GISS-E-generated tropospheric trends were likewise capable. If there is a close structural match between the spatial variance in the model to that in the actual climate then we expect similar results to emerge, namely SAC in the dependent variable but not in the residuals. The next block

14

swaps GISSET for the lower troposphere data, using CRU2v as the dependent variable. Two interesting contrasts emerge. First, the log-likelihood values drop considerably. GISSET has a rather low correlation to the UAH (0.214) or RSS (0.349) tropospheric trends so this is not unexpected. Second, in two of the three specifications, including the one which obtains substantially higher support from the data (in the form of the highest log-likelihood value) the spatial lag term on the residuals is significant. Following from the graphical example above, a possible explanation is the absence of an essential component in GISS-E that replicates the spatial pattern of tropospheric warming, since use of observed data from UAH (and to a somewhat lesser extent RSS) data leaves an independent residual. In other words, something in the observed tropospheric trend field, but missing from the GCM-generated trend field, appears to be essential for explaining the surface trend field in such a way as to leave an uncorrelated residual.

Next we checked the residuals from a regression of CRU2v on, respectively, just RSS plus the geographical variables, the geographical and socioeconomic variables but not RSS, and all variables. The results are in Table 3. In all regressions the no-SAC null is always rejected for the dependent variable so we omit these P values and only show those for the test on the residuals. Under the inverse-square root weights the hypothesis of no SAC in the residuals is not rejected, but the log-likehood value for these weights is the lowest. Using the inverse-linear and inverse-squared weights the no-SAC hypothesis is rejected when either RSS or the socioeconomic variables are omitted, but not rejected when both blocks of variables are included. Hence, only when both the tropospheric trends and the socioeconomic variables are included does SAC disappear from the residuals. This suggests both the tropospheric trend field and the proxy measures of surface data contamination are essential for a well-specified model of the surface trend field. As noted above, use of the GISS-E data as a replacement for tropospheric

15

observations does not suffice to leave an independent residual. The atmospheric data alone (observed or modeled) do not yield a well-specified model—the socioeconomic variables are also necessary. This strengthens the view that there is a causal relationship between the complete set of MM07 variables and the surface trend field.

In the final block of Table 2, GISSES is regressed on GISSET and the other MM07 variables. Here the SAC coefficient is significant in all three weighting schemes on both the dependent variable and the model residuals. The use of modeled data on both the right and left-hand sides fails to leave an uncorrelated residual, in contrast to the case where observational trends are used. This possibly points to a deficiency of the coupling between the surface and troposphere in the GCM. More importantly it indicates that, in the regressions that use modeled data to check for spurious correlations, a correction for SAC is required, and we will show in the next section that this makes a difference.

3 Do GCMS predict the temperature-industrialization correlation pattern? S09 hypothesized that the surface trend field under natural and greenhouse forcing may exhibit a spurious, or fortuitous, match to the pattern of socioeconomic trends, due to natural processes. Since the GCM does not contain a socioeconomic component, if, upon using GISSES and GISSET in place of observations in the MM07 regression model, significant coefficients, of the same approximate size and sign, emerge on the socioeconomic variables, then correlations such as those in MM07 would obviously be coincidental.

16

But this argument runs into two serious problems. First, as shown in Table 2, the regression model using the ensemble mean GISS data yields spatially autocorrelated residuals, requiring an SAC-corrected estimator, which S09 did not apply. When such an estimator is applied, all the socioeconomic coefficients become insignificant. This is shown in the final column of Table 4. While the surface process variables are jointly significant (P = 0.038) the inhomogeneity variables are not (P = 0.991) nor are the socioeconomic variables taken as a complete group (P = 0.075).

Second, even if the coefficients had been significant, the coefficient signs and magnitudes on the modeled data are completely unlike those on the observed data. Table 4 shows that the magnitudes are typically one-tenth those of the observed coefficients. And four of the six socioeconomic coefficients in the GISS-E regression take the opposite sign to those estimated on the observations (denoted with the superscript a.) If the model predicts one pattern and the data exhibit a very different pattern, significance of either set of coefficients increases, rather than decreases, the evidence that the observations contain a pattern at odds with the models. It is worth quoting the argument in S09 to make this point clear.

“There is a relatively easy way to assess whether there is any true significance to these correlations. We can take fully consistent model simulations for the same period and calculate the distribution of the analogous correlations. Those simulations contain no unaccounted-for processes (by definition!) but plenty of internal variability, locally important forcings and spatial correlation. If the distribution encompasses the observed correlations, then the null hypothesis (that there is no contamination) cannot be rejected.

17

(S09, p. 2, emphasis added)

S09 itself shows that the distribution does not encompass the observed correlations. Table 5 shows the 95% confidence interval of the ensemble mean coefficients using the (revised) S09 data with SAC correction applied, and presents the coefficients from the MM07 model estimated using UAH and RSS tropospheric data, which clearly lie outside the range of modeled coefficients. A chi-squared test that there is no difference among the coefficient groups is shown in the bottom row. The test clearly rejects using both the RSS and UAH data sets.

The fact that the distribution does not encompass the observed correlations refutes the stated basis for the conclusion in S09 that the effects are spurious. Indeed the rejection of the hypothesis strengthens the argument that there is a non-climatic pattern in the data, since a model without the contamination processes cannot replicate the observations. The spatial autocorrelation results in the previous section add to the evidence that the GISS model does not capture a structural element in the surface trend field, whereas a model of the surface data contamination processes does.

Further evidence on this point is obtained by repeating the filtering experiment of MM07 on the GISS data. Since the model does not contain any contaminating processes, we should not expect much difference between the raw data from the models and that obtained by applying the MM07 method for removing the socioeconomic effects. Table 6 shows that this is, indeed, the case. In the MM07 case the filtering removed about one-half the surface trend. When RSS data (with CRU3v) are used, the filtering removes about one-third of the mean surface trend. But when GISS data are used, the filtering step does

18

not reduce the mean surface trend at all, indeed, it rises slightly, and ends up midway between the filtered surface values from the models using UAH and RSS data, respectively.

4 4.1

RSS- versus UAH-based results Direct comparison

The results in Table 4 show that the MM07 parameter values remain reasonably stable across the four different combinations of observations, though the individual coefficient results certainly weaken when RSS data are used in place of UAH. This is not so much a test of MM07 as it is a test of similarity between UAH and RSS, since we are not in a position to say whether RSS or UAH is the better data set for this purpose. The two satellite-based trend fields have a correlation of 0.76, which is significant, but surprisingly low considering they come from the same underlying microwave soundings.

The ideal situation would be one in which conclusions did not depend on the RSS-UAH choice. This is the case for the joint significance tests, which remain highly significant with either data set. S09 only reported individual coefficient values, rather than the joint F tests on which the main MM07 conclusions rested. After correcting for the clustering structure, the joint F tests across the parameter groups remain highly significant when using RSS with either CRU2v or CRU3v as the dependent variable. Table 4 lists the results of the different configurations of surface and tropospheric trend fields. The three rows denoted P(I), P(S) and P(all) denote, respectively, the joint significance of the inhomogeneity effects (testing g, e, x all zero), the joint significance of the surface process effects (testing p, m, y, c all zero), and the joint significance of all socioeconomic variables. In every case the hypothesis of no socioeconomic effects is rejected and the prob value is typically of magnitude 10-6 or smaller. In addition, as shown in Table 6,

19

use of the MM07 filtering method and RSS data still yields a substantial drop in the mean trend over land, as does UAH data, whereas using GISS model data does not.

One configuration, however, does provide a case against the results in MM07. Using the CRU data as applied in MM07 as the dependent variable, i.e. prior to both the CRU2v and CRU3v updates, and using RSS data for the troposphere, the robust Lagrange multiplier shown in Table 2 implies a correction for SAC is needed. If it is applied with row standardization in W, the coefficients become jointly insignificant. P(S) rises to 0.358 and P(all) rises to 0.111. This result reverses, however, when row standardization is not applied. In that case (not shown) the robust LM test indicates SAC needs to be corrected, and doing so (with inverse square weights) the coefficients on g, e, m, y and c are individually significant, and the three joint tests strongly reject: P(I) = 5.1× 10 −9 , P(S) = 0.006 and P(all) = 3.2 × 10 −8 . The results also become significant when we reconcile the grid cell dimensions, as discussed in the next section. In that case the SAC disappears from the residuals and the no-contamination hypothesis is once again strongly rejected. Consequently, the one case we found in which the socioeconomic effects could be considered spurious depended on use of an older version of the dependent variable data and was not robust to minor re-specification of the troposphere measure or the SAC weight matrix.

MM07 presented a series of specification tests using the UAH data. We have repeated all of them using CRU3v and RSS, without SAC terms. The Anselin test does not reject the null of no SAC for inverse square root and inverse square weights, though it does for inverse linear weights (the inverse square weights are maximum likelihood, as usual). The results do not depart in any meaningful way from those reported in MM07. Details are available on request, and can be summarized as follows.

20

• The RESET test does not reject a null hypothesis of no un-modeled residual nonlinearity (P = 0.241). • The Hausman test does not reject a null hypothesis of no endogeneity bias (P = 0.963). • The outlier test flags 27 observations as influential. When these are removed the individual and joint

socioeconomic coefficient tests become more significant, yet we do not reject a null hypothesis that the coefficient vectors with and without the outliers are equivalent. • Coefficient results are individually and jointly significant in rich countries but not poor countries,

and in economies with growing but not declining incomes. • After removing a randomly-selected third of the data set and re-estimating the model, the prediction

of the withheld sample scatters along a 45-degree line with the observed values. In 500 repetitions, a regression of the predicted and observed values has a constant of 0.009 and a slope of 0.969, and a test of a perfect fit (constant = 0, slope = 1) obtains an average P value of 0.406, i.e. does not reject on average.

The SAC tests reveal that one of the specification tests in MM07 was done incorrectly. In Section 4.6 of MM07, an alternative estimation is presented in which the surface trends were replaced by the UAHderived tropospheric trends. Had the socioeconomic coefficients retained their size and significance it would provide evidence that the surface results are spurious. As shown in Table 3 of MM07, the socioeconomic coefficients generally lose size and significance, as expected, although an anomalous result emerges whereby the missing variable count in surface data (denoted x) becomes significant. Only 5% of the cells in the sample have at least one missing month. The analysis in MM07 suggests that the

21

occurrence of missing data is possibly acting as a proxy for relatively moist, storm-prone regions, but in any case x is small and insignificant in the main models so its role is unlikely to bias the conclusions.

However, when the tropospheric trends vector is used as the dependent variable, SAC does not disappear in the model residuals, so the regression equation needs to be augmented with a spatial error-correction process. This was not done in MM07. Of the weighting schemes used herein, the inverse-square weights yield much higher likelihood values than the alternatives. Table 7 shows the results of replacing the dependent variable with the tropospheric trends and regressing on the remaining right-hand side variables. The first column reproduces the original (incorrect) results from MM07, and the next two columns show the results after applying the SAC correction.

As in MM07, x acquires significance, indicating it is likely acting as a proxy for some other spatial pattern. In the full regression x is not significant so this effect is spurious, as in MM07. Regarding the other coefficients, there is a plausible attenuation of the surface effects pattern, though one conspicuous anomaly is the e variable when the UAH trends are in the model. The coefficient falls to one quarter its original size and changes sign, but becomes significant. The other coefficients diminish in size and in several cases significance. The population growth measure is smaller but still significant in the RSS record, but insignificant in the UAH record. The joint significance tests, however, do represent evidence against our interpretation of this model, since they yield joint significance for the inhomogeneity and surface process effects taken together. If we are truly measuring socioeconomic influences on surface temperatures, except for some land use-related measures that might be detectable in the troposphere these should be sufficiently attenuated in the lower troposphere as to leave insignificant coefficients.

22

However the lambda coefficients (representing the lag-1 spatial autocorrelation coefficients) indicate that our treatment of the residual dependency is likely inadequate, so judgments about significance in these results should be considered tentative. Each lambda coefficient is greater than 0.99 and has a very large tstatistic. In time-series applications this typically suggests that a lag-1 process is inadequate, and in the present context the results suggest that a spatial analogue to a higher-order autocorrelation process may be needed to obtain unbiased variances.

The results in Table 7 are close enough to our expectations that we will proceed to present our interpretations of the rest of our analysis with respect to the complete model estimations and the comparison of modeled versus observational data. However they do raise the caution that we cannot decisively rule out the influence of some kind of spurious spatial match between climatic variations and socioeconomic development. To provide a decisive test will require development of an updated data base that combines both time series and cross-sectional data on a gridcell basis. This is planned as a subsequent development of this line of inquiry.

4.2

Reconciling gridcell sizes

One potential deficiency of the MM07 data construction is that the satellite data are from 2.5x2.5o grid cells whereas the surface data use 5x5o grid cells. The results in MM07 were based on comparing the topright satellite gridcell to the larger surface gridcell, but for completeness we took the four corresponding satellite grid cells and averaged them up to the matching 5x5o level. The results are in Table 8. Robust

23

LM tests did not show SAC in the residuals. In particular, when using CRU as the dependent variable and UAH as the tropospheric record (on a 5x5o basis) the robust LM statistic (inverse-square weights) has a P value of 0.795, and 0.548 with RSS, thus providing an important contrast with the above-mentioned results on CRU and RSS without reconciling the gridcell size. If the dependent variable is CRU2v the values are 0.880 and 0.743, respectively, and using CRU3v they are 0.821 and 0.691 respectively. Therefore a correction for SAC is not indicated. The coefficient estimates are generally stable. The socioeconomic parameters are somewhat smaller and, as before, the RSS results are weaker than the UAH results. But in both cases the joint parameter estimates are highly significant. For instance, using RSS and CRU2v the hypothesis of no contamination is rejected with a P value of 3.31 × 10 −8 .

5

Conclusions

We have re-examined the question of surface temperature data contamination in light of the arguments of S09 against MM07. S09, similar to Benestad (2004), argued that spatial autocorrelation of the temperature field reduces the number of actual degrees of freedom and biases significance calculations. But neither source estimated an actual SAC model nor tested the regression residuals, as opposed to the dependent variables. We have done so and have shown that in virtually all cases, under a range of spatial weighting rules, SAC is not found to be present in the residuals. Including both the tropospheric trends and the socioeconomic covariates is essential for explaining the structure of the dependent variable and leaving an independent residual, suggesting a causal relationship between the right hand side variables and the surface trend data.

24

S09 also argued that swapping in model-generated data for observations yielded significant results, which are perforce spurious, thereby casting doubts on the conclusions of MM07. However, S09 ignored the fact that when model-generated data are used in place of observations, SAC appears in the residuals. This implies that the GISS-E estimations in S09 should have included a correction for SAC, unlike the observational data, and when this is done the significance of the coefficients disappears, overturning the S09 argument. It also implies that the model-generated data lack the key structural similarity to observations essential for removing the SAC from the error terms. Further evidence of this arises from a comparison of the regression coefficients calculated using the GISS-E data versus those calculated using observations. The distributions do not overlap and we reject at high significance the hypothesis that the pairs of coefficient vectors are identical. Hence we find no support for the view advanced in S09 that the results on observational data are spurious—indeed the data presented in S09 strengthens the argument that they are not.

As noted in S09, substituting RSS for UAH data does weaken the main regression results. We cannot infer from this that the results are therefore compromised, without presuming to decide between alternative data products. Irrespective of choice of data product, however, the joint significance tests remain very significant and in each case uphold the MM07 inference that the surface trend field is strongly affected by industrialization and other forms of socioeconomic changes. Depending on whether the RSS or UAH satellite data are used the filtered data yields a surface trend between one-third and onehalf smaller than the reported average over land.

25

The data set presented in MM07 includes trends up to the end of 2002, and includes coarse resolution of some socioeconomic variables at the national level. Further investigation of the potential surface climatic data problems we have identified herein could involve a reconstruction of the MM07 data base using updated socioeconomic and climatic variables, use of cross-sectional time series (panel) regression rather than trend fields, and use of regional, rather than national, socioeconomic data where available.

26

References Anselin, L, Anil K. Bera, Raymond Florax and Mann J. Yoon (1996). Simple diagnostic tests for spatial dependence. Regional Science and Urban Economics 26: 77—104. Benestad RE (2004) Are temperature trends affected by economic activity? Comment on McKitrick & Michaels (2004). Clim. Res. 27:171–173 Brohan, P., J.J. Kennedy, I. Harris, S.F.B. Tett and P.D. Jones, 2006: Uncertainty estimates in regional and global observed temperature changes: a new dataset from 1850. J. Geophys. Res. 111, D12106, doi:10.1029/2005JD006548 Davidson, R. and J.G. MacKinnon (2004), Econometric Theory and Methods Toronto: Oxford. Jun, Mikyoung, Reto Knutti and Douglas W. Nychka (2008) Spatial analysis to quantify numerical model bias and dependence: How many climate models are there? Journal of the American Statistical Association 108 No. 483 934—947 DOI 10.1198/016214507000001265. McKitrick, R.R. and P. J. Michaels (2004a), A test of corrections for extraneous signals in gridded surface temperature data, Climate Research 26(2) pp. 159-173, Erratum, Clim. Res. 27(3) 265—268. McKitrick, R.R. and P. J. Michaels (2004b), Erratum, Clim. Res. 27(3) 265—268. McKitrick, R.R. and P.J. Michaels (2007), Quantifying the influence of anthropogenic surface processes and inhomogeneities on gridded global climate data, J. Geophys. Res., 112, D24S09, doi:10.1029/2007JD008465. Mears, Carl A., Matthias C. Schabel, and Frank J. Wentz. 2003. A Reanalysis of the MSU Channel 2 Tropospheric Temperature Record. Journal of Climate 16, no. 22 (November 1): 3650-3664 . Pisati, Maurizio (2001) Tools for spatial data analysis. Stata Technical Bulletin STB-60, March 2001, 21—37.

27

Schmidt, Gavin (2009) Spurious correlation between recent warming and indices of local economic activity. International Journal of Climatology 10.1002/joc.1831 Spencer, R.W. and J.C. Christy (1990), Precise monitoring of global temperature trends from satellites, Science 247:1558—1562.

28

Figure 1. Top panel: regression model f(X) fitted to data series y, where f(X) consists of linear trend. Associated residuals e(X) are autocorrelated. Bottom panel: regression model f(Z) fitted to data series y, where f(Z) includes higher-order terms. Associated residuals e(Z) are not autocorrelated.

29

Var CRU

Definition Surface temperature trend from MM07

CRU2v

CRUTEM version 2 trends from S09 CRUTEM version 3 trends from S09 Surface gridcell trend from GISS model (S09) Lower tropospheric gridcell trend from RSS Lower tropospheric gridcell trend from UAH Lower tropospheric gridcell trend from GISS model (S09)

CRU3v GISS-ES RSS UAH GISS-ET

Water Abslat g e x p m y c Rich Grow

Grid cell contains a coast line Absolute latitude 1979 Real National GDP per sq km in millions Literacy +Post-secondary education rates # missing months in grid cell temperature record % growth in population* % growth in real average income* % growth in real national GDP** % growth in coal consumption* 1999 real income > median 1999 real income > 1979 real income

Obs

Mean

Std. Dev.

Min

Max

440 440

0.302 0.296

0.257 0.250

-0.700 -0.699

1.020 1.015

428

0.303

0.253

-0.717

1.042

440

0.141

0.081

-0.091

0.689

434

0.237

0.134

-0.085

0.684

440

0.232

0.184

-0.197

0.683

440

0.160

0.055

-0.016

0.329

440 440

0.6045 40.602

0.4895 17.953

0 2.5

1 82.5

440

0.297

0.600

0.001

3.002

440

106.5

26.20

11.6

144.2

440 440

0.764 0.279

2.552 0.209

0 -0.069

24 1.235

440

0.380

0.614

-0.790

2.147

440

0.771

0.839

-0.669

3.003

440

1.016

4.056

-1

39.333

440

0.493

0.501

0

1

440

0.761

0.427

0

1

30

Table 1: Model Variables. Definitions discussed further in MM07 and S09. *over the interval 1979 to 1999. **Over the interval 1980 to 2000. % Changes should be multiplied by 100, e.g. mean population growth is 27.9%.

31

Estimation Model CRU: UAH + MM07

CRU2v: UAH + MM07

CRU: RSS + MM07

CRU2v: RSS + MM07

CRU2v: GISSET + MM07

GISSES: GISSET + MM07

Weighting Scheme 1/ g ij

Dependent Variable 5.599 (0.018)

Residuals 2.564 (0.109)

LLF 139.250

1/ g ij

7.369 (0.007)

0.032 (0.858)

141.452

1/ g ij2

15.374 (0.000)

0.094 (0.759)

147.844

1/ g ij

7.560 (0.006)

3.411 (0.065)

151.796

1/ g ij

9.318 (0.002)

0.001 (0.974)

153.959

1/ g ij2

17.555 (0.000)

0.493 (0.482)

159.351

1/ g ij

5.831 (0.016)

1.030 (0.310)

132.875

1/ g ij

6.277 (0.012)

4.188 (0.041)

138.335

1/ g ij2

10.794 (0.001)

3.827 (0.050)

160.818

1/ g ij

7.652 (0.006)

1.125 (0.194)

146.139

1/ g ij

7.940 (0.005)

2.983 (0.084)

151.332

1/ g ij2

10.794 (0.001)

2.455 (0.117)

162.541

1/ g ij

14.274 (0.001)

0.922 (0.337)

93.755

1/ g ij

14.468 (0.000)

18.518 (0.000)

106.913

1/ g ij2

15.215 (0.000)

7.237 (0.007)

137.006

1/ g ij

43.108 (0.000)

39.232 (0.000)

587.021

1/ g ij

48.528 (0.000)

84.522 (0.000)

597.260

1/ g ij2

51.478 (0.000)

139.484 (0.000)

785.977

TABLE 2. Spatial autocorrelation tests for 4 regression models. Estimation model described as [surface measure]: [tropospheric measure] + MM07. Surface measure is CRU gridded surface trend vector from McKitrick and Michaels (2007), showing 1979-2002 trend per grid cell, or CRU2v update from Schmidt (2009). Tropospheric measure is either UAH or RSS observational trends, or GISSET, denoting ensemble mean tropospheric trend from Schmidt (2009)). MM07 denotes other dependent variables from McKitrick and Michaels (2007) model, namely SLP through c (see text). Weighting scheme refers to assumed form of spatial dependence. Third and fourth columns: each entry shows estimated autocorrelation parameter and associated prob value of hypothesis that it equals zero. Bold denotes hypothesis rejected at 5% significance. Fifth column: log-likelihood values.

32

Square root weights

P(Robust λ =0) Linear weights

Squared weights

RSS+geog

0.809

0.000

0.001

geog+econ

0.309

0.000

0.006

SAC weights:

RSS+geog+econ 0.194 0.084 0.117 TABLE 3. Spatial autocorrelation tests for CRU2v and RSS configurations. ‘geog’ denotes variable group including mean pressure, absolute latitude, dryness indicator, product of dryness x latitude, and coastline indicator. ‘econ’ denotes variable group g, e, x, p, m, y and c (see Table 1)

33

Variable Trop slp dry dslp water abslat

g e x p m y c Constant

P(I) P(S) P(all)

MM07 0.8631 (8.62) 0.0044 (1.02) 0.5704 (0.10) -0.0005 (-0.09) -0.0289 (-1.37) 0.0006 (0.51)

C2/UAH 0.8323 (8.58) 0.0043 (1.03) 1.6771 (0.33) -0.0016 (-0.31) -0.0293 (-1.37) 0.0010 (0.82)

C2/RSS 0.9872 (14.36) 0.0024 (0.66) 0.1986 (0.04) -0.0001 (-0.03) -0.0169 (-0.68) 0.0028 (2.28)

C3/UAH 0.8146 (8.80) 0.0058 (1.43) 1.5010 (0.30) -0.0014 (-0.28) -0.0240 (-1.06) 0.0014 (1.22)

C3/RSS 0.9627 (13.84) 0.0039 (1.11) 0.2125 (0.04) -0.0001 (-0.03) -0.0117 (-0.46) 0.0033 (2.65)

GISS-E 1.562 (13.43) -0.0042a (-2.55) -0.5264a (-0.31) 0.0005a (0.30) -0.0288 (-3.98) 0.0016 (3.64)

0.0432 (3.36) -0.0027 (-5.14) 0.0041 (1.66) 0.3839 (2.72) 0.4093 (2.39) -0.3047 (-2.22) 0.0062 (3.45) -4.2081 (-0.96)

0.0450 (3.66) -0.0026 (-5.32) 0.0019 (0.73) 0.3665 (2.67) 0.3844 (2.33) -0.2839 (-2.15) 0.0060 (3.46) -4.1522 (-0.97)

0.0444 (3.03) -0.0025 (-4.73) -0.0003 (-0.12) 0.1513 (1.04) 0.2663 (1.52) -0.2160 (-1.55) 0.0076 (3.77) -2.2291 (-0.60)

0.0449 (4.01) -0.0029 (-4.41) 0.0011 (0.42) 0.3524 (2.52) 0.3732 (2.19) -0.2804 (-2.06) 0.0063 (3.38) -5.5546 (-1.36)

0.0446 (3.27) -0.0028 (-4.10) -0.0021 (-0.95) 0.1450 (0.97) 0.2578 (1.42) -0.2139 (-1.48) 0.0079 (3.62) -3.7704 (-1.05)

-0.0009a (-0.13) -0.0000 (0.01) -0.0795a (-1.55) -0.0798a (-1.45) 0.0452a (1.13) 0.0006 (1.07) 4.1237a (2.45)

0.000 0.000 0.000

0.000 0.000 0.000

0.000 0.005 0.000

0.000 0.001 0.000

0.000 0.007 0.000

0.991 0.038 0.075

440 440 434 428 422 440 N 0.53 0.53 0.52 0.53 0.52 0.65 R2 139.22 151.75 145.50 142.64 135.83 640.842 Log-likelihood Table 4 Regression results from a group of data configurations. MM07 shows results from original data set (McKitrick and Michaels 2007). Next four columns show combinations of CRUv2 and CRUv3 surface data sets (denoted C2 and C3 respectively) and the UAH and RSS lower tropospheric data series. Final column shows regression using GISS-E-generated surface and tropospheric data. Numbers in parentheses are t-statistics. For Columns 2-6 these reflect corrections for clustered errors and

34

heteroskedasticity. For Column 7 a correction for spatial autocorrelation (inverse-square weights) is also applied. Bold denotes significant at 95% confidence. a indicates the sign of the GISS-E coefficient is opposite to that in all the observation-based regression. P(I) is prob value of test that all inhomogeneity factors (g–x) are jointly zero; P(S) = prob value of test that all surface process coefficients (p–c) are jointly zero; P(all), prob value of test that g –c are jointly zero. R2 in final column is squared correlation between observed and predicted regression values. Variable x is dropped in GISS-E regression since there are no missing surface values.

35

MM07 MM07 Variable GISS range using UAH using RSS g -0.0136 to -0.0119 0.0432 0.0445 e -0.0005 to 0.0005 -0.0027 -0.0025 p -0.1803 to 0.0213 0.3839 0.1514 m -0.1879 to 0.0282 0.4093 0.2663 y -0.0331 to 0.1235 -0.3048 -0.2160 c -0.0005 to 0.0017 0.0062 0.0076 H0:(results same) 0.0000 0.0000 Table 5. Comparison of coefficient magnitudes between GISS runs and observational model runs. Last two columns show results of MM07 regression using CRU surface data and UAH or RSS satellite data, respectively. The bottom row shows the P-value of a chi-squared test of equivalence between the GISS ensemble mean results and the observational results. Note that in none of the cases do the observed coefficients fall in the corresponding model range and the equivalence test clearly rejects.

36

MM07 method using MM07 method using MM07 method using Mean Trend: CRU & UAH CRU3v & RSS GISS data Surface 0.274 0.277 0.193 Troposphere 0.203 0.228 0.225 Filtered surface 0.128 0.185 0.223 Table 6: Filtering results using UAH/RSS/GISS. Each table entry shows the mean trend in the global sample with gridcells weighted by the cosine of latitude. The second column shows the original MM07 results. The third column shows the results from swapping in CRU3v and RSS data and the fourth column shows the results from swapping in the GISS ensemble mean simulations from S09.

37

ORIG (UAH)

RSS+SAC

UAH+SAC

slp

0.0440 1.02

0.0054 2.14

0.0053 1.85

dry

0.5700 0.1

-4.2987 -1.56

-6.0400 -1.79

dslp

-0.0010 0.09

0.0043 1.57

0.0059 1.79

water

-0.0290 1.37

0.0083 0.89

0.0145 1.33

abslat

0.0006 0.51

0.0027 4.05

0.0048 5.73

g

0.0432 3.36

0.0080 0.72

0.0115 0.87

e

-0.0027 5.14

-0.0001 -0.17

0.0007 2.13

x

0.0040 1.66

-0.0024 -2.05

-0.0048 -2.23

p

0.3839 2.72

0.1969 3.45

0.0997 1.63

m

0.4093 2.39

0.1231 1.76

0.0531 0.68

y

-0.3047 2.22

-0.0752 -1.38

-0.0284 -0.48

c

0.0062 3.45 -4.2080 0.96

0.0006 0.65 -6.1071 -2.23 0.9931 146.28 434 533.58

0.0017 1.35 -6.8253 -2.03 0.9938 161.6 434 449.95

Const Lambda N LL

440 139.22

38

P(I) 0.185 0.003 0.010 P(S) 0.211 0.063 0.000 P(All) 0.010 0.000 0.004 Table 7: Results from replacing the dependent variable with the tropospheric trends and regressing on the remaining variables. First column: from MM07, UAH data, no correction for SAC. 2nd and 3rd columns: RSS and UAH data respectively, SAC correction applied. t statistics underneath coefficients, bold denotes significant at 5%.

39

Variable Trop slp dry dslp water abslat

g e x p m y c Constant

P(I) P(S) P(all)

MM07 0.8631 (8.62) 0.0044 (1.02) 0.5704 (0.10) -0.0005 (-0.09) -0.0289 (-1.37) 0.0006 (0.51)

CRU2/U4 0.9440 (10.75) 0.0056 (1.47) 4.4913 (0.98) -0.0043 (-0.96) -0.0330 (-1.62) 0.0002 (0.20)

CRU2/R4 0.9638 (9.06) 0.0045 (1.22) 4.3373 (0.99) -0.0042 (-0.98) -0.0216 (-0.93) 0.0038 (2.82)

CRU3/U4 0.9330 (11.03) 0.0071 (1.92) 4.5837 (1.01) -0.0044 (-0.99) -0.0278 (-1.29) 0.0005 (0.54)

CRU3/R4 0.9538 (8.77) 0.0059 (1.67) 4.3175 (1.01) -0.0042 (-0.99) -0.0161 (-0.66) 0.0041 (3.00)

0.0432 (3.36) -0.0027 (-5.14) 0.0041 (1.66) 0.3839 (2.72) 0.4093 (2.39) -0.3047 (-2.22) 0.0062 (3.45) -4.2081 (-0.96)

0.0387 (3.27) -0.0025 (-5.30) 0.0020 (0.94) 0.3092 (2.37) 0.2954 (1.86) -0.2245 (-1.76) 0.0056 (3.08) -5.4725 (-1.40)

0.0409 (2.36) -0.0022 (-4.16) 0.0015 (0.56) 0.2206 (1.64) 0.2735 (1.66) -0.2003 (-1.53) 0.0069 (3.74) -4.4279 (-1.19)

0.0383 (3.44) -0.0028 (-4.32) 0.0016 (0.71) 0.2916 (2.18) 0.2810 (1.68) -0.2178 (-1.63) 0.0057 (2.93) -6.9317 (-1.84)

0.0401 (2.58) -0.0025 (-3.73) 0.0005 (0.20) 0.2085 (1.55) 0.2595 (1.50) -0.1944 (-1.42) 0.0071 (3.58) -5.8808 (-1.63)

0.000 0.000 0.000

0.000 0.000 0.000

0.000 0.000 0.000

0.000 0.008 0.000

0.000 0.001 0.000

440 440 440 428 428 N 0.53 0.56 0.54 0.56 0.54 R2 139.22 166.18 155.33 157.10 147.54 Loglikelihood Table 8: Results as for Table 4 after reconciling grid cell size. For notation see notes to Table 4. CRU2, CRU3 denotes CRUTEMv2, CRUTEMv3 surface trends, respectively. U4, R4 denotes UAH and RSS series, respectively, averaged over 4 grid cells in each 5x5 surface grid.

40