Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # 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))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement