##CHEFEN results

   ##LOAD DATA
	url<-"ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean/by_contributor/mann1998/corrigendum2004/figuredata"
	loc<-file.path(url,"fig7-co2.txt")
	co2<-read.table(loc) #1610 1995
	co2<-ts(co2[,2],start=co2[1,1])
	corrs<-read.table(file.path(url,"fig7-corrs.txt"))
	corrs<-ts(corrs[,2:ncol(corrs)],start=corrs[1,1])
	x<-read.table(file.path(url,"fig7-nh.txt"),skip=3)
	nh<-ts(x[,2:ncol(x)],start=x[1,1]) #old series
	x<-read.table(file.path(url,"fig7-solar.txt"),skip=0)
	solar<-ts(x[,2:ncol(x)],start=x[1,1]) #old series
	x<-read.table(file.path(url,"fig7-volcanic.txt"),skip=0)
	volcanic<-ts(x[,2:ncol(x)],start=x[1,1]) #old series
	dimnames(corrs)[2][1:3]<-c("co2","solar","volcanic")
	tsp(corrs)#[1] 1710 1895    1
	
	test<-ts.union(nh, solar,co2,volcanic) #1610 1995
	test<-data.frame(test)

##WINDOW
	N<-100
	K<-(1995-N)-(1609+N)

##PARTIAL CORRELATIONS
	corr.ver.partial<-array(NA,dim=c((1995-N)-(1609+N),3))
	dimnames(corr.ver.partial)[[2]]<-c("solar","co2","volcanic")
	for (k in 1:K){
		index<- k+(1:(2*N)) 
		corry<-cor(test[index,])
		corr.ver.partial[k,]<-corry[2:4,1]
		}


#MULTIPLE CORRELATION
	corr.ver.ols.scaled<-array(NA,dim=c((1995-N)-(1609+N),3))
	dimnames(corr.ver.ols.scaled)[[2]]<-c("solar","co2","volcanic")

	for (k in 1:K){
		index<- k+(1:(2*N)) 
		test.scaled<-data.frame(scale(test[index,]))
		fm<-lm(nh~solar+co2+volcanic,data=test.scaled)
		corr.ver.ols.scaled[k,]<-fm$coef[2:4]
		}


	lwd0<-1
	nf<-layout(array(1:3,dim=c(3,1)),heights=c(1.1,1,1.3))	
	ylim0<-c(-.5,.8)
	#ylim0<-range(c(corrs[,"solar"],corr.ver[,"solar"]))
	par(mar=c(0,3,1,1))
	plot(1710:1895,corrs[,"solar"],type="l",xlim=c(1610,1995),xlab="",ylab="",axes=FALSE,ylim=ylim0,lwd=lwd0)
	axis(side=1,labels=FALSE);axis(side=2,las=1);box()
	lines((1610+N):(1995-N),corr.ver.partial[,"solar"],col="red")
	lines((1610+N):(1995-N),corr.ver.ols.scaled[,"solar"],col="blue")
	text(1610+N,.6,"Solar",font=2,pos=4);abline(h=0)

	par(mar=c(0,3,0,1))
	#ylim0<-range(c(corrs[,"co2"],corr.ver[,"co2"],corr.ver.ols.scaled[,"co2"]))
	plot(1710:1895,corrs[,"co2"],xlim=c(1610,1995),type="l",xlab="",ylab="",axes=FALSE,ylim=ylim0,lwd=lwd0)
	axis(side=1,labels=FALSE);axis(side=2,las=1);box()
	lines((1610+N):(1995-N),corr.ver.partial[,"co2"],col="red")
	lines((1610+N):(1995-N),corr.ver.ols.scaled[,"co2"],col="blue")
	text(1610+N,.6,"CO2",font=2,pos=4);abline(h=0)

	par(mar=c(3,3,0,1))
	#ylim0<-c(-.4,.6)
	plot(1710:1895,corrs[,"volcanic"],type="l",xlim=c(1610,1995),xlab="",ylab="",axes=FALSE,ylim=ylim0,lwd=lwd0)
	axis(side=1);axis(side=2,las=1);box()
	lines((1610+N):(1995-N),corr.ver.partial[,"volcanic"],col="red")
	lines((1610+N):(1995-N),corr.ver.ols.scaled[,"volcanic"],col="blue")
	text(1610+N,.6,"Volcanic",font=2,pos=4);abline(h=0)