Figure 9.32 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?
All of the figures indicate white noise because all of the autocorrelations stay within the blue dashed lines. The differences between the figures are that as the sample gets bigger from 36 to 360 random numbers, the blue dashed lines get closer together and the autocorrelation spikes get smaller. The blue dashed lines get even closer together and the autocorrelation spikes get even smaller once there are 1000 random numbers.
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?
The critical values are at different distances from the mean of zero because as the sample size of random numbers increases, the bounds get tighter. This is due to the fact that the bounds are set at ±1.96/√T where T is the sample size. This means that as the sample increases, the denominator gets bigger as well resulting in a smaller interval.
The autocorrelations are different in each figure when they each refer to white noise because the sample size is different. Since the sample with only 36 numbers is smaller than the sample with 360 or 1000 samples, the autocorrelations will have bigger fluctuations. The same is true in reverse for the bigger samples.
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.
This plot shows that the data is non-stationary because there clearly is a trend present in the data. It trends positively for the most part and then has a negative trend at the end.
Warning: Provided data has an irregular interval, results should be treated
with caution. Computing ACF by observation.
This ACF shows that the data is non-stationary because the ACF does not drop to 0 relatively quickly. Instead, it starts at 1 and then only slowly declines, so slow that at 30 lag, it still looks to be around 0.9.
The original plot shows a positive trend and the acf shows multiple lags where it goes past the blue dashed lines. The pacf also shows a close to 1 value on lag 1. This means that this dataset is non-stationary.
turkey_lambda <- turkey |>features(GDP, features = guerrero) |>pull(lambda_guerrero)turkey_lambda
[1] 0.1572187
I used the Guerrero feature to find the lambda value of 0.1572187.
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()`).
Comparing the plot of the dataset differenced to the original is a night and day difference. Differencing once made the dataset stationary and the acf and pacf plots support that.
b. Accommodation takings in the state of Tasmania from aus_accommodation.
The plot shows a dataset that is seasonal with a slightly positive trend. The acf plot shows multiple lags that go past the blue dashed lines. The pacf plot also shows multiple lags that go past the blue dashed lines.
tasmania_lambda <- tasmania |>features(Takings, features = guerrero) |>pull(lambda_guerrero)tasmania_lambda
[1] 0.001819643
Using the Guerrero method, I got a lambda value of 0.001819643 for the box-cox.
# A tibble: 1 × 3
State nsdiffs ndiffs
<chr> <int> <int>
1 Tasmania 1 1
Using the unitroot_nsdiffs and unitroot_ndiffs features, I got an estimated 1 for the number of differences needed to make the data stationary as well as a 1 for the number of seasonal differences. Since this is quarterly data, seasonal difference would mean a lag of 4.
tasmania |>mutate(d_bc_Takings =difference(box_cox(Takings, tasmania_lambda), lag =4)) |>features(d_bc_Takings, unitroot_ndiffs)
# A tibble: 1 × 2
State ndiffs
<chr> <int>
1 Tasmania 0
After taking the seasonal difference and recalculating with unitroot_ndiffs, it looks like a normal difference is not needed anymore.
tasmania |>mutate(d_bc_Takings =difference(box_cox(Takings, tasmania_lambda), lag =4)) |>gg_tsdisplay(d_bc_Takings, plot_type ="partial")
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()`).
The plot looks to be stationary with the acf dropping rapidly to near zero.
Recalculating for unitroot_nsdiffs and unitroot_ndiffs after taking the seasonal difference shows that no more difference is needed.
souvenirs |>mutate(d_bc_Sales =difference(box_cox(Sales, souvenirs_lambda), lag =12)) |>gg_tsdisplay(d_bc_Sales, plot_type ="partial")
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()`).
Doing a visual test, it shows that the data is no longer trending or seasonal. Also, the acf drops rapidly to 0 and no longer has the spike at lag 12. This shows that the data is now stationary.
Exercise 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.
# A tibble: 1 × 4
State Industry nsdiffs ndiffs
<chr> <chr> <int> <int>
1 New South Wales Other retailing n.e.c. 1 1
Based on the unitroot_nsdiffs and unitroot_ndiffs features, I will use a seasonal difference with lag 12 since its monthly data.
myseries |>mutate(d_bc_Turnover =difference(box_cox(Turnover, retail_lambda), lag =12)) |>features(d_bc_Turnover, list(unitroot_nsdiffs, unitroot_ndiffs))
# A tibble: 1 × 4
State Industry nsdiffs ndiffs
<chr> <chr> <int> <int>
1 New South Wales Other retailing n.e.c. 0 0
After using the seasonal difference of lag 12 and recalculating, it looks like no more difference is required.
myseries |>mutate(d_bc_Turnover =difference(box_cox(Turnover, souvenirs_lambda), lag =12)) |>gg_tsdisplay(d_bc_Turnover, plot_type ="partial")
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()`).
The plot shows that the acf goes to near 0 although it does take until lag 11 to get below the blue dashed line. For my retail data it required just taking the seasonal difference with a lag of 12 to become stationary.
Exercise 9.6
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 <-numeric(100)e <-rnorm(100)for (i in2:100) { y[i] <-0.6* y[i -1] + e[i]}sim <-tsibble(idx =seq_len(100), y = y, index = idx)
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 in2:100) { y[i] <-0* y[i -1] + e[i]}sim <-tsibble(idx =seq_len(100), y = y, index = idx)sim |>autoplot(y)
The plot with a phi value of 0 looks like its white noise.
y <-numeric(100)e <-rnorm(100)for (i in2:100) { y[i] <-0.95* y[i -1] + e[i]}sim <-tsibble(idx =seq_len(100), y = y, index = idx)sim |>autoplot(y)
The plot with a phi value of 0.95 looks like it has a positive trend and a much smoother plot.
y <-numeric(100)e <-rnorm(100)for (i in2:100) { y[i] <--0.5* y[i -1] + e[i]}sim <-tsibble(idx =seq_len(100), y = y, index = idx)sim |>autoplot(y)
The plot with a phi value of -0.5 looks like white noise with a lot of variance from point to point.
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 in2:100) { y[i] <- e[i] +0.6* e[i -1]}sim_ma <-tsibble(idx =seq_len(100), y = y, index = idx)
d. Produce a time plot for the series. How does the plot change as you change θ1?
sim_ma |>autoplot(y)
y <-numeric(100)e <-rnorm(100)for (i in2:100) { y[i] <- e[i] +0.95* e[i -1]}sim_ma <-tsibble(idx =seq_len(100), y = y, index = idx)sim_ma |>autoplot(y)
Increasing the theta value to 0.95 seems to make it smoother with a larger interval between the biggest and smallest y values.
y <-numeric(100)e <-rnorm(100)for (i in2:100) { y[i] <- e[i] +-0.5* e[i -1]}sim_ma <-tsibble(idx =seq_len(100), y = y, index = idx)sim_ma |>autoplot(y)
Using a theta value of -.05 looks to do the opposite by adding even more variance to the plot.
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 in2:100) { y[i] <-0.6* y[i -1] + e[i] +0.6* e[i -1]}sim_arma <-tsibble(idx =seq_len(100), y = y, index = idx)
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 in3: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)
g. Graph the latter two series and compare them.
sim_arma |>gg_tsdisplay(y, plot_type ="partial")
sim_ar2 |>gg_tsdisplay(y, plot_type ="partial")
The plots for the latter 2 series are completely different. The plot for the ARMA(1, 1) looks stationary since there doesn’t appear to be any trend or seasonality, the acf and pacf supports this. The plot for the AR(2) looks like there’s a clear symmetrical pattern and non-stationary since the variance is not constant. The acf plot for AR(2) shows almost every point going past the blue dashed lines.
Exercise 9.7
Consider aus_airpassengers, the total number of passengers (in millions) from Australian air carriers for the period 1970-2011.
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.
ap <- aus_airpassengers |>filter(Year <=2011)
fit_airpassengers <- ap |>model(ARIMA(Passengers))report(fit_airpassengers)
The plot from part a had a much higher projection compared to the plot from part c. This is due to the drift in plot c that takes into account the average change for the whole 1970-2011 period.
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.
fit_d_drift <- ap |>model(ARIMA(Passengers ~1+pdq(2,1,2)))report(fit_d_drift)
fit_d_drift |>forecast(h =10) |>autoplot(ap) +labs(title ="ARIMA(2,1,2) with drift")
The ARIMA(2,1,2) gives a similar forecast to part a and a higher forecast than part c.
fit_d_noconst <- ap |>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(fit_d_noconst)
Series: Passengers
Model: NULL model
NULL model
fit_d_noconst |>forecast(h =10) |>autoplot(ap) +labs(title ="ARIMA(2,1,2) without constant")
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
-Inf
Warning: Removed 10 rows containing missing values or values outside the scale range
(`geom_line()`).
Removing the constant seems to have resulted in a null model.
e. Plot forecasts from an ARIMA(0,2,1) model with a constant. What happens?
fit_e <- ap |>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.
The plot for the ARIMA(0,2,1) model with a constant looks to appear steeper and gives a warning that suggests removing the constant because the model specification induces a quadratic or higher order polynomial trend.
Exercise 9.8
For the United States GDP series (from global_economy):
a. if necessary, find a suitable Box-Cox transformation for the data;
us <- global_economy |>filter(Country =="United States")us |>autoplot(GDP)
us_lambda <- us |>features(GDP, features = guerrero) |>pull(lambda_guerrero)us_lambda
[1] 0.2819443
Using the Guerrero feature, I got a lambda value of 0.2819443.
us |>autoplot(box_cox(GDP, us_lambda))
Using the lambda value, I got the box-cox transformation above.
b. fit a suitable ARIMA model to the transformed data using ARIMA();
fit_us <- us |>model(ARIMA(box_cox(GDP, us_lambda)))report(fit_us)
d. choose what you think is the best model and check the residual diagnostics;
I think the ARIMA(0,1,1) is the best model of the ones I tested since it has the lowest AICc value of 659.18. The plot from the auto ARIMA still has a lower AICc value though at 657.1.
e. produce forecasts of your fitted model. Do the forecasts look reasonable?
fit_us_011 |>forecast(h =10) |>autoplot(us)
The forecasts look reasonable, it follows the upward trend of the previous data at a similar pace.
f. compare the results with what you would obtain using ETS() (with no transformation).
# A tibble: 3 × 11
.model Country .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
<chr> <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 arima_011 United… Test -2.78e11 3.08e11 2.78e11 -1.51 1.51 0.892 0.821 0.419
2 arima_au… United… Test -2.18e11 2.50e11 2.18e11 -1.18 1.18 0.699 0.666 0.400
3 ets United… Test 2.91e11 3.45e11 2.91e11 1.57 1.57 0.935 0.920 0.191
Testing the accuracy of the ETS(), ARIMA(0,1,1), and auto ARIMA() models, it appears that the auto ARIMA() model was the best since it’s RMSSE value was the lowest. The ARIMA(0,1,1) came in second as a middleground with a RMSSE value between that of the auto ARIMA() and ETS() models.
