library(fpp3)
## Warning: package 'fpp3' was built under R version 4.3.3
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──
## ✔ tibble      3.2.1     ✔ tsibble     1.1.5
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.4.1
## ✔ lubridate   1.9.4     ✔ fable       0.4.1
## ✔ ggplot2     3.5.1
## Warning: package 'lubridate' was built under R version 4.3.3
## Warning: package 'tsibble' was built under R version 4.3.3
## Warning: package 'feasts' was built under R version 4.3.3
## Warning: package 'fabletools' was built under R version 4.3.3
## Warning: package 'fable' was built under R version 4.3.3
## ── 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()

Question 9.1

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

Part A

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

The main differences among these figures comes from the different sample sizes. The first shows a small sample size with only 36 observations so we see more variance and fluctuations with a large confidence bounds. Then, the second shows a bigger sample size of 360 observations and we can see how there are smaller fluctuations and the values are closer to zero with a smaller confidence bounds. For the last figure we see an even larger sample size of 1,000 observations and we can see how the values are even closer to zero more stable and has even smaller confidence bounds.

Part B

Why are the critical values at different distances from the mean of zero? Why are the auto correlations different in each figure when they each refer to white noise?

The critical values are at different distances from the mean of zero because the critical values are calculated based on the sample size and for smaller sample sizes the critical values are further away from 0 making the interval wider this is because smaller sample sizes have more variation than larger sample sizes. The auto correlations are different in each figure because even tho they all refer to white noise there is more random variation in smaller sample sizes so as the sample size increases the data will become more stable and the ACF will start to get closer to the ACF that represents white noise.

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

amazon_stock <- gafa_stock %>%
  filter(Symbol == 'AMZN') %>%
  select(Symbol, Close, Date)
amazon_stock %>%
  gg_tsdisplay(Close, plot_type='partial') +
  labs(title = "Amazon Closing Stock Price")
## Warning: Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.
## Provided data has an irregular interval, results should be treated with caution. Computing ACF by observation.

The three plots above for the Amazon stock closing price show that the Amazon stock close price is non-stationary. The first plot of the time series shows that there is a positive trend in the data which right off the bat makes it non stationary. Then the second plot that is the ACF plot shows the values slowly decreasing which shows that the data is non stationary if the data was stationary we would see the values quickly go to zero and they would not start out so large close to one. Finally, the last plot, the pcaf plot we see that the first lag is very close to one this shows again that the data is non-stationary and that it should be differenced.

Question 9.3

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

Part A

Turkish GDP from global_economy.

turk_gdp <- global_economy %>%
  filter(Country == 'Turkey') %>%
  select(Year, Country, GDP)

turk_gdp %>%
  autoplot(GDP) + 
  labs(title = "Turkish GDP")

We can see that there is a clear upward trend so we will need to apply differencing to make it stationary.

turk_gdp_lambda <- turk_gdp %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

turk_gdp %>%
  autoplot(box_cox(GDP,turk_gdp_lambda)) + 
  labs(title = "Turkish GDP transformed")

After applying the box cox we can see there is still trend so we will need to apply differencing to make the data stationary.

turk_gdp %>%
  features(box_cox(GDP,turk_gdp_lambda), unitroot_ndiffs) 
## # A tibble: 1 × 2
##   Country ndiffs
##   <fct>    <int>
## 1 Turkey       1

The KPSS test above calls for one order of differencing to make the data stationary.

turk_gdp <- turk_gdp %>%
  mutate(transformed_gdp = difference(box_cox(GDP, turk_gdp_lambda)))

turk_gdp %>%
  gg_tsdisplay(transformed_gdp, plot_type='partial') +
  labs(title = "Turkish GDP")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

Now we can see that the box cox transformed and differenced GDP is now stationary. The time series plot shows no trend or seasonality, the ACF values drop to zero quickly and the pacf first lag is close to zero and not one.

Part B

Accommodation takings in the state of Tasmania from aus_accommodation.

tasmania_takings <- aus_accommodation %>%
  filter(State == "Tasmania") %>%
  select(Date, State, Takings)

tasmania_takings %>%
  autoplot(Takings)

We can see there is clear trend and seasonality in the data.

tasmania_takings_lambda <- tasmania_takings %>%
  features(Takings, features = guerrero) %>%
  pull(lambda_guerrero)

tasmania_takings %>%
  autoplot(box_cox(Takings,tasmania_takings_lambda)) + 
  labs(title = "Tasmania Takings transformed")

