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

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 memory

lent=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 series

stat = 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=NULL

for(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+msig

x=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.