Data 624: HW 5
Nakesha Fray
2025-03-04
Do exercises 8.1, 8.5, 8.6, 8.7, 8.8, 8.9 in Hyndman. Please submit both the link to your Rpubs and the .pdf file with your run code
8.1 Consider the the number of pigs slaughtered in Victoria, available in the aus_livestock dataset.
data("aus_livestock")
aus_pig <- aus_livestock %>%
filter(Animal == "Pigs",
State == "Victoria")a. Use the ETS() function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of alpha and Initial states: l[0], and generate forecasts for the next four months.
I first filtered the aus_livestock dataset for Victorian pigs. Then used the ETS() function and found the optimal value of alpha to be 0.3221247 and Initial states: l[0] to be 100646.6.
fit <- aus_pig %>%
model(ETS(Count ~ error("A") + trend("N") + season("N")))
report(fit)## Series: Count
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.3221247
##
## Initial states:
## l[0]
## 100646.6
##
## sigma^2: 87480760
##
## AIC AICc BIC
## 13737.10 13737.14 13750.07fc <- fit |> forecast(h = 4)
fc |>
autoplot(aus_pig) +
labs(y="Number of pigs slaughtered", title="Pigs Slaughtered in Victoria, Australia") 
b. Compute a 95% prediction interval for the first forecast using ^y±1.96s where s is the standard deviation of the residuals. Compare your interval with the interval produced by R.
I found the standard deviation to be 9344.66579 and the 95% prediction interval for the first forecast to be (76871.0125, 113502.1024). However, the 95% prediction interval produced by R was (76854.79, 113518.3) for the first forecast. The lower bound that I computed is slightly larger than the lower bound produced by R and the upper bound I computed is slightly smaller than the upper bound produced by R.
#Computed prediction interval
aug_fit <- augment(fit)
s <- sd(aug_fit$.resid)
s## [1] 9344.666lower_bound <- fc$.mean[1] - 1.96 * s
upper_bound <- fc$.mean[1] + 1.96 * s
lower_bound## [1] 76871.01upper_bound## [1] 113502.1#Prediction interval produced by R
fc_pi <- fc %>%
hilo(level = 95)%>%
mutate(lower = `95%`$lower, upper = `95%`$upper)
fc_pi$lower[1]## [1] 76854.79fc_pi$upper[1]## [1] 113518.38.5 Data set global_economy contains the annual Exports from many countries. Select one country to analyze.
uk_economy <- global_economy %>%
filter(Country == "United Kingdom")a. Plot the Exports series and discuss the main features of the data.
The overall trend of the series of UK exports from global_economy is upward with some periods of drastic increasing and decreasing trends throughout. It seems like one of the drastic increases was around the early 1970s and the drastic decrease was during the late 1980s. We cannot determine if there is seasonality in this time series because it is yearly. I also do not see any cyclic behavior in this time series.
uk_economy %>%
autoplot(Exports) +
labs(title = "UK Exports")
b. Use an ETS(A,N,N) model to forecast the series, and plot the forecasts.
fit1 <- uk_economy %>%
model(ETS(Exports ~ error("A") + trend("N") + season("N")))
report(fit1)## Series: Exports
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.9998999
##
## Initial states:
## l[0]
## 20.21233
##
## sigma^2: 1.9003
##
## AIC AICc BIC
## 276.7066 277.1511 282.8880fc1 <- fit1 |> forecast(h = 10)
fc1 |>
autoplot(uk_economy, level=NULL) +
labs(y="Exports of goods and services (% of GDP).", title="Exports of goods and services (% of GDP) in the UK")
fc1 |>
autoplot(uk_economy) +
labs(y="Exports of goods and services (% of GDP).", title="Exports of goods and services (% of GDP) in the UK")
c. Compute the RMSE values for the training data.
The RMSE value for the training data is 1.35.
fit1 |> accuracy()## # A tibble: 1 × 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United Kingdom "ETS(Ex… Trai… 0.178 1.35 1.02 0.564 4.09 0.983 0.991 0.0135d. Compare the results to those from an ETS(A,A,N) model. (Remember that the trended model is using one more parameter than the simpler model.) Discuss the merits of the two forecasting methods for this data set.
When we use the ETS(A,A,N) model, we are accounting for an additive trend in the the time series/dataset compared to an ETS(A,N,N) model which assumes no trend and takes on a constant level. Both of these models do not account for seasonality, which is correct in this case since we cannot look at seasonality in this dataset. However, the overall trend in the UK exports is upward trending with increasing and decreasing trends in the data which should be better for an ETS(A,A,N) model since it assumes trends. We see from the RMSE values that these models predicted the data quite similarly with only a difference of 0.004, with the ETS(A,A,N) model only being slightly better in predicting the values for this dataset.
uk_fit <- uk_economy %>%
model(
ses_ann = ETS(Exports ~ error("A") + trend("N") + season("N")),
holt_aan = ETS(Exports ~ error("A") + trend("A") + season("N")))
tidy(uk_fit)## # A tibble: 6 × 4
## Country .model term estimate
## <fct> <chr> <chr> <dbl>
## 1 United Kingdom ses_ann alpha 1.00
## 2 United Kingdom ses_ann l[0] 20.2
## 3 United Kingdom holt_aan alpha 0.956
## 4 United Kingdom holt_aan beta 0.000100
## 5 United Kingdom holt_aan l[0] 20.9
## 6 United Kingdom holt_aan b[0] 0.116accuracy(uk_fit)## # A tibble: 2 × 11
## Country .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <fct> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 United Kingd… ses_a… Trai… 0.178 1.35 1.02 0.564 4.09 0.983 0.991 0.0135
## 2 United Kingd… holt_… Trai… 0.0494 1.35 1.06 0.00945 4.25 1.01 0.989 0.0559e. Compare the forecasts from both methods. Which do you think is best?
After plotting both the ses and holt’s model, I still think that the holt’s model or the ETS(A,A,N) model is a better model because it accounts for the upward trend in the UK exports dataset. We see that the ETS(A,A,N) model forcast continues at an upward trend considering that the last points of the data was increasing. The ses model or ETS(A,N,N) model did not account for the overall increasing trend in the UK exports dataset. Therefore, we instead see a constant level of predictions after the last few points into the next 10 years.
uk_fc <- uk_fit |>
forecast(h = 15)
uk_fc |>
autoplot(uk_economy, level = NULL) +
labs(title = "Exports of goods and services (% of GDP) in the UK",
y = "Exports of goods and services (% of GDP)") +
guides(colour = guide_legend(title = "Forecast"))
f. Calculate a 95% prediction interval for the first forecast for each model, using the RMSE values and assuming normal errors. Compare your intervals with those produced using R.
The 95% prediction interval for the ses EST(A,A,N) model to be (27.87352, 33.18329) and for the holt’s EST(A,A,N) model to be (27.88065, 33.17616) using the RMSE value instead of the standard deviation. These values are very close together, the ses EST(A,A,N) model has a slightly larger prediction interval margin than holt’s EST(A,A,N) model.
uk_rmse <- accuracy(uk_fit)
uk_ses <- uk_economy %>%
model(ses_ann = ETS(Exports ~ error("A") + trend("N") + season("N")))
uk_holt <- uk_economy %>%
model(holt_aan = ETS(Exports ~ error("A") + trend("A") + season("N")))
ses_fc <- uk_ses |>
forecast(h = 15)
holt_fc <- uk_ses |>
forecast(h = 15)
lower_bound_ses <- ses_fc$.mean[1] - 1.96 * uk_rmse$RMSE[1]
upper_bound_ses <- ses_fc$.mean[1] + 1.96 * uk_rmse$RMSE[1]
lower_bound_ses## [1] 27.87352upper_bound_ses## [1] 33.18329lower_bound_holt <- holt_fc$.mean[1] - 1.96 * uk_rmse$RMSE[2]
upper_bound_holt <- holt_fc$.mean[1] + 1.96 * uk_rmse$RMSE[2]
lower_bound_holt ## [1] 27.88065upper_bound_holt ## [1] 33.176168.6 Forecast the Chinese GDP from the global_economy data set using an ETS model. Experiment with the various options in the ETS() function to see how much the forecasts change with damped trend, or with a Box-Cox transformation. Try to develop an intuition of what each is doing to the forecasts.
[Hint: use a relatively large value of h when forecasting, so you can clearly see the differences between the various options when plotting the forecasts.]
I first plotted China’s GDP and found that there was an overall upward trend, no cyclic behavior, and no way to determine seasonality with the yearly data. Knowing these main features of the dataset, I assumed that either holt’s ETS(A,A,N) or a damped ETS(A,Ad,N) model would perform better than ses. This is because holt’s ETS(A,A,N) and damped ETS(A,Ad,N) models act similarly - they both account for trend but not for seasonality. An ses ETS(A,N,N) model on the other hand, assumes no trend and takes on a constant level. The difference between holt’s ETS(A,A,N) and damped ETS(A,Ad,N) model is that we can expect the damped model to gradually slow down over time to account for flexibility in the way real data may evolve.
The forecast shows exactly what was expected: the ses ETS(A,N,N) model provided a constant rate which does not reflect what we expected the GDP to forecast over the next 40 years, given the increasing upward trend it was moving in. The holt’s ETS(A,A,N) model continued increasing upward with the assumption that the trend would increase at the last data point’s pace over the next 40 years. For the damped ETS(A,Ad,N) model, I created 2 and changed the parameter to see if there was difference. The non-specified damped ETS(A,Ad,N) model while still increasing over time, but eventually begins to slow down and drift away from a very linear upward trend, while the specified damped ETS(A,Ad,N) model of phi=0.9 starts to level off and move towards a constant rate much faster.
I lastly added a box-cox transformation to the ETS() models and found that the ses ETS(A,N,N) model was pretty much the same with a constant level and the holt’s ETS(A,A,N) model was still upward trending but instead appeared more exponential. For the damped ETS(A,Ad,N) model, I again created 2 and changed the parameter to see if there was difference. It was pretty much the same but the specified parameter of .90 for phi was a little lower.
china_economy <- global_economy %>%
filter(Country == "China")
china_economy %>%
autoplot(GDP) +
labs(title = "China's GDP")
fit2 <- china_economy %>%
model(
ses = ETS(GDP ~ error("A") + trend("N") + season("N")),
holt = ETS(GDP ~ error("A") + trend("A") + season("N")),
damped = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
damp_90 = ETS(GDP ~ error("A") + trend("Ad", phi = 0.9) + season("N"))
)
report(fit2)## Warning in report.mdl_df(fit2): 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: 4 × 10
## Country .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 China ses 1.79e23 -1669. 3345. 3345. 3351. 1.73e23 7.22e23 2.13e11
## 2 China holt 3.88e22 -1624. 3258. 3259. 3268. 3.61e22 1.31e23 9.59e10
## 3 China damped 3.96e22 -1624. 3260. 3262. 3273. 3.62e22 1.30e23 9.49e10
## 4 China damp_90 4.17e22 -1626. 3261. 3262. 3271. 3.81e22 1.22e23 9.37e10fc2 <- fit2 |> forecast(h = 40)
fc2 |>
autoplot(china_economy, level = NULL) +
labs(y="GDP ", title="GDP in the China (SES, Holt's, Damped & Damped 90)")
#Box-cox
lambda <- china_economy |>
features(GDP, features = guerrero) |>
pull(lambda_guerrero)
lambda## [1] -0.03446284fit_box <- china_economy %>%
model(
ses = ETS(box_cox(GDP, lambda) ~ error("A") + trend("N") + season("N")),
holt = ETS(box_cox(GDP, lambda) ~ error("A") + trend("A") + season("N")),
damped = ETS(box_cox(GDP, lambda) ~ error("A") + trend("Ad") + season("N")),
damp_90 = ETS(box_cox(GDP, lambda) ~ error("A") + trend("Ad", phi = 0.9) + season("N"))
)
fc_box <- fit_box |> forecast(h = 20)
fc_box |>
autoplot(china_economy, level = NULL) +
labs(y="GDP ", title="GDP in the China (Box-Cox of SES, Holt's, Damped, Damped 90)")
8.7 Find an ETS model for the Gas data from aus_production and forecast the next few years. Why is multiplicative seasonality necessary here? Experiment with making the trend damped. Does it improve the forecasts?
The ETS model recommended is the ETS(M,A,M) which is error(“M”) + trend(“A”) + season(“M”) - multiplicative error, additive trend, and multiplicative seasonality. To experiment with the trend damped, I included a multiplicative ETS(M,A,M) model, a multiplicative ETS(M,Ad,M) model, a multiplicative ETS(M,Ad,M) model with phi=0.9, and a multiplicative ETS(M,Ad,M) model with phi = 0.85. It is hard to tell from the plot which forecast is better, so I decided to look at the RMSE value and found that the multiplicative ETS(M,Ad,M) model with phi = 0.85 performed the best since it has the smallest RMSE value, which means it probably captured the trend and seasonality the best in the data. Multiplicative seasonality is necessary here because the Gas production from aus_production has a seasonality component which actually increases in peaks sizes as the Gas production increases. The multiplicative aspect accounts for this sort of direct proportional relationship between peak size and gas production, compared to the additive.
aus_production %>%
autoplot(Gas)
aus_fit <- aus_production %>%
model(fit = ETS(Gas))
report(aus_fit)## Series: Gas
## Model: ETS(M,A,M)
## Smoothing parameters:
## alpha = 0.6528545
## beta = 0.1441675
## gamma = 0.09784922
##
## Initial states:
## l[0] b[0] s[0] s[-1] s[-2] s[-3]
## 5.945592 0.07062881 0.9309236 1.177883 1.074851 0.8163427
##
## sigma^2: 0.0032
##
## AIC AICc BIC
## 1680.929 1681.794 1711.389aus_fit %>%
forecast(h = '5 years') %>%
autoplot(aus_production)
aus_fit2 <- aus_production |>
model(multiplicative = ETS(Gas ~ error("M") + trend("A") + season("M")),
damp = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
damp_90 = ETS(Gas ~ error("M") + trend("Ad", phi = 0.9) + season("M")),
damp_85 = ETS(Gas ~ error("M") + trend("Ad", phi = 0.85) + season("M")))
accuracy(aus_fit2)## # A tibble: 4 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 multiplicative Training -0.115 4.60 3.02 0.199 4.08 0.542 0.606 -0.0131
## 2 damp Training -0.00439 4.59 3.03 0.326 4.10 0.544 0.606 -0.0217
## 3 damp_90 Training 0.255 4.58 3.04 0.655 4.15 0.545 0.604 -0.00451
## 4 damp_85 Training 0.353 4.57 3.04 0.799 4.18 0.546 0.603 -0.0283aus_fc <- aus_fit2 |> forecast(h = '3 years')
aus_fc |>
autoplot(aus_production, level = NULL) +
labs(y="Gas ", title="Gas Production in Australia (Mult, Damped, Damped 90, Damped 85)")
8.8 Recall your retail time series data (from Exercise 7 in Section 2.10).
a. Why is multiplicative seasonality necessary for this series?
Multiplicative seasonality is necessary here because similar to the last question, the Turnovers from myseries has a seasonality component and variation which actually increases in peaks sizes as the Turover increases. The multiplicative aspect accounts for this sort of direct proportional relationship between peak size/variation and Turnovers, compared to the additive.
set.seed(3458)
myseries <- aus_retail |>
filter(`Series ID` == sample(aus_retail$`Series ID`,1))
myseries_train <- myseries |>
filter(year(Month) < 2011)
myseries_train %>%
autoplot(Turnover)
myseries_train %>%
gg_season(Turnover)
b. Apply Holt-Winters’ multiplicative method to the data. Experiment with making the trend damped.
my_fit <- myseries_train |>
model(multiplicative = ETS(Turnover ~ error("M") + trend("A") + season("M")),
damp = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
my_fc <- my_fit |> forecast(h = "5 years")
my_fc |>
autoplot(myseries_train, level = NULL) +
labs(title="Australian Trade Turnovers (Mult, Damped)",
y="Trade Turnover ($AUD Millions)") 
c. Compare the RMSE of the one-step forecasts from the two methods. Which do you prefer?
Comparing only the multiplicative and multiplicative damped method, I prefer the damped since it has the smaller RMSE value, which means it may have best captured the trend and seasonality better. We can also see all the other damped models with different phi values is still lower than the multiplicative.
accuracy(my_fit)## # A tibble: 2 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Quee… Food re… multi… Trai… 0.244 17.7 12.7 0.0172 1.94 0.227 0.249 0.0521
## 2 Quee… Food re… damp Trai… 1.73 17.4 12.6 0.208 1.93 0.224 0.244 -0.0105d. Check that the residuals from the best method look like white noise.
For myseries, I used the damped ETS(M,Ad,M) model. The histogram of residuals look somewhat unimodal and normally distributed. The ACF of the residuals shows that most of the lags are outside the limit which may indicate autocorrelation. The residual does not look like it is white noise since there does appear to be correlation in the residuals series and there is a pattern/trend.
From the Box-Pierce test and the Ljung-Box test, we see that the results are significant since the p-values are both 0. We can conclude that the residuals are distinguishable from white noise. This means that there is autocorrelation in the residuals and the residuals are not white noise.
damp_myseries <- myseries_train |>
model(damp = ETS(Turnover ~ error("M") + trend("Ad") + season("M")))
damp_myseries |> gg_tsresiduals()
damp_myseries |> forecast() |> autoplot(myseries_train)
augment(damp_myseries) |> features(.innov, box_pierce, lag = 10)## # A tibble: 1 × 5
## State Industry .model bp_stat bp_pvalue
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Queensland Food retailing damp 184. 0augment(damp_myseries) |> features(.innov, ljung_box, lag=10)## # A tibble: 1 × 5
## State Industry .model lb_stat lb_pvalue
## <chr> <chr> <chr> <dbl> <dbl>
## 1 Queensland Food retailing damp 188. 0e. Now find the test set RMSE, while training the model to the end of 2010. Can you beat the seasonal naïve approach from Exercise 7 in Section 5.11?
I beat the seasonal naïve approach from Exercise 7 in Section 5.11. The damped method produced a better forecast with a smaller RMSE values for both training and test compared to the seasonal naïve approach.
fit_series <- myseries_train %>%
model(damp = ETS(Turnover ~ error("M") + trend("Ad") + season("M")),
snaive = SNAIVE(Turnover))
fc_series <- fit_series |>
forecast(new_data = anti_join(myseries, myseries_train))## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fc_series |> autoplot(myseries, level = NULL)
fit_series |> accuracy()## # A tibble: 2 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Queen… Food re… damp Trai… 1.73 17.4 12.6 0.208 1.93 0.224 0.244 -0.0105
## 2 Queen… Food re… snaive Trai… 55.7 71.1 56.2 7.88 7.94 1 1 0.851fc_series |> accuracy(myseries)## # A tibble: 2 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 damp Queensl… Food re… Test 88.4 135. 98.2 3.95 4.50 1.75 1.90 0.858
## 2 snaive Queensl… Food re… Test 326. 362. 326. 15.5 15.5 5.80 5.08 0.9168.9 For the same retail data, try an STL decomposition applied to the Box-Cox transformed series, followed by ETS on the seasonally adjusted data. How does that compare with your best previous forecasts on the test set?
I applied an STL decomposition to the Box-Cox transformed myseries data from us_retail, then used the seasonally adjusted data to replace the turnovers in myseries_train. This was my work around for getting the R to rerun another model on the data and produce results. I found lambda to be 0.09024391 but found the RSME value to be much larger for the test set but lower for the training set. For this reason, I would stick with the damped multiplicative model from the previous forecast.
lambda1 <- myseries_train |>
features(Turnover, features = guerrero) |>
pull(lambda_guerrero)
lambda1## [1] 0.09024391myseries_stl <- myseries_train |>
model(stl = STL((box_cox(Turnover, lambda1)))) %>%
components()
myseries_stl |>
autoplot()
myseries_train$Turnover <- myseries_stl$season_adjust
fit3 <- myseries_train %>%
model(mam = ETS(Turnover ~ error("M") + trend("A") + season("M")))
fc3_series <- fit3 |>
forecast(new_data = anti_join(myseries, myseries_train))## Joining with `by = join_by(State, Industry, `Series ID`, Month, Turnover)`fit3 |> accuracy()## # A tibble: 1 × 12
## State Industry .model .type ME RMSE MAE MPE MAPE MASE RMSSE
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Queen… Food re… mam Trai… -4.92e-4 0.0389 0.0319 -0.00450 0.382 0.216 0.238
## # ℹ 1 more variable: ACF1 <dbl>fc3_series |> accuracy(myseries)## # A tibble: 1 × 12
## .model State Industry .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 mam Queensl… Food re… Test 995. 1205. 995. 97.8 97.8 NaN NaN 0.984