DATA 624 - Homework 6

8.1 - Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

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

The figures seem to show white noise data because the ACF plots do not show a discernable pattern of correlation between lags. The expectation is that at least 95% of autocorrelations are within the critical boundary to be considered white noise¹.

Series X2 does have 2 which nudge slightly above the critical values, which means that slightly less than 95% of the values meet the criteria, however it isn’t egregious enough to say that the data do not resemble 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 are based on the length of the time series. Specifically, they are defined as \(\pm \frac{2}{\sqrt{T}}\) \(^1\). So, the more data (i.e. the longer the time series), the critical values shrink:

tibble(series_len = seq(10,1000,5), crit_width = 2*(2/sqrt(seq(10,1000,5)))) %>% ggplot(aes(x=series_len, y=crit_width)) + geom_line() +
  labs(title="Autocorrelation Critical Values", x="Series Length",
       y="Critical Value Width") + theme_light()

8.2 - A classic example of a non-stationary series is the daily closing IBM stock price series (data set `ibmclose`). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

autoplot(ibmclose) + labs(title="Daily IBM Closing Stock Price", x="Day",
                     y="Price") +
  scale_y_continuous(labels=dollar) + theme_light()

ggAcf(ibmclose)

ggPacf(ibmclose)

The first plot shows that the data has an upward and then downward trend. More importantly, the ACF plot shows a slow decline in the correlations which is indicative of non-stationary data. This data would need to be differenced to make it more stationary.

The Partial-ACF plot doesn’t really tell us anything new about whether the data is stationary or not, but it does show that previous values (beyond \(t-1\)) of the closing price do not really correlate with current values - which is what we’d expect with stock market data.

We can confirm that the series is non-stationary with a KPSS test:

ur.kpss(ibmclose) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 3.6236 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The test statistic is considerably higher than the critical values, confirming that the data are not stationary.

8.3 - For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

a) `usnetelec`

Looking at the time series, we see that there is no apparent seasonality, so we can proceed with transforming and differencing:

autoplot(usnetelec) + labs(title="usnetelec Data") +
  theme_light()

First we find the value of \(\lambda\) for the Box-Cox transform and then check to see how many differences we need to take to get the data to be stationary.

l <- round(BoxCox.lambda(usnetelec),2)

usnetelec_xform <- BoxCox(usnetelec, lambda=l)

diffs <- ndiffs(usnetelec_xform)

The BoxCox.lambda() function shows we should use a \(\lambda\) value of 0.52 and then use 2 differences.

usnetelec_diff <- diff(diff(usnetelec_xform))
subt <- bquote("Box-Cox Transform (" * lambda == .(l) * ") with " * .(diffs) * " Difference(s)")

usnetelec_diff %>% autoplot() +
  labs(title="usnetelec Data",
       subtitle = subt) +
  theme_light()

The transformation and differencing seemed to do ok at making the data stationary. Looking at the ACF plot:

ggAcf(usnetelec_diff)

we see a quick decrease in correlation to zero and, other than the first lag, none of the lags that follow rise above the 95% limits. This seems to support that the data is stationary after the transform and differencing. To be sure, we will check with a KPSS test:

ur.kpss(usnetelec_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 3 lags. 
## 
## Value of test-statistic is: 0.0722 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The small test statistic here confirms the data is stationary.

b) `usgdp`

Looking at the usgdp data we see no seasonality, so we proceed with thransforming and differencing.

autoplot(usgdp) + labs(title="usgdp Data") +
  theme_light()

As before, we find the value of \(\lambda\) for the Box-Cox transform and then check to see how many differences we need to take to get the data to be stationary.

l <- round(BoxCox.lambda(usgdp),2)

usgdp_xform <- BoxCox(usgdp, lambda=l)

diffs <- ndiffs(usgdp_xform)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * ") with " * .(diffs) * " Difference(s)")

The BoxCox.lambda() function shows we should use a \(\lambda\) value of 0.37 and then use 1 differences.

usgdp_diff <- diff(usgdp_xform)

