library(fpp3)
## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──
## ✔ tibble      3.1.8      ✔ tsibble     1.1.3 
## ✔ dplyr       1.0.10     ✔ tsibbledata 0.4.1 
## ✔ tidyr       1.3.0      ✔ feasts      0.3.0 
## ✔ lubridate   1.8.0      ✔ fable       0.3.2 
## ✔ ggplot2     3.4.1      ✔ fabletools  0.3.2
## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()
library(latex2exp)
library(tidyr)
library(ggpubr)

Exercise 9.1

Figure 9.32 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?

We can see that each figure suggests that they are all white noise, because the peaks of each figure are all inside the blue dotted lines. Which means that the correlations in each figure are not significantly different from zero.

  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?

The critical values are at different distances from the mean of zero becuase the length of time, T values are different for each. The formula for the critical values is ±2/T−−√ , so as T increases, the crtical value deceases. The autocorrelations are different in each figure because the T value is increasing thus decreasing the critical value area from left to right.

Exercise 9.2

A classic example of a non-stationary series are stock prices. Plot the daily closing prices for Amazon stock (contained in gafa_stock), along with the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.

Filter for amazon and I plot the time series, acf and pacf:

amzn <- gafa_stock %>% 
  filter(Symbol == "AMZN")

amzn %>% 
  gg_tsdisplay(Close, plot_type = 'partial') +
  labs(title= "Daily Closing Prices", subtitle= "Amazon")

The Amazon Daily Closing Prices plot shows an increasing trend for the most part, but it does not have cyclical behavior or seasonality. The ACF graph is slowly decreasing but it does not have a seasonal pattern. In the PACF graph a first lag of around 1 can be observed, which shows that it is not stationary. It must be differentiated to help stabilize the mean.

Plot the time series, acf, pacf:

amzn %>% 
  gg_tsdisplay(difference(Close), plot_type = 'partial') +
  labs(title= "Differenced Daily Closing Prices",
       subtitle= "Amazon")

The ACF plot of the differenced Amazon Daily Closing prices is not autocorrelation.

Exercise 9.3

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

  1. Turkish GDP from global_economy.
turkey <- global_economy %>% 
  filter(Country=='Turkey') %>% 
  select(Country, GDP)

lambda <- turkey %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

turkey  %>%
  mutate(GDP = box_cox(GDP, lambda)) %>%
  features(GDP, unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1
global_economy %>%
  filter(Country == "Turkey") %>%
  gg_tsdisplay(difference(box_cox(GDP,lambda)), plot_type='partial') +
  labs(title = TeX(paste0("Differenced Turkish GDP with $\\lambda$ = ",
         round(lambda,2))))

For the Turkish GDP from global_economy, we find an appropriate Box-Cox transformation with a lambda of approximately 0.16. With respect to the differencing, using the features function we find the number of differences in order to obtain stationary data as one difference.

  1. Accommodation takings in the state of Tasmania from aus_accommodation.
tasmania <- aus_accommodation %>% 
  filter(State == 'Tasmania')  %>%
  select(State, Takings)

lambda_tasmania <- tasmania %>%
  features(Takings, features = guerrero) %>%
  pull(lambda_guerrero)

tasmania %>%
  mutate(Takings = box_cox(Takings, lambda_tasmania)) %>%
  features(Takings, unitroot_ndiffs) 
## # A tibble: 1 × 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      1
aus_accommodation %>%
  filter(State == "Tasmania") %>%
  gg_tsdisplay(difference(box_cox(Takings,lambda), 4), plot_type='partial') +
  labs(title = TeX(paste0("Differenced Tasmania Accomodation Takings with $\\lambda$ = ",
         round(lambda,3))))

For the Accommodation takings in the state of Tasmania from aus_accommodation, we find an appropriate Box-Cox transformation with a lambda of approximately 0.157. With respect to the differencing, using the features function we find the number of differences in order to obtain stationary data as one difference.

  1. Monthly sales from souvenirs.
lambda_souvenirs <- souvenirs %>% 
  features(Sales, features = guerrero) %>%
  pull(lambda_guerrero)

souvenirs %>%
  mutate(Sales = box_cox(Sales, lambda_souvenirs)) %>%
  features(Sales, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1
souvenirs %>%
  gg_tsdisplay(difference(box_cox(Sales,lambda), 12), plot_type='partial', lag = 36) +
  labs(title = TeX(paste0("Differenced Monthly Souvenir Sales with $\\lambda$ = ",
         round(lambda,3))))

For the Monthly sales from souvenirs, we find an appropriate Box-Cox transformation with a lambda of approximately 0.157. With respect to the differencing, using the features function we find the number of differences in order to obtain stationary data as one difference.

Excercise 9.5

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

set.seed(1975)

myseries <- aus_retail %>%
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

myseries %>%
    gg_tsdisplay(Turnover, plot_type = 'partial', lag_max = 36) +
  labs(title= "Western Australia Turnover",
       subtitle="Industry: Department Stores", y = NULL)

myseries %>% 
  transmute(
    `Turnover` = Turnover,
    `Log Turnover` = log(Turnover),
    `Annual Change Log Turnover` = difference(log(Turnover), 12),
        `Doubly differenced log Turnover` =
                     difference(difference(log(Turnover), 12), 1)
  )%>%
  pivot_longer(-Month, names_to="Type", values_to="Turnover") %>%
  mutate(
    Type = factor(Type, levels = c(
      "Turnover",
      "Log Turnover",
      "Annual Change Log Turnover",
      "Doubly differenced log Turnover"))
  ) %>%
  ggplot(aes(x = Month, y = Turnover)) +
  geom_line() +
  facet_grid(vars(Type), scales = "free_y") +
  labs(title= "Western Australia Turnover",
       subtitle="Industry: Department Stores", y = NULL)

The data appears to present seasonality. There seems to be an increasing trend and a pattern of seasonality. A logarithmic transformation was applied. We can see that a growing trend continues.

Excersice 9.6

Simulate and plot some data from simple ARIMA models.

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

set.seed(1978)

y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

head(sim)
## # A tsibble: 6 x 2 [1]
##     idx      y
##   <int>  <dbl>
## 1     1  0    
## 2     2  0.394
## 3     3 -0.324
## 4     4  2.27 
## 5     5  2.22 
## 6     6  0.284
  1. Produce a time plot for the series. How does the plot change as you change Ï•1?
sim %>% autoplot(y, colour = 'blue') +
  labs(title=  latex2exp::TeX(paste0("AR(1) model with $\\phi_{1}$ = 0.6, $\\sigma^2 = 1$, $y_1=0$")))

Create the data y generate the plot:

for(i in 2:100)
  y[i] <- 0.1*y[i-1] + e[i]
sim_2 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- -1.0*y[i-1] + e[i]
sim_3 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- 1.2*y[i-1] + e[i]
sim_4 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- .9*y[i-1] + e[i]
sim_5 <- tsibble(idx = seq_len(100), y = y, index = idx)

plt1 <- sim_5 %>% autoplot(y) +
  labs(title=  latex2exp::TeX(paste0("AR(1) model with $\\phi_{1}$ = 0.9, $\\sigma^2 = 1$")))
plt2 <- sim_2 %>% autoplot(y) +
  labs(title=  latex2exp::TeX(paste0("AR(1) model with $\\phi_{1}$ = 0.1, $\\sigma^2 = 1$")))
plt3 <- sim_3 %>% autoplot(y) +
  labs(title=  latex2exp::TeX(paste0("AR(1) model with $\\phi_{1}$ = -1.0, $\\sigma^2 = 1$")))
plt4 <- sim_4 %>% autoplot(y) +
  labs(title=  latex2exp::TeX(paste0("AR(1) model with $\\phi_{1}$ = 1.2, $\\sigma^2 = 1$")))

ggarrange(plt1, plt2, plt3, plt4, ncol = 2, nrow = 2)

The wavelength of the time series changes as Ï•1 changes. As Ï•1 increases, the magnitude and wavelength increases and as Ï•1 decreases, so does the magnitude and wavelength for values between 0 and 1. For values less than 0, the wavelength is shorter noting higher frequency with varying amplitudes. Changing the Ï•1 results in different time series patterns. When Ï•1=0, it resemebles white noise, and when Ï•1=1, it resembles a random walk. When it becomes negative, it tends to oscillate around the mean. As Ï•1 decreases, the variation, producing more spikes.

  1. Write your own code to generate data from an MA(1) model with θ1=0.6 and σ2=1.
for(i in 2:100)
  y[i] <- 0.6*e[i-1] + e[i] 

sim_ma <- tsibble(idx = seq_len(100), y = y, index = idx)
sim_ma %>% autoplot(y, colour = 'blue')+
  labs(title=  latex2exp::TeX(paste0("MA(1) model with $\\theta_{1} = 0.6$, $\\sigma^2 = 1$, $y_1=0$")))

for(i in 2:100)
  y[i] <- 0*e[i-1] + e[i] 
sim_ma2 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- .9*e[i-1] + e[i] 
sim_ma3 <- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- 1.5*e[i-1] + e[i] 
sim_ma4<- tsibble(idx = seq_len(100), y = y, index = idx)

for(i in 2:100)
  y[i] <- -1.5*e[i-1] + e[i] 
sim_ma5<- tsibble(idx = seq_len(100), y = y, index = idx)

plt5 <- sim_ma2 %>% autoplot(y)+
  labs(title=  latex2exp::TeX(paste0("MA(1) model with $\\theta_{1} = 0$, $\\sigma^2 = 1$")))

plt6 <- sim_ma3 %>% autoplot(y)+
  labs(title=  latex2exp::TeX(paste0("MA(1) model with $\\theta_{1} = 0.9$, $\\sigma^2 = 1$")))

plt7 <- sim_ma4 %>% autoplot(y)+
  labs(title=  latex2exp::TeX(paste0("MA(1) model with $\\theta_{1} = 1.5$, $\\sigma^2 = 1$")))

plt8 <- sim_ma5 %>% autoplot(y)+
  labs(title=  latex2exp::TeX(paste0("MA(1) model with $\\theta_{1} = -1.5$, $\\sigma^2 = 1$")))

ggarrange(plt5, plt6, plt7, plt8, ncol = 2, nrow = 2)

As θ changes, the plots do not display dramatic changes in shape. All plots show steady variance and suggest the series are stationary. The minimum and maximum amplitude changes with little to no shifts in wavelengths.

  1. Generate data from an ARMA(1,1) model with ϕ1=0.6, θ1=0.6, and σ2=1.
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + 0.6*e[i-1] + e[i]

sim_arma <- tsibble(idx = seq_len(100), y = y, index = idx)
  1. 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.)
for(i in 3:100)
  y[i] <- -0.8*y[i-1] + 0.3*y[i-2] + e[i]

sim_ar2 <- tsibble(idx = seq_len(100), y = y, index = idx)
  1. Graph the latter two series and compare them.
plt9 <- sim_arma %>% autoplot(y)+
  labs(title=  latex2exp::TeX(paste0("ARMA(1,1) model with $\\phi_{1} = 0.6$, $\\theta_1=0.6$, and $\\sigma^2=1$")))

plt10 <- sim_ma5 %>% autoplot(y)+
  labs(title=  latex2exp::TeX(paste0("AR(2) model with $\\phi_1=-0.8$, $\\phi_2=0.3$, and $\\sigma^2=1$")))

ggarrange(plt9, plt10, ncol = 1, nrow = 2)

The ARMA(1,1) model appears to be stationary since it appears to be random, the ACF graph is rapidly decreasing, and the PACF is truncated after the first lag. The wavelength for ARMA(1,1) is slightly wider, the amplitude has minimum and maximum values around ±3, and has slightly steady variance throughout the index range.

The wavelength for AR(2) is narrower, has slightly larger amplitudes with minimum and maximum values around ±5 and has a steady variance. The AR(2) model is not stationary, it oscillates around the mean and increases exponentially in variance as the index increases.

Excersice 9.7

Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.

aus_airpassengers
## # A tsibble: 47 x 2 [1Y]
##     Year Passengers
##    <dbl>      <dbl>
##  1  1970       7.32
##  2  1971       7.33
##  3  1972       7.80
##  4  1973       9.38
##  5  1974      10.7 
##  6  1975      11.1 
##  7  1976      10.9 
##  8  1977      11.3 
##  9  1978      12.1 
## 10  1979      13.0 
## # … with 37 more rows
  1. Use ARIMA() to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
fit <- aus_airpassengers %>%
  model(ARIMA(Passengers))

report(fit)
## Series: Passengers 
## Model: ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8963
## s.e.   0.0594
## 
## sigma^2 estimated as 4.308:  log likelihood=-97.02
## AIC=198.04   AICc=198.32   BIC=201.65

Forecast for 10 periods

fit %>% forecast(h=10) %>%
  autoplot(aus_airpassengers) +
  labs(y = "Millions of Passengers", 
       title = "Australian Air Passengers",
       subtitle = "10 Years")

fit %>% gg_tsresiduals() + 
  labs(title = "Australian Air Passengers",
       subtitle = "10 Year Forecast")

Using the ARIMA() function, the model automatically selected for aus_airpassengers data the was an ARIMA(0,2,1). ARIMA (0,2,1) was selected and the residuals do resemble white noise..

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

yt=−0.8963∗ϵt−1+ϵt

  1. Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.
fit2 <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ pdq(0,1,0)))

