Today, I’m going to try to apply Cochrane-Orcutt a method that could be used to deal with serial autocorrelation in residuals to a linear regression. I hadn’t done this before, so I hunted down some references, and read how it’s supposed to be done.

I have some result– but the problem is I really need someone to look them over and tell me if I’ve made a hash of this. And, even more unfortunately, the only way to do that is to a) publish them and b) do it in a way that makes a real statistician sufficiently interested to get them to read it.

So, I’ll be doing that, and hoping someone tells me what I’ve done right or wrong in applying the method. (Naturally, the actually interesting stuff is at the end– but it’s only interesting if the method is correct.)

Because this is my first attempt to apply Cochrane-Orcutt, I don’t want to state firm conclusions. So, this post organized as follows:

1. Discussion of Least Squares applied to GMST data from 1979-2007 inclusive without correcting for red noise.
2. Discussion of how I corrected for red noise using Cochrane-Orcutt.
3. Discussion of how I think one interpret the uncertainty intervals
4. Questions for real statisticians!
5. Discussion of Least Squares applied to GMST data from 2001-Jan 2008.
6. Results with correction for red noise Cochrane-Orcutt.
7. Repeat questions for real statisticians!

My questions for statisticians will related to:

a) Did I apply Cochrane-Orcutt correctly? I mean conceptually. I’m not asking anyone to plow through my numbers. I’ll check that again, but if the method is wrong, doing the arithmetic correctly can’t fix that and
b) Do the standard errors on transformed variables translate directly to standard errors for the slope for the variables of interest (that is Temperture and time?) Or do I have to do something?

Blind Ordinary Least Squares: Ignore Red Noise.

For the past few blog posts– and in many statistical analysis on the web, we’ve been assuming that GMST anomaly data can be described as obeying a linear trend with time, plus some amount of noise. The noise may be due to weather or measurement error or both. That is, if we have a series of temperature, time pairs (T_i,t_i) if denoted by T_i, the data obeys the following equation:

(1) T_i=( m t_i+ b )+ε_i

where “m” is a trend in temperature over time, b is an intercept at t=0, and ε_i is the ‘noise’ or ‘residuals’ associated with data point ‘i’.

Now, if the values of ε_i — the residuals– are all uncorrelated with each other, one can simply perform an OLS using a LINEST in Excel (or whatever tool you prefer.) The values of the slope, m, and the intercept pop-out. (They may be wrong. But Excel will give you numbers if you ask!)

To the left, I show simple linear regression to monthly data obtained using Excel.

Even better, LINEST provides values of the standard error in the slope, s_m. Knowing that, and the number of data pairs used to create the fit permits the analyst to state 95% uncertainty intervals, or do simple hypothesis tests.

In fact, using OLS results in fits like those illustrated in the figure to the left, where I include for GMST data sets plus the average of the four. As you can see, depending on which of four data sets for temperature measurements I used, the trend in temperature between 1978-2007 ranges from 1.8 C/century to 1.4 C/century.

For the temperature series I obtained by averaging the other four, I get m=1.6C/century with a claimed standard error or s_m=0.1 C/century. (If you download the EXCEL spread sheet, this result is the fit to average temperature and year in columns C and D. LINEST is done in cell C364.)

In principle, if GMST data really follows a linear fit as in (1), with this much data, we expect roughly 2/3 of all repeat ‘experiments’ will have slopes within ±0.1C/century if 1.6C/century. (Of course, we can’t do this. But if somehow running the earth’s climate were like flipping a coin, if we repeated 1979-2007, an infinite number of times, we’d expect this sort of result.)

But of course, “m” and particularly s_m aren’t correct. Why not, some might ask?

The Residuals Are Autocorrelated.

It turns out that if I examine the values of the residuals, ε_i, and calculate something called the “autocorrelation”, the residuals are strongly autocorrelated. (This is because weather trends tend to persist.)

To calculate the residuals, I use the formula:

(2) ε_i = T_i-( m t_i+ b )

