#run using R CMD BATCH amo_regression.r 
#outputs to fig*.png and amo_regression.r.Rout

#options
nsims=1000 #Number of Monte Carlo simulations.
power = 5 #Nonlinearity of human influence, 1 = linear.
recent = 1979 #Calculate trends after this year.
xunits="year"
textsize=1.4
titlesize=1.8
col1="red"
col2="green"
col3="black"
col4="blue"
pch1=20#points
par(bg = "white")
desired_width=800
desired_height=800
desired_res=120

#Define synthetic human, natural influences.
t = 1856:2011
n = length(t)
human = (t-t[1])^power
human = 0.8*human/human[n]
nature = 0.2*cos(2*pi*(t-2000)/70) #70 yr AMO which peaks in year 2000.
#Plot synthetic human, natural influences.
png(file="fig1.png",width=desired_width,height=desired_height,res=desired_res)
title="Synthetic human and natural influences"
yunits="deg. C"
ylims=c(min(c(human,nature)),max(c(human,nature)))
plot(t,human,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col=col1,ylim=ylims)
points(t,nature,type="o",pch=pch1,lwd=2,col=col2)
legend(t[1],ylims[2],c("Human", "Natural"), pch = c(pch1,pch1), col = c(col1,col2), lty = c(1,1))
dev.off
trends = 1:nsims
quadratics = 1:nsims
for(i in 1:nsims){

  #Define synthetic global, N. Atlantic temperatures.
  global = human + nature + rnorm(t,mean=0,sd=0.2)
  n_atlantic = global + rnorm(t,mean=0,sd=0.1)
  #Plot synthetic global, N. Atlantic temperatures.
  if(i==1){
    png(file="fig2.png",width=desired_width,height=desired_height,res=desired_res)
    title="Synthetic global, N. Atlantic temperatures"
    ylims=c(min(c(global,n_atlantic)),max(c(global,n_atlantic)))
    plot(t,global,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col=col3,ylim=ylims)
    points(t,n_atlantic,type="o",pch=pch1,lwd=2,col=col1)
    legend(t[1],ylims[2],c("Global", "N. Atlantic"), pch = c(pch1,pch1), col = c(col3,col1), lty = c(1,1))
    dev.off
  }
  #Define linearly-detrended AMO.
  n_atlantic_trend = lm(n_atlantic~t)
  amo = n_atlantic - coef(summary(n_atlantic_trend))[2,1]*(t-t[1])
  #Plot detrended AMO.
  if(i==1){
    png(file="fig3.png",width=desired_width,height=desired_height,res=desired_res)
    title="Synthetic linearly-detrended AMO index"
    #title="Synthetic ideally-detrended AMO index"
    ylims=c(min(amo),max(amo))
    plot(t,amo,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col=col3,ylim=ylims)
    dev.off
  }
  #Regress global temperatures against human influence and AMO.
  human_p = human - mean(human) #Remove mean so it's handled by the intercept.
  amo_p = amo - mean(amo)
  nature_p = nature - mean(nature)
  regression = lm(global~human_p+amo_p)
  summary(regression)
  est_human = coef(summary(regression))[2,1]*human_p
  est_human = est_human - est_human[1] #So that est_human starts at zero too.
  #Calculate recent trends.
  recent_t = t[which(t>recent)]
  recent_human = human[which(t>recent)]
  recent_human_trend = lm(recent_human~recent_t)
  true_summary = summary(recent_human_trend)
  true_human_trend = coef(true_summary)[2,1]
  #Plot residuals.
  if(i==1){
    png(file="fig4.png",width=desired_width,height=desired_height,res=desired_res)
    title="Residuals"
    residuals = resid(regression)
    plot(t,residuals,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col=col3)
    dev.off
  }
  #Add residual back.
  est_human2 = est_human + residuals
  recent_est_human2 = est_human2[which(t>recent)]
  recent_est_human2_trend = lm(recent_est_human2~recent_t)
  est_summary2 = summary(recent_est_human2_trend)
  est_human2_trend = coef(est_summary2)[2,1]
  true_blurb = sprintf("True human. Post-%d trend = %.2f %s/decade",recent,true_human_trend*10,yunits,xunits) 
  est_blurb2 = sprintf("Est. human. Post-%d trend = %.2f %s/decade",recent,est_human2_trend*10,yunits) 
  #Plot true and estimated human influences.
  if(i==1){
    png(file="fig5.png",width=desired_width,height=desired_height,res=desired_res)
    title="True, estimated human influences"
    ylims=c(min(c(human,est_human2)),max(c(human,est_human2)))
    plot(t,human,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title,col=col3,ylim=ylims)
    points(t,est_human2,type="o",pch=pch1,lwd=2,col=col1)
    recent_trendline = est_human2_trend*(recent_t-recent_t[1])
    recent_trendline = recent_trendline - mean(recent_trendline) +mean(recent_est_human2)
    points(recent_t,recent_trendline,type="o",pch=pch1,lwd=2,col=col4)
    legend(t[1],ylims[2],c(true_blurb, est_blurb2), pch = c(pch1,pch1), col = c(col3,col1), lty = c(1,1))
    dev.off
  }
  trends[i] = est_human2_trend*10
  #Calculate estimated quadratic term.
  est_human2_quadratic = lm(est_human2~t + I(t^2))
  est_quadratic = coef(summary(est_human2_quadratic))[3,1]
  quadratics[i] = est_quadratic
}
trends
png(file="fig6.png",width=desired_width,height=desired_height,res=desired_res)
title = sprintf("Est. human trend since %d",recent) 
xunits = "deg.C/decade" 
hist(trends,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
dev.off
quadratics
png(file="fig7.png",width=desired_width,height=desired_height,res=desired_res)
true_human_quadratic = lm(human~t + I(t^2))
true_quadratic = coef(summary(true_human_quadratic))[3,1]
title = sprintf("Est. human quadratic term since %d",t[1]) 
quadratic_blurb = sprintf("True quadratic = %.3e",true_quadratic) 
quadratic_blurb
xunits = "deg.C/year^2" 
hist(quadratics,density=20,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,main=title)
#text(min(quadratics),1,quadratic_blurb,cex=1,pos=4)
dev.off