library(fpp2)
library(tidyr)
library(dplyr)
library(ggplot2)
library(dse)
library(kableExtra)
library(forecast)
library(smooth)

8.1

a.

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

The differences in the ACF plots can be found in the critical value thresholds and the magnitude of the autocorrelations. As the length of the time series increases in random numbers, the autocorrelations appear to become more narrow and the critical values appear to shrink.

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?

As the observations increase, the large outliers from the mean decreases. A large observation is an outlier with more data. Regarding a white noise series, all autocorrelations converge to zero over long run independent observations in a random process therefore appearing uncorrelated. In the short run there can be fluctuations in autocorrelations for random patterns in a small dataset. This explains why the autocorrelations appear larger with a shorter time series.

8.2

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.

The time series plot displays cyclic patterns and elements of the plot depend on where they are measured. This plot is non-stationary. The ACF plot demonstrates that the time series has autocorrelations for items past 25 days however the autocorrelations decline over the time period. The PACF plot indicates no autocorrelations after a lag of 1 and accounting for lower-lagged autocorrelations.

ggtsdisplay(ibmclose, xlab = "Day")

The ndiffs function can be used to determine the required order of differencing. First-order differencing should be applied to make the time series stationary.

ndiffs(ibmclose) 
## [1] 1

The plot below indicates a stationary plot, after applying first-order differencing.It appears that the time series plot indicates white noise.

ggtsdisplay(diff(ibmclose), xlab = "Day")

8.3

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

a. usnetelec

The below is an annual time series in which no Box-Cox transformation and no seasonal differencing is needed however first-order differencing is needed.

autoplot(usnetelec)

ndiffs(usnetelec)
## [1] 1
 ggtsdisplay(diff(usnetelec))

b. usgdp

The below is a quarterly time series in which the Box-Cox transformationis λ=0.37. No seasonal differencing is needed however first-order differencing is.

autoplot(usgdp)

BoxCox.lambda(usgdp)
## [1] 0.366352
tstransf <- BoxCox(usgdp, lambda = "auto")
nsdiffs(tstransf)
## [1] 0
ndiffs(tstransf)
## [1] 1
ggtsdisplay(diff(tstransf))

c. mcopper

The below is a monthly time series in which the Box-Cox transformation is λ=0.19. No seasonal differencing is needed and first-order differencing is.

autoplot(mcopper)

BoxCox.lambda(mcopper)
## [1] 0.1919047
tstransf <- BoxCox(mcopper, lambda = "auto")
nsdiffs(tstransf)
## [1] 0
ndiffs(tstransf)
## [1] 1
ggtsdisplay(diff(tstransf))

d. enplanements

the below is a monthly time series. The Box-Cox transformation is λ=−0.23. Seasonal differencing occurs with a 12-month lag and there is first-order differencing.

autoplot(enplanements)

BoxCox.lambda(enplanements)
## [1] -0.2269461
tstransf <- BoxCox(enplanements, lambda = "auto")
nsdiffs(tstransf)
## [1] 1
ndiffs(diff(tstransf, lag = 12))
## [1] 1
  ggtsdisplay(diff(lag=12,tstransf))

e. visitors

The below is a monthly time series. The Box-Cox transformation is λ=0.28. Seasonal differencing occurs with a 12-month lag and there is first-order differencing.

autoplot(visitors)

BoxCox.lambda(visitors)
## [1] 0.2775249
tstransf <- BoxCox(visitors, lambda = "auto")
nsdiffs(tstransf)
## [1] 1
ndiffs(diff(tstransf, lag = 12))
## [1] 1
ggtsdisplay(diff(lag = 12, tstransf))

8.5

Find the appropriate order of differencing (after transformation if necessary) to obtain stationary data.

To obtain stationary data, the following steps were taken:

The Box-Cox transformation is λ=0.92. There is seasonal differencing with 12-month lag and there is first-order differencing. Despite the result of the ndiffs function indicating that no further differencing is needed, the plot of the time series suggests first-order differencing.

retaildata <- readxl::read_excel("retail.xlsx", skip=1)
colID <- colnames(retaildata)[8]
myts <- ts(retaildata[ , colID], frequency=12, start=c(1982,4))

autoplot(myts) + labs(x = "Year", y = "")

BoxCox.lambda(myts)
## [1] 0.9165544
tstransf <- BoxCox(myts, lambda = "auto")
nsdiffs(tstransf)
## [1] 1
ndiffs(diff(tstransf, lag = 12))
## [1] 0
  ggtsdisplay(diff(lag=12,tstransf))

8.6

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

set.seed(314)
e <- rnorm(100)

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]