fit2 %>% forecast(h=10) %>%
  autoplot(aus_airpassengers) +
  labs(y = "Millions of Passengers", 
       title = "Australian Air Passengers",
       subtitle = "10 Year")

fit2 %>% gg_tsresiduals() +
  labs(title = "Australian Air Passengers",
       subtitle = "10 Year")

  1. Plot forecasts from an ARIMA(2,1,2) model with drift and compare these to parts a and c. Remove the constant and see what happens.

When the constant is removed, a null model is produced.

fit3 <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ pdq(2,1,2)))
## Warning: It looks like you're trying to fully specify your ARIMA model but have not said if a constant should be included.
## You can include a constant using `ARIMA(y~1)` to the formula or exclude it by adding `ARIMA(y~0)`.
## Warning: 1 error encountered for ARIMA(Passengers ~ pdq(2, 1, 2))
## [1] Could not find an appropriate ARIMA model.
## This is likely because automatic selection does not select models with characteristic roots that may be numerically unstable.
## For more details, refer to https://otexts.com/fpp3/arima-r.html#plotting-the-characteristic-roots
report(fit3)
## Series: Passengers 
## Model: NULL model 
## NULL model
  1. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
fit %>% gg_tsresiduals() + 
  labs(title = "Residuals for ARIMA(0,2,1) with constant")

