# 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))