Data 624 HW6

a. Explain the differences among these figures. Do they all indicate that the data are white noise?

Ans: As sample size increases from 36 to 360 to 1000, we expect the corresponding ACFs to shrink to approach 0. That’s what we observed, i.e., the regions banding the y = 0 are getting narrower. As we don’t see more than 5% of the T values portend beyond the blue dotted lines, I rule that the data are white noise.

b. Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?

The critical values resulted at different distances from the mean of zero as there are random autocorrelations with positive magnitude or negative magnitude. By the sheer definition of white noises, the autocorrelations from 3 independently chosen sets of random numbers should exhibit an expected level of randomness in autocorrelations.

a. Use the following R code to generate data from an AR(1) model with \(\phi_1\)=0.6 and \(\sigma^2\)=1. The process starts with \(y_1\)=0.

y <- ts(numeric(100))
e <- rnorm(100)
for (i in 2:100) y[i] <- 0.6 * y[i - 1] + e[i]

b. Produce a time plot for the series. How does the plot change as you change \(\phi_1\)?

seed = 3098
generate_ts <- function(n, phi) {
    set.seed(seed)
    y <- ts(numeric(n))
    e <- rnorm(n)
    for (i in 2:n) y[i] <- phi * y[i - 1] + e[i]
    return(y)
}


autoplot(generate_ts(n = 100, phi = 0.6), ylim = c(-3, 4), series = "4. Phi=.60", 
    main = "Let's look at the time series generated with phi between 0 and .60") + 
    # autolayer(generate_ts(n=100, phi=0.75), series='Phi=0.75') +
# autolayer(generate_ts(n=100, phi=0.99), series='Phi=0.99') +
autolayer(generate_ts(n = 100, phi = 0.45), series = "3. Phi=0.45") + autolayer(generate_ts(n = 100, 
    phi = 0.2), series = "2. Phi=0.20") + autolayer(generate_ts(n = 100, phi = 0), 
    series = "1. Phi=0")

autoplot(generate_ts(n = 100, phi = 0.6), ylim = c(-7, 3), series = "4. Phi=.60", 
    main = "Let's look at the time series generated with phi between 0.60 and .99") + 
    autolayer(generate_ts(n = 100, phi = 0.75), series = "5. Phi=0.75") + autolayer(generate_ts(n = 100, 
    phi = 0.99), series = "6. Phi=0.99")

# autolayer(generate_ts(n=100, phi=0.45), series='3. Phi=0.45') +
# autolayer(generate_ts(n=100, phi=0.20), series='2. Phi=0.20') +
# autolayer(generate_ts(n=100, phi=0), series='1. Phi=0')

Ans: As \(\phi_1\) changes and within the range of 0 and .75, one can tell the ts oscillates and but trends very similarly to each of the time series generated with n = 100. With the exception of when \(\phi_1 \cong 1\), that’s when the generated ts becomes a random walk where you can’t predict whether T = t+1 is going up or going down. In other words, the time plot itself became unpredictable and dissimilar to the rest of the time plots for any other values of \(\phi_1\)’s.

c. Write your own code to generate data from an MA(1) model with θ₁=0.6 and σ₂=1.

Ans: MA(1) = \(y_t=c+ϵ_t+θ_{1}ϵ_{t−1}\)

seed <- 318
ma.1 <- function(n, theta, c, sigma) {
    set.seed(seed)
    y = ts(numeric(n))
    e = rnorm(n, sigma)
    for (i in 2:n) y[i] = c + e[i] + theta * e[i - 1]
    return(y)
}
autoplot(ma.1(100, 0.6, 0, 1))

d. Produce a time plot for the series. How does the plot change as you change θ₁?

For -1 < θ₁ < 1,

autoplot(ma.1(n = 100, theta = 0.6, c = 0, sigma = 1), ylim = c(-2, 4), series = "4. Theta=.60", 
    main = "Let's look at the time series generated with theta between -.45 and .60") + 
    autolayer(ma.1(n = 100, theta = 0.15, c = 0, sigma = 1), series = "3. Theta=.15") + 
    autolayer(ma.1(n = 100, theta = 0, c = 0, sigma = 1), series = "2. Theta=0") + 
    autolayer(ma.1(n = 100, theta = -0.45, c = 0, sigma = 1), series = "1. Theta=-.45")