usgdp_diff %>% autoplot() +
  labs(title="usgdp Data",
       subtitle = subt) +
  theme_light()

The data seems stationary after the transform and differencing. As before, we will confirm with the ACF plot and a KPSS Root Test:

ggAcf(usgdp_diff)

ur.kpss(usgdp_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 0.2072 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The ACF plot shows some significant correlation in the first two lags, but the rest are under (or nearly under) the critical values. The KPSS Test has a low test statistic which confirms that the data has been made stationary.

c) `mcopper`

This dataset also does not show any seasonality, though there is a significant spike in the last few years.

autoplot(mcopper) + labs(title="mcopper Data") +
  theme_light()

Again, we find the value of \(\lambda\) for the Box-Cox transform and then check to see how many differences we need to take to get the data to be stationary.

l <- round(BoxCox.lambda(mcopper),2)

mcopper_xform <- BoxCox(mcopper, lambda=l)

diffs <- ndiffs(mcopper_xform)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * ") with " * .(diffs) * " Difference(s)")

The BoxCox.lambda() function shows we should use a \(\lambda\) value of 0.19 and then use 1 difference.

mcopper_diff <- diff(mcopper_xform)

mcopper_diff %>% autoplot() +
  labs(title="mcopper Data",
       subtitle = subt) +
  theme_light()

The data seems stationary after the transform and differencing, though there may be a little bit of a trend at the end of the series. To confirm whether the data are stationary, we’ll use an ACF plot and a KPSS Root Test:

ggAcf(mcopper_diff)

ur.kpss(mcopper_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 6 lags. 
## 
## Value of test-statistic is: 0.0571 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The ACF plot shows a quick decline to zero but considerable correlation on lag 1. That may relate to the beginning of a trend that the differencing is not able to remove yet. The KPSS test, however, has a small test statistic, signifying that the data are relatively stationary.

d) `enplanements`

Looking at the time series, we can see that there is seasonality present, so we will need to so seasonal differencing as well.

autoplot(enplanements) + labs(title="enplanements Data") +
  theme_light()

Once again, we find the value of \(\lambda\) for the Box-Cox transform and then check to see how many differences we need to take to get the data to be stationary.

l <- round(BoxCox.lambda(enplanements),2)

plane_xform <- BoxCox(enplanements, lambda=l)

plane_noseas <- diff(plane_xform,lag = 12)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * "), Seasonally Differenced")

autoplot(plane_noseas) + labs(title="enplanements Data",
                              subtitle = subt) +
  theme_light()

The seasonality appears to have been handled, but the data is not stationary yet.

diffs <- ndiffs(plane_noseas)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * ") with " * .(diffs) * " Difference(s) and Seasonally Differenced")

The BoxCox.lambda() function shows we should use a \(\lambda\) value of -0.23 and then use 1 difference.

plane_diff <- diff(plane_noseas)

plane_diff %>% autoplot() +
  labs(title="enplanements Data",
       subtitle = subt) +
  theme_light()

The data seems stationary after the transform and differencing, except for a really large negative spike around 2002 or so. To confirm whether the data are truly stationary, we’ll use an ACF plot and a KPSS Root Test:

ggAcf(plane_diff)

The ACF plot shows a few significant correlations at lags 1 and 12. I’m not sure that is enough to warant further differencing though. We’ll check the KPSS test:

ur.kpss(plane_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0424 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The KPSS test statistic is considerably lower than the significance levels, so the data should be stationary.

e) `visitors`

This data shows definite seasonality so we’ll begin with a seasonal differencing.

autoplot(visitors) + labs(title="visitors Data") +
  theme_light()

l <- round(BoxCox.lambda(visitors),2)

visitors_xform <- BoxCox(visitors, lambda=l)

visitors_noseas <- diff(visitors_xform,lag = 12)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * "), Seasonally Differenced")

autoplot(visitors_noseas) + labs(title="visitors Data",
                              subtitle = subt) +
  theme_light()

The seasonal differencing did well, but a first difference might still be warranted.

