############################
# Author: Joshua M. Tebbs  #
# Date: 20 Dec 2009        #
# Update: 25 Jul 2011      #
# Purpose: STAT 520 R code #
# CHAPTER 4                #
############################

# Example 4.1.
# MA(1) simulation
# Important! R's convention is to use positive thetas for MA models (so we have to negate)
# E.g., ma = 0.9 means theta = -0.9.
# Page 83
ma.sim = arima.sim(list(order = c(0,0,1), ma = 0.9), n = 100)
par(mfrow=c(2,2))
plot(ma.sim,ylab=expression(Y[t]),xlab="Time",type="o",main="MA(1) simulation")
acf(ma.sim,main="Sample ACF")
plot(y=ma.sim,x=zlag(ma.sim,1),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p',main="Lag 1 scatterplot")
plot(y=ma.sim,x=zlag(ma.sim,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p',main="Lag 2 scatterplot")

# Example 4.1.
# MA(1) simulation
# Important! R's convention is to use positive thetas for MA models (so we have to negate)
# E.g., ma = -0.9 means theta = 0.9.
# Page 84
ma.sim = arima.sim(list(order = c(0,0,1), ma = -0.9), n = 100)
par(mfrow=c(2,2))
plot(ma.sim,ylab=expression(Y[t]),xlab="Time",type="o",main="MA(1) simulation")
acf(ma.sim,main="Sample ACF")
plot(y=ma.sim,x=zlag(ma.sim,1),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p',main="Lag 1 scatterplot")
plot(y=ma.sim,x=zlag(ma.sim,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p',main="Lag 2 scatterplot")

# Example 4.2.
# MA(2) simulation
# Important! R's convention is to use positive thetas for MA models (so we have to negate)
# Page 87
ma.sim = arima.sim(list(order = c(0,0,2), ma = c(-0.9,0.7)), n = 100)
par(mfrow=c(2,2))
plot(ma.sim,ylab=expression(Y[t]),xlab="Time",type="o",main="MA(2) simulation")
acf(ma.sim,main="Sample ACF")
plot(y=ma.sim,x=zlag(ma.sim,1),ylab=expression(Y[t]),xlab=expression(Y[t-1]),type='p',main="Lag 1 scatterplot")
plot(y=ma.sim,x=zlag(ma.sim,2),ylab=expression(Y[t]),xlab=expression(Y[t-2]),type='p',main="Lag 2 scatterplot")

# Example 4.3.
# True AR(1) autocorrelation functions (out to k=20 lags)
# Uses ARMAacf function
# ARMAacf function includes the k=0 lag
# Use y = y[2:21] to remove k=0 lag from ARMAacf output
# Page 91
par(mfrow=c(2,2))
y = ARMAacf(ar = c(0.9), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = c(-0.9), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = c(0.5), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = c(-0.5), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)

# Example 4.3.
# AR(1) simulations
# Page 92
par(mfrow=c(2,2))
ar.sim.1 = arima.sim(list(order = c(1,0,0), ar = 0.9), n = 100)
ar.sim.2 = arima.sim(list(order = c(1,0,0), ar = -0.9), n = 100)
ar.sim.3 = arima.sim(list(order = c(1,0,0), ar = 0.5), n = 100)
ar.sim.4 = arima.sim(list(order = c(1,0,0), ar = -0.5), n = 100)
plot(ar.sim.1,ylab=expression(Y[t]),xlab="Time",type="o")
plot(ar.sim.2,ylab=expression(Y[t]),xlab="Time",type="o")
plot(ar.sim.3,ylab=expression(Y[t]),xlab="Time",type="o")
plot(ar.sim.4,ylab=expression(Y[t]),xlab="Time",type="o")

# AR(1) sample ACFs
# Page 93
par(mfrow=c(2,2))
acf(ar.sim.1,main="Sample ACF")
acf(ar.sim.2,main="Sample ACF")
acf(ar.sim.3,main="Sample ACF")
acf(ar.sim.4,main="Sample ACF")

# Figure 4.7.
# AR(2) stationarity region
# Page 97
par(xaxs = "i", yaxs = "i")
plot(x = -2, y = -1, xlim = c(-2,2), ylim = c(-1,1),
    type = "n", frame.plot = FALSE, xlab = expression(phi[1]), ylab = expression(phi[2]))
#abline() draws y = mx + b      
abline(a = 1, b = -1) #Draw line of phi2 = 1 – phi1
abline(a = 1, b =  1) #Draw line of phi2 = 1 + phi1
#Plot the phi1 and phi2 values
points(x = 0, y = 0.5, pch = 1, col = "red")
points(x = 0, y = -0.5, pch = 2, col = "darkgreen")
points(x = -0.2, y = -0.5, pch = 2, col = "darkgreen")
points(x = -1, y = -0.5, pch = 2, col = "darkgreen")
points(x = -1.8, y = -0.9, pch = 2, col = "darkgreen")
points(x = 0.5, y = 0.25, pch = 1, col = "red")
points(x = 1.8, y = 0.9, pch = 3, col = "blue")
points(x = -1.2, y = 0.8, pch = 3, col = "blue")
legend(locator(1), legend = c("Real roots",
      "Complex roots", "Outside stationarity region"), pch =
      c(1,2,3), col = c("red", "darkgreen", "blue"), cex =
      0.75, bty = "n")

# Example 4.4.
# True AR(2) autocorrelation functions (out to k=20 lags)
# Uses ARMAacf function
# ARMAacf function includes the k=0 lag
# Use y = y[2:21] to remove k=0 lag from ARMAacf output
# Page 100
par(mfrow=c(2,2))
y = ARMAacf(ar = c(0.5,-0.5), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = c(1.1,-0.3), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = c(-0.5,0.25), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = c(1,-0.5), lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)

# Example 4.4.
# AR(2) simulations
# Page 101
par(mfrow=c(2,2))
ar2.sim.1 = arima.sim(list(order = c(2,0,0), ar = c(0.5,-0.5)), n = 100)
ar2.sim.2 = arima.sim(list(order = c(2,0,0), ar = c(1.1,-0.3)), n = 100)
ar2.sim.3 = arima.sim(list(order = c(2,0,0), ar = c(-0.5,-0.25)), n = 100)
ar2.sim.4 = arima.sim(list(order = c(2,0,0), ar = c(1,-0.5)), n = 100)
plot(ar2.sim.1,ylab=expression(Y[t]),xlab="Time",type="o")
plot(ar2.sim.2,ylab=expression(Y[t]),xlab="Time",type="o")
plot(ar2.sim.3,ylab=expression(Y[t]),xlab="Time",type="o")
plot(ar2.sim.4,ylab=expression(Y[t]),xlab="Time",type="o")

# AR(2) sample ACFs
# Page 102
par(mfrow=c(2,2))
acf(ar2.sim.1,main="Sample ACF")
acf(ar2.sim.2,main="Sample ACF")
acf(ar2.sim.3,main="Sample ACF")
acf(ar2.sim.4,main="Sample ACF")

# Figure 4.11.
# True ARMA(1,1) autocorrelation functions (out to k=20 lags)
# Uses ARMAacf function
# ARMAacf function includes the k=0 lag
# Use y = y[2:21] to remove k=0 lag from ARMAacf output
# Page 112
par(mfrow=c(2,2))
y = ARMAacf(ar = 0.9, ma = 0.25, lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = -0.9, ma = 0.25, lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = 0.5, ma = 0.25, lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)
y = ARMAacf(ar = -0.5, ma = 0.25, lag.max = 20)
y = y[2:21]
plot(y, x = 1:20, type = "h", ylim = c(-1,1), xlab = "k",
      ylab = "Autocorrelation", main = "Population ACF")
abline(h = 0)