autoplot(ma.1(n = 100, theta = 0.6, c = 0, sigma = 1), ylim = c(-2, 4), series = "4. Theta=.60", 
    main = "Let's look at the time series generated with theta between .60 and .99") + 
    autolayer(ma.1(n = 100, theta = 0.75, c = 0, sigma = 1), series = "3. Theta=.75") + 
    autolayer(ma.1(n = 100, theta = 0.99, c = 0, sigma = 1), series = "2. Theta=0.99")

ma.1(n = 100, theta = 0, c = 0, sigma = 1) %>% ggtsdisplay(main = "theta = 0")

ma.1(n = 100, theta = 0.15, c = 0, sigma = 1) %>% ggtsdisplay(main = "theta = 0.15")

ma.1(n = 100, theta = 0.6, c = 0, sigma = 1) %>% ggtsdisplay(main = "theta = 0.60")

ma.1(n = 100, theta = -0.45, c = 0, sigma = 1) %>% ggtsdisplay(main = "theta = -0.45")

Ans: The effect of theta as displayed in the graphs above is that it tames the autocorrelations (the spikes) within the ACF and PACF plots. The closer the theta is to 0, the smaller the spikes you can see in those 2 plots. Vice versa. On the other hand, I don’t see there are any detectable effects on the time plots themselves basing just on the variations of theta’s.

e. Generate data from an ARMA(1,1) model with ϕ₁=0.6, θ₁=0.6 and σ₂=1.

seed <- 398
ARMA.1.1 <- function(n, phi, theta, c, sigma) {
    set.seed(seed)
    y <- ts(numeric(n))
    e <- rnorm(n, sigma)
    for (i in 2:n) y[i] <- c + phi * y[i - 1] + theta * e[i - 1] + e[i]
    return(y)
}
autoplot(ARMA.1.1(100, 0.6, 0.6, 0, 1), main = "ARMA(1,1)")

f. Generate data from an AR(2) model with ϕ₁=−0.8, ϕ₂=0.3 and σ₂=1. (Note that these parameters will give a non-stationary series.)

set.seed(301)
phi_1 = -0.8
phi_2 = 0.3
sigma = 1
y2 = ts(numeric(100))
e = rnorm(100, sigma)
for (i in 3:100) y2[i] = phi_1 * y2[i - 1] + phi_2 * y2[i - 2] + e[i]

plot_y2 = autoplot(y2) + labs(title = expression(paste("AR(2): ", phi[1], "= -0.8, ", 
    phi[2], "= 0.3")))
acf_y2 = ggAcf(y2)
gridExtra::grid.arrange(plot_y2, acf_y2, ncol = 2)

g. Graph the latter two series and compare them.

ggtsdisplay(ARMA.1.1(100, 0.6, 0.6, 0, 1), main = "ARMA(1,1): phi[1] = 0.6, theta[1] = 0.6")

ggtsdisplay(y2, main = "AR(2): phi[1] = -0.8, phi[2] = 0.3")

Ans:

AR(1,1) is a white noise as its 95% ACF spikes lie within the dotted blue lines.

AR(2) has an oscillating time plot where there is an increase in the variations at time increases. Because of the two phi’s are of opposite signs, the ACF plot does show autocorrelations of alternate signs lag after lag.

Both series have an ACF plot that is decreasing in autocorrelations as lag increases.

a. By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data.

ggtsdisplay(wmurders, main = "Total Murdered Women", ylab = "Unit: 100K")

Ans:

An ARIMA model is defined by,

• p, the order of autoregression (y regressed on previous values)

• d, differencing required for stationarity

• q, the order of the moving average component

First off, do we need any transformations of the time plot? No. We do not. As we don’t see apparent seasonal variations of the time plot, we decide not to pursue any transformations of the series.

To determine the order of differencing, we use ndiffs and nsdiffs().

# Number of differences required for a stationary series
ndiffs(wmurders)

[1] 2

# Number of differences required for a seasonally stationary series
# nsdiffs(wmurders) Non seasonal data.

There is no seasonal differencing needed as this ts is a non-seasonal data.

ur.kpss(diff(wmurders, differences = 2)) %>% summary

####################### 
# KPSS Unit Root Test # 
####################### 

Test is of type: mu with 3 lags. 

Value of test-statistic is: 0.0458 

Critical value for a significance level of: 
                10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

At significance level of .05, we do not reject the null hypothesis so the series is stationary.

To determine the parameter value of p, the autoregressive term, we look at the PACF plot of the differenced series.

ggPacf(diff(wmurders, differences = 2))

It has a significant spike at lag 1, so p = 1.

To determine the parameter value q, the order of moving average, we look at the ACF plot of the differenced series.

ggAcf(diff(wmurders, differences = 2))

The major spike was at lag = 1, so we use q = 1. We can also try q = 2.

Here are the candidates to check:

ARIMA(1,2,1)

ARIMA(1,2,2)

ARIMA(1,2,0)

ARIMA(0,2,1)

ARIMA(0,2,2)

ARIMA(0,2,0)

ARIMA(2,2,1)

ARIMA(2,2,0)

ARIMA(2,2,2)

Arima(wmurders, order = c(1, 2, 1))

Series: wmurders 
ARIMA(1,2,1) 

Coefficients:
          ar1      ma1
      -0.2434  -0.8261
s.e.   0.1553   0.1143

sigma^2 estimated as 0.04632:  log likelihood=6.44
AIC=-6.88   AICc=-6.39   BIC=-0.97

Arima(wmurders, order = c(1, 2, 2))

Series: wmurders 
ARIMA(1,2,2) 

Coefficients:
          ar1      ma1      ma2
      -0.7677  -0.2812  -0.4977
s.e.   0.2351   0.2879   0.2762

sigma^2 estimated as 0.04552:  log likelihood=7.38
AIC=-6.75   AICc=-5.92   BIC=1.13

Arima(wmurders, order = c(1, 2, 0))

Series: wmurders 
ARIMA(1,2,0) 

Coefficients:
          ar1
      -0.6719
s.e.   0.0981

sigma^2 estimated as 0.05471:  log likelihood=2
AIC=0   AICc=0.24   BIC=3.94

Arima(wmurders, order = c(0, 2, 1))

Series: wmurders 
ARIMA(0,2,1) 

Coefficients:
          ma1
      -0.8995
s.e.   0.0669

sigma^2 estimated as 0.04747:  log likelihood=5.24
AIC=-6.48   AICc=-6.24   BIC=-2.54

Arima(wmurders, order = c(0, 2, 2))

Series: wmurders 
ARIMA(0,2,2) 

Coefficients:
          ma1     ma2
      -1.0181  0.1470
s.e.   0.1220  0.1156

sigma^2 estimated as 0.04702:  log likelihood=6.03
AIC=-6.06   AICc=-5.57   BIC=-0.15

Arima(wmurders, order = c(0, 2, 0))

Series: wmurders 
ARIMA(0,2,0) 

sigma^2 estimated as 0.1007:  log likelihood=-14.38
AIC=30.76   AICc=30.84   BIC=32.73

Arima(wmurders, order = c(2, 2, 1))

Series: wmurders 
ARIMA(2,2,1) 

Coefficients:
          ar1     ar2      ma1
      -0.1247  0.2348  -0.9082
s.e.   0.1612  0.1560   0.0866

sigma^2 estimated as 0.04518:  log likelihood=7.56
AIC=-7.12   AICc=-6.28   BIC=0.76

Arima(wmurders, order = c(2, 2, 0))

Series: wmurders 
ARIMA(2,2,0) 

Coefficients:
          ar1      ar2
      -0.8289  -0.2246
s.e.   0.1346   0.1353

sigma^2 estimated as 0.05292:  log likelihood=3.34
AIC=-0.68   AICc=-0.19   BIC=5.23

Arima(wmurders, order = c(2, 2, 2))

Series: wmurders 
ARIMA(2,2,2) 

Coefficients:
          ar1     ar2      ma1      ma2
      -0.3200  0.1947  -0.7027  -0.1776
s.e.   0.5955  0.2161   0.6030   0.5090

sigma^2 estimated as 0.04597:  log likelihood=7.63
AIC=-5.26   AICc=-3.98   BIC=4.59

The final ARIMA model to pick is ARIMA(2,2,0).

b. Should you include a constant in the model? Explain.

Given that d=2, a positive non-zero constant would give forecasts a quadratic trend, which is undesirable in this case.

c. Write this model in terms of the backshift operator.

For ARIMA(2,2,0), p = 2, d = 2, q = 0

according formula 8.2, we have

(1−ϕ₁B-ϕ₂B²)(1−B)²y_t=ϵ_t

d. Fit the model using R and examine the residuals. Is the model satisfactory?