diffs <- ndiffs(visitors_noseas)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * ") with " * .(diffs) * " Difference(s) and Seasonally Differenced")

The BoxCox.lambda() function shows we should use a \(\lambda\) value of 0.28 and then use 1 difference.

visitors_diff <- diff(visitors_noseas)

visitors_diff %>% autoplot() +
  labs(title="visitors Data",
       subtitle = subt) +
  theme_light()

The data looks a bit more like white noise, but there is an increase in variance over time that is concerning. We’ll check the ACF plot and KPSS test:

ggAcf(visitors_diff)

ur.kpss(visitors_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 0.0158 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The ACF plot shows some significant correlations, but the KPSS shows a small enough test statistic to indicate stationary data.

8.5 - For your retail data (from Exercise 3 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

# Load the retail data and select the same column as before
retaildata <- readxl::read_excel("retail.xlsx", skip=1)
retail <- ts(retaildata[,"A3349414R"], frequency = 12, start = c(1982,4))

These data definitely have a seasonality component. So, we’ll start with the seasonal differencing.

autoplot(retail)+
  labs(title="Liquor Retailing Turnover",subtitle="Victoria") +
  xlab("Date") + ylab("Turnover")

l <- round(BoxCox.lambda(retail),2)

retail_xform <- BoxCox(retail, lambda=l)

retail_noseas <- diff(retail_xform,lag = 12)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * "), Seasonally Differenced")

autoplot(retail_noseas) + labs(title="Liquor Retailing Turnover",
                              subtitle = subt) +
  theme_light()

Once we have seasonally differenced the data, we look to see if we need to so first order differencing.

ndiffs(retail_noseas)

## [1] 0

The ndiffs() function indicates that we do not need to difference any further. We’ll confirm this with ACF as before.

ggAcf(retail_noseas)

The ACF plot does not at all support the conclusion that the data is stationary!

So, we will proceed with differencing and see how it works:

retail_diff <- diff(retail_noseas)

subt <- bquote("Box-Cox Transform (" * lambda == .(l) * ") with 1 Difference(s) and Seasonally Differenced")

retail_diff %>% autoplot() +
  labs(title="visitors Data",
       subtitle = subt) +
  theme_light()

Now we check the ACF plot again:

ggAcf(retail_diff)

This ACF plot looks much better than the previous one. Let’s check the KPSS test:

ur.kpss(retail_diff) %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 5 lags. 
## 
## Value of test-statistic is: 0.0235 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The test statistic is considerably smaller than the 1 percent critical value, so that supports the conclusion that the data is now stationary.

8.6 - Use R to simulate and plot some data from simple ARIMA models.

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

set.seed(1234)
  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\)?

We’ll compare \(\phi=0.6\) to \(\phi=0.8\)

y2 <- ts(numeric(100))
y3 <- ts(numeric(100))
y4 <- ts(numeric(100))

for(i in 2:100){
    y2[i] <- 0.8*y2[i-1] + e[i]
    y3[i] <- 0.4*y3[i-1] + e[i]
    y4[i] <- 0.2*y4[i-1] + e[i]
}

autoplot(y, series="0.6 (original)") + autolayer(y2, series = "0.8") +
  scale_color_manual(name=expression(phi),
                     values=c("0.6 (original)" = "black",
                              "0.8" = "brown1"))

Now we compare to a smaller value of \(\phi\):

autoplot(y, series="0.6 (original)") + autolayer(y3, series = "0.4") +
  scale_color_manual(name=expression(phi),
                     values=c("0.6 (original)" = "black",
                              "0.4" = "blue1"))

Now an even smaller value of \(\phi\):

autoplot(y, series="0.6 (original)") + autolayer(y4, series = "0.2") +
  scale_color_manual(name=expression(phi),
                     values=c("0.6 (original)" = "black",
                              "0.2" = "darkorange1"))

The larger value of \(\phi\), the more variability there seems to be in the data.

c) Write your own code to generate data from an MA(1) model with \(\theta_1=0.6\) and \(\sigma^2=1\).

