####################################################################
# Computing theoretical and empirical ACF for simulatd MA(2) process
####################################################################
# MA(2) coefficients
th1 = -0.9
th2 = -0.5

# Theoretical ACF
acf1 = (-th1+th1*th2)/(1+th1^2+th2^2)
acf2 = (-th2)/(1+th1^2+th2^2)

# MA(2) data
set.seed(132235)
n   = 500
sig = 0.1
y   = arima.sim(n=n,model=list(ma=c(-th1,-th2)),sd=sig)

# time series plot and acf plot
par(mfrow=c(1,2))
ts.plot(y)
acf(y)

# Theoretical and empirical ACF
tab = cbind(c(acf1,acf2),acf(y)$acf[2:3])
rownames(tab) = c("lag1","lag2")
colnames(tab) = c("theoretical","empirical")
tab


####################################################################
# Simulating ARMA(2,1) process
####################################################################
n = 10000
y = arima.sim(n=n,model=list(ar=c(0.9,0.05),ma=c(0.9)),sd=sig)
par(mfrow=c(2,2))
ts.plot(y)
acf(y)
pacf(y)

####################################################################
# Random walk and random walk with drift
####################################################################
n = 500
sig=0.0637
mu = 0.0103
y1 = rep(0,n)
y2 = rep(0,n)
for (t in 2:n){
  y1[t] = y1[t-1] + rnorm(1,0,sig)
  y2[t] = mu + y2[t-1] + rnorm(1,0,sig)
}
tt = 1:n
par(mfrow=c(2,2))
ts.plot(y1)
ts.plot(y2)
fit = lm(y2~tt)
abline(fit$coef,col=2)
ts.plot(diff(y1))
ts.plot(y2-fit$fit)

####################################################################
# Simulating ARIMA(0,3,0) and testing for unit roots
####################################################################
n = 500
y = arima.sim(list(order = c(0,3,0)),n = n,sd=sig)
par(mfrow=c(2,4))
ts.plot(y)
ts.plot(diff(y))
ts.plot(diff(diff(y)))
ts.plot(diff(diff(diff(y))))
acf(y)
acf(diff(y))
acf(diff(diff(y)))
acf(diff(diff(diff(y))))

adf.test(y)
adf.test(diff(y))
adf.test(diff(diff(y)))
adf.test(diff(diff(diff(y))))