##calculate all correlation coefficients between instrumental temperature data and proxy data
##plot distribution

#read data, remove "year"-column from proxy data
load("data_allproxy1209.Rdata")
load("data_instrument.Rdata")
t.nh <- instrument$CRU_NH_reform
proxy.inst <- allproxy1209[paste(instrument$Year),-1]

#calculate correlation coefficients and their distribution
cc <- sapply(seq(dim(proxy.inst)[2]),function(z) cor(t.nh,proxy.inst[,z],use="complete.obs"))
dcc <- density(cc)

#now fit an AR(p) model to each of the proxies and build a new matrix of pseudo proxies
#where each "true" proxy is replaced by a realization of its AR(p) fit
build.pseudoproxy <- function(data,ar.order) {
    i.na <- which(is.na(data))
    m <- mean(data[-i.na])
    s <- sd(data[-i.na])
    data <- (data-m)/s
    if (ar.order == 0) {
        pproxy <- rnorm(length(data))
    } else {
        arm <- ar(data,order.max=ar.order,aic=F,na.action=na.omit)
        pproxy <- arima.sim(model=list(order=c(ar.order,0,0),ar=arm$ar),n=length(data))
    }
    as.numeric(pproxy*s+m)
}
pseudo.prox.0 <- sapply( proxy.inst, build.pseudoproxy, ar.order=0 )
pseudo.prox.1 <- sapply( proxy.inst, build.pseudoproxy, ar.order=1 )
pseudo.prox.2 <- sapply( proxy.inst, build.pseudoproxy, ar.order=2 )

#calculate correlation coefficients of observations with pseudo proxies and their distributions
cc0 <- apply(pseudo.prox.0,2,function(z) cor(t.nh,z))
dcc0 <- density(cc0)
cc1 <- apply(pseudo.prox.1,2,function(z) cor(t.nh,z))
dcc1 <- density(cc1)
cc2 <- apply(pseudo.prox.2,2,function(z) cor(t.nh,z))
dcc2 <- density(cc2)

#plot distribution of correlation coefficients
#pdf("cc-dist.pdf",width=6,height=5)
jpeg("cc-dist.jpg",width=600,height=500,quality=90)
plot(NULL,type="l",lwd=2,xlim=c(-1,1),ylim=c(0,5),xlab="corrcoef.",ylab="density")
lines(dcc0$x,dcc0$y,lwd=2,col="orange")
lines(dcc1$x,dcc1$y,lwd=2,col="blue")
lines(dcc2$x,dcc2$y,lwd=2,col="green")
lines(dcc$x,dcc$y,lwd=4,col="black")
legend(x=-1,y=4.5,c("proxies","pseudo AR(0)","pseudo AR(1)","pseudo AR(2)"),lwd=c(4,rep(2,3)),
  lty=rep(1,4),col=c("black","orange","blue","green"))
dev.off()