Exercise 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?

Answer:

ACFs (autocorrelation functions) are plotted for random numbers series of various size (36, 360, and 1000). All ACFs are indicating that the data is white noise since the dashed blue lines indicate whether the correlations are significantly different from zero and for all three plots that is the case.

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?

Answer:

The critical values are getting smaller as the size of our series gets larger. According to chapter 2.9 of the textbook: “For a white noise series, we expect 95% of the spikes in the ACF to lie within ±2/√T where T is the length of the time series.” Basically a larger our T is - the smaller our critical values will be.

Exercise 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)
## Warning in data(ibmclose): data set 'ibmclose' not found
autoplot(ibmclose)

ggAcf(ibmclose)

ggPacf(ibmclose)

The plot of the data has both changing trend and levels which indicate that the series is non-stationary. The ACF plot also confirms this - as ACF of non-stationary data decreases slowly, and the value of r1 is often large and positive. The PACF indicates a strong correlation between the data and the first lag. Differencing can help stabilise the mean of this time series by removing changes in the level of a time series.

Exercise 8.3

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

usnetelec

autoplot(usnetelec)

There is a clear upward trend in the data.

BoxCox.lambda(usnetelec)
## [1] 0.5167714
usnetelec_bc <- BoxCox(usnetelec, lambda = "auto")
autoplot(usnetelec_bc)

A Box-Cox transformation with a lambda = 0.5 reduces variance of the dataset.

ndiffs(usnetelec_bc)
## [1] 2

The number of differences required = 2.

usgdp

autoplot(usgdp)

There is a clear upward trend in the data.

BoxCox.lambda(usgdp)
## [1] 0.366352
usgdp_bc <- BoxCox(usgdp, lambda = "auto")
autoplot(usgdp_bc)

A Box-Cox transformation with a lambda = 0.4 reduces variance of the dataset.

ndiffs(usgdp_bc)
## [1] 1

The number of differences required = 1.

mcopper

autoplot(mcopper)

There is a slight upward trend in the data and a lot of changing variance.

BoxCox.lambda(mcopper)
## [1] 0.1919047
mcopper_bc <- BoxCox(mcopper, lambda = "auto")
autoplot(mcopper_bc)

A Box-Cox transformation with a lambda = 0.2 significantly reduces variance of the dataset.

ndiffs(mcopper_bc)
## [1] 1

The number of differences required = 1.

enplanements

autoplot(enplanements)

There is a clear upward trend in the data and variance seems to increase in later years. There seems to be seaosnality in the data.

BoxCox.lambda(enplanements)
## [1] -0.2269461
enplanements_bc <- BoxCox(enplanements, lambda = "auto")
autoplot(enplanements_bc)

A Box-Cox transformation with a lambda = -0.2 reduces variance of the dataset.

nsdiffs(enplanements_bc)
## [1] 1

Because nsdiffs() returns 1 (indicating one seasonal difference is required), we apply the ndiffs() function to the seasonally differenced data.

enplanements_bc %>% log() %>% diff(lag=12) %>% ndiffs()
## [1] 1

These functions suggest we should do both a seasonal difference and a first difference.

visitors

autoplot(visitors)

There is a clear upward trend in the data and variance seems to increase in later years. There seems to be seaosnality in the data as well.

BoxCox.lambda(visitors)
## [1] 0.2775249
visitors_bc <- BoxCox(visitors, lambda = "auto")
autoplot(visitors_bc)

A Box-Cox transformation with a lambda = 0.3 reduces variance of the dataset.

nsdiffs(visitors_bc)
## [1] 1

Because nsdiffs() returns 1 (indicating one seasonal difference is required), we apply the ndiffs() function to the seasonally differenced data.

visitors_bc %>% log() %>% diff(lag=12) %>% ndiffs()
## [1] 1

These functions suggest we should do both a seasonal difference and a first difference.

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

#setwd("/Users/elinaazrilyan/Documents/Data624/")
retaildata <- readxl::read_excel("retail.xlsx", skip=1)
myts <- ts(retaildata[,"A3349882C"],
  frequency=12, start=c(1982,4))

autoplot(myts)

BoxCox.lambda(myts)
## [1] 0.2324297
myts_bx <- BoxCox(myts, lambda = "auto")
autoplot(myts_bx)

myts_bx %>% log() %>% diff(lag=12) %>% ndiffs()
## [1] 1
myts_bx1 <- myts_bx %>% diff()
autoplot(myts_bx1)

A Box-Cox transformation with a lambda = 0.2 reduces variance of the dataset. The required order of differencing is 1.

Exercise 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 ϕ1=0.6 and σ^2=1.The process starts with y1=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 ϕ1 ?

autoplot(y)

Let try a value of 0.1 for ϕ1.

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

Let try a value of 0.3 for ϕ1.

y2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
  y2[i] <- 0.3*y2[i-1] + e[i]
autoplot(y,series='ϕ1 = 0.6') +
  autolayer(y1, series='ϕ1 = 0.1') +
  autolayer(y2, series='ϕ1 = 0.3') +
  ggtitle("ϕ1 changes illustrated")

It seems that a higher value of ϕ1 results in a greater variance - the fluctuations on the plot are bigger in magnktude, but seem to be less frequent.