(fit = Arima(wmurders, order = c(2, 2, 0)))

Series: wmurders 
ARIMA(2,2,0) 

Coefficients:
          ar1      ar2
      -0.8289  -0.2246
s.e.   0.1346   0.1353

sigma^2 estimated as 0.05292:  log likelihood=3.34
AIC=-0.68   AICc=-0.19   BIC=5.23

checkresiduals(fit)

    Ljung-Box test

data:  Residuals from ARIMA(2,2,0)
Q* = 21.05, df = 8, p-value = 0.007015

Model df: 2.   Total lags used: 10

Ans: The distributions of residuals look okay normal. I think it’s satisfactory.

e. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.

f3 <- forecast(fit, h = 3)
f3

     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
2005       2.381222 2.086401 2.676043 1.930333 2.832112
2006       2.271154 1.817137 2.725171 1.576795 2.965513
2007       2.110176 1.410744 2.809609 1.040487 3.179866

The formula for ARIMA(2,2,0) is:

\(y_{t}^{''} = c + \phi_1 y_{t-1}^{''} + \phi_2 y_{t-2}^{''} + \epsilon_t\)

Or alternatively, using the B notation and knowing that c = 0, we have

$$ \[\begin{multline*} \begin{split} (1 - \phi_1B - \phi_2B^2)(1-B)^2y_t &= \epsilon_t \\ (1 - \phi_1B - \phi_2B^2)(1-2B-B^2)y_t &= \epsilon_t \\ y_t-2y_{t-1}+y_{t-2}-\phi_1y_{t-1}+2\phi_1y_{t-2}-\phi_1y_{t-3}-\phi_2y_{t-2}+2\phi_2y_{t-3}-\phi_2y_{t-4} &=\epsilon_{t} \\ y_t = (2+\phi_1)y_{t-1} + (-1+\phi_2-2\phi_1)y_{t-2}+(\phi_1-2\phi_2)y_{t-3}+\phi_2y_{t-4}+\epsilon_{t} \end{split} \end{multline*}\] $$

So y_{T+1 denoted as yt_h1 in the R code block below.}

y_{T+2 denoted as yt_h2}

y_{T+3 denoted as yt_h3}

res <- resid(fit)
len <- length(res)

phi1 <- coef(fit)["ar1"]
phi2 <- coef(fit)["ar2"]
yt_h1 <- (2 + phi1) * wmurders[length(wmurders)] + (-1 + phi2 - 2 * phi1) * wmurders[length(wmurders) - 
    1] + (phi1 - 2 * phi2) * wmurders[length(wmurders) - 2] + phi2 * wmurders[length(wmurders) - 
    3]
names(yt_h1) <- "yt_h1"
yt_h1

   yt_h1 
2.381222 

yt_h2 <- (2 + phi1) * yt_h1 + (-1 + phi2 - 2 * phi1) * wmurders[length(wmurders)] + 
    (phi1 - 2 * phi2) * wmurders[length(wmurders) - 1] + phi2 * wmurders[length(wmurders) - 
    2]
names(yt_h2) <- "yt_h2"
yt_h2

   yt_h2 
2.271154 

yt_h3 <- (2 + phi1) * yt_h2 + (-1 + phi2 - 2 * phi1) * yt_h1 + (phi1 - 2 * phi2) * 
    wmurders[length(wmurders)] + phi2 * wmurders[length(wmurders) - 1]
names(yt_h3) <- "yt_h3"
yt_h3

   yt_h3 
2.110176 

The results above matched exactly with the ones in the forecast (f3).

f. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.

autoplot(forecast(fit, h = 3), PI = TRUE)

g. Does auto.arima() give the same model you have chosen? If not, which model do you think is better?

(fit2 = auto.arima(wmurders, stepwise = TRUE, seasonal = FALSE))

Series: wmurders 
ARIMA(1,2,1) 

Coefficients:
          ar1      ma1
      -0.2434  -0.8261
s.e.   0.1553   0.1143

sigma^2 estimated as 0.04632:  log likelihood=6.44
AIC=-6.88   AICc=-6.39   BIC=-0.97

Ans: The auto.arima() doesn’t give me the same model as it picked ARIMA(1,2,1) instead of ARIMA(2,2,0) that I had.

When comparing the metric corrected AIC, namely AICc, my handpicked model was better, with

AICc=-0.19

Let me know if you think otherwise.