# load data

# as in all analyses for this course, you may have to change
# file names to reflect full path on your machine
T<-scan("2004GL021750-NOAMER.1400.txt",skip=1)
T<-matrix(T,nrow=581,ncol=72,byrow=TRUE)
tree<-T[1:581,3:72]
tt<- c(1400:1980)
# also load global anomaly data
gl<-matrix(scan('hadcrut3gl.txt'),ncol=27,byrow=T)
gl2<-gl[,c(1,14)]


# MBH- uncentered PCs  GRL Data, Caspar Ammann's code
# =======
m.tree<-apply(tree[503:581,],2,mean)
s.tree<-apply(tree[503:581,],2,sd)
xstd<-scale(tree,center=m.tree,scale=s.tree)
fm<-lm(xstd[503:581,]~c(1:nrow(xstd[503:581,]))) #fit linear model to calibrati$
sdprox<-sd(fm$residuals)                         #get residuals
mxstd<-scale(xstd,center=FALSE,scale=sdprox)     #standardize by residual stddev
w<-svd(mxstd)                                    #SVD
PCs<-w$u[,1:5]                                   #Retain only first 5 PCs for c$
m.pc<-apply(PCs[503:581,],2,mean)
s.pc<-apply(PCs[503:581,],2,sd)
PCs<-scale(PCs,center=m.pc,scale=s.pc)

# set up for 25-year MA calculation
w1<-c(1:13,12:1)
w1<-w1/sum(w1)
wt25<-matrix(0,ncol=581,nrow=557)
for(i in 1:557){wt25[i,i+0:24]<-w1}
wtss<-sum(w1*w1)

# Plot MBH reconstruction with 25-year MA
# original MBH plot
postscript('newfig1.ps',horizo=T)
par(mfrow=c(1,1))
plot(tt,PCs[,1],type='l',col='red',xlab='Year',
ylab='Reconstructed Temperature',main='Original MBH Method')
xnew<-wt25 %*% PCs[,1]
lines(tt[13:569],xnew,col='blue')
dev.off()

# Ordinary PC analysis
PCs.2<- princomp( tree, cor=TRUE)$scores

# plot first centered PC
postscript('newfig2.ps',horizo=T)
par(mfrow=c(1,1))
x0<-(-PCs.2[,1]+mean(PCs.2[503:581,1]))/sqrt(var(PCs.2[503:581,1]))
plot(tt,x0,type='l',col='red',xlab='Year',
ylab='Reconstructed Temperature',main='First PC, Centered to 1902-1980')
xnew<-wt25 %*% x0
lines(tt[13:569],xnew,col='blue')
dev.off()

# calculate PC regression and plot for 6 values of K
# load MASS package first
nk0<-c(1,2,5,10,20,50)
name1<-c('K=1','K=2','K=5','K=10','K=20','K=50')
postscript('newfig3.ps',horizo=T)
par(mfrow=c(3,2))
for(k in 1:6){
nk<-nk0[k]
reg<-lm(gl2[53:131,2]~PCs.2[503:581,1:nk])
xnew<-wt25 %*% PCs.2[,1:nk]
#xnew<-wt15 %*% PCs.2[,1:nk]
if(nk>1){gl3<-reg$coef[1]+PCs.2[,1:nk] %*% reg$coef[2:(nk+1)]
gl4<-reg$coef[1]+xnew %*% reg$coef[2:(nk+1)]}
if(nk==1){gl3<-reg$coef[1]+PCs.2[,1:nk] * reg$coef[2]
gl4<-reg$coef[1]+xnew * reg$coef[2]}
XMAT<-cbind(rep(1,79),PCs.2[503:581,1:nk])
VMAT<-ginv(t(XMAT)%*%XMAT)
predse<-rep(0,557)
for(i in 1:557){predse[i]<-sqrt(wtss+c(1,xnew[i,]) %*% VMAT %*% c(1,xnew[i,]))*summary(reg)$sigma}
ymax<-max(c(gl3,gl4+1.645*predse))
ymin<-min(c(gl3,gl4-1.645*predse))
plot(tt,gl3,type='l',ylim=c(ymin,ymax),col='red',xlab='Year',
ylab='Reconstructed Temperature',main=name1[k])
lines(tt[13:569],gl4,col='blue')
lines(tt[13:569],gl4+1.645*predse,col='green')
lines(tt[13:569],gl4-1.645*predse,col='green')
}
dev.off()