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

*full*Briffa 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.

Briffa continues:

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.

References:

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:

http://www.burtonsys.com/FOIA/0939154709.txt

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:

http://www.duke.edu/~rnau/arimrule.htm

#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.

It looks like I imbedded the wrong link for the data.

ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/n_hem_temp/briffa2001jgr3.txt

I fixed it.

OK, I read it.

“regression”, “stationary”, “residuals”,”stochastic non-stationarities” etc might as well be Martian to me.

I understand the original issue of the “divergence”, but what does this mean:

“…late twentieth century divergence period (for this reconstruction) is an alignment illusion attributable to the deterministic trend difference of origin 1900 rather than 1960.”

Thanks

Google Scholar gets you copies of many papers, free of any paywall. That includes “Low frequency Temperature Variations from a Northern Tree Ring Density Network” (Briffa et al., JGR 2001).

At this point I should note that I have not read Briffa01 which appears to be behind pay wall (abstract here).But somehow that did not prevent you from stating that the decline is “truncated” in Briffa et al. 2001.

You might want to have a quick look at the figures – and perhaps reevaluate the reliability of your sources 🙂

In any case, I now believe that the divergence problem is a much messier issue that I previously thought.Have a look at this for a review.

Matthew, figure 1 shows an image of a reconstruction roughly tracking observations until about 1960. When one subtracts the reconstruction values from the observations, you see an unwavering upward trend in this difference starting in 1900 (IOW the reconstruction values are loosing ground relative to observations ). This trend in the difference starts in 1900 not 1960. Looking closely at figure 1 you will notice that the reconstruction values are higher than the observations. This visually obscures the trend phenomenon until the two series become untangled in the later portion of the graph. This is the illusion.

As for your main point: you might be reading a bit too much in the small 1900-to-1960 variations. How would you distinguish them from just plain noise, especially considering that they do not seem to extend in the pre-1900 period?

This should read: “Looking closely at figure 1 you will notice that the

early part of the centuryreconstruction values are higher than the observations.”Toto

Did CA already do “D’Arrigo, Wilson et al 2007”

http://climateaudit.org/2007/05/07/ghcn-adjusted-data-isnt-good-enough-for-wilson/

#4

Toto, thanks for the link to the paper. I did a search with google scholar and I must have missed it. I figured someone would likely post it on this thread. When I get time I want to study the methods.

The data posted for the Briffa01 reconstruction ends in 1960. Plate 3 of the paper shows a spaghetti graph with the Briffa01 reconstruction truncated at 1960. Figure 6 refers back to Briffa 98

which was processed differently. It is an illustration only as far as I know. Please let us know if you can track down the full data used for that illustration.

#6

There is noise in the difference series from 1900 t0 1994. But it is stationary noise with a constant trend. As I noted in the post there appears to be two different trends in these differences which change from – to + slope circa 1900. I chose the post 1900 segment to analyze because it is this segment which contains the divergence period. It is a fair question to ask whether or not the pre 1900 segment has similar noise characteristics (but different slope) as post 1900. I don’t know why the upward trend starts when it does but I am confident that it is not due to random noise.

Very interesting post. I wonder if the divergence is a symptom of simple nonlinearity in the response of the trees in the series Briffa studied. Suppose for example the true response to rising temps is proportional to the square root of the rise over the minimum temperature that allows survival. That would show up as an almost linear deviation, at least over modest temperature changes. This sort of thing should be testable as well: the same species of trees should show lower response to changing temperatures the further away from the tree line; that change in response with distance from the tree line would then describe the true temperature response if the local temperature were known at each site.

#10

Thanks Steve. I’ll make the point now that I am open to all causal explanations. If there is a real life explanation for this, and it can be modelled, then it looks like mxd is a lot better proxy than a lot of people would have guessed.

Two graphs that tend to confirm a divergence appearing in the early twentieth century for a regional case :

(With a divergence which also appears to affect the glaciers).

And a hypothesis on the origin of the divergence:

I think toto believes that the only series that suffer from a divergence are north of 55°N. D’Arrigo gets pretty rounded criticized in the mails, I’m not sure if it’s because what she says is nonsense (like her “you have to pick cherries if you want to make cherry pie”) or what. In any case, the idea that the divergence is only in the ≥ 55°N series is flatly not true.

Anyway, my admonition is to be careful of who you read and always read it critically.

Layman brings up a pretty interesting question for me, which is the reliability of the early temperature record. We know that prior to 1950 the uncertainty gets pretty large. See e.g. this

Shouldn’t they include the uncertainty in the independent variable (temperature) when they regress against it?

Phi, I just tweaked the code to adjust the observations down by the trend in the differences so that it matched the reconstruction. The adjusted R^2 drops to 0.5018 compared to 0.6366 when you adjust the reconstruction up to match observations.

Phi: The MXD series extend way back in the past, right? Why not use the Luterbacher reconstruction (which go back to ~1500AD and are generally regarded as reliable) for a longer comparison and calibration?

Layman Lurker,

“The adjusted R^2 drops to 0.5018 compared to 0.6366 when you adjust the reconstruction up to match observations.”

It’s possible but you must also correct glaciers and TLT, this is all very strange.

Toto,

For regional data, I use only the raw ones and for HN, the data from Schweingruber before the stranges corrections made for the period after 1960.

It is also a bit odd but I do not like the data which were too sorted and manipulated.

Ooops – I linked to the wrong Luterbacher paper, I was thinking of this one.

#15

If I understand Phi’s hypothesis correctly (correct me if I’m wrong), he is suggesting that it might be observational record which contains the source of the trend (some type of bias presumably) in the differences, and that if we match observations trend to fit the Schwiengruber mxd reconstruction that this would also bring surface observations into line with the hypothesized tropospheric amplification,

Insteresting hypothesis. The goodness of fit is not quite as good when looking at that angle, but I’m not in any position to make causal inferences with this analysis.

Toto, have you read Jeff’s stuff on Luterbacher?

Layman Lurker,

Yes, basically that’s what I think. Instrumental temperatures pose a number of issues especially for the continents of the Northern Hemisphere:

– TLT divergence of 1 ° C per century since 1979,

– MXD divergence of 1 ° C per century since the beginning of the twentieth century,

– Divergence of melting anomalies of 1.5 ° C per century in the Alps,

– Unexplained bias from homogenization of 0.5 ° C per century,

– Surprising behavior T individual stations compared to regional averages.

All these points can be simply explained by the perturbations of the stations and there are simple physical explanation: urban drainage and energy consumption.

Have you done the analysis for the MXD only? It would be interesting to know whether the same result is obtained by correcting the proxy or T.

Toto,

Thank you for the links but I have not seen MXD series in your references.

Phi, I will say this much. A little over a week ago I would have argued against your hypothesis because of my understanding of the divergence period. After doing this analysis I think I would still be inclined to argue against your hypothesis, but I have now changed my views on divergence to the point that it would no longer support my argument (also keeping in mind this is only one case).

Layman Lurker,

Oops, I said something stupid, Briffa 2001 is pure MXD. So I have the answer.

There is something I do not understand well, how can you have such a difference (r^2) by correcting one or the other curve?

On the other hand, what happens if we change the origin of the correction (for global HN, it seemed to me that the divergence appeared in 1920 rather than 1900)?

Phi, I would rather refer to my exercise as “adjusting” rather than “correcting” because the cause is not clear.

I think that if the source of the trend is some type of bias in the reconstruction (my first inclination would be to look at possible standardization artifacts based on what I remember with Yamal and RCS) then I think that a justified upward adjustment of the reconstruction to equate with observations (assuming that observations don’t have trend error) would optimize the fit. Don’t you think? I know that R^2 is sensitive to a change in trend slope (all else being equal).

For my analysis, I doubt that shifting the start date forward would change things much. In the situation I have looked at here, I would say that 1920 would be to great of a shift.

“If I understand Phi’s hypothesis correctly (correct me if I’m wrong), he is suggesting that it might be observational record which contains the source of the trend (some type of bias presumably) in the differences, and that if we match observations trend to fit the Schwiengruber mxd reconstruction that this would also bring surface observations into line with the hypothesized tropospheric amplification,”

There is a paper ( esper I believe) who looks at the issue of divergence from the standpoint of errors in the temperature record.

basically the grid cells in that area have gone through a warming adjustment.

cant find the reference. Ive discussed this on CA

lemme look

its not this

http://www.springerlink.com/content/873486u687j56246/

try here

http://www.geo.uni-mainz.de/1358.php

look at this below.

120. Esper J, Frank DC, Büntgen U, Verstege A, Hantemirov RM, Kirdyanov AV (2010) Trends and uncertainties in Siberian indicators of 20th century warming. Global Change Biology 16, 386-398. (Esper_2010_GCB.pdf) (Esper_2010_GCB_sup.pdf)

Tree-ring series are adjusted for age-bias, i. e. that tree-growth slows down as a tree grows older. This is a somewhat controversial matter since this de-trending may also remove real low-frequency signals. It would seem that if such de-trending is applied a little too strongly the result might be a spurious trend like in Figure #4.

Tty, yes that very phenomenon was much discussed a couple of years ago here and at CA regarding RCS standardization and Yamal. There is an assumption of population homogeneity when mxd series are standardized by the same growth function. The potential problems if the homogeneity assumption doesn’t hold is just as you describe (discussed by Esper in the literature), and was suspected in the case of Yamal (search the tAV posts and you will find a great post by Jeff on this). I will have to study the methods used to see if similar problems could potentially arise.

