{"id":173,"date":"2008-03-07T20:20:33","date_gmt":"2008-03-08T04:20:33","guid":{"rendered":"http:\/\/rankexploits.com\/musings\/2008\/correcting-for-serial-autocorrelation-cochrane-orcutt\/"},"modified":"2008-03-07T20:34:37","modified_gmt":"2008-03-08T04:34:37","slug":"correcting-for-serial-autocorrelation-cochrane-orcutt","status":"publish","type":"post","link":"https:\/\/rankexploits.com\/musings\/2008\/correcting-for-serial-autocorrelation-cochrane-orcutt\/","title":{"rendered":"Correcting for Serial Autocorrelation: Cochrane Orcutt"},"content":{"rendered":"

Today, I’m going to try to apply Cochrane-Orcutt<\/a> 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.<\/p>\n

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. <\/p>\n

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

So, this post organized as follows:<\/p>\n

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

    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.<\/p>\n

    So, here goes the boring post! <\/p>\n

    Ordinary Least Squares: Ignore Red Noise.<\/h3>\n

    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 (Ti<\/sub>,ti<\/sub>) if denoted by Ti<\/sub>, the data obeys the following equation: <\/p>\n

    (1)<\/eqnumber> Ti<\/sub>=( m ti<\/sub>+ b )+\u00ce\u00b5i<\/sub><\/equation><\/p>\n

    where “m” is a trend in temperature over time, b is an intercept at t=0, and \u00ce\u00b5i<\/sub> is the ‘noise’ or ‘residuals’ associated with data point ‘i’. <\/p>\n

    Now, if the values of \u00ce\u00b5i<\/sub> — 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!)<\/p>\n

    These are anomalies with 20 year trends shown.<\/a>To the left, I show simple linear regression to monthly data obtained using Excel.<\/p>\n

    Even better, LINEST provides values of the standard error in the slope, sm<\/sub>. 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. <\/p>\n

    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. <\/p>\n

    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 sm<\/sub>=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.)<\/p>\n

    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.)<\/p>\n

    Using this “t”, I found confidence bands by taking the product product (t sm<\/sub>= \u00c2\u00b1(1.967) (1 C\/century). This results in: 1.5 C\/century < m< 1.8 C\/century. \n\nBut there is a problem: The trend, \"m\" and particularly<\/em> sm<\/sub> aren’t correct. Why not, some might ask?<\/p>\n

    The Residuals Are Autocorrelated.<\/h3>\n

    Residuals from Least Squares Fit to Global Mean Temperature Anomalies.<\/a>It turns out that if I examine the values of the residuals, \u00ce\u00b5i<\/sub>, and calculate something called the “autocorrelation”, the residuals are strongly autocorrelated. (This is because weather trends tend to persist.)<\/p>\n

    To calculate the residuals, I use the formula:<\/p>\n

    (2)<\/eqnumber> \u00ce\u00b5i<\/sub> = Ti<\/sub>-( m ti<\/sub>+ b )<\/equation><\/p>\n

    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 ( Ti<\/sub>, ti<\/sub>), enter the formula that matches (2) and create a column of residuals.<\/p>\n

    I then look examine the autocorrelation of adjacent residuals by taking the ratio of two sums:
    \n(3)<\/eqnumber> \u00cf\u0081 = \u00ce\u00a3\u00ce\u00b5i<\/sub>\u00ce\u00b5i-1<\/sub>\/ \u00ce\u00a3\u00ce\u00b5i<\/sub>\u00ce\u00b5i<\/sub>)<\/equation><\/p>\n

    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 \u00cf\u0081=0.78. (See cell E10, shown in green.)<\/p>\n

    Now, it happens that there is a method to diagnose whether \u00cf\u0081=0.78 is statistically significant– but I haven’t learned it yet. Plus, for a correlation coefficient this<\/i> large, I don’t need to. With 30 years worth of monthly data, I know<\/em> \u00cf\u0081=0.78 is statistically significant. <\/p>\n

    Cochrane-Orcutt: Correction for Red Noise<\/h3>\n

    I read about Cochrane-Orcutt. It appears to be very easy. <\/p>\n

    The first step is to assume the errors in the residual are AR(1) noise with a serial autocorrelation of \u00cf\u0081. If we know the ‘true’ correlation for residuals, \u00cf\u0081, 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<\/i> variables. If we had really known the true correlation, \u00cf\u0081, we would now be done. The trend, mnew<\/sub>, for the transformed variables is the same as for the original variables. We can also obtain a new intercept, bnew<\/sub><\/p>\n

    However, if we didn’t know the true \u00cf\u0081, we do something a bit more involved. It requires iteration, involving a number of relatively simple steps. <\/p>\n

    The iteration: <\/p>\n

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

      What I do involves a few steps:
      \nComputational Step 1<\/strong>: Guess \u00cf\u0081old<\/sub>. On my first try, I guessed \u00cf\u0081old<\/sub>=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.)<\/p>\n

      Computational Step 2<\/strong>: Define and calculate new transformed variables (X,Y) using these two formulas:<\/p>\n

      (4a&b)<\/eqnumber>
      \nYi<\/sub>= Ti<\/sub>-\u00cf\u0081old<\/sub>Ti-1<\/sub> and Xi<\/sub>= ti<\/sub>-\u00cf\u0081old<\/sub>ti-1<\/sub><\/equation><\/p>\n

      To plug and chug, I just created columns for ( Yi<\/sub>, Xi<\/sub>) adjacent to those for (Ti<\/sub>, ti<\/sub>) 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. )<\/p>\n

      Computational Step 3:<\/strong> Perform a linear regression on (Y,X) to get new values for “m” and “\u00ce\u00b2”. (I just use LINEST on (Y,X); m is the slope, \u00ce\u00b2 is the intercept.)<\/p>\n

      That is, I assume Y and X obey:<\/p>\n

      (5)<\/eqnumber>Yi<\/sub>=( m Xi<\/sub>+ \u00ce\u00b2 )+\u00ce\u00b5i<\/sub><\/equation><\/p>\n

      I obtained the slope and intercept for Y vs. X using LINEST. The formula is in H364.<\/p>\n

      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 “mnew<\/sub>“.<\/p>\n

      However, the new intercept is not “bnew<\/sub>“. The intercept from this linest is \u00ce\u00b2; this is related to the intercept in (1) as \u00ce\u00b2= b*(1-\u00cf\u0081). So, I to get a new improved ‘bnew<\/sub>‘ . So, it can be calculated as bnew<\/sub>=\u00ce\u00b2new<\/sub>\/(1-\u00cf\u0081old<\/sub>).
      \n(This isn’t actually in a cell. I use it algebraically in step 5.)
      \nNow, in principle, I have “new improved” values of “m” and “b”.<\/p>\n

      Computational Step 5<\/strong>: Now, using the new<\/em> values of “mnew<\/sub>” and “bnew<\/sub>” calculate use equation (2) and calculate new<\/em> residuals, \u00ce\u00b5i<\/sub>, using the original (T, t) pairs. (I turned this into a new column, rather than replacing the old ones. It’s column F. )<\/p>\n

      Computational Step 6<\/strong>: Now, equation (3) to re-compute a new estimate for \u00cf\u0081new<\/sub>. (I used the “CORREL()” in Excel; see cell F10, in yellow.) <\/p>\n

      Compare the new value of \u00cf\u0081 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.<\/p>\n

      If the \u00cf\u0081 differ too much, guess a new value for \u00cf\u0081old<\/sub>. Your new guess should, of course be, \u00cf\u0081new<\/sub>! <\/p>\n

      Repeat computational steps 1-5 but starting with the new value of \u00cf\u0081. (I substituted the new \u00cf\u0081new<\/sub> for \u00cf\u0081old<\/sub> by copying the new value, and then pasting it into E9 using “paste special, pasting only values.). <\/p>\n

      It happens that this computation converged in one iteration. (I set the cells to read to 4 decimal places.)<\/p>\n

      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.<\/p>\n

      Results so far<\/h3>\n

      Having done this, I now found m=1.5 C\/century and sm<\/sub>=0.3 C\/century for JAN, 1979-now.<\/p>\n

      Now, I think<\/i> that standard error is<\/i> an honest to goodness standard error on the slope, “m”. I think<\/i> 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.<\/p>\n

      Question for Real Statisticians<\/h3>\n

      Can I do hypothesis tests using sm<\/sub>=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.)<\/p>\n

      But you aren’t really<\/i> done.<\/h3>\n

      Ok, that seems easy enough. Tedious, but easy.<\/p>\n

      Unfortunately, there are a few more steps. Once you have the value of \u00cf\u0081, you should check the residuals on the transformed<\/i> variables, X and Y, and see if those residuals are correlated. Mine had a correlation coefficient of -0.21. <\/p>\n

      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. <\/p>\n

      I skipped the test to see if \u00cf\u0081=-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.)<\/p>\n

      When I finished this, my new results were:<\/p>\n\n\n\n\n\n\n\n
      Temperature Trend and Standard Errors for “mean” GMST anomaly.<\/caption>\n
      <\/td>\nNo Correction<\/td>\nAR(1) correction.<\/td>\nAR(2) correction.<\/td>\n<\/tr>\n
      m<\/td>\n1.6 C\/century<\/td>\n1.5 C\/century<\/td>\n1.5 C\/centuy<\/td>\n<\/tr>\n
      sm<\/sub><\/td>\n\u00c2\u00b10.1 C\/century<\/td>\n\u00c2\u00b10.3 C\/century<\/td>\n\u00c2\u00b10.2 C\/century<\/td>\n<\/tr>\n
      95% confidence intervals.<\/td>\n1.5 C\/century < m< 1.8 C\/century<\/td>\n1.0 C\/century < m< 2.1 C\/century<\/td>\n1.1 C\/century < m< 2.0 C\/century<\/td>\n<\/tr>\n
      ρ of residuals.<\/td>\n0.78<\/td>\n-0.21 <\/td>\n0.03<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

      Conclusions<\/h3>\n

      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.<\/p>\n

      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 ρ) <\/p>\n

      So, if there are any real statisticians reading:<\/p>\n

        \n
      1. Did I apply this technique correctly.<\/li>\n
      2. Do I need to modify my uncertainty intervals? Of so, how do I do this? <\/li>\n<\/ol>\n

        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<\/i> correctly.<\/p>\n

        —–
        \nAttachment: Excel SpreadSheet Cochrane Orcutt\u00e2\u20ac\u201c since 1979<\/a><\/p>\n

        Question I asked at a blog: Standard Errors.<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"

        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 … Continue reading Correcting for Serial Autocorrelation: Cochrane Orcutt<\/span> →<\/span><\/a><\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[5,3],"tags":[],"class_list":["post-173","post","type-post","status-publish","format-standard","hentry","category-global-climate-change","category-statistics"],"_links":{"self":[{"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/posts\/173","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/comments?post=173"}],"version-history":[{"count":0,"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/posts\/173\/revisions"}],"wp:attachment":[{"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/media?parent=173"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/categories?post=173"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/rankexploits.com\/musings\/wp-json\/wp\/v2\/tags?post=173"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}