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 an analysis, even if it might be wrong.

So, I’ll be doing that, and hoping someone tells me what I’ve done right or wrong in applying the method.

So, this post organized as follows:

1. Discussion of Least Squares applied to GMST data from 1979-Jan 2008 without correcting for red noise.
2. Discussion of how I corrected for red noise using Cochrane-Orcutt.
3. Discussion of how I think one interprets the uncertainty intervals
4. Questions for real statisticians!

I’m asking questions because the issue of calculating the uncertainty intervals is key to validating or falsifying hypothesis. Today, I’m limiting my post to examining a long run of data, because the conclusions will be no surprise: Temperature rose since the late 70s. This lets me illustrate what I have done, and get answers to questions without worrying dealing with any potentially controversial issues.

So, here goes the boring post!

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 my test case, I averaged the temperature measurements from all four souces; I get a 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), I can look up the “t” for 5% confidence bands, and I find t=1.967. (I used TINV with probability 0.05 and 346 degrees of freedom.)

Using this “t”, I found confidence bands by taking the product product (t s_m= ±(1.967) (1 C/century). This results in: 1.5 C/century < m< 1.8 C/century. But there is a problem: The trend, "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.)

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 is statistically significant, 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
95% confidence intervals.	1.5 C/century < m< 1.8 C/century	1.0 C/century < m< 2.1 C/century	1.1 C/century < m< 2.0 C/century
ρ of residuals.	0.78	-0.21	0.03

Conclusions

I think I have applied Cochrane-Orcutt correctly to the monthly data from 1979-Jan. 2008. My goal to day is simply to request feed back from statistician in the know: Did I do this correctly.

I have done some reading, and understand that uncertainty intervals, while improved by using Cochrane-Orcut still tend to be too narrow. (I suspect this is because of the uncertainty in determining ρ)

So, if there are any real statisticians reading:

1. Did I apply this technique correctly.
2. Do I need to modify my uncertainty intervals? Of so, how do I do this?

To those wondering: Yes, I have done this for another more interesting data set. I have also attached the Excel spread sheet. Obviously, you can modify to look at other data sets. But, before you fiddle, let’s wait to hear back from the econometricians to see if I understood the method correctly.

—–
Attachment: Excel SpreadSheet Cochrane Orcutt– since 1979

Question I asked at a blog: Standard Errors.