The slope become steeper.

Excersice 9.8

For the United States GDP series (from global_economy):

  1. if necessary, find a suitable Box-Cox transformation for the data;

The variation doesn’t seem to increase or decrease with the level of the series, so a transformation isn’t necessary.

us_gdp <- global_economy %>% 
  filter(Country=="United States")%>%
  select(Country, GDP)

us_gdp %>% autoplot(GDP, colour = 'blue') +
  labs(title = "United States", subtitle = "GDP")

  1. fit a suitable ARIMA model to the transformed data using ARIMA();
fit <- us_gdp %>%
  model(
    arima = ARIMA(GDP, stepwise = FALSE, approx = FALSE))

report(fit)
## Series: GDP 
## Model: ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.4206  -0.3048
## s.e.   0.1197   0.1078
## 
## sigma^2 estimated as 2.615e+22:  log likelihood=-1524.08
## AIC=3054.15   AICc=3054.61   BIC=3060.23

Using the ARIMA() function, ARIMA(0,2,2) was found to be the best fit.

  1. try some other plausible models by experimenting with the orders chosen;
us_gdp %>%
  features(GDP, unitroot_ndiffs)
## # A tibble: 1 × 2
##   Country       ndiffs
##   <fct>          <int>
## 1 United States      2
us_gdp %>%
  gg_tsdisplay(GDP, plot_type = 'partial') +
  labs(title = "United States GDP")

