Questions

Question 8.1

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

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

All three figures indicate the data is white noise. When using an ACF plot to determine whether the data resembles white noise you should look for two main things. First, we do not want values to pass the critical value threshold. It is not a problem is a small number of values surpase this threshold as long they appear randomly. The critical value threshold represents a 95% confidence interval. Thus if we were to plot 20 random values we would expect 1 value to surpass the critical value threshold. The second aspect we want to examine is whether there is any apparent pattern to the data. None of the plots exhibit this feature.

  1. 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?

ACF measures the correlation between values in the data. Since we know the dataset is randomly generated there should be no correlation between values besides what occurs by chance. By the law of large numbers, as the number of observations increases, the sample will begin to more closely represent the population. Thus, as the number of observations increases to 1000 it is incredibly unlikely that the data will be correlated in any meaningful way. This is reflected in the smaller bounds for the critical values (which are calculated based on the square root of the number of observations).

The autocorrelations are different because the data is randomly generated. Each sample consists of different unique values that were generated from the same underlying formula.

Question 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.

data("ibmclose")
ibmclose %>%
  ggtsdisplay()

The ibmclose data is non-stationary. The plot of the data shows a trend and changing levels over time. In the ACF plot, a non-stationary dataset will slowly drop towards 0 while a stationary dataset will drop to 0 quickly. The size of lag 1 in the PACF plot indicates that most of the correlation seen in the later lags of the ACF plot are due to the high autocorrelation in the first lag. This sizeable spike can be adressed by taking the first difference. This would be another indication that the data is non-stationary.

Question 8.3

usnetelec

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

  1. usnetelec

The plot of the data shows no need for a transformation. The Pacf plot shows a large spike at lag 1 which indicates that differencing once should make the data stationary.

data(usnetelec)
usnetelec %>%
  ggtsdisplay()

The ACF and PACF plots indicate the data is stationary. The single value in above the critical value in the ACF plot is not problematic as the rest of the data appears to be white noise.

usnetelec %>%
  diff() %>%
  ggtsdisplay()

usgdp

  1. usgdp

The trend line in the original data indicates that a boxcox transformation is likelye necessary.

data(usgdp)
usgdp %>%
  ggtsdisplay()

After box cox transformation and differencing the data appears to still be non-stationary as the ACF plot does not quickly tend towards 0. Taking the second difference creates a plot that appears stationary both in the ACF plot and in the plot of the data.

usgdp %>%
  BoxCox(lambda=0.366352) %>%
  diff() %>%
  diff() %>%
  ggtsdisplay()

mcopper

  1. mcopper

The changing levels in the original data with the highly autocorrelated data indicates that it is not stationary. The variance also appears to be increasing.

data(mcopper)
mcopper %>%
  ggtsdisplay()

Using a BoxCox transformation and taking the difference once appears to lead to stationary data. The ACF plot generally trends towards 0 but there are a few values outside the critical value bound. A second differencing does not appear to address the issue and in fact makes it worse.

mcopper %>%
  BoxCox(lambda=0.1919047) %>%
  diff() %>%
  ggtsdisplay()

enplanements

  1. enplanements

The data is in clear need of transformation due to the increasing variance.

data(enplanements)
enplanements %>%
  ggtsdisplay()

The data appears to have yearly seasonality which will need to be addressed. Afterward a standard difference can be taken to make the data stationary. The big dip in the data is due to the gigantic dip in late 2001. This can likely be attributed to the terrorist attack of September 11th.

enplanements %>%
  BoxCox(lambda=-.2269461) %>%
  diff(lag=12) %>%
  diff() %>%
  ggtsdisplay()

visitors

  1. visitors

The data appears to have growing variance and yearly seasonality. This will be addressed with a BoxCox transformation and a seasonal difference.

data(visitors)
visitors %>%
  ggtsdisplay()

After completing the above transformations the data was still not stationary and thus was differenced again. The resulting data is stationary but the spikes in the PACF at 12 and 24 may indicate a seasonal AR(2) term.

