9.1

9.1.1

All appear to indicate that the data is white noise as none of the autocorrelations rise to the level of significance, although I will admit the plot of 1000 random numbers is difficult to determine with the naked eye due to the size of the spikes and low critical value.

9.1.2

The critical values decrease with the number of observations because the critical value is defined as \(±1.96/\sqrt{T}\) where T is the length of the time series. This means that as the time series grows longer, the critical values get smaller. The autocorrelations differ due to simple random variation in the series. As long as there is not a statistically significant number of autocorrelations exceeding the critical values or autocorrelations significantly exceeding the critical values the series is indistinguishable from white noise.

9.2

Beginning with the plot of the daily closing price we can see that there is an obvious upward trend across the series with an inversion near the end of the trading period. The presence of a trend in a series means it is non-stationary. Moving on to the ACF plot we can see that there is a significant autocorrelation that lasts for over 30 lags. A stationary series would have an ACF plot that drops to zero relatively quickly. Additionally the value of the first lag is large and positive which is another indicator on non-stationary data. Finally a white noise series would have no significance in the PACF plot but there is a clear highly significant PACF at the first lag. All these indicators lead us to conclude the data is non-stationary and must be differenced to achieve stationarity.

amzn_stock <- gafa_stock |> 
  filter(Symbol == "AMZN") |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)

amzn_stock |>
  gg_tsdisplay(Close, plot_type = "partial")

9.3

9.3.1

Looking at the initial Turkish GDP data we can see the data exhibits an exponential trend, indicating it would benefit from a Box-Cox transformation.

turk_gdp <- global_economy |>
  filter(Country == "Turkey") |>
  select(GDP)

turk_gdp |>
  autoplot(GDP) +
  labs(title = "Untransformed Turkish GDP Data")

After applying the Box-Cox transformation we have a much smoother trend and decreased variance.

lambda <- turk_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

transformed_gdp <- turk_gdp |>
  mutate(trans_gdp = box_cox(GDP, lambda)) |>
  select(!GDP)

transformed_gdp |>
  autoplot(trans_gdp) +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed Turkish GDP with $\\lambda$ = ",
         round(lambda,2))))

Performing a KPSS test lets us know how many differences we need to take to achieve stationary data. In this case the transformed data needs only one differencing. Since we can see from the previous plots there is no seasonality in the data we will not bother checking if seasonal differencing is needed.

transformed_gdp |>
  features(trans_gdp, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1

Upon plotting the transformed, differenced data we can see the mean is constant and the variance and ACF do not change over time indicating the series is now stationary

differenced_gdp <- transformed_gdp |>
  mutate(dif_gdp = difference(trans_gdp))

differenced_gdp |>
  gg_tsdisplay(dif_gdp, plot_type = "partial")

9.3.2

The Tasmanian accomodation takings show a positive trend with strong seasonality that increases with time which indicates we need to apply a transformation before differencing.

tas_take <- aus_accommodation |>
  filter(State == "Tasmania") |>
  select(Takings)

tas_take |>
  autoplot(Takings) +
  labs(title = "Untransformed Tasmanian Accomodation Takings Data")

After applying a Box-Cox transformation of 0 which is equivalent to a log transformation we can see the seasonality is no longer increasing with time.

lambda <- tas_take |>
  features(Takings, features = guerrero) |>
  pull(lambda_guerrero)

trans_tas_take <- tas_take |>
  mutate(trans_take = box_cox(Takings, lambda)) |>
  select(!Takings)

trans_tas_take |>
  autoplot(trans_take) +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed Tasmanian Takings with $\\lambda$ = ",
         round(lambda,2))))

Given the obvious seasonality we will start with determining what seasonal differencing is needed and apply the single seasonal difference indicated. We can see from the plots that the data is relatively stationary, and performing a non-seasonal KPSS test we can see that additional differncing would not improve the stationarity.

