the Air Vent

Because the world needs another opinion

Building a ‘Hokey’ stick

Posted by Jeff Id on October 12, 2008

If this is your first time here recently see this article Will the Real Hockey Stick Please Stand Up?

The first graph of this post is the final result after the 1209 series were magnified and averaged together. Vertical scale is SD units calibrated from 1850-1995. This is my first post which is based on R software.

The point of this post is just to see the types of data going into the latest hokey hockey stick and what shape the groups provide. If you look at all the graphs after the first one, you might see components which add up to the shape of the final graph. Also, you will notice shapes which you can’t find. The fact that you cannot find certain valleys and peaks means very clearly that individual proxies had substantial influence in creating historic trends.

Luterbacher instrument proxies with measured temp data.

Tree ring composite data possibly already HS sorted. Not too many series.

3001 – Conglomeration of marine sediments

4000 -Lake sediments

4001 Rainfall and isotope

5000 – Moisture and grape ripening records – for god’s sakes I bet grape ripening is linear with temperature.

5001-Historic documents!!

6000 – Speleothum – Cave stalactite isotope data.

More cave speleothum data

7000 Isotope ratios in coral

7001 Coral growth (I think) two proxies only

Schweingruber latewood MXD data.  I bet $30 this is our little ice age.  What tweaks me about this one is the fact that the last 38 years are painted in upslope.  (Fake data!)  Without it these series don’t get accepted and I bet no little ice age.

Ice core isotope and precipitation data

Ice core isotope records

9000 – Tree ring width data. 76.7% of the 1209 series

Did you catch the proxy which was created from ripening grapes? Imagine converting that plot to temperature.