visitors %>%
  BoxCox(lambda=0.2775249) %>%
  diff(lag=12) %>%
  diff() %>%
  ggtsdisplay()

Question 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.

I reloaded my old data and it is clear that the variance needs to be addressed with a BoxCox transformation. The data also appears to have yearly seasonality that will need to be addressed as well.

retaildata <- readxl::read_excel('./retail.xlsx', skip=1)
myts <- retaildata %>%
  select(A3349335T) %>%
  ts(frequency=12, start=c(1982, 4))
myts %>%
  ggtsdisplay()

After performing the Box Cox transformation and addressing the seasonal difference, I was still left with non-stationary data. Exploring the nsdiffs and ndiffs functions indicated that I would need to take the difference one more time in order to achieve stationary data.

myts %>%
  BoxCox(0.193853) %>%
  diff(lag=12) %>%
  diff() %>%
  ggtsdisplay()

Question 8.6

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

  1. 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\).
AR.Model <- function(phi, e){
  y <- ts(numeric(100))
  for(i in 2:100)
    y[i] <- phi*y[i-1] + e[i]
  return(y)
}
  1. Produce a time plot for the series. How does the plot change as you change \(\phi_1\)?
e <- rnorm(100)
AR.Model(0.6, e) %>% 
  autoplot()

AR.Model(0.05, e) %>%
  autoplot()

AR.Model(0.9, e) %>%
  autoplot()

As \(\phi\) grows closer to 1 there is a greater amount of smoothing of the data. This is because a larger proportion of the value of \(y[i]\) is coming from \(y[i-1]\).

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

Following the set up provided in part A, I substituted the function for the MA(1) generation. I moved theta to a parameter to aid in future exploration.

MA.Model <- function(theta, e){
  y <- ts(numeric(100))
  for(i in 2:100){
    y[i] <- e[i] + theta*e[i-1]
  }
  return(y)
}
  1. Produce a time plot for the series. How does the plot change as you change \(\theta_1\)?

The increase in theta determines to what effect the previous error term affects the subsiquent term in the function. The larger the value of theta, the larger the proportion of y’s value comes from the previous error. At a value of 1, the current value of y is even balanced between the two previous errors. This has the effect of smoothing out sharp deviations in errors and “pulling” y’s values closer to the previous error.

MA.Model(0.2, e) %>%
  autoplot()

MA.Model(0.6, e) %>%
  autoplot()

MA.Model(1, e) %>%
  autoplot()

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

I can create this data by combining the two previous functions from part d and part c. I will create a new third function for clarity.

ARMA.Model <- function(phi, theta, e){
  y <- ts(numeric(100))
  for(i in 2:100){
    y[i] <- phi*y[i-1] + e[i] + theta*e[i-1]
  }
  return(y)
}
ARMA.Model(0.6, 0.6, e) %>%
  autoplot()

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

Following the pattern of the previous parts of the question, I created a function that will consider AR(2). This required restricting the range of i values for the first for loop as not to go out of bounds. As indicated by the textbook, a value of \(\phi_1 < 0\) results in data that oscilates between positive and negative values. Furthermore, we violate contraint \(\phi_2 - \phi_1 < 1\) which helps explain the odd resulting graph. Changing the values to fit the requirements results in a more accurate plot.

AR2.Model <- function(phi1, phi2, e){
  y <- ts(numeric(100))
  for(i in 3:100)
    y[i] <- phi1*y[i-1] + phi2*y[i-2] + e[i]
  return(y)
}
AR2.Model(-0.8, 0.3, e) %>%
  autoplot()

Graph the latter two series and compare them.

I’m a bit confused by the selection of \(\phi_1\) and \(\phi_2\) in the previous question. As noted in that problem, the selection of those values means that we will not have a workable model produced by the function. For this problem I have select more reasonable parameters.

The two plots are very similar. As noted, the AR2 model is smoother and has fewer extreme valleys and peaks as the ARMA model. This is due to the additional smoothing brought about by incorporating the y[i-2] value.

autoplot(ARMA.Model(0.6, 0.6, e)) +
  autolayer(AR2.Model(0.5, 0.3, e))

