# A script that compares a 20+60-yr fit to 100-yr sequences against a # 100,000yr time series. # Test of Loehle & Scafetta paper. # R.E. Benestad # 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); N <- length(x) # 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 <- ( A[i,1]*cos(i*w) + A[i,2]*sin(i*w))/sqrt(i) x <- x + y } x <- ma.filt(x,12) # Define harmonics for entire series: 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) cal <- data.frame(y=x,x1=c60,x2=s60,x3=c20,x4=s20) # Regression analysis & prediction fitm <- lm(y ~ x1 + x2 + x3 + x4,data=cal) s0 <- sd(predict(fitm)) # Define harmonics for the sequences: n <- 1200 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) s <- rep(NA,100) ii <- 0 for (i in seq(1,10000*12,b=n)) { CAL <- data.frame(y=x[i:(i+n-1)],x1=c60,x2=s60,x3=c20,x4=s20) # Regression analysis & prediction FITm <- lm(y ~ x1 + x2 + x3 + x4,data=CAL) FIT <- predict(FITm) ii <- ii + 1 s[ii] <- sd(FIT)/s0 } hist(s); print(summary(s))