trans_tas_take |>
  features(trans_take, unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
s_diff_tas_take <- trans_tas_take |>
  mutate(s_dif_take = difference(trans_take, 4))

s_diff_tas_take |>
  gg_tsdisplay(s_dif_take, plot_type = "partial")

s_diff_tas_take |>
  features(s_dif_take, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0

9.3.3

The monthly sales data shows some very strong and increasing seasonality, with a relatively small positive trend. The rest of this answer will be slightly abridged due to it being largely the same process as the previous two questions.

souvenirs |>
  autoplot(Sales) +
  labs(title = "Untransformed Sales Data")

Once again we will be using a log transformation to stabilize the variance, and move on to differencing. Normally we perform seasonal differencing first, then ordinary differencing but we can see the PKSS indicates no ordinary differencing is required although our data doesn’t look extremely stationary.

lambda <- souvenirs |>
  features(Sales, features = guerrero) |>
  pull(lambda_guerrero)

trans_sales <- souvenirs |>
  mutate(trans_sales = box_cox(Sales, lambda)) |>
  select(!Sales)

trans_sales |>
  features(trans_sales, unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
s_diff_trans_sales <- trans_sales |>
  mutate(s_dif_sales = difference(trans_sales, 12))

s_diff_trans_sales |>
  features(s_dif_sales, unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      0
s_diff_trans_sales |>
  gg_tsdisplay(s_dif_sales, plot_type = "partial")

If we instead do a first order difference, then a seasonal difference we get a better looking distribution. There are no longer large postive ACF and PACF values at the first lag, and the mean is more constant without the apparent trends in the previous transformation. Whether the previous set of differencing would be preferable due to a lower degree of complexity depends on the accuracy difference in the models produced.

dif_trans_sales <- trans_sales |>
  mutate(dif_sales = difference(trans_sales, 1))

dif_trans_sales |>
  features(dif_sales, unitroot_nsdiffs)
## # A tibble: 1 × 1
##   nsdiffs
##     <int>
## 1       1
dif_trans_sales |>
  mutate(s_dif_sales = difference(dif_sales, 12)) |>
  gg_tsdisplay(s_dif_sales, plot_type = "partial")

9.5

The set of data we will be looking at is the footwear and personal accessory turnover data from the aus_retail dataset. Using the guerrero feature we are able to determine the optimal Box-Cox transformation of -0.22 which ends up producing a series with a bit of a logarithmic trend line but standardizes the seasonal variance.

set.seed(38472)
myseries <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

myseries |>
  autoplot(Turnover) +
  labs(title = "Queensland Footwear and Personal Accessory Turnover")

lambda <- myseries |>
  features(Turnover, features = guerrero) |>
  pull(lambda_guerrero)

trans_retail <- myseries |>
  mutate(trans_turnover = box_cox(Turnover, lambda)) |>
  select(!Turnover)

trans_retail |>
  autoplot(trans_turnover) +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed Retail Turnover with $\\lambda$ = ",
         round(lambda,2))))

Beginning with a round of seasonal differencing and another round of ordinary differencing we get the dataset below. The data appears largely normal looking at the data plot, but we do notice both the ACF and PACF plots have multiple observations rising above the level of significance although the correlation at the first lag is not positive. Overall we can conclude the data is fairly stationary, but there are enough indicators present to warrant caution when developing models relying on startionary data.

trans_retail |>
  features(trans_turnover, unitroot_nsdiffs)
## # A tibble: 1 × 3
##   State      Industry                                        nsdiffs
##   <chr>      <chr>                                             <int>
## 1 Queensland Footwear and other personal accessory retailing       1
s_dif_retail <- trans_retail |>
  mutate(s_dif_turnover = difference(trans_turnover, 12))

s_dif_retail |>
  features(s_dif_turnover, unitroot_ndiffs)
## # A tibble: 1 × 3
##   State      Industry                                        ndiffs
##   <chr>      <chr>                                            <int>
## 1 Queensland Footwear and other personal accessory retailing      1
s_dif_retail |>
  mutate(dif_turnover = difference(s_dif_turnover, 1)) |>
  gg_tsdisplay(dif_turnover, plot_type = "partial")

9.6

9.6.1

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)

