#	To simulate two AR(N) series, we use the built-in arima functions.  
# This simulation will show how autocorrelated but independent 
# time-series will often show spurious correlations.

n.simulations <- 100
n.ar.coefficients <- 100
n.points <- 100

#	Note this is a three-dimensional array.
#	The third dimension lists metadata: 
#  R Squareds, Adjusted R Squareds, correlations (absolute value).
#  the fourth spot in the third dimension is reserved for p-values.
metaData <- array(0, dim = c(n.ar.coefficients, n.simulations, 4))

ar.coefficients <- seq(
	from = -0.999, 
	to = 0.999, 
	length.out = n.ar.coefficients)

#	This may take some time to run.  It is going through 10,000 
# regressions (with n.ar.coefficients and n.simulations both at 100).
for (i in 1:n.ar.coefficients) {
	for (j in 1:n.simulations) {
		coef_i <- list(ar = c(ar.coefficients[i]))
		seriesA <- arima.sim(model = coef_i, n = n.points)
		seriesB <- arima.sim(model = coef_i, n = n.points)

		model <- lm(seriesB ~ seriesA)
		metaData[i, j, 1] <- summary(model)$r.squared
		metaData[i, j, 2] <- summary(model)$adj.r.squared
		metaData[i, j, 3] <- sqrt(summary(model)$r.squared)
		#metaData[i, j, 4] <- getpvalue(model)
	}
}


#	Now set up the pdf and some graphical parameters.
png("independentARseries.png")
point.size = 0.3
point.type = 19
line.width = 3


plot(x = ar.coefficients, 
	y = metaData[,1,1],
	type = "points",
	pch = point.type,
	cex = point.size,
	xlab = "Autocorrelation Coefficient",
	ylab = "R-Squared",
	xlim = c(-1,1),
	ylim = c(0,1))

for (i in 2:n.simulations)
{
	points(x = ar.coefficients, y = metaData[,i,1], pch = point.type,
		cex = point.size)
}

#	Now we add lines at the 5 and 1% significance levels.  
significances <- c(0.95, 0.99)
f.values <- sapply(significances, function(x) {qf(x, 1, 98)})
denominator.dof <- n.points - 2
R2s <- (f.values / denominator.dof) / (1 + f.values / denominator.dof)
abline(a = R2s[1], 0, lty = 1, lwd = line.width)
abline(a = R2s[2], 0, lty = 2, lwd = line.width)

legend("center", c("p = 0.05", "p = 0.01"), 
                  lty = c(1,2), lwd = line.width)

#	Close the graphics device.
dev.off()