############################

# 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)