A MA(1) model will have the following formula: \(y_t=c+\epsilon_t+\theta_1\epsilon_{t-1}\).

set.seed(3456)
ma1 <- ts(numeric(100))
e <- rnorm(100,0,1)
c <- 1

for(i in 2:100){
  ma1[i] <- c + e[i] + 0.6*e[i-1]
}

autoplot(ma1) + labs(title="MA(1) Model", subtitle=bquote(theta[1] == 0.6 ~ "   " ~ sigma[epsilon]^2 == 1))

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

First we compare to a really small value of \(\theta_1\):

ma2 <- ts(numeric(100))
ma3 <- ts(numeric(100))
ma4 <- ts(numeric(100))

for(i in 2:100){
  ma2[i] <- c + + e[i] + 0.1*e[i-1]
}

for(i in 2:100){
  ma3[i] <- c + + e[i] + 0.4*e[i-1]
}

for(i in 2:100){
  ma4[i] <- c + + e[i] + 0.9*e[i-1]
}

autoplot(ma1, series="0.6 (original)") + autolayer(ma2, series="0.1") +
  scale_color_manual(name=bquote(theta[1]),
                     values=c("0.6 (original)" = "black",
                              "0.1" = "darkorange1"))

Now we’ll compare to a modest value of \(\theta_1\):

autoplot(ma1, series="0.6 (original)") + autolayer(ma3, series="0.4") +
  scale_color_manual(name=bquote(theta[1]),
                     values=c("0.6 (original)" = "black",
                              "0.4" = "yellow3"))

Finally we’ll compare to a really high \(\theta_1\):

autoplot(ma1, series="0.6 (original)") + autolayer(ma4, series="0.9") +
  scale_color_manual(name=bquote(theta[1]),
                     values=c("0.6 (original)" = "black",
                              "0.9" = "brown1"))

A MA(1) model takes into consideration the previous error (\(\epsilon_{t-1}\)) multiplied by \(\theta_1\). So, as the value of that \(\theta_1\) increases, the more that error value (at \(t-1\)) influences the next value.

If \(\theta_1\) approaches 0, the equation \(y_t=c+\epsilon_t+\theta_1\epsilon_{t-1}\) approaches \(y_t=c+\epsilon_t\) all but ignoring the previous prediction error.

e) Generate data from an ARMA(1,1) model with \(\phi_1=0.6\), \(\theta_1=0.6\) and \(\sigma^2=1\).

set.seed(9876)
arma1 <- ts(numeric(100))
e <- rnorm(100,0,1)
c <- 1

for(i in 2:100){
  y[i] <- 0.6*y[i-1] + e[i]
  arma1[i] <- c + 0.6*arma1[i] + 0.6*e[i-1] + e[i] 
}

autoplot(arma1) + labs(title="ARMA(1,1) Model", subtitle=bquote(phi[1] == 0.6 ~ "  " ~ theta[1] == 0.6 ~ "   " ~ sigma[epsilon]^2 == 1))

f) Generate data from an AR(2) model with \(\phi_1 = -.08\), \(\phi_2=0.3\), and \(\sigma^2=1\). (Note that these parameters will give a non-stationary series).

  set.seed(9876)
  y <- ts(numeric(100))
  e <- rnorm(100)
  c <- 1
  
  for(i in 3:100){
    y[i] <- c + -0.8*y[i-1] + 0.3*y[i-2] + e[i]
  }
autoplot(y) + labs(title="AR(2) Model", subtitle=bquote(phi[1] == -0.8 ~ "  " ~ phi[2] == 0.3 ~ "   " ~ sigma[epsilon]^2 == 1))

g) Graph the latter two series and compare them

cbind(`ARMA(1,1)` = arma1, `AR(2)` = y) %>%
autoplot(facets=TRUE)

Interestingly, the AR(2) model oscillates about the value of c and it’s variance quickly increases. This is because the value of \(\phi_1\) is a negative value, flipping the direction of the new point. Then \(\phi_2\) is added, and that point at \(y_{t-2}\) is already in the opposite direction. This creates the oscillation. The error term \(\epsilon\) is quickly drowned out.