If I’ve found “m” and “b” with excel, this is easy to do. I just put that in a column next to the (temperature, time) data ( T_i, t_i), enter the formula that matches (2) and create a column of residuals.

I then look examine the autocorrelation of adjacent residuals by taking the ratio of two sums:
(3) ρ = Σ (ε_iε_i-1/ ε_iε_i)

In Excel, is most easily done using the CORREL function. I did this for the ‘average’ monthly GMST measured using four standard methods, and got ρ=0.78. (See cell E10, shown in green.)

Now, it happens that there is a method to diagnose whether ρ=0.78 is statistically significant– but I haven’t learned it yet. Plus, for a correlation coefficient this large, I don’t need to. With 30 years worth of monthly data, I know ρ=0.78 is statistically significant.

Cochrane-Orcutt: Correction for Red Noise

I read about Cochrane-Orcutt. It appears to be very easy.

The first step is to assume the errors in the residual are AR(1) noise with a serial autocorrelation of ρ. If we know the ‘true’ correlation for residuals, ρ, we can linearly transform the variables of interest (i.e. temperature and time) to new variables using some simple algebra. Once we do this, we apply OLS to the transformed variables. If we had really known the true correlation, ρ, we would now be done. The trend, m_new, for the transformed variables is the same as for the original variables. We can also obtain a new intercept, b_new

However, if we didn’t know the true ρ, we do something a bit more involved. It requires iteration, involving a number of relatively simple steps.

The iteration:

1. We guess for the ‘true’ serial correlation for residuals, ρ; I’ll call the current guess ρ_old. (It makes sense to pick the current guess from the OLS we already did. So, I pick ρ_old=0.78 from above.)
2. Using ρ_old, transform the data set (T,t) into a new data set (X,Y). (Equation given later.)
3. After doing this, we obtain new estimates for slope and intercept (m_new, b_new).
4. We then compute new residuals, using the updated slope and intercept. Calculate ρ_new based on the new residuals.
5. If ρ_old=ρ_new to some limit, we just decide that’s ρ, and use these current values of “m” and “b” as correct. Otherwise, we

What I do involves a few steps:
Computational Step 1: Guess ρ_old. On my first try, I guessed ρ_old=0.78 because that’s what I got for the OLS on T and t. (I put this in cell E9 and named the variable Correl so I could later proof-read.)

Computational Step 2: Define and calculate new transformed variables (X,Y) using these two formulas:

(4a&b)
Y_i= T_i-ρ_oldT_i-1 and X_i= t_i-ρ_oldt_i-1

To plug and chug, I just created columns for ( Y_i, X_i) adjacent to those for (T_i, t_i) entered the formula and filled down. (Naturally, I have 1 less set of (Y,X) compared to (T,t). You can find Y and X in columns H and I. )

Computational Step 3: Perform a linear regression on (Y,X) to get new values for “m” and “β”. (I just use LINEST on (Y,X); m is the slope, β is the intercept.)

That is, I assume Y and X obey:

(5) Y_i=( m X_i+ β )+ε_i

I obtained the slope and intercept for Y vs. X using LINEST. The formula is in H364.

It turns out that the best fit slope for equation (5) is expected to return the same value as equation (1) for (T,t) pairs. So, “m” from this new fit is the new improved “m_new“.

However, the new intercept is not “b_new“. The intercept from this linest is β; this is related to the intercept in (1) as β= b*(1-ρ). So, I to get a new improved ‘b_new‘ . So, it can be calculated as b_new=β_new/(1-ρ_old).
(This isn’t actually in a cell. I use it algebraically in step 5.)
Now, in principle, I have “new improved” values of “m” and “b”.

Computational Step 5: Now, using the new values of “m_new” and “b_new” calculate use equation (2) and calculate new residuals, ε_i, using the original (T, t) pairs. (I turned this into a new column, rather than replacing the old ones. It’s column F. )

Computational Step 6: Now, equation (3) to re-compute a new estimate for ρ_new. (I used the “CORREL()” in Excel; see cell F10, in yellow.)

