### Simulate means and covariances of processes of an AR(1) process

sim.ar <- function(input, coef) {
	filter(input, coef, method="recursive")
}

cov.ar <- function(alpha, len) {
	alpha^{0:(len-1)}/(1-alpha^2)
}

corr.ar <- function(alpha, len) {
	alpha^{0:(len-1)}
}

generalized.difference <- function(xt, alpha) {
	n<-length(xt)
	c(xt[1], xt[-1]-alpha*xt[-n])
	}

	
mean.na <- function(x) { mean(x,na.rm = TRUE) }
sd.na   <- function(x) { sd  (x,na.rm = TRUE) }


# ---
# Estimate an AR(1) and then estimate the mean and correlation
#

n     <- 200;
mu    <- 100;
alpha <- 0.7;
xt <- mu + sim.ar(rnorm(n),alpha)

plot(xt); abline(h=mu, col="gray")


# Fit the regression of xt on lag
plot(xt[-n],xt[-1])
regr <- lm(xt[-1]~xt[-n]); summary(regr)
abline(regr, col="red")


# Estimated correlation not quite the same
acf.hat <- acf(xt); alpha.hat <- acf.hat$acf[2]
alpha.hat


# Drop the first case since lack the predecessor
mu.hat <- c(mean(xt[-1]),  # ols
            mean(generalized.difference(xt,alpha)[-1])/(1-alpha), # utopian gls estimator
            mean(generalized.difference(xt,alpha.hat)[-1])/(1-alpha.hat)
         	)

# ---
# Simulate the estimators
# Names attached to the columns make it easier to keep track of content

n.reps <- 500;

n     <-  400;
mu    <-  -5;
alpha <- 0.862;
results <- matrix(0, nrow=n.reps, ncol=4)
colnames(results) <- c("alpha.hat", "ols", "gls","est.gls")
for (i in 1:n.reps) {
	xt <- mu + sim.ar(rnorm(n),alpha);
	acf.hat <- acf(xt, lag.max=2, plot=FALSE);
	alpha.hat <- acf.hat$acf[2];
	results[i,1] <- alpha.hat;
	results[i,2] <- mean(xt[-1]);
   	results[i,3] <- mean(generalized.difference(xt,alpha)[-1])/(1-alpha);
   	results[i,4] <- mean(generalized.difference(xt,alpha.hat)[-1])/(1-alpha.hat);
}


# Plot a histogram and kernel density of the results for alpha.hat
hist(results[,1], freq=FALSE, main="Estimates of Correlation"); 
abline(v=alpha, col="red")  # appears to be a nontrivial amount of bias
lines(density(results[,1]))

results[1,]

# Collect into a data frame (which is a list) or get just one boxplot
boxplot(as.data.frame(results[,2:4]))

apply(results, 2, mean)
apply(results, 2, sd)