In contrast, the ARMA(1,1) model only takes the previous value and the previous error (six tenths of each) and uses that to compute the new value. There is no oscillation created as a result.

8.7 Consider `wmurders`, the number of women murdered each year (per 100,000 standard population) in the United States.

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

We’ll start with a plot of the raw data:

autoplot(wmurders) +
  labs(title="Total Murdered Women", subtitle="Per 100,000 population",
       y="Victims",x="Year") + theme_light()

The first thing we see is that there is a definite trend. In fact two trends. One is an upward trend from about 1955 through 1973 and then a downward trend from 1993 until the data ends.

There also appears to be no seasonality to speak of (the data is also an annual frequency).

We’ll start by checking if a transformation is needed:

l <- BoxCox.lambda(wmurders)
l

## [1] -0.09529835

xform <- BoxCox(wmurders,lambda=l)

autoplot(xform) +
  labs(title="Total Murdered Women", subtitle=bquote("Per 100,000 population; Box-Cox transformed with " ~ lambda == ~ .(round(l,3))),
       y="Victims",x="Year") + theme_light()

Now we look to see about differencing:

diff1 <- diff(xform)

autoplot(diff1) +
  labs(title="Total Murdered Women", subtitle=bquote("Per 100,000 population; Box-Cox transformed with " ~ lambda == ~ .(round(l,3)) ~ "; 1 Difference"),
       y="Victims",x="Year") + theme_light()

There appears to still be an element of a trend here. So we will try a second differencing:

diff2 <- diff(diff1)

autoplot(diff2) +
  labs(title="Total Murdered Women", subtitle=bquote("Per 100,000 population; Box-Cox transformed with " ~ lambda == ~ .(round(l,3)) ~ "; 2 Differences"),
       y="Victims",x="Year") + theme_light()

Now the data looks a lot more stationary.

Next we will use ACF and PACF plots to try and reason what values our ARMIA model should use:

ggAcf(diff2)

ggPacf(diff2)

The sinusoidal nature of the ACF plot and the significant spike in the PACF at lag 1, leads us to believe that this might follow an ARIMA(1,2,0) model.

We’ll try this model and measure AICc:

arima1 <- Arima(wmurders, order=c(1,2,0))

autoplot(wmurders, series="Original Data") +
  autolayer(arima1$fitted, series="ARIMA(1,2,0)", lty=2) +
  annotate("text", x=1975,y=3,label=paste0("AICc = ",round(arima1$aicc,3))) +
  scale_color_manual(name="Series",
                     values=c("Original Data" = "black", "ARIMA(1,2,0)"="red")) +
  labs(title="Total Murdered Women", subtitle="Per 100,000 population",
       y="Victims",x="Year") + theme_light()

The fit isn’t too bad, but we’ll try another model to see if it fits better by changing \(q\) from 0 to 1:

arima2 <- Arima(wmurders, order=c(1,2,1))

autoplot(wmurders, series="Original Data") +
  autolayer(arima2$fitted, series="ARIMA(1,2,1)", lty=2) +
  annotate("text", x=1975,y=3,label=paste0("AICc = ",round(arima2$aicc,3))) +
  scale_color_manual(name="Series",
                     values=c("Original Data" = "black", "ARIMA(1,2,1)"="red")) +
  labs(title="Total Murdered Women", subtitle="Per 100,000 population",
       y="Victims",x="Year") + theme_light()

This new ARIMA(1,2,1) model fits a lot better visually and has a much lower value for AICc than the first model.

Let’s try adjusting the \(p\) parameter instead:

arima3 <- Arima(wmurders, order=c(2,2,0))

autoplot(wmurders, series="Original Data") +
  autolayer(arima3$fitted, series="ARIMA(2,2,0)", lty=2) +
  annotate("text", x=1975,y=3,label=paste0("AICc = ",round(arima3$aicc,3))) +
  scale_color_manual(name="Series",
                     values=c("Original Data" = "black", "ARIMA(2,2,0)"="red")) +
  labs(title="Total Murdered Women", subtitle="Per 100,000 population",
       y="Victims",x="Year") + theme_light()