9.6.2

As we increase the value of \(\phi_1\) we can see an increasing prevalence of the trend as well as an increase in the variance which can be seen best in the differing y scale of the first and last graphs.

for(i in 2:100)
  y[i] <- 0*y[i-1] + e[i]

plot_0 <- tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("AR(1) model with  $\\phi_1$ = 0")))

for(i in 2:100)
  y[i] <- 0.3*y[i-1] + e[i]

plot_03 <- tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("AR(1) model with  $\\phi_1$ = 0.3")))

plot_06 <- sim |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("AR(1) model with  $\\phi_1$ = 0.6")))

for(i in 2:100)
  y[i] <- 1*y[i-1] + e[i]

plot_1 <- tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("AR(1) model with  $\\phi_1$ = 1")))

grid.arrange(plot_0, plot_03, plot_06, plot_1, nrow = 2)

9.6.3

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

sim_ma <- tsibble(idx = seq_len(100), y = y, index = idx)

9.6.4

Plotting our MA(1) series and changing the value of \(\phi_1\) we can see that as \(\phi_1\) increases the variance increases but the trend stays the same. Additionally, as \(\phi_1\) increases the observation to observation variance smooths out leading to a less “spiky” plot.

for(i in 2:100)
   y[i] <- e[i] + 0.0 * e[i-1]

plot_0 <- tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("MA(1) model with  $\\phi_1$ = 0")))

for(i in 2:100)
   y[i] <- e[i] + 0.3 * e[i-1]

plot_03 <- tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("MA(1) model with  $\\phi_1$ = 0.3")))

plot_06 <- sim_ma |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("MA(1) model with  $\\phi_1$ = 0.6")))

for(i in 2:100)
   y[i] <- e[i] + 1 * e[i-1]

plot_1 <- tsibble(idx = seq_len(100), y = y, index = idx) |>
  autoplot(y) +
  labs(title = latex2exp::TeX(paste0("MA(1) model with  $\\phi_1$ = 1")))

grid.arrange(plot_0, plot_03, plot_06, plot_1, nrow = 2)

9.6.5

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

sim_arma <- tsibble(idx = seq_len(100), y = y, index = idx)

9.6.6

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)

9.6.7

Looking at the two series we can see neither is stationary. The ARMA(1,1) model has a gradual decline in the ACF plot although the PACF does only have two observations with significance indicating that the first lag is a strong predictor of the current value. The plot of the data does not have a high degree of variance but does show some obvious trends. The AR(2) plot shows continuously increasing variance around a stable mean. The alternating and increasing values are reflected in the ACF and PACF as one would expect.

sim_arma |>
  gg_tsdisplay(y, plot_type = "partial") +
  labs(title = latex2exp::TeX(paste0("ARMA(1,1) model with  $\\phi_1 = 0.6$, $\\theta_1 = 0.6$, and $\\sigma^2 = 1$")))

sim_ar2 |>
  gg_tsdisplay(y, plot_type = "partial") +
  labs(title = latex2exp::TeX(paste0("AR(2) model with  $\\phi_1 = -0.8$, $\\phi_2 = 0.6$, and $\\sigma^2 = 1$")))

9.7

9.7.1

Looking at the raw passenger data we see that there is a positive trend over the time period and no evidence of seasonality. The slight bowing of the trend indicates a Box-Cox transformation would be beneficial, but following the instructions of this question no transformation will be applied. The fable ARIMA model chose an ARIMA(0,2,1) model.

aus_airpassengers |>
  autoplot(Passengers) +
  labs(title = "Passengers  on Australian Air Carriers")

arima <- aus_airpassengers |>
  model(ARIMA(Passengers))

