the Air Vent

Because the world needs another opinion

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


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 )

Google+ photo

You are commenting using your Google+ 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 )

Connecting to %s

 
%d bloggers like this: