# The code below compares three approaches to autocorrelated data
# 1. Just take the mean of the raw data
# 2. Adjust the raw data to remove autocorrelation
# 3. Trim the raw data to take only every second data point

constant <- 2 # Set constant to whatever

# Make the values below divisible by 2 or it won't work
max_iterations <- 5000
data_set_size <- 100

rawmeanz <- vector(length=max_iterations)
rawsdz <- vector(length=max_iterations)
adjmeanz <- vector(length=max_iterations)
adjsdz <- vector(length=max_iterations)
trimmeanz <- vector(length=max_iterations/2)
trimsdz <- vector(length=max_iterations/2)

for(i in 1:max_iterations)
{
  x <- filter(rnorm(data_set_size), filter=rep(1,3), circular=TRUE)+constant

  rawmeanz[i] <- mean(x)
  rawsdz[i] <- sd(x)
  trimmeanz[i] <- mean(x[seq(2, data_set_size, 2)])
  trimsdz[i] <- sd(x[seq(2, data_set_size, 2)])

  model <- arima(x,order=c(1,0,0))

  adjres <- residuals(model) + coef(model)[2]

  # coef(model)[2] is the unconditional mean, although
  # inexplicably the output says 'intercept'.
  # See http://www.stat.pitt.edu/stoffer/tsa3/Rissues.htm Issue 1
  # for a partial explanation of how this happened.

  adjmeanz[i] <- mean(adjres)
  adjsdz[i] <- sd(adjres)

}

sprintf("Mean error of raw data is %f",mean(abs(constant-rawmeanz)))
sprintf("Mean error of adjusted data is %f",mean(abs(constant-adjmeanz)))
sprintf("Mean error of trimmed data is %f",mean(abs(constant-trimmeanz)))

sprintf("Mean SD of raw data is %f",mean(rawsdz))
sprintf("Mean SD of adjusted data is %f",mean(adjsdz))
sprintf("Mean SD of trimmed data is %f",mean(trimsdz))