tasmania_takings %>%
  features(box_cox(Takings,tasmania_takings_lambda), unitroot_ndiffs) 
## # A tibble: 1 × 2
##   State    ndiffs
##   <chr>     <int>
## 1 Tasmania      1

The KPSS test calls for one differencing. Since we can see strong seasonality in the data we will do a seasonal differencing.

tasmania_takings <- tasmania_takings %>%
  mutate(transformed_takings = difference(box_cox(Takings, tasmania_takings_lambda), 4))

tasmania_takings %>%
  gg_tsdisplay(transformed_takings, plot_type='partial') +
  labs(title = "Tasmania Takings")
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

Looking at the plots above the transformed and differenced data seems to be stationary. The time series plot shows no trend and no seasonality, the ACF plot values go to zero fast.

Part C

Monthly sales from souvenirs

souvenirs %>%
  autoplot(Sales)

We can see that there is some seasonality and trend and the variance increases over time in the data making it non-stationary.

sales_lambda <- souvenirs %>%
  features(Sales, features = guerrero) %>%
  pull(lambda_guerrero)

souvenirs %>%
  autoplot(box_cox(Sales,sales_lambda)) + 
  labs(title = "Monthly Sales transformed")

There is still some seasonality so lets check to see if a difference is needed

souvenirs %>%
  features(box_cox(Sales,sales_lambda), unitroot_ndiffs) 
## # A tibble: 1 × 1
##   ndiffs
##    <int>
## 1      1
souvenirs <- souvenirs %>%
  mutate(transformed_sales = difference(box_cox(Sales, sales_lambda), 12))

souvenirs %>%
  gg_tsdisplay(transformed_sales, plot_type='partial') +
  labs(title = "Monthly Sales")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

Question 9.5

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

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

Looking at the turnover data we can see that it is not stationary it has a clear trend and strong seasonality.

retail_lambda <- myseries %>%
  features(Turnover, features = guerrero) %>%
  pull(lambda_guerrero)

myseries %>%
  autoplot(box_cox(Turnover,retail_lambda)) + 
  labs(title = "Retail Turnover transformed")

myseries %>%
  features(box_cox(Turnover,retail_lambda), unitroot_ndiffs) 
## # A tibble: 1 × 3
##   State           Industry                                         ndiffs
##   <chr>           <chr>                                             <int>
## 1 South Australia Hardware, building and garden supplies retailing      1
myseries <- myseries %>%
  mutate(transformed_turnover = difference(box_cox(Turnover, retail_lambda), 12))

myseries %>%
  gg_tsdisplay(transformed_turnover, plot_type='partial') +
  labs(title = "Monthly Retail Turnover")
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).

After applying seasonal difference we can see that the data is now stationary as theres no trend or seasonality in the time series, the ACF values quickly go to zero.

Question 9.6

Simulate and plot some data from simple ARIMA models.

Part 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 <- 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)

Part B

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

sim %>%
  autoplot(y)

y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 1.2*y[i-1] + e[i]
sim2 <- tsibble(idx = seq_len(100), y = y, index = idx)
sim2 %>%
  autoplot(y)

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

Part C

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

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

Part D

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

simC %>%
  autoplot(y)

y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- e[i] + 1.2 * e[i-1]
simC2 <- tsibble(idx = seq_len(100), y = y, index = idx)
simC2 %>%
  autoplot(y)

y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- e[i] + -0.3 * e[i-1]
simC3 <- tsibble(idx = seq_len(100), y = y, index = idx)
simC3 %>%
  autoplot(y)

As theta changes we can see that the range of values changes as theta decreases we can see the range of values decrease. We can also see that as theta increases the spikes become less frequent while when its smaller we see more frequent spikes and the variance seems more consistent than with larger values of theta.

Part E

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

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

Part 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.)

y <- numeric(100)
e <- rnorm(100)
for(i in 3:100)
  y[i] <- -0.8 * y[i-1] + 0.3 * y[i-2] + e[i]
simF <- tsibble(idx = seq_len(100), y = y, index = idx)
simF %>%
  autoplot(y)

Part G

Graph the latter two series and compare them.

simE %>%
  gg_tsdisplay(y, plot_type='partial') +
  labs(title = "ARMA(1,1) Model")

simF %>%
  gg_tsdisplay(y, plot_type='partial') +
  labs(title = "AR(2) Model")