c)

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

yt = εt + θ1εt − 1 .

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

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

autoplot(k)

Let try a value of 0.1 for θ1.

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

Let try a value of 0.3 for θ1.

k2 <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
  k2[i] <- 0.3*e[i-1] + e[i]
autoplot(k,series='θ1 = 0.6') +
  autolayer(k1, series='θ1 = 0.1') +
  autolayer(k2, series='θ1 = 0.3') +
  ggtitle("θ1 changes illustrated")

It seems that a higher value of θ1 results in fluctuations on the plot that are bigger in magnitude, but seem to be less frequent.

e)

Generate data from an ARMA(1,1) model with ϕ1=0.6 , θ1=0.6 and σ^2=1 .

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

f)

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

z <- ts(numeric(100))
e <- rnorm(100)
for(i in 3:100)
  z[i] <- (-0.8)*z[i-1] + 0.3*z[i-2] + e[i]

g)

Graph the latter two series and compare them.

autoplot(x,series='ARMA(1,1)') +
  ggtitle("ARMA(1,1)")

autoplot(z, series='AR(2)') +
  ggtitle("AR(2)")

ARMA(1,1) and AR(2) graphs are extremely different. ARMA(1,1) has a lot of peaks but seems to stay near 0 while AR(2) seems to increase significantly as time increases.

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

autoplot(wmurders)

Let’s see what transformation are required to make the data stationary:

Log Transformation:

wmurders1 <- wmurders %>% log() 
autoplot(wmurders1)

Order of differencing:

ndiffs(wmurders1)
## [1] 2
wmurders2 <- diff(diff(wmurders1))
autoplot(wmurders2)

ggAcf(wmurders2)

ggPacf(wmurders2)

Now that the data is stationary, and we can study the graphs above to make an assumtion about the best ARIMA model.

The data may follow an ARIMA(0, d , q ) model if the ACF and PACF plots of the differenced data show the following patterns:

  • the PACF is exponentially decaying or sinusoidal;

  • there is a significant spike at lag q in the ACF, but none beyond lag q .

Based on the above - I am going to try (0, 1, 2)

fit2 <- Arima(wmurders, order=c(0,1,2))
fit2
## Series: wmurders 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1     ma2
##       -0.0660  0.3712
## s.e.   0.1263  0.1640
## 
## sigma^2 estimated as 0.0422:  log likelihood=9.71
## AIC=-13.43   AICc=-12.95   BIC=-7.46

b)

Should you include a constant in the model? Explain.

The inclusion of a constant in a non-stationary ARIMA model is equivalent to inducing a polynomial trend of order d in the forecast function. Let’s take a look what the impact of including a contact is:

fit3 <- Arima(wmurders, order=c(0,1,2), include.constant = TRUE)
fit3
## Series: wmurders 
## ARIMA(0,1,2) with drift 
## 
## Coefficients:
##           ma1     ma2   drift
##       -0.0659  0.3711  0.0007
## s.e.   0.1263  0.1641  0.0355
## 
## sigma^2 estimated as 0.04302:  log likelihood=9.71
## AIC=-11.43   AICc=-10.61   BIC=-3.47

There is no improvement in AIC when a constant is included so we won’t include it.

c)

Write this model in terms of the backshift operator.

  1. (1−B)yt = c+(1+θ1B+θ2B^2)et

d)

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

checkresiduals(fit2)

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

Resuduals are resembling white noise, the only concern is that they appear to be slightly not normal (left-skewed) but overall - they indicate a satisfactory model.

e)

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

f3 <- forecast(fit2, h=3)
f3
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005       2.458450 2.195194 2.721707 2.055834 2.861066
## 2006       2.477101 2.116875 2.837327 1.926183 3.028018
## 2007       2.477101 1.979272 2.974929 1.715738 3.238464
  1. (1−B)yt = c+(1+θ1B+θ2B^2)et yt = c+(1+θ1B+θ1B^2)et / (1-B)

yt = c+(yt-1+01(yt-1)+01(yt-2)*et) / (1 - et-1)

tail(wmurders)
## Time Series:
## Start = 1999 
## End = 2004 
## Frequency = 1 
## [1] 3.034263 2.805041 3.304467 2.797697 2.662227 2.589383
tail(residuals(fit2))
## Time Series:
## Start = 1999 
## End = 2004 
## Frequency = 1 
## [1] -0.01715313 -0.22005936  0.49127148 -0.39265501 -0.34376126  0.05023862

ϕ1= -0.0660 and θ1=0.3712

(0 + (2.589383+0.3712*2.589383)+0.3712*-2.662227*0.05023862)/(1+0.34376126)
## [1] 2.60531

My manual calculation is close but I can’t get it to match the forecast exactly, so I must have an error.

f)

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

autoplot(forecast(fit2, h=3))

g)

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

auto.arima(wmurders)
## 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

auto.arima() found a model ARIMA(1,2,1) - the AIC of this model is higher than ARIMA(0,1,2) so I am going to conclude that my original model is better.

autoplot(forecast(auto.arima(wmurders), h=3))

From the graph above it seems that the forecast is showing a continued downward trend.