Provenance of the Decline, a Forensic Analysis
Posted by Jeff Id on November 27, 2011
There is much about paleoclimatology we don’t know. Key among the questions is, “What exactly was done to that innocent, unsuspecting data?” The stealthy nature of these statistical oddities has led to the practice of an informal new field of blog science which could accurately be called forensiclimatology. Long time reader Layman Lurker has discovered an interesting characteristic in Briffa’s famous hide-the-decline series. - Jeff
Exploring the Divergence Problem in the Briffa01 Timeseires
Guest Post By – Layman Lurker
Both before and after climategate 2 broke out, there has been ongoing discussion at tAV, Lucia’s and CA about divergence and particularly about whether the 1960 data truncation we saw in the TAR and other Briffa / Osborne publications was justified. It occurred to me during these discussions that an OLS regression of the full, non-truncated Briffa series in question and annual northern hemisphere temperature observations might show evidence of the divergence problem in the residuals. We are fortunate to have a non-truncated version thanks to a climategate (#1) email which Tim Osborne sent to Michael Mann on October 5, 1999 in the lead up to the IPCC TAR. Later on in another email, Osborne sent the data to Mann again but truncated at 1960, explaining that data after this point was unreliable due to the divergence and an apparent loss of temperature sensitivity. This reconstruction was unpublished at the time of the TAR (and was actually tagged as Briffa ’98 in the spaghetti graph), but was later published (truncated) in the 2001 paper: “Low frequency Temperature Variations from a Northern Tree Ring Density Network” in the Jounal of Geophysical Research. This paper also utilized a time series of NH annual temperature observations prepared by Phil Jones which I used for regression analysis in this post. At this point I should note that I have not read Briffa01 which appears to be behind pay wall (abstract here). However the time series data for both the truncated reconstruction, and the annual NH observations, is posted here.
As always, one should start by looking at the data. In this case we will plot the full, non-truncated (climategate) version of the Briffa reconstruction along side the Jones annual NH observations for comparison.If the truncated Briffa01 climate reconstruction is a good proxy for NH temperature as the paper implies, then ideally the ols regression residuals should fluctuate randomly around a stationary mean. In contrast the regression residuals of observations vs the fullBriffa should at least show a hint of the divergence issue. As the above figure #3 shows, the divergence period is not a distinguishable feature of the residuals using the full Briffa reconstruction. In addition, the overall regression fit is poor and the residuals are not stationary. Figure #2 residuals using the truncated version of Briffa01 confirms that in spite of being a better fit, that truncation does not eliminate the issue of non-stationary residuals and lack of cointegration.
Anyone familiar with R and who looks through my code will notice that I spent some time with the individual time series (both reconstruction and observations) to look at their stationarity. This is somewhat peripheral to the core of the analysis so I won’t get into the details here, but I found that observations and both Briffa series are not stationary. No real surprise there, however the book says some order of differencing is usually required to address the issue and render the series stationary. There are some simple tests one can do to confirm the level of differencing needed to achieve a stationary series. For observations it went according to script and diagnostics showed that the series was I(1) stationary (1 order of differencing). But when it came to the truncated and full versions of the Briffa reconstruction, every test showed that first order differencing resulted in indications of over differencing. A bit puzzling, but the reason for this should become clear as you continue to read.
At this point I wanted to look at the data from a bit of a different angle so I subtracted the full untruncated Briffa time series from observations and plotted.I found this image striking. The difference between observations and the full non-truncated Briffa series seems to come down to two linear trends (plus what looks like random unbiased noise) which changes from negative to positive slope circa 1900 and shows no hint of a divergence period starting in 1960. In other words, from 1900 and on, the difference between NH observations and the untruncated Briffa series appears to be characterized by a “trend stationary”process. If this is true, then there is a simple diagnostic indicator which should confirm this. If one differences a time series with a stationary trend, then the fitted MA coefficient should approach a unit root. When I fit an arima(0,1,1) model to the difference series, the fitted MA = -1.0 therefore the 1900 to 1994 difference between observations and the full Briffa reconstruction is very likely a “deterministic” trend and not a random realization of noise that just happens to look like a trend.
We would expect now that if the trend in the difference series accounted for the lack of cointegration between the Briffa series and observations, that detrending or a trend adjustment (rather than differencing) should cointegrate the two variables. This in turn should enable a valid regression with stationary residuals. This comes back to the puzzle of why the non-stationary Briffa reconstruction can’t be rendered stationary through differencing. My view is that a deterministic trend embedded within this series dominates the stochastic processes and throws off the effect of differencing. If this is true, then if we eliminate the effects of the trend difference we should be left with a stochastic, non-stationary time series. I’ll come back to this point later.
To adjust the Briffa series to account for the trend, I fitted a trend line to the 1900 to 1994 segment of the difference series from figure #2 then added the fitted values of the trend to the unadjusted Briffa reconstruction over the same time span. Now for the good part. The plot of the “trend adjusted” Briffa series along side the NH observations is shown below.
When I first looked at this image my jaw dropped. The “goodness of fit” R^2 statistic is much stronger for this than any of the other scenarios looked at. Besides the great fit, the residuals are stationary and consistent with expectations of a meaningful regression. The autocorrelation function for the residuals does not show significant autocorrelation at any time lag. Nor does it show any sign of over differencing (lag 1 autocorrelation with negative value greater than -0.5). These residuals are shown below in figure #6.
So what happens to the non-stationary character of the Briffa series now after the trend adjustment? The series is still non stationary, but adjusting the series to account for the deterministic trend leaves only stochastic non-stationarities behind. Simple diagnostics show that the adjusted series would now be rendered stationary after first order differencing.
If this analysis is correct, I believe a very compelling case is made that the so called late twentieth century divergence period (for this reconstruction) is an alignment illusion attributable to the deterministic trend difference of origin 1900 rather than 1960.
For me at least, this raises many questions: 1. What accounts for the trend? Is it real and natural? artifact of standardization? a bodge? 2. Is it possible that a simple trend adjustment to a proxy reconstruction could lead to such a remarkably good fit to observations?
In any event, the conventional divergence hypothesis does not hold up for this reconstruction.
There is follow up analysis that could be done with this. Comparisons at the grid level would be interesting. Also, after actually reading the paper and familiarizing myself with the methods (I am not familiar with either Hugershoff standardization or age banded decomposition) one might be able to track down artifacts of methodology as a possible cause.
There is an interesting quote from Briffa in the climategate 2 email #1055 commenting to colleagues about a presentation by Rosanne D’Arrigo on divergence at an NRC workshop which caused a bit of a stir. There seem to be strong hints from Briffa that at least part of the divergence problem could be an artifact of methodology. The following quote from Briffa is from email text file 1055 dated March 11, 2006:
First, it is important to note that the phenomena is complicated because it is not clearly identifiable as a ubiquitous problem. Rather it is a mix of possible regionally distinct indications, a possible mix of phenomena that is almost certainly in part due to the methodological aspects of the way tree-ring series are produced. This applies to my own work, but also very likely to other work. The reason we did not discuss it in detail in our Chapter is that the evidence is at present, in part contradictory, and not well defined. Giving it too much space would be tantamount to building a straw man at this stage; one that would then necessitate complicated details and space to knock down.
As theses are our data, I am able to say that initial unpublished work will show that the problem can be mitigated with the use of new, and again unpublished, chronology construction methods. In the case of the work by Rosanne and colleagues, I offer my educated opinion that the phenomenon they describe is likely also, at least in part, a chronology construction issue. I am not saying that this is a full explanation, and certainly there is the possibility of increased moisture stress on these trees, but at present the issue is still being defined and explored. As the issue needs more work, this is only an opinion, and until there is peer-reviewed and published evidence as to the degree of methodological uncertainty , it is not appropriate to criticize this or other work.
Could Briffa be discussing the effects we see in this post? Perhaps. In any case, I now believe that the divergence problem is a much messier issue that I previously thought.
I want to close by saying that I am open to discussion of any flaws. Undoubtedly many readers will know much more than I do about this type of time series analysis. Critiques are welcome. I will read and respond as I am able.
Abstract; Reconstruction and Observational data:
Briffa, K.R., T.J. Osborn, F.H. Schweingruber,
I.C. Harris, P.D. Jones, S.G. Shiyatov, and E.A. Vaganov, 2001,
Low-frequency temperature variations from a northern tree ring density network,
Journal of Geophysical Research, 106 D3 (16-Feb-2001) 2929-2941.
Email time series data version including post 1960 values:
On stationary and non-stationary processes of time series modelling the “Decision 411
Forecasting ” class notes from Duke University were referenced – as an example here is a
summary page on ARIMA modeling:
#Refer to Briffa 01 data (ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/briffa2001jgr3.txt) #for reconstruction data as published and for the Jones temperature observations. For the full, untruncated Briffa 01 #data refer to the climategate email from Tim osborne to Michael Mann dated Oct 5, 1999 which shows the full data set ending #at 1994. instbrif #objects: brif.trun brif.div brif.1900 inst.trun inst.div inst.1900 diff.1871 diff.1900 trend adj.brif.1900 year.div year.trun year.1900 #data plots of individual original series from 1871 to 1994 plot(inst.div~year.div,type="l") plot(brif.div~year.div,type="l") ts.plot(ts(inst.div,start=1871),ts(brif.div,start=1871),col=1:2,lwd=2,main="Jones99 NH annual temp (black) and full non-truncated Briffa01 (red) 1871 to 1994") #regress annual instrumental vs both the original truncated Briffa01 and the full untruncated (climategate) version of Briffa01 summary(lm(inst.trun~brif.trun))#adjusted R^2 0.4036 summary(lm(inst.div~brif.div))#adjusted R^2 0.1189 #The nature of the regression residuals were such that one should subtract one time series from the other for another visual perspective on the data. ts.plot(ts(inst.div-brif.div,start=1871),ylab="anomaly",type="l",lwd=2,main="Annual NH observations minus Full Non-Truncated Briffa01") #The nature of the 1900 to 1994 difference between observations and the full Briffa reconstruction appears visually from this graph to be "trend stationary" #Are the time series variations to consider stationary? #for brif.trun (1871 to 1960) acf(brif.trun)#significant autocorrelation to lag 9 brif.trun.diff plot(brif.trun.diff,type="l") acf(brif.trun.diff)#lag 1 negative autocorrelation with first order differencing is greater negative than a -0.5 threshold (likely overdifferenced) sd(brif.trun)#0.1772024 sd(brif.trun.diff)#0.2151509 differencing increases the sd suggesting the series should not be differenced #this series is not stationary but should not be differenced #for brif.div (1871 to 1994) acf(brif.div)#significant autocorrelation to lag 13 brif.div.diff plot(brif.div.diff,type="l") acf(brif.div.diff)#lag 1 negative autocorrelation exceeds -0.5 (overdifferenced) sd(brif.div)#0.1888438 sd(brif.div.diff)#0.2151509 differencing increases the sd #this series is not stationary but should not be differenced #for inst.trun (1871 to 1960) acf(inst.trun)#significant autocorrelation to lag 11 inst.trun.diff plot(inst.trun.diff,type="l") acf(inst.trun.diff)#lag 1 autocorrelation does not exceed -0.5 (not overdifferenced) sd(inst.trun)#0.2568897 sd(inst.trun.diff)#0.2359495 sd minimized in 1st difference #The series appears to be I(1) stationary #for inst.div (1871 to 1994) acf(inst.div)significant autocorrelation to lag 16 inst.div.diff plot(inst.div.diff,type="l") acf(inst.div.diff)#lag 1 autocorrelation does not exceed -0.5 (not overdifferenced) sd(inst.div)#0.2676663 sd(inst.div.diff)#0.2436191 sd minimized in 1st difference #The series appears to be I(1) stationary #At this point we see that none of the observational or reconstruction time series we wish to consider are stationary. To further complicate the situation, #the instrumental observations *can* be rendered stationary through 1st order differencing but *not* Briffa01 (truncated or full). #Why is differencing not effective to stationarize the Briffa series? #Are the time series variations to be considered cointegrated (stationary regression residuals)? rez ts.plot(ts(rez,start=1871),lwd=2,ylab="Anomaly",main="OLS Regression residuals of Observations and 1960 Truncated Briffa01") acf(rez,main="ACF of Briffa / Instrumental ols Residuals 1871 to 1960") #regression residuals are not stationary with significant autocorrelation extending to approx lag 10 therefore no cointegration rez2 ts.plot(ts(rez2,start=1871),lwd=2,ylab="Anomaly",main="OLS Regression residuals of Observations and Full Non-Truncated Briffa01") acf(rez2,main="ACF of Briffa / Instrumental ols Residuals 1900 to 1994") #regression residuals are not stationary with significant autocorrelation extending to approx lag 20 therefore no cointegration #Now the issue becomes whether or not one or both series can be transformed in such a way as to cointegrate the two series and allow for a meaningful regression. #Are the post 1900 differences between Briffa and observations stochastic or deterministic trend? #One diagnositic is to difference a time series and check the fitted MA coefficient for a unit root. #If a series with a deterministic trend is differenced, the fitted MA coefficient will approach a unit root. arima(diff.1900,order=c(0,1,1))#MA is -1.0 a perfect unit root acf(diff.1900) #According to this diagnostic the Briffa and observations series need to be trend equated or detrended (not #differenced) to achieve cointegration. ###plot comparison of trend adjusted Briffa to original briffa (starting at 1900) ts.plot(ts(brif.1900,start=1900),ts(adj.brif.1900,start=1900),col=1:2,lwd=2,main="Original Briffa Series (black) vs Trend Adjusted Briffa Series (red) 1900 to 1994") #summary statistics of ols regression with 1960 truncated Briffa01 (1871 to 1960) and analysis of residuals #plot of data: ts.plot(ts(inst.trun,start=1871),ts(brif.trun,start=1871),col=1:2,lwd=2,main="Original Scenario: NH Annual Observations (black) with 1960 truncated Briffa01 (red)") #summary stats and residuals: summary(lm(inst.trun~brif.trun))#adj R^2=0.4036; p-value: 1.055e-11 rez2 ts.plot(ts(rez2,start=1871),main="OLS Residuals of Observations vs truncated Briffa01 1871 to 1960") acf(rez2)#significant autocorrelation extending out first 6 lags with scattered significance up to 15 lags #these residuals are not stationary therefore observations and truncated Briffa01 are not cointegrated and regression results can be misleading #summary statistics of ols regression of NH annual observations vs trend adjusted Briffa01 and analysis of residuals #plot of data: ts.plot(ts(inst.1900,start=1900),ts(adj.brif.1900,start=1900),col=1:2,lwd=2,main="Jones99 NH annual temp (black) vs trend adjusted Briffa01 (red) 1900 to 1994 (adj. R^2=0.6366)") #summary stats and residuals: summary(lm(inst.1900~(brif.1900+trend)))#adjusted R^2 = 0.6366; p-value: < 2.2e-16 rez ts.plot(ts(rez,start=1900),ylab="anomaly",lwd=2,main="OLS Residuals of Observations vs trend Adjusted Briffa01 1900 to 1994") acf(rez,main="ACF of Observations vs trend adjusted Briffa01 OLS regression residuals ")#no significant correlations at any lag #residuals appear to be stationary therefore the trend equated regression succesfully cointegrated the two time series.