Looking at the plots above we can see that the ARIMA(1,1) model appears to be stationary. The time plot shows no trend or seasonality, the ACF plot shows the values quickly going to zero and the PACF plot has a significant spike at lag 1. The ARIMA(2) model appears to not be stationary, the time plot shows values oscillating and the variance increases as time goes on. The ACF plot shows values not going to zero but oscillating between significant negative and positive values. The PACF plot shows one significant spike at lag one and the rest are 0. This could be due to the stationary condition of ϕ2 - ϕ1 > 1 not being met.

Question 9.7

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

Part A

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.

aus_airpassengers <- aus_airpassengers %>%
  filter(Year < 2012)
aus_airpassengers %>%
  autoplot(Passengers)

arima_partA_fit <- aus_airpassengers %>%
  model(ARIMA(Passengers))

report(arima_partA_fit)
## Series: Passengers 
## Model: ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -0.8756
## s.e.   0.0722
## 
## sigma^2 estimated as 4.671:  log likelihood=-87.8
## AIC=179.61   AICc=179.93   BIC=182.99

From the output above we can see that the model that was selected was a ARIMA(0, 2, 1) model.

arima_partA_fit  %>%
  gg_tsresiduals()

Looking at the ACF plot above none of the points are significant so the residuals do look like they are white noise.

arima_partA_fit %>% 
  forecast(h=10) %>%
  autoplot(aus_airpassengers) +
  labs(title = "Australian Aircraft Passengers ARIMA(0,2,1) Model Forecasts", y = "Passengers (in millions)")

Part B

Write the model in terms of the backshift operator.

We have a ARIMA(0,2,1) model so we have no auto regressive term, second order differencing and one moving average.

For the second order differencing: \[ (1 -B)^2 y_t \]

Then for the moving average part: \[ (1 - 0.8963B) e_t \]

putting it all together we get: \[ (1 -B)^2 y_t = (1 - 0.8963B) e_t \]

Part C

Plot forecasts from an ARIMA(0,1,0) model with drift and compare these to part a.

arima_partC_fit <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ pdq(0,1,0)))
arima_partC_fit  %>%
  gg_tsresiduals()

arima_partC_fit %>%
  forecast(h=10) %>%
  autoplot(aus_airpassengers)

report(arima_partC_fit)
## Series: Passengers 
## Model: ARIMA(0,1,0) w/ drift 
## 
## Coefficients:
##       constant
##         1.3669
## s.e.    0.3319
## 
## sigma^2 estimated as 4.629:  log likelihood=-89.08
## AIC=182.17   AICc=182.48   BIC=185.59

Looking at the results above the ARIMA(0,1,0) model has a slightly higher AICc of 182.48 compared to the ARIMA(0,2,1) model which had a AICc of 179.93 . The forecasts look pretty similar looking at the plots. Also, the residuals still look like white noise.

Part D

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.

arima_partD_fit <- aus_airpassengers %>%
  model(ARIMA(Passengers ~ pdq(2,1,2)))
arima_partD_fit  %>%
  gg_tsresiduals()

arima_partD_fit %>%
  forecast(h=10) %>%
  autoplot(aus_airpassengers)

report(arima_partD_fit)
## Series: Passengers 
## Model: ARIMA(2,1,2) w/ drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  constant
##       1.4694  -0.5103  -1.5736  0.6780    0.0650
## s.e.  0.3780   0.3558   0.3081  0.2688    0.0294
## 
## sigma^2 estimated as 4.748:  log likelihood=-87.74
## AIC=187.47   AICc=189.94   BIC=197.75

From the output above we can see that the ARIMA(2,1,2) model gave a AICc of 189.94 which was higher than the models in parts A and C which gave AICcs of 179.93 and 182.48 meaning this model was the worst of the three. We can also see that the residuals are still white noise.

arima_partD_no_constant_fit <- 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
report(arima_partD_no_constant_fit)
## Series: Passengers 
## Model: NULL model 
## NULL model

By removing the constant we can see that the model becomes the null model.

Part D

Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?

arima_partE_fit <- 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_partE_fit  %>%
  gg_tsresiduals()

arima_partE_fit %>%
  forecast(h=10) %>%
  autoplot(aus_airpassengers)

report(arima_partE_fit)
## Series: Passengers 
## Model: ARIMA(0,2,1) w/ poly 
## 
## Coefficients:
##           ma1  constant
##       -1.0000    0.0721
## s.e.   0.0709    0.0260
## 
## sigma^2 estimated as 4.086:  log likelihood=-85.74
## AIC=177.48   AICc=178.15   BIC=182.55