12 Responses to “Building a ‘Hokey’ stick”

  1. Andy said

    Erm….if they all measure “temperature” why are they all so different?

    As they are all different, which are truthful regarding “temperature” and which are false and how do you tell?

    Or is it just those that are Hocky stick shaped measure temperature?

  2. Jeff Id said

    Andy,

    You are exactly where I was not that long ago. I have a link which answers some of your questions.

    http://noconsensus.wordpress.com/2008/09/11/ten-things-everyone-should-know-about-the-global-warming-hockey-stick/

    I hope you like it.

  3. Demesure said

    Hi Jeff
    Can you post your R script please ? (don’t forget to enclose it with )
    I would like to change the graphs' title to "Icecores : x / 1029 proxies", "Tree ring width: y / 1029 proxies" ...

  4. Demesure said

    @3 I meant “don’t forget to enclose it with the tag CODE “.

  5. Jeff Id said

    Demesure,

    I tried to upload it today and apparently I mad a mistake and didn’t save the final version of my code. I have something close though so I will fix it and post tonight after work. That’s the main reason I am learning R, so people can review my work easily. C++ would be huge to go through.

    Perhaps you can help me, I needed to make an indexed array of time series so I could join them in a ts.union. The function becomes very slow when the union grows large. I gave up after a couple of hours but what I needed was to address a bunch of time series like this x[i]. I ended up just typing in 4 series like xa, xb, xc, xd. It was a pain but it got the job done.

    You’ll see what I mean later.

  6. Gary said

    Jeff, You are doing amazing stuff. Keep it up. Is there some reason that we can’t just be shown a graph of the actual proxies results without CPS, EIV or any other selection process? I did read your original blog on how the data gets chosen. But I would really like to just see all the data in it’s original form as a graph. Thanks for your efforts.

  7. Jeff Id said

    The original data has values as high as 3000 (tree ring) and others around 1 (pre-calibrated isotope). The only way to look at it together is to standardize. With all 1209 proxy series on a single graph it is just a mess.

    There is an answer though. On Climate Audit, (Link on the right) Steve McIntyre made gif based movies of the proxies which play one after another on the screen in standardized units, he even had a way to download the movie file.

    The posts were around the beginning of September08 and there were several of them. It was a bit jaw dropping for me.

    I hope that helps.

  8. Jeff Id said

    R script used to generate above graph. If anyone knows the correct way to replace chronA,B,C,D into an array like chron[i], it would be appreciated.

    ———————
    library(gplots)
    source(“http://www.climateaudit.org/scripts/utilities.txt”)
    download.file(“http://www.climateaudit.org/data/mann.2008/mann.tab”,”temp.dat”,mode=”wb”);load(“temp.dat”)
    download.file(“http://www.climateaudit.org/data/mann.2008/details.tab”,”temp.dat”,mode=”wb”);load(“temp.dat”)

    codelist=NULL
    for(j in 1:10000)
    {
    temp=(details$code==j)
    K=sum(temp)
    if(K>0)
    {
    codelist[length(codelist)+1]=j
    }
    }

    for(j in 1:length(codelist))
    {
    temp=header$Data.type.code==codelist[j]
    K=sum (temp)
    index=(1:1209)[temp]

    chronA=NULL
    chronB=NULL
    chronC=NULL
    chronD=NULL
    for ( i in 1:K)
    {
    X=g(mann[[ index[i] ]]) ##convert to time series
    temp2=(time(X)>=1850)&(time(X)<=1995) ##collect mann proxy calibration into temp
    m0=mean(X[temp2],na.rm=T)
    sd0=sd(X[temp2],na.rm=T) ##Calc mean and standard deviation
    Xscaled=(scale(X[],center=m0,scale=sd0)) ##Linear scaling of temp
    Xscaled=ts(Xscaled,start=time(X)[1]) ##Cram it back into time series
    if(i<250)
    {
    chronA=ts.union(chronA,Xscaled ) ##stick everything into time series array
    }
    else
    {
    if(i<500)
    {
    chronB=ts.union(chronB,Xscaled ) ##stick everything into time series array
    }
    else
    {
    if(i<750)
    {
    chronC=ts.union(chronC,Xscaled ) ##stick everything into time series array
    }
    else
    {
    chronD=ts.union(chronD,Xscaled ) ##stick everything into time series array
    }
    }
    }
    }

    if(length(chronB)!=0)
    {
    chronA=ts.union(chronA,chronB)
    }
    if(length(chronC)!=0)
    {
    chronA=ts.union(chronA,chronC)
    }
    if(length(chronD)!=0)
    {
    chronA=ts.union(chronA,chronD)
    }

    mannavg=ts(apply(chronA,1,mean,na.rm=T),start=tsp(chronA)[1]) ##average time series

    Y=mannavg;##//ts.union(annual,mannavg) ##cram time series into Y
    f=function(x) filter.combine.pad(x,truncated.gauss.weights(11) )[,2] ##filter functionh
    plot(time(Y),f(Y),col=2,type=”l”,lwd=2,las=1,ylab=””,xlab=”Year”,xlim=c(0,2000), ylim = c(-5,5))
    title(paste(“Average of Proxy Type =”,codelist[j]))
    ##savePlot(paste(“c:/agw/r/Average of Proxy Type =”,codelist[j]),type=”jpg”)
    }

    It should run as is.

  9. Jeff Id said

    When I pasted it the tabs disappeared. Makes it tough to read.

  10. Andy said

    Thanks Jeff, if only I could be where you are, I’m afraid my non-existant computer/maths skills are dwarfed by those here and on other blogs.

    Keep up the good work.

  11. Demesure said

    Jeff,
    I can’t get your code running (“header$Data.type.code”, “g(mann[[ index[i] ]])” and “filter.combine.pad” don’t run) so I tweaked it my way. I don’t know why your chronA..D so don’t know where I can help.

    Now things are much clearer to me : (see my graphs here
    http://pichuile.free.fr/ecad/mann2008 )

    You see strange things when data are not filtered in code 5000 and 8001. And of course, that 927 proxies of 1209 are tree rings.


    library(gplots)
    source("http://www.climateaudit.org/scripts/utilities.txt")
    download.file("http://www.climateaudit.org/data/mann.2008/mann.tab","temp.dat",mode="wb");load("temp.dat")
    download.file("http://www.climateaudit.org/data/mann.2008/details.tab","temp.dat",mode="wb");load("temp.dat")

    proxiname=c("Miscellaneous","Composite","Marine sediments","Miscellaneous","Miscellaneous",
    "Miscellaneous","Miscellaneous","Miscellaneous" ,"Miscellaneous" ,"Miscellaneous" ,"Miscellaneous" ,
    "Tree latewood density","Oxygen isotopes","Miscellaneous","Tree width")

    codelist=as.numeric(levels(details$code))

    for(j in 1:length(codelist))
    {
    y=mann[details$code== codelist[j]]

    chronA=NULL
    for ( dat in y)
    {
    X=ts(dat[,2], start=dat[1,1]) ##convert to time series
    temp2=(time(X)>=1850)&(time(X)<=1995) ##collect mann proxy calibration into temp
    m0=mean(X[temp2],na.rm=T)
    sd0=sd(X[temp2],na.rm=T) ##Calc mean and standard deviation
    Xscaled=(scale(X[],center=m0,scale=sd0)) ##Linear scaling of temp
    Xscaled=ts(Xscaled,start=time(X)[1]) ##Cram it back into time series
    chronA=ts.union(chronA,Xscaled ) ##stick everything into time series array
    }

    Y=ts(apply(chronA,1,mean,na.rm=T),start=tsp(chronA)[1]) ##average time series
    plot(Y,col=2,type="l",las=1,xlab="Year",ylab="values scaled to 1 sigma mean",xlim=c(0,2000), ylim = c(-3,3))
    title(paste(proxiname[j]," (code=",codelist[j], ") / ",length(y)," proxies",sep=""))
    abline(h=seq(-3,3,1), v=seq(0,2000,500),col="gray", lty="dotted", add=T)
    lines(Y,col=2)
    savePlot(paste("c:/agw/Average of Proxy Type =",codelist[j]),type="jpg")
    }
    #

  12. Jeff Id said

    I finally got a chance to play around with your code. Much cleaner. The chronA..D is a method for increasing the speed of the union for the tree ring proxies. It cuts the time down approximately 2^4 or 16 times. I just figured out how to dim the array beforehand. Works great.

    Thanks for the example.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

 
Follow

Get every new post delivered to your Inbox.

Join 148 other followers

%d bloggers like this: