## Build your own HS, Using the Same Data and Math as Michael Mann

Posted by Jeff Id on October 19, 2008

Ok, sorry about the delay in posting. I have been trying to recreate my correlation C++ software in R and it takes about 3X longer while I’m learning.

This post lets everyone see how the hockey stick CPS graphs are made. You can put your own patterns in and see how the actual highly random data which is not temperature can make any graph you want and can fool people into thinking it might be temperature.

If you don’t want to do it yourself you can see this link below.

## Will the Real Hockey Stick Please Stand Up?

I made an R script which generates plots like the one below. R is a open source FREE programming software which is easily obtained through a google search for “R”

I have commented the script quite a bit. I can’t fix the damn quote marks. It copies them into italicized quotes. You need to copy and paste the code into an R script and then correct all the quote marks. Sorry, I can’t get wordpress to cooperate.

After you copy and fix the quotes, you can modify the section with ######### below using your own functions or commenting out different ones provided, change all the values in this section and make your own hockey sticks using M08 data.

I have examples of positive and negative linear slopes and sin functions. You can adjust the end time and length for each pattern you want to look for.

##TEMPERATURE CORRELATION

#this calculates statlist with 6 items in list#Mann PRoxy Data and Info

method=”online”

if (method==”offline”) {load(“d:/climate/data/mann08/mann.tab”)

load(“d:/climate/data/mann08/mannorig.tab”)

load(“d:/climate/data/mann08/idorig.tab”)

names(mannorig)=idorig$idorig

load(“d:/climate/data/mann08/details.tab”)

source(“d:/climate/scripts/mann.2008/instrumental.txt”)} else {

url=”http://www.climateaudit.org/data/mann.2008″

download.file(file.path(url,”mann.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)

download.file(file.path(url,”mannorig.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)

download.file(file.path(url,”details.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)

download.file(file.path(url,”idorig.tab”),”temp.dat”,mode=”wb”);load(“temp.dat”)

names(mannorig)=idorig$idorig

source(“http://www.climateaudit.org/scripts/mann.2008/instrumental.txt”)}#notes: details.tab here incorporates lat-longs from rtable.

##FUNCTIONS

g=function(X) ts(X[,2],start=X[1,1])

use0=”pairwise.complete.obs”

absmax= function(x) if(sum(!is.na(x))>0) unique(x[(abs(x)==max(abs(x),na.rm=T))&!is.na(abs(x))])[1] else NA

library(gplots)

source(“http://www.climateaudit.org/scripts/utilities.txt”)

f=function(x) filter.combine.pad(x,truncated.gauss.weights(11) )[,2]########### mannorig is 1357 infilled series,,, mann is the 1209 series after infilling ##########

mano=mannorig #copy non-infilled series to different memory

#mano=mann #copy infilled series to different memorylent=length(mano) ##lent = number of series

#filter all series using gaussian CA filter

for(i in 1:lent)

{

mano[[i]][,2]=f(mano[[i]][,2])

}#apply functions to look for change values in this section

#change values here and run from here down only

#############################################################################################

#############################################################################################

#############################################################################################

################## function generation variables

endyr=1995 #end year for function correlation 100-2000

leng=200 #lenght of function to correlate

cycles=4 #cycles of sine wave

#################################### functions – choose by commenting out all but one

m=c(1:leng)/-leng ##0 to -1 linear slope

#m=c(1:leng)/leng ##0 to -1 linear slope

#m=sin( (1:leng)/(leng/cycles)*pi) #sin wave with # cycles

#m=-sin( (1:leng)/(leng/cycles)*pi) #-sin wave with # cycles

##################r=.3 #sorting correlation criteria

####################

####################

#############################################################################################

#############################################################################################

#############################################################################################leng=length(m) ## recalc leng

m=ts(m,endyr-leng) ## convert to time seriesstat = dim(lent)

for(i in 1:lent)

{

y=g(mano[[i]])

y=ts.union(m,y)

tim = (time(y)>=endyr-leng & time(y)<=endyr)

stat[i]=cor(y[tim,],use=use0)[1,2]

}

pass= stat>r # build logical r array for pass fail#######calculate SD for sorting function m — used in CPS to match standard

#######deviation and offset of all proxies to the function sd and offset

msig=mean (m,na.rm=T)

ssig=sd(m,na.rm=T)#clear values

y=NULL;

x=NULL;

Z=NULLfor(j in 1:lent)

{

if(pass[j]==TRUE)

{

y=(g(mano[[j]])) #make y a time series

tim = (time(y)>=endyr-leng&(time(y)<=endyr)) #use only calibration range

m0=mean(y[tim],na.rm=T)

sd0=sd(y[tim],na.rm=T)

################# CPS method for rescaling according to function calibration area #######################

yscale=y-m0

yscale=yscale/sd0

yscale=yscale*ssig

yscale=yscale+msigx=ts.union(x,ts (yscale,time(y)[1])) # mash everything into single time series array

# plot(time(y)[1]) #debug plot

}

}

clip = (time(x)>=0 & time(x)<=1995) #clip years to 0 to 1995 for plot

Z=apply(x,1,mean,na.rm=T) #average all ts

Z=ts(Z,time(x)[1]) #make Z a ts###plot

plot(Z[clip],xlim=c(0,2000),ylim=c(-2,2),type=”l”,xlab=”Year”, ylab=”Degrees C”)

lines(m,col=”red”)

smartlegend(x=”left”,y=”bottom”,inset=0,fill=1:2,legend=c(“Mann 2008 Proxy Average”,”Sort Pattern”) ,cex=.7)

title(paste(“M08 Data Sort r >”,r, “\nPercent Prox Used =”,round(sum(pass)/lent*100,2),”%”))

savePlot(“c:/agw/Find Correlation”,type=”jpg”)

Have fun.

## Leave a Reply