Compare the new value of ρ to the value you had before you began step 1. If the new and old values agree withing some criteria you set, you are done! The most recent values of “m” and “b” are the converged values, and apply to this data sample.

If the ρ differ too much, guess a new value for ρ_old. Your new guess should, of course be, ρ_new!

Repeat computational steps 1-5 but starting with the new value of ρ. (I substituted the new ρ_new for ρ_old by copying the new value, and then pasting it into E9 using “paste special, pasting only values.).

It happens that this computation converged in one iteration. (I set the cells to read to 4 decimal places.)

If a real statistician confirms this is correct, I’ll probably write a macro to simplify this. (I couldn’t find a free one for Excel.) But, I’m pretty sure that’s what all the various articles I found are telling me to do.

Results so far

Having done this, I now found m=1.5 C/century and s_m=0.3 C/century for JAN, 1979-now.

Now, I think that standard error is an honest to goodness standard error on the slope, “m”. I think I have 346 degrees of freedom for any t tests. The transformations (4 a&b) were supposed to give me a set of independent uncorrelated variables. So, I think I could now go ahead and do t-test with this.

Question for Real Statisticians

Can I do hypothesis tests using s_m=0.3 C/century as the standard error to get confidence intervals? Are the number of degrees of freedom driven 346? (That’s three less than the 349 (T,t) data pairs I began with.)

That’s my major question, and I’ll get back to it. (The reason it’s important has to do with later tests.)

But you aren’t really done.

Ok, that seems easy enough. Tedious, but easy.

Unfortunately, there are a few more steps. Once you have the value of ρ, you should check the residuals on the transformed variables, X and Y, and see if those residuals are correlated. Mine had a correlation coefficient of -0.21.

If this second set of residuals are correlated to a level that is is statistically significant, you repeat to basically pretend your (X,Y) variables are “original” variables, and redo the fit.

I skipped the test to see if ρ=-0.21 and just assumed it was. (I will teach myself how to test later on. But for now, I figured, just assume it is.)

When I finished this, my new results were:

Temperature Trend and Standard Errors for “mean” GMST anomaly.

	No Correction	AR(1) correction.	AR(2) correction.
m	1.6 C/century	1.5 C/century	1.5 C/centuy
sm	±0.1 C/century	±0.3 C/century	±0.2 C/century
ρ	0.78	-0.21	0.03
Note: These are one sigma, or roughly 67% confidence intervals for two-tailed tests. A little higher for 1 tailed tests.

Results for Data Starting Jan 1, 2001.

Here are the OLS curve fits look like to monthly data collected since 2001; they are not corrected for red noise.

Note that the trends range from m=+0.5 C/century to -0.8 C/century, spanning zero.

But, of course, this doesn’t account for serial autocorrelation. It mattered for the 30 year string of data. The same physics apply to the climate in more recent times so we expect the serial autocorrelation in residuals to remain constant.

So, to correct, I assumed I could use the value of ρ obtained from the previous fit. (That is, I used the 0.78 value etc. in the bottom row in the table above. Using these values consistent with assuming warming continues, and any pause is due to random weather noise. It also means I don’t recalculate it ρ based on the smaller data set. If I find the trend does not apply, the values of ρ calculated above would be dubious, so I would need to re-peat the analysis. You’ll see I didn’t do that yet.)

Other than my assumption for ρ I repeated the analysis above but applying it to the shorter set of data. This results in new values for the slope and the intercept.

What do I get?

Temperature Trend and Standard Errors for “mean” GMST anomaly: Since 2001.

	No Correction	AR(1) correction.	AR(2) correction.
m	-0.1 C/century	-2.4 C/century	-2.2 C/centuy
sm	±0.5 C/century	±2.0 C/century	±1.6 C/century
Note: These are one sigma, or roughly 67% confidence intervals for two-tailed tests. A little higher for 1 tailed tests. More important note: I’m actually asking if these are even the 1-sigma error bars.

Hmmm…. Hmmmmm…..

If you are a real statistician, and you see the magnitude of standard errors listed in the table above, you likely now understand why I am asking if these standard errors actually apply to “m”, and whether the number of degrees of freedom at each step is (N-2-number of red noise shifts.)

Because if that is correct…. Well, what would you say falsification of the IPCC’s upper range of projected values for the warming trend for the current, and upcoming decade.

Naturally, some will suggest this is meaningless because I must have cherry picked. In fact, I did not cherry pick a start date of 2001. The IPCC projections are supposed to apply starting 2001; they are based on a document published in 2000. So, that is why I did start this analysis in 2001!

(Although, in all honesty, if I were cherry picking, and I were a denialists, it turns out that 2001 is the year to start! That said, many denialists will be happy starting with 2000, 1998, 1997 and not too displeased starting with 1996. I couldn’t help checking out other choices of years after I’d done this.

As it happens, I believe AGW is real, and if the statisticians tell me I’ve done this more or less correctly, tomorrow I’ll explain why denialists can’t expect the whole world to throw IPCC predictions out the window… unless the recent cool weather spell persists or deepens. The drop over the previous year does affect the results quite a bit, and I want to be careful about that.)

Shortly, I’ll be talking a bit about why cherry picking intentionally is really really bad. Cherry picking unintentionally is also bad. And, I’ll also talk about why, potentially, monthly data gives cherry-pickers, so much to choose from. When I get confirmation on my interpretation of the meaning of thse standard errors, I’ll also explain why, if you are a denialist eager to decree global warming has ended, you need to wait a little longer. On the other hand, if you just want to suggest the IPCC’s most recent projections err alittle on the high side, it’a serioualy looking like they do.)

Conclusions

Honestly, I’m afraid to make conclusions at this point.
Clearly, I need to go back to the blackboard and find what I’ve done wrong. Maybe I’ve:

1. Filled some cells incorrectly in EXCEL? I’ve been known to do that!
2. Mis-understood how to handle red noise with Cochrane-Orcutt?
3. The standard errors, s_m need to be multiplied by something? To make them way bigger? Otherwise, I’m going to have to conclude “The data that has trickled in since Jan 1, 2001 is inconsistent with predictions of warming greater than 2.0C/century to a pretty decent degree of confidence!”

But…. no this is not cherry picked. The start date was selected as the earliest plausible start date to test the most recent IPCC projections; those were first published in 2000 and later included in the 2007 report!

————-
Notes:

1. EXCEL Spread Sheet include Cochrane Orcutt results form 2/6/2008

2. Reference on Cochrane Orcutt PDF read page 4 for iterative procedure.

3. Blog reference to Cochrane Orcutt http://economistsview.typepad.com/economics421/2006/02/the_cochraneorc.html

4. Hat tip to Basil at Watts Up With That for mentioning the technique by name so I could look it up!

5. Thejll & Schmidtt 2005, JGR Vol 110.. Paragraph 30 on page D18103 suggests a way to test whether two regressions differ statistically. So, if I can’t get an answer on the “s_m” question, I should be able to use that method. It’s straightforward.

6. Trends in climatological series.

7. Applied Statistics for Political Scientists. Explains things alternatives to Cochrane-Orcutt. But, reading the reasons CO is out of favor, I suspect these are reasons why I should use it. I think it’s telling me that, C-O has problems if the data model shifts over time, such that the correlation residuals changes with time . But in the case of climate, governed by physical laws like the Navier-Stokes, continuity, conservation of energy, the spectral properties should be stationary if the process is stationary. There are problems with the fact that I’m looking at non-stationarity, but the theory of AGW suggests we are choosing between weakly non-stationary and stationary possibilities. So, forcing the stastitical model to stay constant may be an advantage. (But, I need to make sure I really understand what the statistical papers are saying.)

8. Practical Considerations when estimating the presence of Autocorrelation. Page 4 of 13 tells me that I can use the estimates of the slope and it’s standard error going forward. Note: If I use the EXCEL method I described, I need to adjust the R².