pigs <- aus_livestock |> filter(Animal == 'Pigs', State == 'Victoria')
autoplot(pigs, Count)
data_624_hw05
8.1
a.
A plot of the pigs in the state of Victoria, Australia.
The code below fits a simple exponential smoothing model to the pig data. Forecasts 4 years, and extracts the alpha value and initial level for display.
pig_fit <- pigs |>
model(ses = ETS(Count ~ error('A') + trend('N') + season('N')))
pig_forecast <- pig_fit |>
forecast(h = 4)
alpha <- pig_fit |>
tidy() |>
filter(term == 'alpha') |>
pull(estimate)
initial_level <- pig_fit |>
tidy() |>
filter(term == 'l[0]') |>
pull(estimate)
cat('The optimal alpha value was calculated as:',alpha, '\n\n')The optimal alpha value was calculated as: 0.3221247 cat('The optimal inital level was calculates as:', initial_level, '\n\n')The optimal inital level was calculates as: 100646.6 report(pig_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.07 A plot of the pig data and the corresponding forecast.
pig_forecast |>
autoplot(pigs) +
geom_line(aes(y = .fitted), col = 'orange',
data = augment(pig_fit)) +
guides(colour = 'none')
b.
The prediction interval is calculated in the code below.
first <- pig_forecast |>
as_tibble() |>
slice(1) |>
pull(.mean)
s <- pig_fit |>
augment() |>
pull(.resid) |>
sd(na.rm = TRUE)
lower <- first - 1.96 * s
upper <- first + 1.96 * sThe following code uses r to calculate the prediction interval and then compares it to the manually calculated one from the code just above.
interval <- pig_forecast |>
hilo(level = 95) |>
as_tibble() |>
mutate(
lower95 = `95%`$lower,
upper95 = `95%`$upper) |>
slice(1)
lower_r <- interval$lower95
upper_r <- interval$upper95
cat('The manually calculated lower bound of the 95% interval:', lower, '\n')The manually calculated lower bound of the 95% interval: 76871.01 cat('The model calculated lower bound of the 95% interval:', lower_r, '\n')The model calculated lower bound of the 95% interval: 76854.79 cat('A difference of', lower - lower_r, '\n\n')A difference of 16.22359 cat('The manually calculated upper bound of the 95% interval:', upper, '\n')The manually calculated upper bound of the 95% interval: 113502.1 cat('The model calculated upper bound of the 95% interval:', upper_r, '\n')The model calculated upper bound of the 95% interval: 113518.3 cat('A difference of', upper_r - upper, '\n\n')A difference of 16.22359 As shown above the manually calculated and model calculated are close but differ by a value of 16.22359 for each bound.
8.5
a.
A plot of Mexico’s exports as a percentage of GDP over time.
mexports <- global_economy |>
filter(Country == 'Mexico') |> select(Country, Year, Exports)
mexports |> autoplot(Exports)
The plot of Mexico’s exports by year since 1960 shows an upward trend over time starting about 1975. There is ongoing variation throughout this trend with a notable decline and subsequent recovery in the mid 1990s.
b.
The code fits a simple exponential smoothing model to the Mexican exports data and forecasts the following ten years.
mexports_fit <- mexports |>
model(ses = ETS(Exports ~ error('A') + trend('N') + season('N')))
mexports_forecast <- mexports_fit |>
forecast(h = 10)
mexports_forecast |> autoplot(mexports) +
geom_line(aes(y = .fitted), col = 'orange',
data = augment(mexports_fit)) +
guides(colour = 'none')
c.
The RMSE for the forecast is calculated in the following code and is then compared with the r calculated RMSE, which is calculated in the code block following the one below.
rmse_simple <- mexports_fit |> augment() |>
as_tibble() |>
group_by(Country, .model) |>
summarise(
rmse = sqrt(mean(.resid^2)),
.groups ='drop'
) |> pull(rmse)
cat(rmse_simple)2.154425The RMSE calculated above matches the one returned by the accuracy() function below.
mexports_fit |> 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 Mexico ses Training 0.506 2.15 1.38 1.83 7.78 0.983 0.991 0.203d.
The Mexican export data is then fit to a model with a trend component.
mexports_fit_trend <- mexports |>
model(trend = ETS(Exports ~ error('A') + trend('A') + season('N')))
mexports_forecast_trend <- mexports_fit_trend |>
forecast(h = 10)
mexports_forecast_trend |> autoplot(mexports) +
geom_line(aes(y = .fitted), col = 'orange',
data = augment(mexports_fit_trend)) +
guides(colour = 'none')
The model, with the additional parameter for trend, has a slightly improved RMSE over the model without trend, from a 2.154425 to 2.093999.
rmse_trend <- mexports_fit_trend |> augment() |>
as_tibble() |>
group_by(Country, .model) |>
summarise(
rmse = sqrt(mean(.resid^2)),
.groups ='drop'
) |>
pull(rmse)
cat(rmse_trend)2.093999mexports_fit_trend |> 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 Mexico trend Training -0.00816 2.09 1.41 -1.97 8.68 1.01 0.964 0.203Ultimately, when considering the data for the exports, there is a clear trend upwards. With this in mind, utilizing the model that includes a trend component will most likely be beneficial when forecasting this data.
e.
The code below builds a table that calculates the difference between the two models forecasts.
mexports_forecast |>
as_tibble() |>
left_join(
mexports_forecast_trend |>
as_tibble() |>
transmute(Country, Year, mean_trend = .mean),
by = c('Country', 'Year')
) |>
mutate(
trend_mean_difference = mean_trend - .mean
)# A tibble: 10 × 7
Country .model Year
<fct> <chr> <dbl>
1 Mexico ses 2018
2 Mexico ses 2019
3 Mexico ses 2020
4 Mexico ses 2021
5 Mexico ses 2022
6 Mexico ses 2023
7 Mexico ses 2024
8 Mexico ses 2025
9 Mexico ses 2026
10 Mexico ses 2027
# ℹ 4 more variables: Exports <dist>, .mean <dbl>, mean_trend <dbl>,
# trend_mean_difference <dbl>Inspecting the forecasts, the SES forecast is similar to a Naive forecast with each forecast equal to the last estimated level, whereas the trend forecast increases over time. Considering the noticeable trend in the observed data, the trended forecast is most likely the one that would be best to use in this instance.
f.
The code below manually calculates the prediction interval for the simple exponential smoothing model (without trend).
mean_mex_simple <- mexports_forecast |>
as_tibble() |>
slice(1) |>
pull(.mean)
mex_simple_lower <- mean_mex_simple - 1.96 * rmse_simple
mex_simple_upper <- mean_mex_simple + 1.96 * rmse_simple
cat(mex_simple_lower, '\n')33.64367 cat(mex_simple_upper)42.08901Then the prediction interval is manually calculated for the model with a trend component.
mean_mex_trend <- mexports_forecast_trend |>
as_tibble() |>
slice(1) |>
pull(.mean)
mex_trend_lower <- mean_mex_trend - 1.96 * rmse_trend
mex_trend_upper <- mean_mex_trend + 1.96 * rmse_trend
cat(mex_trend_lower, '\n')34.278 cat(mex_trend_upper)42.48647The prediction interval for the simple model is calculated using the R functions below.
mexport_interval <- mexports_forecast |>
hilo(level = 95) |>
as_tibble() |>
mutate(
lower95 = `95%`$lower,
upper95 = `95%`$upper) |>
slice(1)
mexports_lower_r <- mexport_interval$lower95
mexports_upper_r <- mexport_interval$upper95
cat(mexports_lower_r, '\n')33.569 cat(mexports_upper_r)42.16368The prediction interval for the trended model is calculated using the R functions below.
mexport_trend_interval <- mexports_forecast_trend |>
hilo(level = 95) |>
as_tibble() |>
mutate(
lower95 = `95%`$lower,
upper95 = `95%`$upper) |>
slice(1)
mexports_trend_lower_r <- mexport_trend_interval$lower95
mexports_trend_upper_r <- mexport_trend_interval$upper95
cat(mexports_trend_lower_r, '\n')34.12878 cat(mexports_trend_upper_r)42.63569The table below combines all the prediction intervals for comparison.
interval_table <- tibble(
model = c('Simple exponential smoothing', 'Trend model'),
manual_lower = c(mex_simple_lower, mex_trend_lower),
manual_upper = c(mex_simple_upper, mex_trend_upper),
r_lower = c(mexports_lower_r, mexports_trend_lower_r),
r_upper = c(mexports_upper_r, mexports_trend_upper_r)
)
interval_table# A tibble: 2 × 5
model manual_lower manual_upper r_lower r_upper
<chr> <dbl> <dbl> <dbl> <dbl>
1 Simple exponential smoothing 33.6 42.1 33.6 42.2
2 Trend model 34.3 42.5 34.1 42.6Again, as in the previous problem, the ranges of the prediction intervals are very similar for the manually calculated and r-calculated ranges, however they are not identical.
8.6
Below is a plot for the GDP of China’s economy over time from 1960 to 2018.
china <- global_economy |>
filter(Country == 'China') |>
select(Country, Year, GDP)autoplot(china, GDP)
China’s GDP shows a strong growth trend.
Below, China’s GDP data is fit to four different models and plotted.
lambda <- china |>
features(GDP, guerrero) |>
pull(lambda_guerrero)
cat('lambda:', lambda)lambda: -0.03446284china_fit <- china |>
model(
ses = ETS(GDP ~ error('A') + trend('N') + season('N')),
trend = ETS(GDP ~ error('A') + trend('A') + season('N')),
damped = ETS(GDP ~ error('A') + trend('Ad') + season('N')),
box_cox = ETS(box_cox(GDP, lambda))
)
china_forecast <- china_fit |>
forecast(h = 20)
china_forecast |>
autoplot(china)
Below is table of the different models’ forecasts.
china_forecast |>
as_tibble() |>
select(.model, Year, .mean) |>
pivot_wider(
names_from = .model,
values_from = .mean
) |> head()# A tibble: 6 × 5
Year ses trend damped box_cox
<dbl> <dbl> <dbl> <dbl> <dbl>
1 2018 1.22e13 1.30e13 1.30e13 1.39e13
2 2019 1.22e13 1.38e13 1.37e13 1.58e13
3 2020 1.22e13 1.45e13 1.44e13 1.81e13
4 2021 1.22e13 1.53e13 1.51e13 2.06e13
5 2022 1.22e13 1.60e13 1.58e13 2.36e13
6 2023 1.22e13 1.68e13 1.65e13 2.71e13The models show the different behaviors of their forecasts. Simple exponential smoothing behaves as Naive forecasting, continuing to forecast the final estimated value, while a model including a trend component increases progressively over time. The damped trend also increases but the rate of increases is damped over time. As seen above the trend forecasts increase at a faster rate than the damped forecasts. The Box-Cox transformation forecast increases the forecast for this model at the fastest rate. The Box-Cox model changes the scale on which ETS models the growth and variation and then the transformed forecasts are returned to the GDP scale. The strong increasing trend of the observed China GDP data, supports the idea of the trend models all having continued growth rates.
8.7
Gas production data for Australia is plotted below.
gas <- aus_production |>
select(Quarter, Gas)
gas |> autoplot(Gas)
Looking at the plot of the Australian gas production shows quarterly gas production with a strong upward trend and a seasonal variation component.
The data is fit to seven different models and then each is plotted below.
lambda <- gas |>
features(Gas, guerrero) |>
pull(lambda_guerrero)
gas_fit <- gas |>
model(
additive = ETS(Gas ~ error('A') + trend('A') + season('A')),
multiplicative = ETS(Gas ~ error('M') + trend('M') + season('M')),
additive_damped = ETS(Gas ~ error('A') + trend('Ad') + season('A')),
multiplicative_damped = ETS(Gas ~ error('M') + trend('Md') + season('M')),
multiplicative_additive_damped = ETS(Gas ~ error('M') + trend('Ad') + season('M')),
box_cox = ETS(box_cox(Gas, lambda)),
automatic = ETS((Gas))
)
gas_fit |>
select(additive) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(additive) |>
report()Series: Gas
Model: ETS(A,A,A)
Smoothing parameters:
alpha = 0.5209159
beta = 0.0001000024
gamma = 0.4180453
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
10.06526 0.9870514 -5.890283 15.61932 -2.811314 -6.917728
sigma^2: 23.5601
AIC AICc BIC
1872.452 1873.318 1902.913 gas_fit |>
select(multiplicative) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(multiplicative) |>
report()Series: Gas
Model: ETS(M,M,M)
Smoothing parameters:
alpha = 0.6947632
beta = 0.1182963
gamma = 0.08784581
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
5.913798 1.010604 0.9320972 1.180305 1.070716 0.8168816
sigma^2: 0.0032
AIC AICc BIC
1680.654 1681.519 1711.115 gas_fit |>
select(additive_damped) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(additive_damped) |>
report()Series: Gas
Model: ETS(A,Ad,A)
Smoothing parameters:
alpha = 0.5759917
beta = 0.04513765
gamma = 0.3538354
phi = 0.9558689
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
16.69437 0.2260366 -0.7505321 2.699016 2.663774 -4.612259
sigma^2: 21.8606
AIC AICc BIC
1857.090 1858.153 1890.935 gas_fit |>
select(multiplicative_damped) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(multiplicative_damped) |>
report()Series: Gas
Model: ETS(M,Md,M)
Smoothing parameters:
alpha = 0.6418229
beta = 0.1296012
gamma = 0.0939101
phi = 0.9615711
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
5.853198 1.013851 0.9292935 1.183223 1.069365 0.818119
sigma^2: 0.0032
AIC AICc BIC
1679.596 1680.658 1713.440 gas_fit |>
select(multiplicative_additive_damped) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(multiplicative_damped) |>
report()Series: Gas
Model: ETS(M,Md,M)
Smoothing parameters:
alpha = 0.6418229
beta = 0.1296012
gamma = 0.0939101
phi = 0.9615711
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
5.853198 1.013851 0.9292935 1.183223 1.069365 0.818119
sigma^2: 0.0032
AIC AICc BIC
1679.596 1680.658 1713.440 gas_fit |>
select(box_cox) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(box_cox) |>
report()Series: Gas
Model: ETS(A,A,A)
Transformation: box_cox(Gas, lambda)
Smoothing parameters:
alpha = 0.81796
beta = 0.1216283
gamma = 0.0001000577
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3]
1.888944 0.03198818 -0.1065298 0.2508721 0.1003452 -0.2446875
sigma^2: 0.0059
AIC AICc BIC
63.28945 64.15484 93.74991 gas_fit |>
select(automatic) |>
forecast(h = 40) |>
autoplot(gas, level = NULL)
gas_fit |>
select(automatic) |>
report()Series: Gas
Model: ETS(M,A,M)
Transformation: (Gas)
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.389 This table shows the report data for each of the seven models plotted above.
glance(gas_fit) |>
arrange(AICc)# A tibble: 7 × 9
.model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 box_cox 5.86e-3 -22.6 63.3 64.2 93.7 5.65e-3 0.0102 0.0570
2 multiplicative_da… 3.21e-3 -830. 1680. 1681. 1713. 2.08e+1 31.5 0.0410
3 multiplicative 3.21e-3 -831. 1681. 1682. 1711. 2.20e+1 33.9 0.0413
4 automatic 3.24e-3 -831. 1681. 1682. 1711. 2.11e+1 32.2 0.0413
5 multiplicative_ad… 3.29e-3 -832. 1684. 1685. 1718. 2.11e+1 32.0 0.0417
6 additive_damped 2.19e+1 -919. 1857. 1858. 1891. 2.10e+1 29.6 3.16
7 additive 2.36e+1 -927. 1872. 1873. 1903. 2.27e+1 29.7 3.35 Inspecting the individual models and their corresponding AICc scores, it appears that a multiplicative-damped model ETS(M, Md, M) provides the preferred model, followed closely by the fully multiplicative, and then an ETS(M, A, M) model. The small difference provides little evidence to prefer one over the other. The Box-Cox model cannot be directly compared as it underwent transformation and is fitted on a transformed scale.
The seasonality component is necessarily multiplicative because the size of the seasonal fluctuations increases as the level of gas production rises. It models the seasonal effects proportionally to the level. Additive seasonality assumes fluctuations are roughly constant in size.
The damped trend model improves the additive and multiplicative model’s AICc score. While there is a relatively small difference in AICc values, results do support damping the multiplicative trend.
8.8
The retail data is plotted below.
# set the seed
set.seed(47)
# select the retail data
myseries <- aus_retail |> filter(`Series ID` == sample(aus_retail$`Series ID`, 1))
# plot the training data and the entire data
autoplot(myseries, Turnover)
# create a training dataset of observations before 2010
myseries_train <- myseries |> filter(year(Month) <= 2010)
# plot the training data and the entire data
autoplot(myseries, Turnover) +
autolayer(myseries_train, Turnover, colour = "red")
a.
Inspecting the plot of the retail data series, multiplicative seasonality is necessary because the magnitude of the seasonality increases with the level of the series. A multiplicative seasonal component allows the forecast to adjust proportional to the level of the data series.
b.
Two models are fit to the retail data below and then both are plotted.
lambda <- myseries |>
features(Turnover, guerrero) |>
pull(lambda_guerrero)
myseries_fit <- myseries |>
model(
ets_mult = ETS(Turnover ~ error('M') + trend('A') + season('M')),
ets_damped = ETS(Turnover ~ error('M') + trend('Ad') + season('M'))
)
myseries_fit |>
select(State, Industry, ets_mult) |>
forecast(h = 36) |>
autoplot(myseries, level = NULL)
myseries_fit |>
select(ets_mult) |>
report()Series: Turnover
Model: ETS(M,A,M)
Smoothing parameters:
alpha = 0.3891239
beta = 0.0001001378
gamma = 0.1833412
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
45.40407 0.7237567 0.8641859 0.7649106 0.9353169 1.450089 1.047977 1.036834
s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
1.027067 0.9674972 0.9716836 0.9911705 1.011353 0.9319146
sigma^2: 0.0025
AIC AICc BIC
4597.666 4599.112 4667.179 myseries_fit |>
select(State, Industry, ets_damped) |>
forecast(h = 36) |>
autoplot(myseries, level = NULL)
myseries_fit |>
select(ets_damped) |>
report()Series: Turnover
Model: ETS(M,Ad,M)
Smoothing parameters:
alpha = 0.4848494
beta = 0.01207443
gamma = 0.1545561
phi = 0.9799992
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
44.99204 0.338861 0.8642028 0.7674952 0.9394509 1.453263 1.064383 1.044351
s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
1.014907 0.9744656 0.9820932 0.9772068 1.008291 0.9098907
sigma^2: 0.0026
AIC AICc BIC
4605.481 4607.102 4679.084 report(myseries_fit)# A tibble: 2 × 11
State Industry .model sigma2 log_lik AIC AICc BIC MSE AMSE MAE
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Queensla… Clothin… ets_m… 0.00248 -2282. 4598. 4599. 4667. 123. 151. 0.0378
2 Queensla… Clothin… ets_d… 0.00255 -2285. 4605. 4607. 4679. 128. 163. 0.0380Inspecting the results of the models, they are all very similar. The non-damped version is the preferred version based on the AICc. Ultimately, the difference is small enough that there is not a compelling reason to select one model over the other.
c.
myseries_fit |>
accuracy() |>
select(.model, RMSE)# A tibble: 2 × 2
.model RMSE
<chr> <dbl>
1 ets_mult 11.1
2 ets_damped 11.3The RMSE for both models are also relatively similar. The Holt-Winters undamped model has a slightly lower value make it preferred in this instance but the difference to the damped version is minimal.
d.
myseries_fit |>
select(ets_mult) |>
augment() |>
autoplot(.resid)
The residuals are plotted above. The residual variation does not appear exclusively as white noise. The magnitude of the variation increases along with the increase in level of the series data but there is no discernible pattern within the residuals.
myseries_fit |>
select(ets_mult) |>
gg_tsresiduals()
The innovation residuals fluctuate around zero with roughly constant variance. Some of the spikes exceed the dashed bounds, especially at a 12 month lag, indicating that the model has not captured all of the seasonality and showing substantial autocorrelation. While they have a bell shaped distribution, the residuals do not fully resemble white noise.
e.
The data is split into a training and test set in the code below. Then an ETS(M, A, M) model is fit to the training data and forecasted to the test data.
# create a training dataset of observations before 2011 and a test set from 2011 on
myseries_train <- myseries |> filter(year(Month) <= 2010)
myseries_test <- myseries |> filter(year(Month) > 2010)
train_fit <- myseries_train |>
model(
ets_mult_train = ETS(Turnover ~ error('M') + trend('A') + season('M')))
train_forecast <- train_fit |>
forecast(new_data = myseries_test)
train_forecast |>
autoplot(myseries_train, level = NULL)
Below the RMSE for the is calculated below.
train_forecast |>
as_tibble() |>
select(State, Industry, .model, Month, .mean) |>
left_join(
myseries_test |>
as_tibble() |>
select(State, Industry, Month, Turnover),
by = c('State', 'Industry','Month')
) |> group_by(State, Industry, .model) |>
summarise(
RMSE = sqrt(mean((Turnover - .mean)^2)),
.groups = 'drop'
)# A tibble: 1 × 4
State Industry .model RMSE
<chr> <chr> <chr> <dbl>
1 Queensland Clothing, footwear and personal accessory retailing ets_mult… 28.3A seasonal naive model is fitted to the training data below and then a forecast is made and the RMSE calculated.
# fit a seasonal naive model to the training data
fit <- myseries_train |> model(SNAIVE(Turnover))
forecast <- fit |>
forecast(myseries_test)
forecast |>
as_tibble() |>
select(State, Industry, .model, Month, .mean) |>
left_join(
myseries_test |>
as_tibble() |>
select(State, Industry, Month, Turnover),
by = c('State', 'Industry','Month')
) |> group_by(State, Industry, .model) |>
summarise(
RMSE = sqrt(mean((Turnover - .mean)^2)),
.groups = 'drop'
) |> select(RMSE)# A tibble: 1 × 1
RMSE
<dbl>
1 52.3The RMSE from the Holt-Winters model has an RMSE of 28.25067. The seasonal naive model has an RMSE of 52.28899. The Holt-Winters model reduced the RMSE substantially from the seasonal naive model.
8.9
The retail data is plotted again below.
myseries# A tsibble: 441 x 5 [1M]
# Key: State, Industry [1]
State Industry `Series ID` Month Turnover
<chr> <chr> <chr> <mth> <dbl>
1 Queensland Clothing, footwear and personal acc… A3349884J 1982 Apr 42.8
2 Queensland Clothing, footwear and personal acc… A3349884J 1982 May 45.3
3 Queensland Clothing, footwear and personal acc… A3349884J 1982 Jun 45.1
4 Queensland Clothing, footwear and personal acc… A3349884J 1982 Jul 43.1
5 Queensland Clothing, footwear and personal acc… A3349884J 1982 Aug 41.5
6 Queensland Clothing, footwear and personal acc… A3349884J 1982 Sep 44.5
7 Queensland Clothing, footwear and personal acc… A3349884J 1982 Oct 48.6
8 Queensland Clothing, footwear and personal acc… A3349884J 1982 Nov 51.4
9 Queensland Clothing, footwear and personal acc… A3349884J 1982 Dec 70.9
10 Queensland Clothing, footwear and personal acc… A3349884J 1983 Jan 44.6
# ℹ 431 more rowsautoplot(myseries, Turnover) +
autolayer(myseries_train, Turnover, colour = "red")
A model is fit to the data below consisting of applying an ETS model to an STL decomposition of the Box-Cox transformation of the retail data.
lambda <- myseries_train |>
features(Turnover, guerrero) |>
pull(lambda_guerrero)
stl_fit <-myseries_train |>
model(
stl_ets = decomposition_model(
STL(
box_cox(Turnover, lambda) ~ season(window = 'periodic')
),
ETS(season_adjust ~ season('N'))
)
)
stl_forecast <- stl_fit |>
forecast(new_data = myseries_test)
stl_forecast |>
accuracy(myseries_test) |>
select(.model, RMSE)# A tibble: 1 × 2
.model RMSE
<chr> <dbl>
1 stl_ets 83.0The model created from by applying an ETS model to a STL decomposition of the Box-Cox transformation of the data results in an RMSE of 82.97691. When compared with the 28.25 for the Holt-Winters model and the 52.29 seasonal naive, the ETS/STL/Box-Cox model performs worse than both of the alternative on this test set. From these three models, the undamped Holt-Winters is preferred.