“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.”

LL, I am not at all sure what a stationary trend is (existence of a deterministic trend implies the series is not stationary and must be differenced) and obtaining an MA coefficient near positive or negative 1 is a sign that your model is not fit properly and is probably over differenced. You should look at other models including an ARFIMA model with fractional integration (00.5) . Certainly a deterministic trend in a model would require at least one differencing to make it stationary. I am not certain whether what you end up with here is a second difference because you differenced an already differenced series.

It is rather well established that a temperature series over the instrumental period should contain a deterministic trend, but, with a proxy response, that is not necessarily true, and thus, I think you complicate things by subtracting the two series which could be the result of very different models.

I would strongly suggest that what you see here may well be something unique to the given proxy/reconstruction and that other individual proxies would show different divergence characteristics. If what you see here is something general you should see it when using other dendro proxies.

I’m glad there has been some interest in LL’s post. Oddly enough, my take on this is different from many of the comments. The fit LL achieved from a single trend makes me question whether the data after 1900 is actually temp. There are several peaks with absolutely perfect alignment to the temp series pre-1960. The fact that a single linear trend achieved such perfect residuals is telling. I wonder if we used Schweingruber local temps, whether the match would even be better?

We don’t really know what these scribbles are. We do know that the boys are prone to replacing “bad” data with good a the drop of a hat. It seems unlikely that the tree data would have the local peak correlations so perfect and have the correct amount of variance without having been influenced by a creative mind. I’ve been wrong before though.

Steven Mosher and LL, below is a paper to which I linked in another thread here at tav. It lists most of the current conjectures for explaining divergence. Perhaps LL you should like to add to the list.

Click to access DArrigo_GPC2007.pdf

Also I am curious to see who will write the review of conjectures for the non dendro divergence.

re 29.

Jeff dont you think the right approach is to

1. Calculate a regional temp for the proxy area. you specify a radius, you collect the stations, you cacluate a regional average

( using romans or taminos method) Forget using the damn CRU grids.

2. Then you have to calibrate that region to the hemisphere.

or am I stupid and can 1&2 actually be done in one step.Done is 2 steps you are going to have 2 different types of residuals

one that is due to the fit of the proxy to the regional temp and one that is due to the regional fit to the hemisphere.

err. needs an equation

Steven Mosher at 31, if you can posterior select the months, the trees and the temperature region,and the max, mean or minimum temperature and adjusted or unadjusted temperature and TRW or MXD or some combinatination why not go all the way and tune the factors that give the best calibration relationship. I don’t think you can be a little bit pregnant when it comes to in- and out-of-sample modeling and model calibration.

Steven,

I think your approach makes sense considering your software sophistication with temp series to which there is no current match. My eyeball guess here is that there is little or no proxy info since 1900 in this proxy series, it may be since 1960 but it is a very good match pre-1960. The boys would have chosen CRU gridcells if that is what this data actually is so if someone gets time, a processing of the CRU grids might prove to be an exact (or near exact) match.

The big point though if the data is not MXD after 1900 rather than 1960, is that nobody can claim ‘hide the decline’ was known in publication. Instead of malfeasance, it could be true temp signal in the MXD data due to my less recent familiarity with the curve but several of the extreme peaks align very well pre-1960.

“As these 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.”

LL, he is talking here about the use of RCS chronology development where the tree ring growth is adjusted for age. And it reminds me that I forgot another factor that can be used for obtaining the “correct” answer.

Kenneth, a deterministic trend implies stationary noise around a constant trend. Therefore you render the series stationary by detrending, not differencing. If you difference a series with a deterministic trend, you will induce a unit root in the fitted MA coefficient so then yes it would be overdifferenced. Sorry if I wasn’t clear on that, but that is what I showed in the post .

arima(0,1,1) on the difference series is a first order difference model which generates an MA coefficient with a unit root. This means that the difference series should not be differenced and serves to confirm the deterministic trend.

Steven Mosher:

As long as you really calibrate it (that is allow for a complex valued transfer function to account for the difference lags associated with the effects long-period versus short-period climate variability on dendro proxies), this has a chance of working. Even then, the nonlinearity of the response to the environmental forcings could still kill you.

“extreme peaks align very well pre-1960.”

Extreme peaks usually line up well in these dendro series such as the date of a major volcanic eruption and a blip on the TR series, but that says nothing about the magnitude of these bilps which as I recall do not necessarily scale with the expected response.

#29

Good point Jeff. This fit looks too good (to me anyway). I don’t have your feel for proxy data though. Given what we have always seen in the reconstructions I sure would have expected a lot more noise even if we are confident it has been properly fitted to observations.

“Kenneth, a deterministic trend implies stationary noise around a constant trend. Therefore you render the series stationary by detrending, not differencing. If you difference a series with a deterministic trend, you will induce a unit root in the fitted MA coefficient so then yes it would be overdifferenced. Sorry if I wasn’t clear on that, but that is what I showed in the post ”

I believe you render a series with a deterministic trend stationary by using the residuals or by first differencing. I have not seen the case where first differencing a series with a trend will result in overdifferencing and a unit AR or MA coefficient. Could you link me to a source that explains this relationship? I should obtain a unit root by first differencing the instrumental.record. .

LL, a trivial point, but when I model the CRU temperatures the best fit (of the few that I tried) is for arima(0,1,1) and I obtain an MA coefficient of -0.735 and an AIC of -39. If I difference again, i.e. arima(0,2,1) I obtain an MA coefficient of -1.00 and an AIC of 12.

Kenneth, if you consult the Duke stuff that I referenced, you should find your answere somewhere but I don’t think you’ll need it if you think through a simple simulation. White noise plus a trend.

It is a non-stationary series (non constant mean) but what do you have if you detrend it? Stationary white noise. What would you get if you difference it? I added white noise to the calculated trend from my code, then ran arima(0,1,1) again. Here is the call:

Call:

arima(x = (trend + sim.noise), order = c(0, 1, 1))

Coefficients:

ma1

-0.9382

s.e. 0.0320

Kenneth, actually, if you look at the Wiki link that I put up in the post, you’ll see the confirmation that you need in the first sentence of the article. “Trend stationary” and “Deterministic trend” are interchangeable terms. http://en.wikipedia.org/wiki/Trend_stationary

I think your example with CRU does not consider the possibility of stochastic non stationary procesess being significant – this does not rule out a deterministic trend being part of the mix.

LL, I graphed the Briffa proxies, S1 through S6 and the CRU temperature data S7 and linked the graphs below.

LL, I agree with you that the optimum way of making a series with a deterministic trend stationary is to use the residuals from a time regression. What I do not agree with is that first differencing a series with deterministic trend will necessarily give an AR or MA coefficient of one, but that a second differencing often will. It is not an easy task to determine whether a deterministic trend or a stochastic one are present in the series and I do not believe the test would be to first difference the series and look at a MA or AR coefficients.

I looked at the CRU residuals (detrended) and the series still needed to be first differenced to get the lower AIC value and a second differencing gave a -1.00 for MA..

Click to access Briffa.2001.jgr.pdf

Here is my marble graph for the Mann 2008 data. Blue is negative slope, red is positive.

Fans of tree-ring proxies will be cheered to learn that 664 of 927 tree-ring series had positive trends for this period most of them in the US.

An irony here is that D’Arrigo claims that the divergence in the US (which mostly doesn’t happen except in mountain species which ironically are supposed to be temperature-limited) is due to “global dimming”. Apparently that claim was made without reference to the basic facts, which is that most of the US since 1960 has seen an increase in sunlight (probably due to improved pollution controls) in rural areas, with the observed dimming in the US concentrated near large cities (i.e., not close to where the measurements were actually taken).

(I have data on this too, from nearly 300 sites around the US measured between 1960 and 1990. This is part of a publicly available data set from the DOE.)

D’Arrigo and others (Cook 2004 I think toto said?) have made a big deal about how the divergence is limited to species about 55°N. Nothing could be further from the truth. In the subsample above 55°N, 39 of 58 series show a positive trend over this period (a full 2/3s of the data), which compares to 71% for the total sample. Without running statistics, they look “the same”.

In fact, the Southern Hemisphere seems to suffer a larger divergence problem, and again it seems to be the “super special mountain species” that have the biggest problem. I’m wishing at this point, Mann would have included the elevations that the measurements were taken at, because it would be interesting to split it out that way.

Finally a comment on correlation: Correlation is only useful in comparing series X to Y when Y is responding with the same amplitude and phase as X in the frequency range being compared (the transfer function between Y and X is “flat”). Otherwise, you have to use a more sophisticated method, or carefully limit yourself to regions where agreement is good.

For another way of looking at these proxies, Here’s a histogram of tree-ring trends, and basic statistical information: mean = 3.451365382794 for 927 points, std. dev. = 6.8329260275118 (the error of the mean was 0.224327, so yes the mean slope is significant).

My opinion: You can make almost any claim you want from plots of spaghetti curves or looking at individual series because these plots say virtually nothing and are mostly just noise. You have to go beyond that to more rigorous testing if you are going to firm up any conclusions.

I think this is appropriate, especially in light of Craig’s comment “like seeing cartoons in the clouds” in reference to the black art of plotting stupid pictures without providing any accompanying quantitative analysis.

The work that Briffa was talking about is his work with Melvin on chronologies.