This ARIMA(2,2,0) fits much better than our original ARIMA(1,2,0) but it is not as good as the ARIMA(1,2,1) model.

We’ll try one more just to be sure. This time we will adjust both \(p\) and \(q\) by 1:

arima4 <- Arima(wmurders, order=c(2,2,1))

autoplot(wmurders, series="Original Data") +
  autolayer(arima4$fitted, series="ARIMA(2,2,1)", lty=2) +
  annotate("text", x=1975,y=3,label=paste0("AICc = ",round(arima4$aicc,3))) +
  scale_color_manual(name="Series",
                     values=c("Original Data" = "black", "ARIMA(2,2,1)"="red")) +
  labs(title="Total Murdered Women", subtitle="Per 100,000 population",
       y="Victims",x="Year") + theme_light()

Interestingly, this latest attempt did get better but it is ever-so-slightly worse than the ARIMA(1,2,1) model; so we will stick with that model.

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

No, I would not include a constant. The original data does not seem to have an overall trend (drift). Plus, as the text\(^1\) states, the auto.arima() function does not do so when \(d=2\), which it is here.

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

\[ (1-\phi_1B) (1-B)^2y_t = (1+\theta_1B)\epsilon_t \]

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

checkresiduals(arima2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,2,1)
## Q* = 12.419, df = 8, p-value = 0.1335
## 
## Model df: 2.   Total lags used: 10

The redisuals don’t have any discernable pattern and look like white noise. The (relatively) high p-value of the Ljung-Box test confirms this.

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

pred <- forecast(arima2,h=3)
pred

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005       2.470660 2.194836 2.746484 2.048824 2.892496
## 2006       2.363106 1.986351 2.739862 1.786908 2.939304
## 2007       2.252833 1.765391 2.740276 1.507354 2.998313

Manually, we can calculate the predictions like this:

\[ (1-\phi_1B) (1-B)^2y_t = (1+\theta_1B)\epsilon_t \] or…

\[ (1-\phi_1B) (1-2B+B^2)y_t = (1+\theta_1B)\epsilon_t \] First, we distribute:

\[ [(1-\phi_1B)-(2B)(1-\phi_1B)+(B^2)(1-\phi_1B)]y_t = \epsilon_t (1+\theta_1B) \]

\[ [1-\phi_1B-2B+2B\phi_1B+B^2-B^2\phi_1B]y_t = \epsilon_t(1+\theta_1B) \]

\[ [1-\phi_1B-2B+2\phi_1B^2+B^2-\phi_1B^3]y_t = \epsilon_t(1+\theta_1B) \] Then apply the backshift:

\[ y_t - \phi_1y_{t-1}-2y_{t-1}+2\phi_1y_{t-2}+y_{t-2}-\phi_1y_{t-3} = \epsilon_t +\theta_1\epsilon_{t-1} \]

Finally, we isolate \(y_t\) on one side of the equation:

\[ y_t = \phi_1y_{t-1}+2y_{t-1}-2\phi_1y_{t-2}-y_{t-2}+\phi_1y_{t-3} + \epsilon_t +\theta_1\epsilon_{t-1} \]

Now we can use this formula to predict \(h=1\):

\[ y_{t+1} = \phi_1y_{t}+2y_{t}-2\phi_1y_{t-1}-y_{t-1}+\phi_1y_{t-2} + \epsilon_{t+1} +\theta_1\epsilon_{t} \]

phi <- arima2$coef["ar1"]
theta <- arima2$coef["ma1"]

# Predict t+1
p1 <- phi*wmurders[55]+2*wmurders[55]-2*phi*wmurders[54]-wmurders[54]+phi*wmurders[53]+0+theta*arima2$residuals[55]

# Predict t+2
p2 <- phi*p1+2*p1-2*phi*wmurders[55]-wmurders[55]+phi*wmurders[54]+0

