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

#	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,
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()```