### 19.R

### Random walk or trend

n      <- 100
time <- 1:n - (n+1)/2
ssx  <- time %*% time

n.reps <- 2000
ests   <-matrix(0,nrow=n.reps,ncol=3)
colnames(ests) <- c("a","b","t")
for(rep in 1:n.reps) {
	xt <- cumsum(rnorm(n))
	a <- mean(xt)
	b <- (xt %*% time)/ssx
	e <- xt - (a+b*time)
	rmse <- sqrt((e%*%e)/(n-2))
	se <- rmse/sqrt(ssx)
	ests[rep,]<- c(a,b,b/se)
}

# plot data for last realization
plot(time,xt)
regr <- lm(xt~time); summary(regr); abline(regr,col="red")

# histogram
hist(ests[,"t"], probability=TRUE, main="t-stats of linear fit, n=100", xlab="t-stat")
d <- density(ests[,"t"],bw=5)
lines(d$x, d$y, col="red")
lines(d$x, dt(d$x, df=n-2), col="green")

# shape is right, but scale is wrong
qqnorm(ests[,"t"])
abline(0,sd(ests[,"t"]))

length(which(abs(ests[,"t"]) > 2))/n.reps




### Code for simulating AR(1)

sim.ar1 <- function(phi,errors) {
	n <- length(errors)
	result <- errors
	for(t in 2:n)  # x_0 = 0
		result[t] <- result[t] + phi*result[t-1]
	result
	}

est.ar1 <- function(xt) { # assume mean is 0
	n <- length(xt)
	(xt[-n]%*%xt[-1])/(xt[-n]%*%xt[-n])
	}


# -- simulate
n.reps  <- 1000
n.flush <- 50         # burn-in
n       <- 200

phi.vec<- c(0.90, 0.95, 0.99, 1)

ests <- matrix(0, nrow=n.reps, ncol=length(phi.vec))
colnames(ests) <- format(phi.vec)

for(rep in 1:n.reps) {
	wt <- rnorm(n+n.flush)
	for(j in 1:length(phi.vec)) {
		 xt <- sim.ar1(phi.vec[j],wt)[-(1:n.flush)]
		 ests[rep,j] <- est.ar1(xt) - phi.vec[j]  # deviation from true value
	}
}

boxplot(as.data.frame(ests))

par(mfrow=c(1,4))
	for(j in 1:4) {
		hist(ests[,j], probability = TRUE,
					main=paste("phi=",phi.vec[j]), xlim=c(-0.21,0.06), ylim=c(0,30),xlab="phi^")
		d <- density(ests[,j], bw=0.01);
		lines(d$x,d$y, col="red")
		lines(d$x, dnorm(d$x,0,sqrt((1-phi.vec[j]^2)/n)), col="green")
	}
par(mfrow=c(1,1))


###  Dickey-Fuller tests

source("http://www.stat.pitt.edu/stoffer/tsa2/Rcode/itall.R")

# global temperature (n=142)
temperature <- scan("http://www.stat.pitt.edu/stoffer/tsa2/data/globtemp.dat")

plot(temperature)

# fit manually
n      <- length(temperature)
n.lags <- 3
time   <- (3+2):n; time <- time - mean(time)

xt     <- temperature[(n.lags+2):n]
xtm1   <- temperature[(n.lags+1):(n-1)]
xtm2   <- temperature[(n.lags  ):(n-2)]
xtm3   <- temperature[(n.lags-1):(n-3)]
xtm4   <- temperature[(n.lags-2):(n-4)]
xt.d1  <- xt   - xtm1  # dep var
xt.d2  <- xtm1 - xtm2  # lagged differences
xt.d3  <- xtm2 - xtm3
xt.d4  <- xtm3 - xtm4

regr <- lm(xt.d1 ~ time + xtm1 + xt.d2 + xt.d3 + xt.d4); summary(regr)

acf2(residuals(regr))

#
# urca package
# library(urca)
#
# need to peek inside this function to parse the results  

ur.df

r1 <- ur.df(temperature, type="trend", lags=3)  # Test statistics: -2.7181 3.1843 4.1428 
summary(r1)


r2 <- ur.df(temperature, type="drift", lags=3)
r3 <- ur.df(temperature, type="none",  lags=3)