Climate Audit

Steve McIntyre

Here is implementation of UC’s code in R:

##LOAD MBH data
#load MBHsmooths1.txt # http://www.climateaudit.info/wp-content/uploads/2009/11/mbhsmooths1.txt
## MBHsmooths1 = [ Year AnnualRecon Instrumental MB98SmoothTrick MBH98SmoothNoTrick MBH99SmoothTrick MBH99SmoothNoTrick ]

MBH=MBHsmooths1=read.table(“http://www.climateaudit.info/wp-content/uploads/2009/11/mbhsmooths1.txt”)
names(MBH)=names(MBHsmooths1)=c(“Year”, “AnnualRecon”, “Instrumental” ,”MBH98SmoothTrick”, “MBH98SmoothNoTrick”,
“MBH99SmoothTrick”, “MBH99SmoothNoTrick”)
dim(MBH) #999 7
MBH[1,]

#[B98,A98]=butter(10,2/50); # ’50 year lowpass’

library(signal)
bf=butter(10,2/50); # ’50 year lowpass’
B98=bf$b;A98=bf$a
B98 # [1] 6.500312e-13 6.500312e-12 2.925140e-11

#in98=MBHsmooths1(401:981,2); # Annual Recon MBH98,
##source:
## http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt ,
## col 2

in98= MBHsmooths1[401:981,”AnnualRecon”];
length(in98) #581

#Trick=MBHsmooths1(982:996,3) #; # pad with instrumental (1981..1995) ..
## source:
## http://www.nature.com/nature/journal/v430/n6995/extref/FigureData/nhmean.txt
## col 3

Trick=MBHsmooths1[982:996,”Instrumental”]#; # pad with instrumental (1981..1995) ..
length(Trick) #15

##function flipud
#tfilt=conv(E,flipud(E))/M;
#Oh, I forgot to explain the flipud-thing:
#Suppose E=[1;2;3;4;5] (a column vector)
# Then flipud simply reverses the values (each column), i.e. flipud(E)=[5;4;3;2;1] and

flipud=function(x) rev(x)

# conv(E,flipud(E)) is
# 5 14 26 40 55 40 26 14 5

E=1:5
conv(E,flipud(E))

#They should be convolved
#tfilt=conv(E,flipud(E))/M;
#Reading the point 145 in the document I referred that would be in R something like tfilt=convolve(E,E,type=’open’)/M

convolve(E,E,type=”open”)
# 5 14 26 40 55 40 26 14 5

#out98t=flipud(filter(B98,A98,flipud( [filter(B98,A98, [in98;Trick;zeros(100,1)] )] ))); #..then smooth
#% Note that filter.m initializes with zeros
#out98t=out98t(1:596);
#out98tr=out98t(26:576); % smoothed reconstruction with the trick
#flipud(filter(B98,A98,flipud( [filter(B98,A98, c(in98,Trick, rep(0,100) ) )] )));
x= c( in98,Trick )
#flipud(filter(B98,A98,flipud( [ filter(B98,A98, c(in98,Trick, rep(0,100) ) ) ] )));
# y= filter(B98,A98, x )
out98t= y= rev (filter(B98,A98,rev ( filter(B98,A98,c(rep(0,100),x,rep(0,100) ) ) ) ))

out98t=flipud(filter(B98,A98, flipud( c(filter(B98,A98, c(in98,Trick,rep(0,100)) )) ))); #..then smooth
length(out98t)
out98t=ts(out98t[1:596],start=1400)
out98tr=ts( out98t[26:576],start=1425) # % smoothed reconstruction with the trick

plot(1400:1980,in98,type=”l”,col=”grey50″,ylim=c(-.8,.5))
lines(out98tr,lty=2,lwd=3)

cbind(out98tr ,MBH[(1425:1975)-999,”MBH98SmoothTrick”]) #1425 to 1975
range( out98tr -MBH[(1425:1975)-999,”MBH98SmoothTrick”],na.rm=TRUE)
# -0.000729800 0.001258341

# no trick
#out98nt=flipud(filter(B98,A98,flipud([filter(B98,A98,[in98;zeros(115,1)])]))); %..then smooth
#out98nt=out98nt(1:596);
#out98ntr=out98nt(26:576); % smoothed reconstruction without the trick

out98nt=flipud(filter(B98,A98, flipud( c(filter(B98,A98, c(in98,rep(0,115)) )) ))); #..then smooth
out98nt=ts(out98nt[1:596],start=1400)
out98ntr=ts( out98nt[26:576],start=1425) # % smoothed reconstruction with the trick

range( out98ntr -MBH[(1425:1975)-999,”MBH98SmoothNoTrick”],na.rm=TRUE)
# -5.796738e-05 1.489949e-04

#% read image, http://uc00.files.wordpress.com/2010/05/mbh985b.png
#A=imread(‘mbh985b.png’); % from MBH98, copied for scientific purposes
#XS=1.1691; % x-scale
#YS=-365.7848; % y-scale
#XB=-1.5749e+003; % x-bias
#YB=281.8316; % y-bias
#Year=(1400:1995)’;
#close all
#image(A)

library(png)
download.file(“http://uc00.files.wordpress.com/2010/05/mbh985b.png”,”temp.png”,mode=”wb”)
img=readPNG(“temp.png”)
dim(img)
#610 781 3

XS=1.1691; # % x-scale
YS=-365.7848; #% y-scale
XB=-1.5749e+003;# % x-bias
YB=281.8316;# % y-bias
Year=(1400:1995);

library(grDevices)
par(mar=c(0,0,0,0))
par(bg=”lightblue”)
plot(MBH[,1],xlim=c(1300,2025),ylim=c(-.8,.5),type=”n”,xaxs=”i”,axes=FALSE)
rasterImage(img , 1320, -.65, 2025, .5)
lines(1275+(Year-1300)*1.05,out98nt*.7-.03,col=3,lwd=1)
lines(1275+(Year-1300)*1.05,out98t*.7-.03,col=2,lwd=1)

#temp=Year>=1980
#lines(1275+(Year[temp]-1300)*1.05,out98t[temp]*.7-.03,col=2,lwd=2)
temp=Year>=1940
lines(1275+(Year[temp]-1300)*1.05,out98nt[temp]*.7-.03,col=”yellow”,lwd=2)