# A reproduction of Lohle & Scafetta 2011 # Read and organise the data from cRU: CRU <- read.table("CRU-T2m.mon.txt") yymm <- sort(rep(CRU[,1],12)) + (rep(1:12,length(CRU[,1]))-0.5)/12 # It seems they osed annual data: t2m <- ma.filt(c(t(CRU[,2:13])),12); N <- length(t2m) # Select the data before 1950 i1950 <- yymm < 1950 n <- sum(i1950) # Construct similar harmonics w20 <- pi/(10*12)*seq(1,sum(i1950),by=1) w60 <- pi/(30*12)*seq(1,sum(i1950),by=1) c20 <- cos(w20); s20 <- sin(w20) c60 <- cos(w60); s60 <- sin(w60) W20 <- pi/(10*12)*seq(1,N,by=1) W60 <- pi/(30*12)*seq(1,N,by=1) C20 <- cos(W20); S20 <- sin(W20) C60 <- cos(W60); S60 <- sin(W60) cal60 <- data.frame(y=t2m[i1950],x1=c60,x2=s60,x3=c20,x4=s20) pre60 <- data.frame(y=t2m,x1=C60,x2=S60,x3=C20,x4=S20) # Regression analysis & prediction fit60m <- lm(y ~ x1 + x2 + x3 + x4,data=cal60) fit60 <- predict(fit60m,newdata=pre60) # Plot the data plot(yymm,t2m,type="l",col="grey") lines(yymm[i1950],t2m[i1950]) lines(yymm,fit60,col="red",lwd=2) dev2bitmap("loehlesscafetta2011redux.png") # Explore other frequencies too, keeping the best-fit amplitudes: dev.new() plot(yymm,t2m,type="l",col="grey") lines(yymm[i1950],t2m[i1950]) A <- rep(NA,50); A1 <- A; A2 <- A for (i in 1:50) { wX <- pi/(i*12)*seq(1,sum(i1950),by=1) cX <- cos(wX); sX <- sin(wX) calX <- data.frame(y=t2m[i1950],x1=cX,x2=sX) fitXm <- lm(y ~ x1 + x2,data=calX) fitXs <- summary(fitXm) fitX <- predict(fitXm,newdata=calX) lines(yymm[i1950],fitX,col="blue",lwd=2) A[i] <- sqrt(fitXs\$coefficients[2]^2 + fitXs\$coefficients[3]^2) A1[i] <- sqrt((fitXs\$coefficients[2]-fitXs\$coefficients[5])^2 + (fitXs\$coefficients[3]-fitXs\$coefficients[6])^2) A2[i] <- sqrt((fitXs\$coefficients[2]+fitXs\$coefficients[5])^2 + (fitXs\$coefficients[3]+fitXs\$coefficients[6])^2) } lines(yymm,fit60,col="red",lwd=2) dev2bitmap("loehlesscafetta2011-manycycles.png") # Plot the amplitudes: dev.new() plot(c(2,100),range(c(A1,A2)),type="n",lwd=2, main="Amplitude of best-fit cycle", xlab="Periodicity (Year)",ylab="Amplitude (C)") polygon(c(seq(2,100,by=2),reverse(seq(2,100,by=2))), c(A1,reverse(A2)),col="grey",border="grey") lines(seq(2,100,by=2),A,lwd=3) grid() dev2bitmap("loehlesscafetta2011amplitudes.png") #----------------------- # Construct a synthetic time series: 10.000 year long # Set up data matrix representing random amplitudes, and split between # cosine and sine: A <- rnorm(1000); dim(A) <- c(500,2) # Set up matrix and vector holding the results; x holds the random curve and # y the different cosine+sinus components. x <- rep(0,10000*12); # y <- matrix(rep(x,5000),5000,10000*12) # Define the basis frequency of the cosine and sine curves: the largest # frequency corresponds to one cycle for the whole curve: w <- pi*seq(0,2,length=10000*12) # Carry out the construction of the random curves: for (i in 1:500) { # y[i,] <- ( A[i,1]*cos(i*w) + A[i,2]*sin(i*w))/sqrt(i) # x <- x + y[i,] y <- ( A[i,1]*cos(i*w) + A[i,2]*sin(i*w))/sqrt(i) x <- x + y } x <- ma.filt(x,12) dev.new() plot(x,type="l",lwd=2, main="10.000yr synthetic series") dev2bitmap("loehlesscafetta2011syntheticseries.png",res=150) dev.new() plot(c(1,n),range(x,na.rm=TRUE),type="n", main=paste("Looking at",n/12,"year seqments"), xlab="month",ylab="") for (i in seq(1,10000*12,b=n)) { CAL60 <- data.frame(y=x[i:(i+n-1)],x1=c60,x2=s60,x3=c20,x4=s20) # Regression analysis & prediction FIT60m <- lm(y ~ x1 + x2 + x3 + x4,data=CAL60) FIT60 <- predict(FIT60m,newdata=CAL60) lines(x[i:(i+n-1)],lwd=2,col="grey") lines(FIT60,col="red",lwd=2) } dev2bitmap("loehlescafetta2011sequences.png",res=150)