report(arima)
## 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

The residuals for the model look very good, there are a couple large spikes in the innovation residuals but overall they look like white noise. The residuals are fairly normally distributed with the exception of a single outlier and the ACF plot shows no significant correlations at any lag. Thus we can conclude the residuals look like white noise and can have confidence in the accuracy of our 10 period forecast.

arima |> 
  forecast(h=10) |>
  autoplot(aus_airpassengers) + 
  labs(title = "Modeling Passenger Numbers using ARIMA(0,2,1)")

arima |>
  gg_tsresiduals() +
  labs(title = "Residuals for ARIMA(0,2,1)")

9.7.2

The ARIMA(0,2,1) model can be written in terms of the backshift operator as:

\((1-B)^2y_t = -0.8963\varepsilon_{t-1} + \varepsilon_t\)

9.7.3

The ARIMA(0,1,0) with drift model has a very similar forecast compared to the ARIMA(0,2,1) model. The new model has a less steep trend line to its forecasts and tighter confidence intervals than the autoselected model.

arima_dif <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 1 + pdq(0,1,0)))

arima_dif |> 
  forecast(h=10) |>
  autoplot(aus_airpassengers) + 
  labs(title = "Modeling Passenger Numbers using ARIMA(0,1,0) with drift")

9.7.4

The ARIMA(2,1,2) has a more interesting behavior than the prior models. The slope of the forecast is no longer constant and appears to match the behavior of the series better, although this does not necessarily make for a better model

arima_d <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 1 + pdq(2,1,2)))

arima_d |> 
  forecast(h=10) |>
  autoplot(aus_airpassengers) + 
  labs(title = "Modeling Passenger Numbers using ARIMA(2,1,2)")

Removing the constant appears to make the data unforecastable due to the non-stationary AR component.

arima_nc <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 0 + pdq(2,1,2)))
## Warning: 1 error encountered for ARIMA(Passengers ~ 0 + pdq(2, 1, 2))
## [1] non-stationary AR part from CSS
arima_nc |> 
  forecast(h=10)
## # A fable: 10 x 4 [1Y]
## # Key:     .model [1]
##    .model                                Year Passengers .mean
##    <chr>                                <dbl>     <dist> <dbl>
##  1 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2017         NA    NA
##  2 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2018         NA    NA
##  3 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2019         NA    NA
##  4 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2020         NA    NA
##  5 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2021         NA    NA
##  6 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2022         NA    NA
##  7 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2023         NA    NA
##  8 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2024         NA    NA
##  9 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2025         NA    NA
## 10 ARIMA(Passengers ~ 0 + pdq(2, 1, 2))  2026         NA    NA

9.7.5

Adding a constant to the ARIMA(0,2,1) model causes the point forecast to become quadratic, which is generally an undesirable trait in a model.

arima <- aus_airpassengers |>
  model(ARIMA(Passengers ~ 1 + pdq(0,2,1)))
## Warning: Model specification induces a quadratic or higher order polynomial trend. 
## This is generally discouraged, consider removing the constant or reducing the number of differences.
arima |> 
  forecast(h=10) |>
  autoplot(aus_airpassengers) + 
  labs(title = "Modeling Passenger Numbers using ARIMA(0,2,1) with a constant")

9.8

9.8.1

Looking at the US GDP data we see a slight exponential trend which we can address with a Box-Cox transformation of 0.28.

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

us_gdp |> gg_tsdisplay(GDP, plot_type = "partial") + labs(title = "Raw US GDP Data")

lambda <- us_gdp |>
  features(GDP, features = guerrero) |>
  pull(lambda_guerrero)

us_gdp |>
  gg_tsdisplay(box_cox(GDP, lambda), plot_type = "partial") +
  labs(y = "",
       title = latex2exp::TeX(paste0(
         "Transformed US GDP with $\\lambda$ = ",
         round(lambda,2))))

9.8.2