The famous “Briffa bodge”, interestingly,was a forcible linear adjustment of data to remove the decline – originally done on Tornetrask in Briffa et al 1992 Clim Dyn).

the argument that divergence doesn’t affect species below 55N relies on the very data that has been most questioned. The Cook article that is the source of this meme was based on a small sample of sites, with strip bark (foxtails) contributing much of the overall increase.

Changing the interval to 1968-1998 doesn’t dramatically affect the results, the variance just goes up in a predictable fashion, and the median (but not mean) trend is increased slightly, as would be expected.

Shortening the interval further with relatively few species and a pretty scare global coverage isn’t likely to a sensible act.

Doing some sort of gridded calibration of the form described by Steven Mosher is probably doable, with the caveat that we need to look at coherence between temperature and tree ring proxies. If we are restricting ourselves to the same grid cell for temperature and tree ring proxies, “errors in variables” or TLS probably isn’t necessary but it’s something that should be explored.

Steve McIntyre:

And as is typical of that field, the memes have nothing to do with reality—they are totally orthogonal to the truth (not even just wrong, they are complete poppycock).

This does illustrate a problem with my little foray into tree-ring proxies: I really need better meta data than I have so far, so I can split the trends out between different species and growing conditions, including especially elevation effect.

Layman Lurker #24, Tty #26,

The only systematic bias that I know for MXD is that of the increase in ring density with age (optimization of the resistance of trunk by a greater strength in the periphery, analogy of hollow sections). For the particular case of the twentieth century, this potential bias uncorrected may lead to an overestimate of the warming. I would add that for the regional series (see # 12), I used the raw data that I have normalized but there is no weighting and no density correction.

I would also be suspicious that the fit is too good to be true for a tree-based reconstruction.

Craig Loehle,

MXD are quite remarkable. In my opinion, there is no reason to doubt of the raw data. If there had been manipulation at this level, the authors would not had fun adding a linear divergence.

I am starting to think the divergence problem may not be with Briffa’s trees.

What if they are a fairly accurate gauge of temperature, and it is in fact the thermometers that are diverging? Maybe the biases (UHI, LHI, siting, bucket/intake) and adjustments to the temperature trends are actually the things causing the divergence.

Phi,

I would agree with you that there is sometimes a correlation of tree data with temp but with so much processing going on, how do we know what we are looking at?

Would you agree that the correlation does not demonstrate a long term signal is present?

Are you aware of the uneven expression of variance which occurs in any regression of hockeystick data on temperatures. Mann07 used unrealistically low autocorrelation values to show the effect doesn’t exist – I confirmed his results and then used better values. My estimates show that the variance difference from the calibration to historic range can be near 50%. It doesn’t matter which regression method a person used either.

Layman Lurker,

Maybe acidic rain (sulfate aerosols) which limits nutrient availability in certain forest environments, explains the divergence.

See graph for total sulfate emissions in http://www.pnl.gov/main/publications/external/technical_reports/PNNL-14537.pdf

I fit three alternative arima models to the the regression residuals (observations vs Briffa01 + trend). The model orders were: (0,0,0), (1,0,0), (0,0,1). The aic score was minimized with white noise model.

Therefore the Briffa01 series after 1900 could accurately be described as follows: Briffa01 = NH instrumental observations + linear trend + white noise

Truthfully, my gut reaction to this is: “You gotta be shitten’ me.”

D’Arrigo and others (Cook 2004 I think toto said?) have made a big deal about how the divergence is limited to species about 55°N. Nothing could be further from the truth. In the subsample above 55°N, 39 of 58 series show a positive trend over this period (a full 2/3s of the data), which compares to 71% for the total sample. Without running statistics, they look “the same”.In a world where “post-1960 divergence” would somehow be identical to “negative trend 1968-1998”, that would be a point.

Regarding Cook 2004, you might want to RTFriendlyP.

the argument that divergence doesn’t affect species below 55N relies on the very data that has been most questioned. The Cook article that is the source of this meme was based on a small sample of sites, with strip bark (foxtails) contributing much of the overall increase.That’s very possible, but :

1- What happens when you include reconstructions from other extra-tropical NH locales (while still ensuring that you’re only looking at locations that have indeed warmed post-1960)?

2- The point is not that the South trees are the one true arbiter of temperature truth (as Cook mentions, they may well have an inverse, positive bias in the recent period – apparently the observation is from Briffa himself – I suppose that’s what you’re hinting at). The point is that while North and South diverge strongly post 1960, they remain highly similar pre-1960, all the way to the Middle Ages. As do the comparisons between diverging and non-diverging that I know of. If you have some counter-examples I’d be glad to have a look at them.

Carrick @ 46:

I like, of course, your showing results of individual proxies and how, in doing so, we get a different and more detailed picture of a reconstruction. What it does not show is the magnitude of the trends and further what we do not know is how much posterior selection was carried out that might have weeded out some negative trending proxies. What I would find even more informative would be to use an a prior selection criteria for tree candidates followed by a complete publication of the tree ring responses to local temperatures and noting the portion of trees that show negative, positive and no trends. That result could be compared to some model to estimate those portions that could occur by chance. I would think such an excercise would be a logical analysis for publication, but I have never found such a paper.

To me the situation with divergence and other problematic issues in these matters seems always to be approached with the idea that the proxies are basically valid thermometers, both past and present, and any contradictions or anomalies need only be analyzed and explained (away), but never by looking squarely at the issue that the proxies might not be valid thermometers.

Jeff Id,

I do not have great knowledge in statistics and my approach of the issue is necessarily limited. I have no particular opinion on the possibility of using MXD to measure temperatures of a thousand years ago, I assume that adaptation and natural selection can distort the results. I am mainly interested in temperatures of the twentieth century and I try to see what can be drawn from MXD for this limited period. In this context, I found an excellent high-frequency correlation with the T but a divergence similar to that presented by glaciers in a regional case and similar to the divergence of TLT for NH land (and some other stuff). But, from Mann cooking, I only know what I read especially on your blog.

Toto:

Whether it is or isn’t depends on how reconstruction is done of course. If you start out with 2/3s of the series responding positively from 1968-1998, and your method flips the sign of these, I’d suspect there is a problem in the reconstruction not in the original data.

I’ve had a look at Cook. Incidentally this is the unbroken link. Since he restricts himself to northern hemisphere, and has fewer proxies, I don’t think there is any fundamental disagreement between him and Mann.

Here’s the relevant figure.

He only had 8 specimens north of 55°N, and obviously you can’t draw very many conclusions from that. Using his method, my suspicion is there isn’t statistically any difference in the performance of the ≥ 55°N series compared to the sub 55°N series.

Cook does do a coherence study by the way that I really liked. It’s moving things in the right direction, IMO.

LL @ 42:

The value that you obtain for the MA coefficient after first differencing a series with a linear trend depends on the amount of noise that the linear trend contains. A linear trend with 0 noise will indeed yield an MA=-1.00, but one with moderate noise will yield MA coefficients closer to zero and will yield an MA of -1.00 only when you do a second differencing.

Kenneth Fritsch,

Posterior selection is something I don’t have control for. It does appears that Mann 08 used complete series from different sources, so he at least was aware of the bias issues associated with posterior selection.

In terms of the magnitude of the trends, here is a csv file of all of the reconstructions if you want to look at what I’ve found in more detail.

ID=9000 is the one for tree-rings.

I am surmising that some reading here are surprised by the coherence of the various proxies/reconstructions. The proxies can respond synchronously to some variable or variables where temperature is only a minor contributor or where the level of response varies from proxy to proxy. What is also obvious from viewing these series is that the amplitudes of these synchronous responses are significantly variable amongst the proxies and reconstructions. In fact when one looks at individual proxy time series it is rather obvious that the variation over time varies from proxy to proxy, i.e. the standard deviations of the residuals of the series vary greatly. I think that looking at reconstructions for these features might be misleading since, I think, the methods used for reconstruction could affect these variations.

“Posterior selection is something I don’t have control for. It does appears that Mann 08 used complete series from different sources, so he at least was aware of the bias issues associated with posterior selection.”

Mann (08) does, of course, use a posterior screening process, that selected the instrumental record (as proxies) and the Schweingruber MXD series that was truncated at 1960 and added back with something not well explained.and the Tijlander proxies used upside down.

I have the Mann (08) proxies and have in fact done at good deal of analysis on them. It was these proxies and my anlysis that really reawakened me to the problems that these proxies/reconstructions have – even though others, such as Jeff ID and Steve M, had blogged on these issues previously.

#63 Kenneth Fritsch

Kenneth, I don’t disagree with that other than to say that white noise is not 0 noise. And you can increase the amplitude of the white noise case and still almost always get unit roots in a fitted MA of a first differenced series with a trend. So consider this: if you have a noisy time series which is non-stationary, do you not agree now that you almost certainly can confirm a deterministic trend as the *only* non-stationary property of the series if you difference the series and find a unit root in the fitted MA? I concede that as you increase the influence of autocorrelation in the errors, that you get a corresponding decrease in the incidence of a unit root.

Better choice of words for #67 would be “autocorrelation in the noise” rather than “autocorrelation in the errors”.

Kenneth, here is a short script for a 1000 simulation monte carlo. The simulation(s) will construct a time series consisting of a linear trend plus varying levels of 1st order autocorrelated noise. The MA coefficient will be extracted and plotted in a frequency histogram. Here is the script:

n=95

ar=0.3

p=1000

{

a<-matrix(NA,n,p);

for (j in 1:p ) {

b<-trend+arima.sim(model=list(ar=ar),n)

a[,j]<-arima(b,order=c(0,1,1))$coef

}

}

hist(a[1,])

It is not quite turnkey as it uses the "trend" object from my script from this post. So you will have to calculate trend. It is set up to tweak the AR(1) coefficient which is currently set at 0.3.

The histogram using ar = 0.3 is here.

It shows pretty much as I expected (I forget what that type of distribution is called) but you will see the shape change starting at ar=0.4 and by ar=0.5 the distribution changes to more of a normal shape with the peak frequency gradually moving away from a unit root with increasing AR(1). I believe this type of change makes sense as the AR coefficient itself moves toward a unit root.

Kenneth Fritsch:

I’m not sure what you mean by “selected the instrumental record … as proxies”. Wouldn’t you select the temperature record that is geographically in the same region as your measurements? That’s not post-hoc, as long as you select it based on an objective criterion that doesn’t depend on what the data look like.

I don’t like or trust the Schweingruber MXD series, and won’t use them, so there’s that.

If there are say 80 different measured series by a given observer and you use the entire set, I’m not really sure you can accuse Mann or anybody else of anything nefarious here.

The judgement is out there for me still whether or not there is low-frequency signal in the data. Cook set out to answer it with his coherence study, but used too coarse of a bin size. Right method wrong execution.

“Kenneth, I don’t disagree with that other than to say that white noise is not 0 noise. And you can increase the amplitude of the white noise case and still almost always get unit roots in a fitted MA of a first differenced series with a trend. So consider this: if you have a noisy time series which is non-stationary, do you not agree now that you almost certainly can confirm a deterministic trend as the *only* non-stationary property of the series if you difference the series and find a unit root in the fitted MA? I concede that as you increase the influence of autocorrelation in the errors, that you get a corresponding decrease in the incidence of a unit root.”

LL, I have added different amounts of white noise (random/normal) to a trend in my analysis and the MA difference from zero does vary with that noise.

A unit root in MA can occur if a series not requiring it (to make it stationary) is first differenced. That is why a second differencing often leads to AR and MA unit roots. I have to reread and analyze what you did with the instrumental record and the Briffa reconstruction. I was under the impression that an implied differencing was carried out and that by differencing you were perhaps actually doing a second differencing.

I am assuming that you subtracted the instrumental data from the reconstruction data and then worked with the residuals of that series resulting from the initial subtraction and then regressing that subtracted difference versus time. But it was subtracted series that gave you a unit root for MA on differencing – and not the residuals of the subtracted series. Is that correct?

LL, I ran a regression for all 6 of the proxies and on the CRU instrumental data during the period 1871 to 1992. All seven series had statistically significant trends and when I ran an arima(0,1,1) on all series the proxies were not that different than the temperature series and in fact the temperature series was near the middle of the range of the proxy series. Range for MA was -0.48 to -0.83 and the temperature MA was -0.68.

I then did a difference of all 6 proxy series with the temperature series. Of the 6 difference series 2 had a statistically significant trends, while 4 did not. The arima(0,1,1) on the 6 difference series did not produce any unit MAs. The MAs values were -0.66, -0.70, -0.81, -0.81, -0.79 and -0.77. A second difference gave MAs of -1.00, i.e. for arima(0,2,1).

I have to reread and analyze what you did with the instrumental record and the Briffa reconstruction. I was under the impression that an implied differencing was carried out and that by differencing you were perhaps actually doing a second differencing.

Were you perhaps confused by me refering to the “difference series”? This is the time series calculated when Briffa01 was subtracted from observations, but it was never “differenced”. I definitely should have used a different handle when refering to that. To be clear – only 1 order of differencing on that series with a fitted MA = -1.0

Regarding your #72: are you making a distinction between a “series with a trend” and a series with a “deterministic trend”? Two different things when one is talking about stationarizing a time series. I believe I have demonstrated that the time series ~ “observations minus Briffa01” (formerly refered to by me as the “difference series” 🙂 ) is a non stationary series characterized by a deterministic trend. When I do *first order* differencing of “OMB” I get an MA unit root. When I leave the “OMB” series undifferenced and instead detrend, I get stationary white noise.

This type of result is only possible in a series characterized soley by a deterministic trend. As soon as I start throwing in factors (on top of the deterministic trend) plus white noise.

The histograms generated by my script above demonstrate this. I’ll try and post up the trend component to make it turn key a little later

“Were you perhaps confused by me refering to the “difference series”?”

That had me for a bit. I read your referenced links which helped with the terminology and concepts you were expressing. The code helped me work things out.

#66 Ken: [my bold]

Jeff quoting the Osborne email (#2346)in his head post, Paleoclimate – Rotten to the core :

Not previously “well explained”, but looking at the Mann 2008 Schweingruber data, it seems pretty obvious now.

As well, there was an email quoted in the Readme file in the emails which caught my attention:

One of the arguments for a genuine “decline” is that the data closely matches temperatures in the pre-1961 time period. The 1998 paper quotes only statistics for the post-1880 time period and there is NO mention that temperature data is available before that. However, the above quote makes it apparent that there may have been some undisclosed “cherrypicking” on the other end in order to boost “validation statistics and provide justification that the reconstruction had

anyvalue whatsoever.#76 Rats! Firefox acting up and deleting my name in the previous comment.

LL at Post #73, I agree with what you are saying about detrending and trends in general. Unfortunately, a trend in a time series in combination with the expected noise in most time series and where only part of the time series is analyzed, one could have a deterministic trend or even deterministic segmented trends that would appear statistically indistinguishable from a trend that arose an ARFIMA or ARIMA model. In other words, if a model with a deterministic trend plus noise (white and/or red) can yield essentially the same data as that of an ARFIMA model, I do not think there would be a method that would allow you to determine the model used.

I need to state also that I am not anything close to an expert on these models or analyzing, but that only I have worked with modeling data and that is why I was surprised that you obtained a -1.00 MA on a first differencing of a series that likely contained a deterministic trend (or two) and had the expected noise.

My previous post says that I do not get MA = -1.00 when I do for the individual proxies that you did for the combined proxies. I will now attempt to do duplicate your arima(0,1,1) on the combined proxy result minus the temperature data.

Hi Kenneth. I actually ran the code again to check and am now getting MA=- -0.7947. I have another file of code that I did most of my work in. The file that I posted here was a cleaned up version. I’ll look at the original file to see if I got crossed up somewhere but in the meantime it sure looks like I made a mistake (allmost looks like I mistakenly did a second order diff).

LL, I just made the arima(0,1,1) calculation on the combined series (of the 6 proxies series) minus the temperature series and obtained MA = -0.74. I used the period 1871 to 1991. When I second differenced (arima(0,2,1), I obtain MA=-1.00. I am therefore not duplicating your calculation and one or both of us must be doing something wrong/different.

RomanM at Post #76:

Mann (08) did some cutting to remove a part of a proxy that was “under responding” to recent temperatures. What might be surprising to some is that in Mann(98) the authors actually “adjusted” data in a proxy that was showing too much response to temperature. In Mann (98) the North American TR series PCs were adjusted for showing unrealistic growth in the latter parts of the series. The adjustment was based on another TR series. The rationale for the adjustment was that the North America TR PC showed growth that was conjectured to be caused by CO2 fertilization in trees at high altitudes. That growth leveled off in the latest part of the series and this was attributed by the authors to a “saturation” effect. The authors ended up stating that regardless of the exact mechanism the growth was unnatural and with anthropogenic causes. Would said you cannot have it both ways – or any ways?

Kenneth, try using the composite, non-truncated series from the email here.

Kenneth, I think we cross posted in #79 and #80.

Kenneth Fritsch:

Why we really need the raw data sets. You just throw out obviously screwy data (and explain why you know it’s “screwy”), you don’t tamper with it to try and get it to look like the rest of the data.

LL, we did indeed cross-post. If you are saying that using the linked composite and the one I derived provides the differences in MAs that you noted -.079 and -0.74 I will take your word for it. Neither of these values is -1.00. Have you tracked down how you obtained -1.00? I would continue to guess that it was a second difference.

The above should be -0.79 and -0.74.

#85 Kenneth

Right now I am confused. My original code and the call which I just ran again is shown below along with my comment on the same line (diff.1900 is the “O minus B” time series):

As you can see I differenced the series first prior to arima. Then I fitted arima for MA(1) without differencing. In the version of the code I used for this post, I transcribed to what I thought was a more simple and straight forward version of the same process, removing the separate differencing operation and then including it in the arima model. The line of code and the call as follows:

AFAIK the above two operations should be identical. WTF?

This change of operation dramatically affects the shape of the histogram which my monte carlo code generates with white noise.Here is a modified monte carlo script from #75 reflecting the change of operation described in #87.

The code now first order differences the simulated series seperately, then does an MA(1) fit. A white noise histogram from the modified script is here.

The original monte carlo script, which performed first order differencing within the arima model, produced this histogram.

Aha, the issue I encountered is documented here as issue #2. The correct way to do this (which was a pure fluke on my part) is the original way I coded it, IOW MA = -1.0. This is now confirmed by the histogram from the modified script in #88 which shows that first order differencing a series with deterministic trend + white noise is virtually certain to induce a unit root in a fitted MA(1). Even introducing progressively higher autocorrelation changes a bit (particulary with values up to +/- 0.4) from the script in #75.

LL, I did the calculation both ways, i.e. arima(x, order=c(0,1,1)) and arima(diff(x),order=c(0,0,1)) and I get the same answer MA=-0.74. If I do arima(diff(x),order=c(0,1,1)) I get MA= -1.00. We still have a “dfifference” here.

I did correlations of all the pair-wise combinations of the 6 Briffa proxies linked by LL and obtained the following results with combinations noted by column number, numbered left to right, of proxy used:

The combinations are: B12, B13, B14, B15, B16, B23, B24, B25, B26, B34, B35, B36, B45, B46 and B56.

The correlations are: 0.50, 0.62, 0.43, 0.72, 0.54, 0.37, 0.35, 0.61, 0.55, 0.56, 0.49, 0.51, 0.68, 0.49 and 0.84.

I did correlation of the proxies in the first 2 columns of the Briffa reconstruction(B12) in series created by dividing the time series into 100 year periods. The correlations were starting with the earliest 100 year period: 0.32, 0.22, 0.36, 0.35, 0.11, 0.17, 0.46, 0.44, 0.28 and 0.65.

In summary these results would indicate that over the entire time period the proxies responses correlated fairly well, but over any 100 year period the correlations were considerably poorer. That the proxy series show some synchronous blips seems to belie the deterioration of proxy correlations over shorter time periods.

Based on the correlation of proxy series and synchronous blips in most series the critical issue is how do the proxies respond to temperature in the instrumental period. Below I show in the link the six Briffa proxies in the instrumental period with the temperature series. One can see an alignment of some of the proxy peaks and valleys with temperature but the pattern is rather hit and miss in my view. Obviously six proxies is merely anecdotal evidence given the thousand of dendro proxies used in reconstructions. Showing these data does, however, allow me to discuss my point in all this which is the synchronous behavior of proxies is something that I think Carrick has pointed to as being important in attempting to understand to what exactly the proxies are responding. I believe that I have seen proxies that are more and less synchronous as these 6 in the Briffa reconstruction. I do not recall, however, seeing a published analysis of this synchronous feature of proxies that goes into much detail.

Firstly, if proxies show reasonably good synchronous responses throughout the period of interest it says that proxies are reacting similarly to some variable (probably climate related), but not temperature necessarily. Secondly if proxies show reasonably good synchronous responses with temperature in the instrumental period that would indicate that the synchronous blips we see in the entire proxy series might well be reactions to temperature. Thirdly there remains the problem of these synchronous peaks and valleys responding with the same amplitude with temperature or with each other for that matter. Without a reasonably close response in amplitude the use of the proxy as a good thermometer breaks down because we do not know how to select the proxies that might be responding correctly over the entire time period. I think these amplitude differences are the most easily observed weakness in the use of proxies as thermometers.

My alternate view of all this is that the proxies follow a low frequency model response with a long term persistence but with an overlay of higher frequency that results from reactions to a climate variable(s) with an unknown portion of which is temperature. The divergence problem is trying to tell us something, but unfortunately I suspect that the explanations depend a good deal on the models we all have in mind about the proxies generally. I think divergence could merely be the result of many of the original proxies being selected posterior (cherry picked) to fit the instrumental temperature and the divergence is a result more or less of an out-of-sample test of the goodness of the proxy responses to temperature. That Mann (08) mentions that even non dendro proxies show some divergence would seem to reinforce this view. I know from looking at lots of ARIMA and ARFIMA models that one can obtain series excursions up and down at the end of the series that look much like those we see in the selected proxies for reconstructions. What I do not think we can obtain is as much synchronous peaks and valleys as we see in the real proxy series.

LL, I finally noted that you were using the 1900-1994 data (I used 1871 to 1991 data) with the Briffa reconstruction minus the temperatures. Now when I use arima(x, order=c(0,1,1)) I obtain MA= -0.79, but when I use arima(diff(x), order=c(0,0,1)) I obtain MA=-1.00. I further checked by first differencing the Briffa minus temperature series in Excel abd then doing arima(0,0,1) in R. From this I obtained MA=-1.00. Now I agree with your results but I still do not know why the arima(0,1,1) does not give MA=-1.00. And why the methods give the same result when using 1871 to 1991 data.

LL, I went back to the Briffa Con -Temperature series (DiffX) and in the R library(forecast) I used the auto.arima function to find the best fit for DiffX (limited to p and q max of 3 and obtained the following results. Notice that auto.arima selects an arima model with arima(2,1,1) and that does not give an MA=-1.00 in that result or when I do arima(2,1,1) direct. If however I difference first and then use arima(2,0,1) I obtain MA=-1.00.

auto.arima(DiffX, step=TRUE, allowdrift=TRUE, max.p=3, max.q=3)

Series: DiffX

ARIMA(2,1,1)

Coefficients:

ar1 ar2 ma1

-0.0227 -0.1843 -0.7480

s.e. 0.1320 0.1238 0.0819

sigma^2 estimated as 0.03081: log likelihood = 29.61

AIC = -51.23 AICc = -50.78 BIC = -41.05

> arima(DiffX,order=c(2,1,1))

Series: DiffX

ARIMA(2,1,1)

Coefficients:

ar1 ar2 ma1

-0.0227 -0.1843 -0.7480

s.e. 0.1320 0.1238 0.0819

sigma^2 estimated as 0.03081: log likelihood = 29.61

AIC = -51.23 AICc = -50.78 BIC = -41.05

arima(diff(DiffX),order=c(2,0,1))

Series: diff(DiffX)

ARIMA(2,0,1) with non-zero mean

Coefficients:

ar1 ar2 ma1 intercept

0.0664 -0.1364 -1.0000 -0.0083

s.e. 0.1084 0.1109 0.0362 0.0006

sigma^2 estimated as 0.02639: log likelihood = 35.09

AIC = -60.18 AICc = -59.5 BIC = -47.47

Kenneth, I had spotted that too and I’m glad we got that reconciled. It was a great discussion, without which I never would have uncovered this issue with R. Crap, I have a couple of other scripts which I was playing around with arima modelling and always using arima do the differencing. I have some tweaking to do now.

When I get time I will modify my monte carlo script again to include a second trend segement similar to pre 1900 of Briffa01 with it’s opposite slope and adding in the same noise. I think that the introduction of a second source of non-stationarity just muddies the “unit root” waters. I have been doing some arima modelling following an iterative process as laid out in the Duke material I referenced. Absolutely fascinating. On one particular time series I did 14 iterations or something when I found myself looped back to a duplicate configuration from a previous iteration so it was never going to converge to an optimized solution. We’ll see where this goes. If it looks interesting I may send it to Jeff to see if he wants to put it up.

#93 Kenneth.

Isn’t that interesting. So auto.arima has the same issue (at least it looks that way). In this case it looks like it actually prevents auto.arima from converging on the right solution. Using the Duke (manual) process in your second example, having run your arima(diff(DiffX),order=c(2,0,1)) configuration and obtaining a unit root, you would then reduce 1 order of differencing and 1 order of MA and then run the next iteration which would be arima((DiffX),order=c(2,0,0)).

Of course if the only source of non stationarity is a linear trend and not a stochastic process then I don’t think auto.arima would work anyway because it can’t detrend.

LL, it is interesting that what we are discussing here, the unit-root versus trend-stationary, is part of a major debate in economics about the GDP growth.

http://practicalquant.blogspot.com/2009/03/unit-root-or-trend-stationary.html

http://en.wikipedia.org/wiki/Unit_root

The proper (or perhaps best effort) test is evidently the Dickey Fuller Test which is described here :

http://en.wikipedia.org/wiki/Dickey-Fuller_test

And for use with R here:

http://rss.acs.unt.edu/Rdoc/library/uroot/html/ADF.test.html

LL @ 96:

“Of course if the only source of non stationarity is a linear trend and not a stochastic process then I don’t think auto.arima would work anyway because it can’t detrend”

I think it can, but I have not experimented with it. The option, xreg, I think is what one would use.

LL, you might want to try this R function: adf.test {tseries}. I could not find the library {uroot}.

LL, the Augmented Dickey-Fuller test run on R below indicates that the series Briffa reconstruction – NH temperature (DiffX) is stationary using p <= 0.05 for statistical significance. I will now shut up and let you proceed.

library(tseries)

TS=ts(DiffX,start=1900)

adf.test(TS)

Augmented Dickey-Fuller Test

data: TS

Dickey-Fuller = -3.4968, Lag order = 4, p-value = 0.04639

alternative hypothesis: stationary

Kenneth, I differenced the full “OMB” series going back to 1871 which of course appears to have two different and opposite slope trends. Then I fit MA(1). Then I did the same excercise for the 1871 to 1899 (which I haven’t done previously). The fitted MA on the full series was -0.8131. The 1871 to 1899 segment was -0.9994. So when the series is broken and analyzed like this as two separate segments, the MA(1) unit root occurs in each segment.

I duplicated the same analysis as posted for the 1871 to 1899 segment. It is also a trend + stationary white noise.

I then detrended each segment of the Briffa01 series by their corresponding slopes and now plotted the full combined adjusted series along side NH observations from 1871 to 1994 here. Residuals here.

LL, I was thinking that looking at the earlier occuring line segment in your difference would be the next step – and primarily to determine if that segment could be adjusted to fit the Biffra series as does the later occurring segment. Now you have a two segment linear trend differences between the Briffa reconstruction and the temperture series that have an apparent breakpoint near the year 1900. That evidence could imply that the basis for the divergence is in the entire instrumental period and that around 1900 we see an abrupt change in the divergence going from negative to positive – as opposed to neutral to positive. I would think that neutral to positive would be easier to explain than negative to positive.

I would think of further interest would be to look at the 6 individual proxies that make up the Briffa reconstruction the same way you have for the proxy composite (reconstruction).

Also I wanted here to reiterate what I said earlier about series generated from ARIMA and ARFIMA models that can have apparent trends in 100 year or longer time segments of the series that would be difficult to differentiate from a series with a deterministic trend, i.e. the segments would pass the Augmented Dickey-Fuller test. I may generate some series and test my supposition.

Also, LL, I have noted that both R functions arima and auto.arima have a provision for xreg which is there in order to regress the series and obtain the residuals from which to model. As it then turns out if you suspect that you have a deterministic trend throughout the series you are working with you have to use that option or you will get wrong answers with these functions. I have learned something from this excercise – besides learning that I know less about these things than I thought I did.

As it stands right now the two detrending functions do not intersect at 1900 and are therefore offset. When I get time, I may shift the breakpoint of the two segments in order to eliminate the offset. Then I can see if this results in a better “goodness of fit”.

LL, I did a breakpoint calculation in R and came up with 1 breakdate for the Briffa Recon – temperature series at 1901.

TS_Briff_Temp=ts(DiffTX,start=1871)

library(strucchange)

Year=time(TS_Briff_Temp)

breakpoints(TS_Briff_Temp~Year)

Optimal 2-segment partition:

Call:

breakpoints.formula(formula = TS_Briff_Temp ~ Year)

Breakpoints at observation number:

31

Corresponding to breakdates:

1901

I find it interesting that when you subtract the temperature series from the Briffa reconstruction and then compare that subtracted series back to the temperature series (after detrending the subtracted series) you obtain reasonably good coherence of the 2 series. If the original Briffa Recon and temperature series had good coherence a subtracted series would not necessarily show coherence but something more like flat plateaus of various heights. I am attempting to put my head around what I am looking at here and I thought the reason we do not see plateaus or something approximating them is because the temperature series might be showing more high frequency variation than the reconstruction – but that is not obvious in looking at your graphs.

Just to make sure we are on the same page: Regression of the original Briffa series from 1871 to 1994 resulted in an R^2 fit of 0.1189 to go with non stationary residuals proving that the two series are not cointegrated. I used the “Observation minus Briffa” series to determine that the difference between the two series (1900 to 94) is a deterministic trend. I used the calculated trend to apply detrending on

Briffa01(assuming the observations to be correct) which cointegrated Briffa01 and Observations. This was confirmed by regression of the observations with the detrended Briffa01 showing residuals consisting of stationary white noise.LL @107:

That would explain why I do not see plateaus and in fact I should be looking at the same coherence as in the orginal Briffa01 versus temperature – just without the trend difference.

Kenneth the image of original Briffa01 vs “detrended” Briffa01 1900 – 1994 is here.

LL, with the R code below I did an arfima model of the first proxy listed in the Briffa reconstruction (column 1) and then used those coefficients to obtain an arfima model and selected approximate century long trends from parts of that series and tested those trends for being trend stationary using the augmented Dickey-Fuller test. The 1000 year long arfima generated series is shown in the link below.

My point here is merely to show that a trend indistinguishable from a deterministic trend can be generated by an arfima model (and I assume also by an arima model) provided that I am allowed to select a portion of a long series. This exercise does not preclude that you found a deterministic trend that might have a basis in the physics of the matter but rather only is a precautionary note.

library(forecast)

TSB1=ts(B1,start=1000)

B1_Arfima=arfima(TSB1)

Coefficients:

d ma.ma1

0.4040910 0.1016734

FS_B1=fracdiff.sim(1000, ma = 0.102, d=0.404)

TS_FS_B1=ts(unlist(FS_B1), start=1000)

plot(TS_FS_B1)

Wind=window(TS_FS_B1,start=1590, end=1700)

adf.test(Wind)

data: Wind

Dickey-Fuller = -4.2594, Lag order = 4, p-value = 0.01

alternative hypothesis: stationary

Wind2=window(TS_FS_B1,start=1240, end=1360)

adf.test(Wind2)

Augmented Dickey-Fuller Test

data: Wind2

Dickey-Fuller = -4.4755, Lag order = 4, p-value = 0.01

alternative hypothesis: stationary

Wind3=window(TS_FS_B1,start=1280, end=1360)

adf.test(Wind3)

data: Wind3

Dickey-Fuller = -5.171, Lag order = 4, p-value = 0.01

alternative hypothesis: stationary

Interesting Kenneth. Column 1 is the Jones et al climate reconstruction. It wouldn’t surprise me if segments from Briffa01 or any other reconstructions could be shown to be deterministic trends. I should point out that neither of the 1900 to 1994 segments of either the observations or Briffa01 are characterized by a deterministic trend. It is the

differencebetween the two series which is deterministic.LL, I looked at the 6 individual Briffa proxies from the Briffa reconstruction in the same manner as you and I did for the composite reconstruction by subtracting the temperature from each proxy series for the period 1900-1993. I recorded the MA coefficient from the arima(diff(x),order=c(0,0,1) )and the probability from the Augmented Dickey_Fuller test (tseries) for each of proxy subtracted series.

I obtained the following for the subtracted proxies proceeding from B1 through B6:

MA coefficients: -.082, -0.88, -1.00, -0.88, -0.70 and -0.70.

.

adf.test probabilities: 0.258, 0.1012, 0.047, 0.053, 0.620 and 0.520.

From these results and using a rejection level from the adf.test of 0.05 that a series with a unit root could occur by chance that looks like the series in question we see that we could accept the alternative hypothesis that the series is trend stationary for B3 and marginally for B4 and perhaps even B2, but not for B1, B5 and B6.

Note that the individual proxy series do not all extend to 1993 with the shortest being to 1961, but nonetheless if there is linear trend present in the proxy if should be manifested in the shorter segment.

#3 and #4 are Briffa01 and Briffa2000. #2 is Mann99

It is rather late in the day for me to realize this but what I thought were proxies for a Briffa reconstruction are actually reconstructions in their own right and the one in questions is Briffa (01) which appears as a truncated reconstruction that I have called B3. The untruncated Briffa (01) reconstruction is what you linked separately. LL, do you know if the proxies that make up the Briffa (01) reconstruction are available online?

I don’t think the data as used has been posted. Tim Osborne’s web page links to this table which provides documentation on proxy names and locations. He then refers to the ITRDB which has a search engine using key words to identify type of proxy and the investigator (Schweingruber). Ex: Arrowsmith Mountain is one of the sites. The raw measurement files (width, density, etc) are posted for each site. The ITRDB page links to a page which explains the file format to identify the various raw data and site chronology files for each sample site. There are “Latewood Density” files and “Maximum Density” files but not “Maximum Latewood Density”.

I have done a first pass read of Briffa01 to review the methods. A script would need to be written to download the data into R. A full replication would be an ambitious undertaking to say the least, but it would still be interesting to even look at the data from one site, then follow the described methods to begin the “Age Banding” process.

I found some reconstructions (including the ones I misrepresented as proxies from Briffa (01)) that appear to have been used in Briffa(01) in the links below. Even if they are not the reconstructions used in Briffa (01) they show a difinite divergence and thus I’ll download the data and have a look at how they perform after subtracting the temperature data.

http://www.ncdc.noaa.gov/paleo/metadata/noaa-recon-6274.html

http://www.ncdc.noaa.gov/paleo/pubs/briffa2001/briffa2001.html

LL, I looked at those 9 reconstructions (from the link in a previous post that I think are at least part of the Briffa (01) Reconstruction) by subtracting the temperature series from the reconstructions for the period 1900-1994 and then doing the arima(diff(x), order=c(0,0,1)) and the ADF test. The reconstructions in order of the results below were: NEUR, SEUR, NSIB, ESIB, CAS, TIBP, WNA, NWNA and ECCA.

The MA coefficients calculated were: -1.00, -1.00, -1.00, -0.85, -1.00, -0.78, -1.00, -1.00 and -1.00.

The adf.test probabilities were: <0.01, 0.01, 0.01, 0.48, <0.01, 0.57, <0.01, <0.01, and <0.01.

We can therefore see that 7 of the 9 reconstruction subtraction series show strong evidence of being trend stationary while 2 do not for the period starting from 1900 and ending at various dates for reconstructions on or before 1994.

The subtracted series are shown in the plots linked below with the order of appearance being from left to right and then down the same as given in the first paragraph above.

Layman Lurker #115,

The data I used for the first two charts that I linked in # 12 are from the source you provide.

Unless mistake on my part, Latewood density and Maximum density mean the same thing. The summer wood is simply more dense than spring wood.

I noticed something else, in terms of comparisons global NH, it seems that the correlation is better with annual temperatures than the one you use (probably from April to September ?). This is especially true for the period before 1920, I don’t know why.

Layman Lurker #115,

The data I used for the first two charts that I linked in # 12 are from the source you provide.

Unless mistake on my part, Latewood density and Maximum density mean the same thing. The summer wood is simply more dense than spring wood.

I noticed something else, in terms of comparisons global NH, it seems that the correlation is better with annual temperatures than the one you use (probably from April to September ?). This is especially true for the period before 1920, I don’t know why.

Response to comment / question by Matt Skaggs at Climate Audit. http://climateaudit.org/2011/12/01/hide-the-decline-plus/#comment-314618

Matt, I belive that the simple difference time series (figure #4) between observations and Briffa01 is much more revealing than the regression residuals of the same series in figure #3. Hard to read more into residuals other than to show that the fit is not good and the series are not cointegrated. After the post I looked at the the pre 1900 segment of this series and showed that it is also characterized by an deterministic linear trend and stationary white noise (but opposite sign). See the Briffa01 detrended plot alongside the observations (link in comment #101) to see the fit. If you have not already, I would suggest reading the methods section of Briffa01 regarding the “age band decomposition” method of standardizing. I think that Briffa is right that the method is susceptible to bias in that it could affect the quality of the age adjustment. Whether it accounts for a relationship whereby temperatures differ from a reconstruction by two linear trends with a sharp knuckle at 1900 is a different story.

Reposting a couple climategate E-mails from CA:

from 0493, anonymous paper reviewer:

“It is well known that individual standardisation of tree-ring series does not allow reconstructions for frequency bands greater than the mean segment length of the single series. Additionally this standardization produces even end effects that might be at least partly responsible for the a divergence problem in recent years.”

and:

Briffa 4450:

“The divergence is only apparent (at the NH average scale) in the smoothed i.e. lower-frequency domain. You need to smooth the tree-ring records and the temperature to see it. However, the divergence is largely an artifact of using curve fitting (i.e. based on least-squares fitted regression lines or functions ) to estimate the unwanted (biological) growth trend in the tree-ring data. These fits are influenced by climate warming signals in the recent data, and this signal is inadvertently removed in the standardising process. When non-curve fitting methods are used (such as RCS) this problem is largely removed.”

The first quote implies that the problem is somehow related to the lower curve in Figure 2 of Briffa 01. The second quote seems to imply that older trees are being interpreted as younger trees by the ABD algorithm (and so the ABD correction is swamping the temperature signal), but that just makes no sense if longer series are underrepresented.

LL, I have the Schweingruber MXD series from Mann (08). There are 105 individual series that are cut-off at 1960 due the divergence problem. In the final series used in the Mann (08) reconstruction the 1960 forward portion was replaced with some not explained data. SteveM at CA and others have attempted to find the portions cut-off without success to my recollections. I was pondering whether the reconstructions used in Briffa (01) and that I have analyzed in a previous post on this thread used at least some of these Schweingruber proxies – and reported the post 1960 data. The dates correspond well and the series are, I believe, both MXD and not TRW.

Regardless of this use, we have a MXD series with an admitted case of divergence with the Schweingruber series. You have shown in the Briffa (01) reconstruction (and we see nearly the same in Briffa (00)) that the divergence could well start in 1900 and thus even though we have a shortened Schweingruber series we have the 1900 to 1960 portion to test for a trend stationary as you did for Briffa (01). We could even do the 1871-1899 portion. As I recall the anomalies for the Briffa data and the NH temperature series was based on the time period 1961-1990 and that period would have to adjusted for this comparison.

Do you think this analysis could provide useful information?

Yes Kenneth I think that might be a handy way to check to see if these proxies can replicate the regional chronologies. If they can then we know that the trend was introduced at the proxy level. One would have to check Mann08 to see how these proxies were utilized.

LL, I did the following with the Schweingruber MXD proxy series (SGS ) with recognized divergence from the period 1960 to 2000 for the periods 1900-1960, 1871-1899 and 1871-1960:

I re-baselined the 1871-1960 SGS and the NH temperature series to the period 1920-1950 and obtained 105 anomaly series from the SGS and the one for the NH temperature based on the new baseline. I subtracted the NH temperature from each of the 105 SGS and obtained 105 series I’ll call DSGS. I then did the following for each of the time periods noted above, i.e. 1900-1960, 1871-1899 and 1871-1960:

Calculated the MA coefficients from the R function arima(diff(x), order=c(0,0,1)) for each series, the adf.test p.value in R for accepting the alternate hypothesis of a stationary series and the trend slope and p.values for zero slope for each series. Included with these results are the latitude and longitude of each series location and a proxy identifier.

I show in links below the tabled results for the 1871-1960 period and will summarize the results for all periods here.

The number of DSGS series that were indicated as stationary using the MA coefficient as an indication were 65 , 87 and 64 for the periods 1900-1960, 1871-1899 and 1871-1960, respectively. The Augmented Dickey-Fuller test in R indicated stationary series for the time periods in the order above of 80, 9 and 82. All tests used a rejection level of the null hypothesis of 0.05. You can see that these two different tests have some disagreements on what is stationary and particularly so for the 1871-1899 period.

I then looked at the trend slopes for all these series in all 3 time periods and the probability of rejecting the null hypothesis that the slope was zero. The trend results show a large variation in slope for all time periods that does not appear to be related to a geographic location. The number of DSGS series with slopes greater than zero were 61, 23, 72 for the time periods 1900-1960, 1871-1899 and 1871-1960, respectively. All the slopes, regardless of the p.value indicating zero slope, all were negative for all series for the periods 1900-1960 and 1871-1960, while for the period 1871-1899 25 DSGS series had negative slopes.

The 105 DSGS series were plotted and a visually representative sample of 48 are linked below. Observing these series and considering the resulted listed above indicates to me that the trend difference seen in the Briffa (01) reconstruction follows through to the Schweingruber proxies and a case can be made that the divergence could start at least as far back as 1871. The divergence varies significantly from proxy to proxy even in those from the same approximate location where one might assume that local temperatures and variations are nearly the same. The DSGS series plots shows some proxies have a continuation of the negative trend back to 1871 from 1960 others appear to be flatter in the 1871-1899 period, while some show what the Briffa (01) reconstruction showed, i.e. a positive slope in the 1871-1899 period.

The link below shows a statistical summary table of series created by subtracting the NH Temperature from the Schweingruber proxies. The table specifically shows augmented Dickey_Fuller test probabilities of the null hypothesis that the series has unit root, slopes of series trends and probabilities of the null hypothesis that the series slopes are zero along with latitude and longitude of each series and its Proxy ID number.

The four linked images below are plots of the Schweingruber minus the NH temperature series for the period 1871-1960 showing the results for the Proxy IDs 1 through 48

Kenneth,

I’ve followed your discussion all the way. This is an interesting result that seems to confirm the possibility that something unexplained was done to the data. .You have put a lot of time into this, is it terribly difficult to look at CRU for local temps in these series? It seems to me that there is a chance of a near perfect match. I suppose the assumption is that the data has been altered by trend but it might simply be replaced.

#124

Great work Kenneth. Do you know how the proxies you were working with were standardized (I am assuming that they were not raw)? One would need to confirm that the standardization was the same or similar to Briffa01. Also, Are you able to cross reference your proxy data with the ITRDB files? I would be interesting to look at the raw versions in comparison.

Intuitively, I think the variability in the individual proxies is to be expected. Combining the proxies should cancel out much of the offsetting noise.

Some points I would like to clarify:

1. When you did arima(diff(x), order=c(0,0,1)) when you speak of “stationary” am I correct in interpreting this as an MA(1) unit root which indicates “trend stationarity”. I say this because a DSGS series with a known trend would not be “mean stationary” unless you detrend it first.

2. You state:

Did you mean to say “while for the period 1871-1899 25 DSGS series had

positiveslopes.”LL to answer your questions and reply to your comments in order:

I think that the standardization used may have been different, but showing that will take some additional digging. The Schweingruber series as used in and downloaded from Mann (08) was standardized.

I do not buy, at this point in time, that one can take individual proxies with greatly varying responses and “average” them together to obtain a true picture of the past climate and specifically the temperature. Even if that were a legitimate method, I would think that the CIs resulting from acknowledging all that proxy noise would be from floor to ceiling.

Answers to points 1 and 2 are yes and no, respectively. Remember that the DSGS series are derived from subtracting the NH temperature series from the Schweingruber proxies. What you found in the Briffa (01) reconstruction would have given a negative slope in the 1900-1994 period (and thus the 1900-1960 period) while the 1871-1899 period would have yielded a positive slope. What I found in the 1871 to 1899 period was 25 DSGS series with negative trend slopes that were not necessarily different than zero on statistical testing while the remaining 80 series had positive slopes that were not necessarily different than zero on statistical testing. The number of DSGS series in the 1871-1899 period that had statistically significant slopes were 23 and of those 23, 21 had positive slopes. I suspect when I average all series together I will see the trends breakpoint at or near 1900 as you did with Briffa(01). I plan to look at the breakpoint data in the averaged series.

I was going to reply with a post with the comments below and also note that I would like to present some probabilities of trends in pair-wise difference series from the individual DSGS series to give a feel for what would be significantly different and what would not. Also I want to average together the individual DSGS series to see what the composite shows with regards to what you found with Briffa (01).

I should try to summarize my findings, but I do think that at this point the analyses are far from complete. After finally understanding what Layman Lurker (LL) found and looking at individual reconstructions and proxies I would say that MXD proxies do indeed suffer from a divergence as LL found and going further back in time than is currently suggested by climate scientists. In other words, I think LL is onto something here and has presented an analysis that deserves attention.

I do see differences in trends between the 105 individual MXD Schweingruber series minus NH temperature in the 1871-1960 and 1900-1960 periods, but the series are noisy and I am finding there are not differences between a majority of the series but indeed that there are statistically significant differences between a number of the individual series in pair-wise comparisons.

I can also link all the DSGS plots and tables for all time periods if you are interested.

LL, I did not address the following question in my post above. “Also, Are you able to cross reference your proxy data with the ITRDB files? I would be interesting to look at the raw versions in comparison.”

That was a task I attempted when I was hoping to crack case of the disappearance of the 1960 to present cut-off data from the Schweingruber MXD series as shown in Mann(08). I have been able to approximately connect the time and length of the series of some of the 105 series from an index in what I recall is called The North America Tree Ring Index or some such name. The Index presents the raw data and a choice of three standardizations as I recall. I thought if one of those standardizations fit the data for the period it was available in Mann(08) I might be onto something. I think the hang up for me was that the index deals with individual trees while a proxy fits these data together. I need to go back and determine what actually stymied me. Perhaps it was that I thought others, including SteveM, wanted to find that missing information also and if they had not found it I was probably not either.

The link to the tree ring indexes starts here in the link below and you link to a specific researcher, like Schweingruber, from that linked page. I obviously failed to point to the latitude and longitude available from Mann (08) and from the index is information that can used for matching. Schweingruber and Schweingruber/Briffa have 492 chronologies listed.

http://www.ncdc.noaa.gov/paleo/treering.html

Kenneth Fritsch,

I follow with interest your dialogues but I have not the level in English nor in statistics to well understand what you are doing. But just a note, I’m not sure that comparing individual series with global temperature is the right method. Much of what you see comes from the spatial variability. A grid of 5 ° by 5 ° seems to be a good level of analysis.

Another thing, I had an exchange on a parallel thread with Blog Lurker, it is possible that the observed inverse divergence in the nineteenth century comes from a problem with thermometers exposure.

Phi and Kenneth: I found Mann08 proxies (including Schweingruber mxd) which I had downloaded some time ago. I believe these are the same ones used by Kenneth in 124. An example of one of the file names is :schweingruber_mxdabd_grid1″ which suggests that these are abd (age band decomposition) standardized chronologies. The grid referal I think shows that the age band standardized time series were combined into chronologies on a grid scale. Still not sure whether this is an apples to apples comparison with Briffa01, however Kenneth’s results in #124 suggest it must be similar at the very least. A next step is to download grid scale time series of observations and analyse against chronologies with corresponding coordinates.

LL, I obtained a composite of the 105 Schweingruber Proxy minus NH Temperature series by averaging the individual series. The plot is linked below.

I initially did a breakpoint analysis and found breakpoints at 1883 and 1944 which are seen in observing the linked plot. We have three line segments for this composite compared to the two observed in Briffa (01) with one major negative trend (T1) in the middle and 2 positive trends at either end of the series, T3 at the beginning and T2 at the end of the series. I looked at the three segments with the arima(diff(x), order=c(0,0,1) ) and all three had MA=-1.00 indicating trend stationary. The adf.test function in R gave a marginal result for rejecting the null hypothesis of a unit root for T1 and did not reject the null for T2 and T3. I think that the augmented Dickey_Fuller test is conservative and will not reject with small sample sizes. The trends for T1 had slope = -0.015 degrees C per year, P.value=10^-9, AR1=-0.17; for T2 slope = 0.036 degrees C per year, P.value=0.018, AR1=0.17; for T3 slope = 0.023 degrees C per year, P.value=0.24, AR1=0.14.

I also calculated breakpoints for the Briffa (01) subtracted series for the period 1871 to 1960 in order to determine whether the breakpoint initially found for that series at 1901 and using the 18971 to 1994 time period was sensitive to the end date of the series. It was not as I continue to obtain 1901 as the breakpoint.

“I follow with interest your dialogues but I have not the level in English nor in statistics to well understand what you are doing. But just a note, I’m not sure that comparing individual series with global temperature is the right method. Much of what you see comes from the spatial variability. A grid of 5 ° by 5 ° seems to be a good level of analysis.”

You make a good point and that is why I showed the latitude and longitude for the proxy locations. A quick look analysis of the locations and the propensity to see differing trends did not appear to indicate that the differences are related to location and thus local temperatures. Surely the general observation of a trend in the amount of divergence is not dependent on local temperature variations. I think we started here with the difference to NH temperatures because that is what these reconstructions are compared and in some case calibrated. I would agree whole-heartedly that local temperatures are the key, but I do not see those variations changing much of what we have analyzed here.

Kenneth I think you forgot to include a link in #132.

LL at Post # 131:

LL, I think what you say about combining the TR series on a grid level for forming proxies is probably what threw me off track in my search for the missing Scweingruber data. If I resurrect my search I’ll take that into account – thanks for the lead. Is abd (age band decomposition) standardized chronologies what Briffa(01) used?

Sorry about that, LL. The link is below:

I flipped my series around by subtracting observations from Briffa01 and truncated at 1960 to get a proper comparison to Kenneth’s in #136. Similar but definitely not apples and apples.

#135

Yes. I’ll re-read the methods when I get time but I believe they were first combined at the site rather than the grid.

I have not studied this article in detail, but it points to differing chronology methods dealing with TR can dramatically affect the end result in a reconstruction. It specifically talks about Briffa’s handling of the Yamal proxy and how the use of a different methodology by Briffa changed the original results (no hockey stick).

http://lustiag.pp.fi/ClimateFromTreeRings_gb.htm

It would appear that the threads which have been listed in chronological order have been shuffled. It took me awhile to find this thread.

The article in the link above points to another potential parameter for tuning proxy data in order to obtain the “correct” answer – as the excerpts below indicate.

“Climatic interpretation of the Yamal chronology has triggered a lot of discussion for its hockey stick shaped curve. There was no hockey stick in the earlier climatic presentation of the Yamal chronology by Hantemirov et al. 2002 (Abstract). The hockey stick was reported later on, in 2009 Briffa’s paper. Why did Briffa get different results compared to Hantemirov et al.? Some possible reasons can be listed: 1) Different methodologies: Hantemirov applied the so called Corridor Method, but Briffa his fine-tuned RCS method. 2) Cambium age distribution over time (calendar year range) in the Yamal data is unbalanced because the oldest tree-rings appear only during the last couple of hundred years. A good data for building basic RCS ring width – cambium age models requests an evenly distributed data, both over the cambium age and the calendar year. Accepting all the available ring-width data for model processing inevitably results an unbalanced time-related age distribution, which may lead to strange results as seen in the Yamal case.”

“This important view of RCS modelling should be noted: RCS is highly dependent on successful ring width – ring age modelling. Poor data (like in the case of the Yamal modelling problems) may lead to biased models and misleading conclusions. Ignoring the role of ring age in RCS modelling in fact seems to be a major problem in Dendrochronology. The Yamal chronology data is a typical example of a dataset with distorted age distribution over time (calendar year) causing “funny” results.

As seen in the Yamal case (Age Banding – Scatter), tree-ring age distribution over the calendar year axis does not fulfil the age restriction of 70-110 years. All the tree-ring cambium ages in the 1900s exceed 110 years, because the tree-rings of the same living trees since the 1750s were used for filling the 250-yr gap between the youngest subfossils and the present. Considering the use of BFM approach in the Yamal case, it is not possible, because there is not enough data available for passing the Age Band limitation. Arguing on this fact, my strong opinion is: the Yamal data represents a typical tree-ring dataset that never should be applied in RCS modelling! Compare some age-banded models in the poster.”

“The maximum trend or cycle length in the Yamal case would thus be some 130 years. That is not enough for identifying periods like MWP, LIA and MWP. Comparison of the RCS and STS chronologies, however, show that RCS produces a hockey stick shaped curve but the STS method not. Check the poster.”

“Five important principles

Heavy criticism focused on the dendroclimatic use of the Yamal chronology can be considered as a school example of how science, close to important political decision making, should not be performed. Although the Yamal chronology technically is almost flawless and it is perfect for dating purposes, its validity and reliability considering climatic modelling is poor. There are at least five important questions to be answered:

1) What is the validity of this kind of data? What do the trees actually represent? The risk they are exposing something else than climatic responses is obvious.

2) What is the reliability of the sampled data set? Knowing well that variation in tree growth at the Finnish timberline areas, expressed in terms of the coefficient of variation (CV), is 40-50 % (timberline pine in Finnish Lapland over 40 % and Larix in the Yamal peninsula over 45 %), a lot of careful investigation should be paid to sampling. And this is not enough, we definitely have also to also to consider, what kind of predefined ruling to adapt in order to maximise the climatic signal.

3) What research method is most suitable for producing proper conclusions? Considering the two climatic signal exposing methods, STS and RCS, it is crucial to make difference that STS exposes only local annual and decadal variation while RCS can expose also long-term variation. Age Banding should always play an important role in applying the both methods. Because RCS is a very data sensitive method, considerable attention should be paid for analysis of all the relevant data characteristics. If the data set does not fulfil the predefined information requested by the BFM, RCS should not be used for climatic analyses at all. This is actually the case also with the Yamal data set.

4) How to generalize results? Small data sets, as used e.g. in the Yamal case, represent only local climatic conditions, provided that the validity and reliability of the data is high enough. If not, climatic conclusions may be even locally badly misleading.

5) One more thing to be considered: the whole research process, the results and the conclusions should be exposed to public criticism before publishing. As seen every now and then, peer-reviewing does not always guarantee the proper scientific conclusions! Unfortunately. “