By generating a ARIMA(0,2,1) model with a constant we get a warning that the model has a quadratic or higher order polynomial and it discourages it. When we plot the forecasts we can see that they are higher values and the forecasts line becomes much steeper. However, we have a good AICc of 178.15 which is actually better than all other models that we created.

Question 9.8

For the United States GDP series (from global_economy):

Part A

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

usa_gdp <- global_economy %>%
  filter(Country == "United States") %>%
  select(Country, GDP, Year)
 
usa_gdp %>% 
  gg_tsdisplay(GDP, plot_type='partial')

We can see that the data is clearly not stationary as it has a positive trend. So, apply a box-cox transformation to the data

usa_gdp_lambda <- usa_gdp %>%
  features(GDP, features = guerrero) %>%
  pull(lambda_guerrero)

usa_gdp <- usa_gdp %>%
  mutate(transformed_gdp = box_cox(GDP, usa_gdp_lambda))

usa_gdp %>%
  autoplot(transformed_gdp)

Part B

fit a suitable ARIMA model to the transformed data using ARIMA();

usa_arima <- usa_gdp %>%
  model(ARIMA(transformed_gdp))

report(usa_arima)
## Series: transformed_gdp 
## Model: ARIMA(1,1,0) w/ drift 
## 
## 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

Part C

try some other plausible models by experimenting with the orders chosen;

usa_gdp %>%
  gg_tsdisplay(transformed_gdp, plot_type = "partial")

Looking at the box-cox transformed data we can see that there is still a clear trend and the ACF plot steadily goes to zero over multiple lags so it is non stationary and differencing will be needed at least to the first order. We can also see that the PACF plot has a huge spike at the first lag meaning it could be an AR(1) model. I will try a ARIMA(1,1,0) model which is the one that was suggested and I will also try a ARIMA(1,2,0) and ARIMA(1,1,1) model to see what happens.

usa_arima_models <- usa_gdp %>%
  model(
    arima1 = ARIMA(transformed_gdp ~ pdq(1,1,0)),
    arima2 = ARIMA(transformed_gdp ~ pdq(1,2,0)),
    arima3 = ARIMA(transformed_gdp ~ pdq(1,1,1))
    )

report(usa_arima_models)
## Warning in report.mdl_df(usa_arima_models): Model reporting is only supported
## for individual models, so a glance will be shown. To see the report for a
## specific model, use `select()` and `filter()` to identify a single model.
## # A tibble: 3 × 9
##   Country       .model sigma2 log_lik   AIC  AICc   BIC ar_roots  ma_roots 
##   <fct>         <chr>   <dbl>   <dbl> <dbl> <dbl> <dbl> <list>    <list>   
## 1 United States arima1  5479.   -325.  657.  657.  663. <cpl [1]> <cpl [0]>
## 2 United States arima2  6780.   -326.  656.  656.  660. <cpl [1]> <cpl [0]>
## 3 United States arima3  5580.   -325.  659.  659.  667. <cpl [1]> <cpl [1]>

Looking at the results from the models created it looks like the ARIMA(1,1,0) and ARIMA(1,2,0) models performed the best as they gave the best AICcs of 657.1005 and 656.2160 since the ARIMA(1,2,0) model gave a slighty better AICc I will choose this as the “best” model.

Part D

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

usa_arima_models %>%
  select(arima2) %>%
  gg_tsresiduals()

Looking at the residuals for the ARIMA(1,2,0) model they appear to be white noise. There is only one value outside the confidence interval for the ACF plot.

Part E

produce forecasts of your fitted model. Do the forecasts look reasonable?

usa_graph <- usa_gdp %>%
  select(Year, transformed_gdp)
usa_arima_models %>%
  select(arima2) %>%
  forecast(h=25) %>%
  autoplot(usa_graph)

The above plot shows the plotted forecasts with the original transformed GDP data the forecasts appear to be reasonable.

Part E

compare the results with what you would obtain using ETS() (with no transformation).

ets_model <- usa_gdp %>%
  model(ETS(transformed_gdp))

report(ets_model)
## Series: transformed_gdp 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.9998985 
##     beta  = 0.3731217 
## 
##   Initial states:
##      l[0]     b[0]
##  7045.398 198.9275
## 
##   sigma^2:  0
## 
##      AIC     AICc      BIC 
## 744.1762 745.3300 754.4784
ets_model %>%
  forecast(h=25) %>%
  select(-Country) %>%
  autoplot(usa_graph)

Looking at the output for the ETS model it attained a AICc of 745.3300 which was significantly higher than the AICc of 656.2160 that was attained by the ARIMA(1,2,0) model. The forecasts appear reasonable.