Data for the AR(1) model is generated with the parameter Ï• ranging from -1 to 1

phi <- seq(-1, 1, 0.2)
y <- ts(matrix(0, nrow = 100, ncol = length(phi)))
colnames(y) <- paste0("phi = ", phi)
for (j in 1:length(phi)) {for (i in 2:100) {y[i, j] <- phi[j] * y[i-1, j] + e[i]}}
phi
##  [1] -1.0 -0.8 -0.6 -0.4 -0.2  0.0  0.2  0.4  0.6  0.8  1.0

b.

Produce a time plot for the series. How does the plot change as you change Ï•1?

The plot below displays the time series generated from the AR(1) model for ϕ∈, while the sample errors ϵt are held constant across scenarios. For ϕ=−1, the time series grows in magnitude over time. As ϕ increases from -0.8 to -0.2, the time series gets smaller and the cycle appear to weaken. For ϕ=0, the time series displays white noise. As ϕ increases from 0.2 to 0.8, there are fewer changes of direction and the magnitude grows. For ϕ=1, the time series is a random walk which shows long periods of positive and negatve trends.

autoplot(y, facets=TRUE)

c.

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

theta <- seq(-1, 1, 0.2)
y <- ts(matrix(0, nrow = 100, ncol = length(theta)))
colnames(y) <- paste0("theta=", theta)
for (j in 1:length(theta)) {for (i in 2:100) {y[i, j] <- e[i] + theta[j] * e[i-1]}}
theta
##  [1] -1.0 -0.8 -0.6 -0.4 -0.2  0.0  0.2  0.4  0.6  0.8  1.0

d.

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

The plot below displays the time series generated from the MA(1) model for θ∈{−1,−0.8,−0.6,⋯,0.6,0.8,1}, while holding the sample errors ϵt constant across scenarios. As θ increases from -1 to 0, the curve appears to smoothen out as sudden changes of direction seem to weaken. For θ=0, the time series is white noise and is identical to the AR(1) model above.As θ increases from 0 to 1, the magnitude of the time series increases and trend lines are more pronounced while changes of direction appear to weaken.

autoplot(y, facets=TRUE)

## e.

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

phi <- 0.6
theta <- 0.6
y1 <- ts(numeric(100))
for (i in 2:100) {y1[i] <- phi * y1[i-1] + theta * 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.)

phi1 <- -0.8
phi2 <- 0.3
y2 <- ts(numeric(100))
for (i in 3:100) {y2[i] <- phi1 * y2[i-1] + phi2 * y2[i-2] + e[i]}

g.

Graph the latter two series and compare them.

The plot below displays the time series generated from the ARMA(1,1) (for ϕ1=0.6 and θ1=0.6) and the AR(2) model (for ϕ1=−0.8 and ϕ2=0.3)and uses same sample errors ϵt in both models (ϵt∼N(μ=0,σ2=1)). In the ARMA(1,1) there are periods where cyclical patterns are evident, but trends tend to revert to the mean over time. In the AR(2) time series rapid cycles occur and the amplitude grows over time.

cbind("ARMA(1,1)" = y1, "AR(2)" = y2) %>% 
  autoplot(facets = TRUE)

8.7

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

autoplot(wmurders)

In the plot the Box-Cox transformation is λ=−0.10. There appears to be Second-order differencing The resulting time series appears stationary.

a.

By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data. Based on the ACF plot at lags 1 and 2 and in the PACF plot at lags 1 and 2, the following ARIMA models are considered: ARIMA(2,2,0), ARIMA(0,2,2), ARIMA(1,2,1).

BoxCox.lambda(wmurders)
## [1] -0.09529835
tstransf <- BoxCox(wmurders, lambda = "auto")
ndiffs(tstransf)
## [1] 2
  ggtsdisplay(diff(tstransf))

b.

Should you include a constant in the model? Explain.

Adding a constant to the model will induce a straight line trend for d=1 and a quadratic trend for d=2 and since the original time series does not show a trend and therefore the constant can be ommited.

c.

Write this model in terms of the backshift operator.

Using backshift notation, the general non-seasonal ARIMA model ARIMA(p,d,q) can be written as the following:(1−ϕ1B−ϕ2B2)(1−B)2yt=ϵt, (1−B)2yt=(1+θ1B+θ2B2)ϵt, (1−ϕ1B)(1−B)2yt=(1+θ1B)ϵt

d.

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

The three candidate models and their residuals are examined. Also considered are two models using the auto.arima function.

In Fit 1, ARIMA(2,2,0) the Ljung-Box test indicates they don’t resemble white noise.

(fit1 <- Arima(wmurders, lambda = "auto", order = c(2,2,0), include.constant = FALSE))
## Series: wmurders 
## ARIMA(2,2,0) 
## Box Cox transformation: lambda= -0.09529835 
## 
## Coefficients:
##           ar1      ar2
##       -0.8786  -0.2906
## s.e.   0.1330   0.1341
## 
## sigma^2 estimated as 0.003242:  log likelihood=77.31
## AIC=-148.61   AICc=-148.12   BIC=-142.7
checkresiduals(fit1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,2,0)
## Q* = 17.975, df = 8, p-value = 0.02141
## 
## Model df: 2.   Total lags used: 10
(fit2 <- Arima(wmurders, lambda = "auto", order = c(0,2,2), include.constant = FALSE))
## Series: wmurders 
## ARIMA(0,2,2) 
## Box Cox transformation: lambda= -0.09529835 
## 
## Coefficients:
##           ma1     ma2
##       -1.0375  0.1950
## s.e.   0.1246  0.1196
## 
## sigma^2 estimated as 0.002908:  log likelihood=79.85
## AIC=-153.69   AICc=-153.2   BIC=-147.78
checkresiduals(fit2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,2,2)
## Q* = 10.119, df = 8, p-value = 0.2568
## 
## Model df: 2.   Total lags used: 10
(fit3 <- Arima(wmurders, lambda = "auto", order = c(1,2,1), include.constant = FALSE))
## 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
checkresiduals(fit3)

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

Auto fit method

(fit4 <- auto.arima(wmurders, lambda = "auto", seasonal = FALSE, allowdrift = FALSE, allowmean = FALSE))
## 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

Auto method with arguments, stepwise and approximation set to false.

(fit5 <- auto.arima(wmurders, lambda = "auto", seasonal = FALSE, allowdrift = FALSE, allowmean = FALSE, stepwise = FALSE, approximation = FALSE))
## 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

RSE, AIC, AICc, and BIC comparison for our candidate models. Based on the below, fit3 is the best model. The residual plots given above for fit3 appear satisfactory.

fitlist <- list(fit1, fit2, fit3)
metrics <- NULL 
for (j in 1:length(fitlist)) {temp <- fitlist[[j]]
  metrics <- rbind(metrics, c(sqrt(temp$sigma2),temp$aic, temp$aicc,temp$bic))
}
row.names(metrics) <- paste0("fit", 1:length(fitlist)) 
metrics %>% kable(digits = c(5, 2, 2, 2), col.names=c("Sigma^2", "AIC", "AICc", "BIC"))
Sigma^2 AIC AICc BIC
fit1 0.05694 -148.61 -148.12 -142.70
fit2 0.05392 -153.69 -153.20 -147.78
fit3 0.05339 -154.74 -154.25 -148.83

e.

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

Checking the forecasts by hand, the model without constant could be written as follows

(1−B)2yt=(1+θ1B)ϵt(1−ϕ1B)(1−2B+B2)yt= (1+θ1B)ϵt(1−(ϕ1+2)B+(2ϕ1+1)B2−ϕ1B3)yt=(1+θ1B)ϵt yt=(ϕ1+2)yt−1−(2ϕ1+1)yt−2+ϕ1yt−3+ϵt+θ1ϵt−1

For future periods, the formula could be written as yT+h=(ϕ1+2)yT+h−1−(2ϕ1+1)yT+h−2+ϕ1yT+h−3+ϵT+h+θ1ϵT+h−1

Future observations are replaced with their forecasts and future errors are set to zero and past errors with model residuals are replaced. Three point forecasts are given by:

y^T+1 = (ϕ1+2)yT−(2ϕ1+1)yT−1+ϕ1yT−2+θ1eT y^T+2 = (ϕ1+2)yT+1−(2ϕ1+1)yT+ϕ1yT−1) yT+3 = (ϕ1+2)yT+2−(2ϕ1+1)yT+1+ϕ1yT

The same point forecasts using our manual calculation matches. The observations and forecasts inverse-transform are transformed.

fit3 %>% forecast(h = 3) 
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2005       2.487322 2.309061 2.680766 2.220398 2.789794
## 2006       2.398638 2.169841 2.654125 2.058488 2.801324
## 2007       2.310805 2.027049 2.638650 1.892470 2.832570

f.

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

fit3 %>% forecast(h = 3) %>% autoplot() + labs(x = "Year", y = "") + geom_vline(xintercept = 2005, lty = 3, lwd = 1)

## g.

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

The auto.arima function does gives the same model as ARIMA(1,2,1) [fit3]. Whether using the default search method or using the long search method (Stepwise & approximation set to FALSE), the same model is arrived at.