arima_gdp <- us_gdp |>
  model(ARIMA(box_cox(GDP, lambda)))

report(arima_gdp)
## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: box_cox(GDP, lambda) 
## 
## Coefficients:
##          ar1  constant
##       0.4586  118.1822
## s.e.  0.1198    9.5047
## 
## sigma^2 estimated as 5479:  log likelihood=-325.32
## AIC=656.65   AICc=657.1   BIC=662.78

9.8.3

Testing how many rounds of differentiating are needed we also arrive at 1. Looking at the ACF and PACF of the transformed data made earlier both plots do indicate an AR(1) component would be appropriate, but for the sake of comparison we can compare to several models that include an MA component with and without drift. Testing the models we arrive at the conclusion the ARIMA(0,2,1) without drift model is actually the best performing.

us_gdp |>
  features(box_cox(GDP, lambda), unitroot_ndiffs)
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1
custom_model <- us_gdp |>
  model(arima120nd = ARIMA(box_cox(GDP, lambda) ~ pdq(1,2,0)),
        arima120d = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(1,2,0)),
        auto = ARIMA(box_cox(GDP, lambda)),
        arima110nd = ARIMA(box_cox(GDP, lambda) ~ pdq(1,1,0)),
        arima021nd = ARIMA(box_cox(GDP, lambda) ~ pdq(0,2,1)),
        arima021d = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0,2,1)),
        arima011nd = ARIMA(box_cox(GDP, lambda) ~ pdq(0,1,1)),
        arima011d = ARIMA(box_cox(GDP, lambda) ~ 1 + pdq(0,1,1)))

glance(custom_model) |> arrange(AICc) |> select(.model:BIC)
## # A tibble: 8 × 6
##   .model     sigma2 log_lik   AIC  AICc   BIC
##   <chr>       <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 arima021nd  6020.   -323.  650.  650.  654.
## 2 arima021d   6115.   -323.  652.  652.  658.
## 3 arima120nd  6780.   -326.  656.  656.  660.
## 4 auto        5479.   -325.  657.  657.  663.
## 5 arima110nd  5479.   -325.  657.  657.  663.
## 6 arima120d   6898.   -326.  658.  658.  664.
## 7 arima011nd  5689.   -326.  659.  659.  665.
## 8 arima011d   5689.   -326.  659.  659.  665.

9.8.4

Looking at the residuals for the ARIMA(0,2,1) no drift model we can see they look relatively good. The innovation residuals largely match what we would expect from white noise although there are a couple large negative spikes. There are no significant autocorrelations, and while the plot of the residuals looks fairly left skewed this is more an artifact of the bin width of the plot and a few large outliers and the actual distribution of residuals is more normal than the plot indicates, although still slightly skewed left.

custom_model |>
  select(arima021nd) |>
  gg_tsresiduals() +
  labs(title = "ARIMA(0,2,1) no drift")

9.8.5

The forecasts produced by our model look very good.

custom_model |>
  forecast(h=10) |>
  filter(.model == 'arima021nd') |>
  autoplot(us_gdp)

9.8.6

While the numbers are a tad large for the measures in the same scale as the data we can see that the ARIMA model outperformed the ETS model with no transformation by all measures.

gdp <- us_gdp |>
  slice(-n()) |>
  stretch_tsibble(.init = 10, .step = 1)

forecast_compare <- gdp |>
  model(arima110nd = ARIMA(box_cox(GDP, lambda) ~ pdq(1,1,0)),
        ets = ETS(GDP)) |>
  forecast(h = 1)

forecast_compare |>
  accuracy(us_gdp) |>
  select(.model, RMSE:MAPE)
## # A tibble: 2 × 5
##   .model              RMSE           MAE   MPE  MAPE
##   <chr>              <dbl>         <dbl> <dbl> <dbl>
## 1 arima110nd 173955083878. 109850422846. 0.294  1.61
## 2 ets        190604101250. 126812323851. 0.594  1.73