Question 8.7

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

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

Following the steps laid out in section 8.7, I first performed a BoxCox transformation and then differenced the data twice in order to make it stationary. The ndiffs function suggested 2 differences would be necessary.

data("wmurders")
wmurders %>%
  autoplot()

Now that the data is stationary, I can estimate the values of p, d and q. The value of d is 1 and the value of p is 2. This can be determined by the spike in lag2 in the PACF plot. Thus my initial arima assumption is (2, 1, 0).

wmurders %>%
  BoxCox(lambda=0) %>%
  diff() %>%
  ggtsdisplay()

From here, I would have to explore nearby values to determine whether there is a model that is a better fit.

wmurders %>%
  Arima(order=c(2, 1, 0))
## Series: . 
## ARIMA(2,1,0) 
## 
## Coefficients:
##           ar1     ar2
##       -0.0572  0.2967
## s.e.   0.1277  0.1275
## 
## sigma^2 estimated as 0.04265:  log likelihood=9.48
## AIC=-12.96   AICc=-12.48   BIC=-6.99
  1. Should you include a constant in the model? Explain.

A model with \(c\neq 0\) and \(d=1\) results in a long-term forecast that follows a straight line whereas when \(c=0\) the long-term forecast will trend towards a non-zero constant. In this case I find no compelling argument for adding a trend to the forecast. As a result I would not include a constant in the model.

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

There are two different methods of writing in backshift notation according to the text book. The first results in:

\[(1-\phi_1B-\phi_2B^2)(1-B)y_t=c\] and the second (which is the one that R uses is):

\[(1-\phi_1B-\phi_2B^2)((1-B)y_t-\mu)=0\]

  1. Fit the model using R and examine the residuals. Is the model satisfactory?
wmurders %>%
  Arima(order=c(2, 1, 0)) %>%
  checkresiduals()

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

The model is satisfactory. The lag plot shows no values outside the critical range and the portmanteau test has a large p-value.

  1. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.
wmurders %>%
  Arima(order=c(2, 1, 0)) %>%
  forecast(3)
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005       2.553355 2.288694 2.818015 2.148591 2.958118
## 2006       2.533802 2.170056 2.897548 1.977501 3.090104
## 2007       2.524231 2.033825 3.014637 1.774220 3.274242

Starting with the equation identified in part b:

\[(1-\phi_1B-\phi_2B^2)(1-B)y_t=0\] Some algebraic regrouding results in:

\[[1-B(1+\phi_1)+B^2(\phi_1-\phi_2)+\phi_2B^3]y_t=0\]

Multiplying out and applying the backshift operator:

\[y-y_{t-1}(1+\phi_1)+y_{t-2}(\phi_1-\phi_2)+y_{t-3}\phi_2=0\]

As per the steps outlined in the book, we replace \(t\) with \(T+1\) and isolate the term we wish to solve for:

\[y_{T+1}=y_T(1+\phi_1)-y_{T-1}(\phi_1-\phi_2)-y_{T-2}\phi_2\]

All the value on the right are known and can be substituted in:

\[y_{T+1}=2.589382(1+-.0572)-2.662227(-.0572-.2967)-2.797697(.2967)\]

Finally:

2.589382*(1+-.0572)-2.662227*(-.0572-.2967)-2.797697*.2967
## [1] 2.553355

\[y_{T+1}=2.553355\]

This process can be repeated for as many forecasts as desired.

  1. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
wmurders %>%
  Arima(order=c(2, 1, 0)) %>%
  forecast(3) %>%
  autoplot()

  1. Does auto.arima() give the same model you have chosen? If not, which model do you think is better?
wmurders %>%
  auto.arima()
## Series: . 
## 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

Auto arima does not give the same select model. The AIC is worse for the auto model than for my selected model. Comparing the below plot to mine, I prefer mine. The steep decrease in future values to rates below anything seen in the data seem unlikely. The general trend of leveling out seen in my plot appears to follow the expontentially decaying rates from the moore recent years.

wmurders %>%
  auto.arima() %>%
  forecast(3) %>%
  autoplot()