# Predict t+3
p3 <- phi*p2+2*p2-2*phi*p1-p1+phi*wmurders[55]+0

pred.manual <-  rbind(p1,p2,p3)
rownames(pred.manual) <- c("h+1","h+2","h+3")
colnames(pred.manual) <- "Point Forecast"

pred.manual

##     Point Forecast
## h+1       2.470660
## h+2       2.363106
## h+3       2.252833

The manual calculations match the automatic ones.

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

autoplot(wmurders, series="Actual") +
  autolayer(pred, series=paste("Forecast","ARIMA(1,2,1)",sep="\n")) +
  labs(title="Total Murdered Women", subtitle="Per 100,000 population",
       y="Victims",x="Year") + theme_light() +
  scale_color_manual(name="",
                     values=c("Actual"="black","Forecast\nARIMA(1,2,1)"="blue"))

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

auto.arima(wmurders,lambda="auto")

## Series: wmurders 
## ARIMA(1,2,1) 
## Box Cox transformation: lambda= -0.09529835 
## 
## Coefficients:
##           ar1     ma1
##       -0.3006  -0.786
## s.e.   0.1529   0.119
## 
## sigma^2 estimated as 0.002851:  log likelihood=80.37
## AIC=-154.74   AICc=-154.25   BIC=-148.83

Yes, the auto.arima function gave the same model, ARIMA(1,2,1) that was chosen manually.

Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. OTexts.com/fpp2.↩︎

DATA 624 - Homework 6

Adam Douglas

10/18/2020

8.1 - Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

a) Explain the differences among these figures. Do they all indicate 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?

8.2 - A classic example of a non-stationary series is the daily closing IBM stock price series (data set `ibmclose`). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

8.3 - For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

a) `usnetelec`

b) `usgdp`

c) `mcopper`

d) `enplanements`

e) `visitors`

8.5 - For your retail data (from Exercise 3 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

8.6 - Use R to simulate and plot some data from simple ARIMA models.

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

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

c) Write your own code to generate data from an MA(1) model with \(\theta_1=0.6\) and \(\sigma^2=1\).

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

e) Generate data from an ARMA(1,1) model with \(\phi_1=0.6\), \(\theta_1=0.6\) and \(\sigma^2=1\).

f) Generate data from an AR(2) model with \(\phi_1 = -.08\), \(\phi_2=0.3\), and \(\sigma^2=1\). (Note that these parameters will give a non-stationary series).

g) Graph the latter two series and compare them

8.7 Consider `wmurders`, the number of women murdered each year (per 100,000 standard population) in the United States.

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

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

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

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

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

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

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

DATA 624 - Homework 6

Adam Douglas

10/18/2020

8.1 - Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.

a) Explain the differences among these figures. Do they all indicate 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?

8.2 - A classic example of a non-stationary series is the daily closing IBM stock price series (data set ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

8.3 - For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.

a) usnetelec

b) usgdp

c) mcopper

d) enplanements

e) visitors

8.5 - For your retail data (from Exercise 3 in Section 2.10), find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

8.6 - Use R to simulate and plot some data from simple ARIMA models.

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

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

c) Write your own code to generate data from an MA(1) model with \(\theta_1=0.6\) and \(\sigma^2=1\).

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

e) Generate data from an ARMA(1,1) model with \(\phi_1=0.6\), \(\theta_1=0.6\) and \(\sigma^2=1\).

f) Generate data from an AR(2) model with \(\phi_1 = -.08\), \(\phi_2=0.3\), and \(\sigma^2=1\). (Note that these parameters will give a non-stationary series).

g) Graph the latter two series and compare them

8.7 Consider wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.

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

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

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

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

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

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

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

8.2 - A classic example of a non-stationary series is the daily closing IBM stock price series (data set `ibmclose`). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

a) `usnetelec`

b) `usgdp`

c) `mcopper`

d) `enplanements`

e) `visitors`

8.7 Consider `wmurders`, the number of women murdered each year (per 100,000 standard population) in the United States.

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