Using part b. and the output of the features function, the p and q values for the ARIMA() model are altered below.

P=2:

fit_p2 <- us_gdp %>%
  model(ARIMA(GDP ~ pdq(2,2,2)))

report(fit_p2)
## Series: GDP 
## Model: ARIMA(2,2,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       1.3764  -0.4780  -1.9659  1.0000
## s.e.  0.1216   0.1354   0.0723  0.0719
## 
## sigma^2 estimated as 2.283e+22:  log likelihood=-1521.14
## AIC=3052.27   AICc=3053.47   BIC=3062.4

P=1:

fit_p1 <- us_gdp %>%
  model(ARIMA(GDP ~ pdq(1,2,2)))

report(fit_p1)
## Series: GDP 
## Model: ARIMA(1,2,2) 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.2053  -0.5912  -0.1928
## s.e.  0.3008   0.2886   0.2102
## 
## sigma^2 estimated as 2.646e+22:  log likelihood=-1523.86
## AIC=3055.72   AICc=3056.51   BIC=3063.82

Changing the p value from 2 to 1 increases each metric, AIC, AICc, and the BIC.

Q=0:

fit_q0 <- us_gdp %>%
  model(ARIMA(GDP ~ pdq(2,2,0)))

report(fit_q0)
## Series: GDP 
## Model: ARIMA(2,2,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.2088  -0.2059
## s.e.   0.1320   0.1317
## 
## sigma^2 estimated as 2.984e+22:  log likelihood=-1527.51
## AIC=3061.02   AICc=3061.48   BIC=3067.09

Q=1:

fit_q1 <- us_gdp %>%
  model(ARIMA(GDP ~ pdq(2,2,1)))

report(fit_q1)
## Series: GDP 
## Model: ARIMA(2,2,1) 
## 
## Coefficients:
##          ar1      ar2      ma1
##       0.4321  -0.1605  -0.8028
## s.e.  0.1537   0.1405   0.0908
## 
## sigma^2 estimated as 2.619e+22:  log likelihood=-1523.61
## AIC=3055.22   AICc=3056   BIC=3063.32

Changing the q value from 2 to 1 increased each AIC, AICc, and BIC values; decreasing from 1 to 0 further decreased the metrics as well. This leaves the 2,2,2 as the optimal value of the combinations atempted.

  1. choose what you think is the best model and check the residual diagnostics;

Based on the AICc value, the best performing model is that of ARIMA(2,2,2) with AICc=3053.47.

fit_p2 %>% gg_tsresiduals() +
  labs(title = "USA GDP 10 Year Forecast",
       subtitle = "ARIMA(2,2,2)")

The residual plot above shows a left skew distribution and the ACF plot resembles white noise.

  1. produce forecasts of your fitted model. Do the forecasts look reasonable?
fit_p2 %>% forecast(h=10) %>%
  autoplot(us_gdp) +
  labs(title = "USA GDP 10 Year Forecast",
       subtitle = "ARIMA(2,2,2)")

Yes, the forecast looks reasonable considering the previous data and the

  1. compare the results with what you would obtain using ETS() (with no transformation).
fit_ets <- us_gdp %>% 
  model(ETS(GDP))

report(fit_ets)
## Series: GDP 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9990876 
##     beta  = 0.5011949 
## 
##   Initial states:
##          l[0]        b[0]
##  448093333334 64917355687
## 
##   sigma^2:  7e-04
## 
##      AIC     AICc      BIC 
## 3190.787 3191.941 3201.089
fit_ets %>% forecast(h = 10) %>%
  autoplot(us_gdp) +
  labs(title = "USA GDP 10 Year Forecast",
       subtitle = "ETS(M,A,N)")

fit_ets %>% gg_tsresiduals()  +
  labs(title = "USA GDP 10 Years Forecast",
       subtitle = "ETS(M,A,N)")

ARIMA(2,2,2) model: AICc = 3053.47 ETS(M,A,N) model : AICc = 3191.941

The ACF plots for each have a lack of correlation, suggesting that both forecasts are good. ARIMA(2,2,2) is the model with the best performance, Comparing the values of AICc.