Data624

library(fpp3)

## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──

## ✔ tibble      3.2.1     ✔ tsibble     1.1.4
## ✔ dplyr       1.1.3     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.0     ✔ feasts      0.3.1
## ✔ lubridate   1.9.3     ✔ fable       0.3.3
## ✔ ggplot2     3.4.4     ✔ fabletools  0.3.4

## ── 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()

library(caret)

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following objects are masked from 'package:fabletools':
## 
##     MAE, RMSE

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ readr   2.1.4
## ✔ purrr   1.0.2     ✔ stringr 1.5.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ gridExtra::combine() masks dplyr::combine()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ purrr::lift()        masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tsibble)

Do the exercises 9.1, 9.2, 9.3, 9.5, 9.6, 9.7, 9.8 in Hyndman.

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

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

Yes, if the alpha level for any of these generated series is at 5% than out of the 20 lagged correlations there would typically be about 1 significant lag in white noise data. All three of the ACF plots have less significant correlation values than that level which would appear to indicate they are equivalent to random white noise data. The significance line for each of the plots is different and at a alternative critical value. These differencees make sense given that they are all sourced from different input data.

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?

\(\frac{2}{\sqrt{T}}\)

The calculation for the 95% confidence interval is dependent on the number of time periods or items in the sample as identified in the formula above and so the intervals for large time ranges are going to lower the thresholds allowable for residuals to be considered white noise. The relationships among lagged values will be different among each group of random numbers and white noise has many varieties that are “random”. The variance among the different series is also likely to be different, but it’s hard to say the exact relationship between each set’s variabililty; however, on paper the expectation is that smaller datasets would tend to have more spread than larger ones on average.

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.

gafa_stock |> filter(Symbol=='AMZN') |> gg_tsdisplay(Close,
               plot_type='partial', lag=10) +
  labs(title="Amazon Closing Stock Price: Non-Stationary", y="")

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

It is immediately apparent from the initial line plot of the closing price that taking different cuts of data along the time frame would result in completely variable results and distributions. This time series by definition breaks the assumptions needed to classify data as stationary and also shows very strong lagged correlations in the ACF due to the strong cyclic trends that are present. The partial lag plot if we were evaluating the data would indicate a potential MA term of 5 based on the last significant correlation, but otherwise does not indicate that much in this case.

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

Turkish GDP from global_economy.

turk_gdp <- global_economy |> filter(Country=='Turkey') |> dplyr::select(GDP)

turk_gdp |> autoplot(GDP)

turk_bc <- caret::preProcess(x=as.data.frame(turk_gdp),method='BoxCox')
turk_bc$bc$GDP$lambda

## [1] 0

Based on the maximized log likelihood lambda returned it appears that a log transformation is the best option for this data.

autoplot(turk_gdp,log(GDP))

Given that this is annual data there isn’t seasonal differencing that can applied to prepare the time series for stationarity.

turk_gdp |>
  gg_tsdisplay(difference(log(GDP)),
               plot_type='partial', lag=36) +
  labs(title="First Level differenced GDP", y="")

## Warning: Removed 1 row containing missing values (`geom_line()`).

## Warning: Removed 1 rows containing missing values (`geom_point()`).

Aside from an initial change after the first recorded GDP, the log transformed data does appear to be approximately stationary after differencing one time.

Accommodation takings in the state of Tasmania from aus_accommodation.

aus_accommodation |> filter(grepl('Tasmania',State)) |> autoplot(Takings)

At first glance the variance within the seasonality of the data appears to be growing after 2005 and will likely require a transformation.

tas_tak <- aus_accommodation |> filter(State=='Tasmania') |> dplyr::select(Takings)
tas_bc <- caret::preProcess(as.data.frame(tas_tak),method='BoxCox')
tas_bc$bc$Takings$lambda

## [1] 0.3

autoplot(tas_tak,Takings^tas_bc$bc$Takings$lambda)

This plot appears to have much more consistent variance after applying the cube root transformation.

tas_tak |>
  gg_tsdisplay(difference(Takings^tas_bc$bc$Takings$lambda,4),
               plot_type='partial', lag=12) +
  labs(title="Seasonal differenced Takings", y="")

## Warning: Removed 4 rows containing missing values (`geom_line()`).

## Warning: Removed 4 rows containing missing values (`geom_point()`).

It would appear that taking different cuts of time periods will result in potentially different data values and distributions. Therefore let’s attempt one more differencing step to see how it will impact the dataset.

tas_tak |>
  gg_tsdisplay(difference(Takings^tas_bc$bc$Takings$lambda,4) |> difference(),
               plot_type='partial', lag=12)  +
  labs(title="Seasonal and second level differenced Takings", y="")

## Warning: Removed 5 rows containing missing values (`geom_line()`).

## Warning: Removed 5 rows containing missing values (`geom_point()`).

Taking a second difference appears to do a better job at creating a time series that is more stationary that would better satisfy the requirements of an ARIMA model.

Monthly sales from souvenirs.

souvenirs |> autoplot(Sales)

The souvenir sales are clearly exhibiting seasonal patterns and need to be seasonally differenced at a minimum to satisfy stationarity conditions. Let’s determine if any transformation makes sense as there does appear to be a larger seasonal spike towards the end of the dataset.

souv_bc <- preProcess(as.data.frame(souvenirs),method='BoxCox')
sales_lambda <- souv_bc$bc$Sales$lambda
sales_lambda

## [1] -0.2

Let’s review what the power transformed data would look like given that -0.2 is not the most interpretable transformation.

autoplot(souvenirs,Sales^sales_lambda)

The seasonal patterns appears to have inverted which makes some sense given the negative power transform. From reviewing this graphic it does appear that that seasonal variability is a bit more consistent.

souvenirs |> gg_tsdisplay(difference(Sales^sales_lambda,lag=12),
                          plot_type='partial',lag=36) +
    labs(title = 'Seasonally Differenced Souvenir Sales')

## Warning: Removed 12 rows containing missing values (`geom_line()`).

## Warning: Removed 12 rows containing missing values (`geom_point()`).

There is still some trend exhibited in the seasonally differenced data and it would likely require one additional differencing step to try to make the data stationary.

souvenirs |> gg_tsdisplay(difference(Sales^sales_lambda,lag=12) |> difference(),
                          plot_type='partial',lag=36) +
    labs(title = 'Double Differenced Souvenir Sales')

## Warning: Removed 13 rows containing missing values (`geom_line()`).

## Warning: Removed 13 rows containing missing values (`geom_point()`).

The time series does appear to be mostly stationary by adding an additional layer of differencing after first applying seasonal step.

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(12345678)

sample_aus_retail <- aus_retail |>
  filter(`Series ID` == sample(aus_retail$`Series ID`,1))

sample_aus_retail |> autoplot(Turnover) + labs(title='Clothing, footwear turnover')

There is considerable change in the seasonal variance that exists in this model and it very likely needs a power transformation to further optimize downstream forecasting.

aus_ret_bc <- caret::preProcess(as.data.frame(sample_aus_retail),method='BoxCox')
aus_bc <- aus_ret_bc$bc$Turnover$lambda

aus_bc

## [1] 0.5

The recommended power transformation is the square root of the data.

autoplot(sample_aus_retail,Turnover^aus_bc)

It is still rather difficult to tell if the transformed series has elinimated the change in variance across the time frame but it is less noticeable than without the transformation.

sample_aus_retail |> gg_tsdisplay(difference(Turnover^aus_bc,12),
                                  plot_type='partial',lag=36) +
    labs(title='Seasonally Differenced Turnover Data')

## Warning: Removed 12 rows containing missing values (`geom_line()`).

## Warning: Removed 12 rows containing missing values (`geom_point()`).

This seasonally differenced data may be stationary although there are very strong autocorrelations that remain as evidenced by the ACF and PACF. Let’s review an additional difference step to see if that improves the results.

sample_aus_retail |> gg_tsdisplay(difference(Turnover^aus_bc,12) |> difference(),
                                  plot_type='partial',lag=36) +
    labs(title='Double Differenced Turnover Data')

## Warning: Removed 13 rows containing missing values (`geom_line()`).

## Warning: Removed 13 rows containing missing values (`geom_point()`).

From reviewing the time series after second level of differencing it is very likely stationary after this second step; however, there are still a good amount of autocorrelation that exists that could potentially be used for parameterization.

Simulate and plot some data from simple ARIMA models.

Use the following R code to generate data from an AR(1) model with \(\phi=0.6\) and \(\sigma^2=1\) The process starts with \(y_1=0\).

y <- numeric(100)
e <- rnorm(100)
phis <- seq(0.6,0.8,by=0.05)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

create_arma<-function(phi,phi2,theta,theta2) {
    y <- numeric(100)
    e <- rnorm(100)
    for(i in 3:100)
        y[i] <- phi*y[i-1]+phi2*y[i-2] + e[i] + theta*e[i-1]+theta2*e[i-2]
    sim <- tsibble(idx = seq_len(100), y = y, index = idx)
    return(sim)
}

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

autoplot(sim,y)

phi_07 <- create_arma(phi=0.7,phi2=0,theta=1,theta2=0)
phi_075 <- create_arma(phi=0.75,phi2=0,theta=1,theta2=0)
phi_08 <- create_arma(phi=0.8,phi2=0,theta=1,theta2=0)
phi_085 <- create_arma(phi=0.85,phi2=0,theta=1,theta2=0)
phi_09 <- create_arma(phi=0.9,phi2=0,theta=1,theta2=0)
cowplot::plot_grid(autoplot(sim,y)+labs(title='Phi=0.6'),autoplot(phi_07,y) +labs(title='Phi = 0.7'),autoplot(phi_075,y)+labs(title='Phi = 0.75'),autoplot(phi_08,y)+labs(title='Phi = 0.8'),autoplot(phi_085,y)+labs(title='Phi = 0.85'),autoplot(phi_09,y)+labs(title='Phi = 0.9'),nrow=3,ncol=2)

As \(\phi_1\) changes in value it affects the time series adjustment to account for the prior lags values. The higher values are impacting the localized auto regression but it’s somewhat hard to tell exactly what is occurring with the increase.

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

Repurposing the same function for AR models with MA

ma_1_06 <- create_arma(phi=0,phi2=0,theta=0.6,theta2=0)
ma_1_07 <- create_arma(phi=0,phi2=0,theta=0.7,theta2=0)
ma_1_075 <- create_arma(phi=0,phi2=0,theta=0.75,theta2=0)
ma_1_08 <- create_arma(phi=0,phi2=0,theta=0.8,theta2=0)
ma_1_085 <- create_arma(phi=0,phi2=0,theta=0.85,theta2=0)
ma_1_09 <- create_arma(phi=0,phi2=0,theta=0.9,theta2=0)

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

cowplot::plot_grid(autoplot(ma_1_06,y)+labs(title='Theta=0.6'),autoplot(ma_1_07,y) +labs(title='Theta = 0.7'),autoplot(ma_1_075,y)+labs(title='Theta = 0.75'),autoplot(ma_1_08,y)+labs(title='Theta = 0.8'),autoplot(ma_1_085,y)+labs(title='Theta = 0.85'),autoplot(ma_1_09,y)+labs(title='Theta = 0.9'),nrow=3,ncol=2)

As you change the parameter impacting the first error term it changes the ability of the model to be able to adjust for the prior lagged residual. It is somewhat hard to say that the parameter consistently impacts the rest of the time series.

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

arma_101 <- create_arma(phi=0.6,phi2=0,theta=0.6,theta2=0)

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

arma_200 <- create_arma(phi=-0.8,phi2=0.3,theta=0,theta2=0)

Graph the latter two series and compare them.

cowplot::plot_grid(autoplot(arma_101,y)+labs(title='ARMA(1,0,1)'),autoplot(arma_200,y)+labs(title='AR(2)'),nrow=2)

The ARMA(1,0,1) appears to show some minor cycle that exists and does not seem like it stationary while the AR(2) is gradually increasing the seasonality with every few iterations as it goes through the sequence.

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

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 |> autoplot(Passengers)

From a preliminary review of the passengers time series it does not appear that there is any seasonality to address and it’s not immediately clear if any transformation is necessary as the relationship appears to be somewhat linear. There will need to be differencing done in order to get rid of the trend and create stationarity.

Let’s run a Box Cox to cover our bases:

air_bc <- caret::preProcess(as.data.frame(aus_airpassengers),method='BoxCox')
air_bc$bc$Passengers$lambda

## [1] 0

The proposed transformation recommends a log transform and the plot below shows the time plot after applying that power transformation:

autoplot(aus_airpassengers,log(Passengers))

Let’s take a quick pass at the distributions of each series:

orig_pass <- aus_airpassengers |>
    ggplot(aes(x=Passengers)) +
    geom_histogram() +
    labs(title='Australian Air Passengers')

log_pass <- aus_airpassengers |>
    ggplot(aes(x=log(Passengers))) +
    geom_histogram() +
    labs(title=' Log Australian Air Passengers')

cowplot::plot_grid(orig_pass,log_pass)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

The regular non-transformed data set is moderately skewed to the right,but the log power transform still appears to show most of the same skew. Therefore, the ARIMA model will be run on the standard dataset.

pass_fit <- aus_airpassengers |>
    model(auto=ARIMA(Passengers))

report(pass_fit)

## 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 default automatically selected ARIMA model was (0,2,1) which means that two differences were taken to get to stationarity and 1 \(\theta\) backshift operator was determined adequate to account for the error in the model.

gg_tsresiduals(pass_fit |> dplyr::select(auto))

The residuals do appear to satisfy optimal assumptions for an ARIMA model as there are not significant autocorrelations, the spread of the data is fairly Gaussian shaped (with a few large outliers that might need to be reviewed), and homoskedastic nature of the innovations residuals.

augment(pass_fit) |>
  filter(.model == "auto") |>
  features(.innov, ljung_box, lag=24, dof=1)

## # A tibble: 1 × 3
##   .model lb_stat lb_pvalue
##   <chr>    <dbl>     <dbl>
## 1 auto      24.0     0.405

The Ljung Box test would indicate that we fail to reject that the residuals appear to be white noise.

pass_fit |>
    forecast(h=10) |>
    autoplot(aus_airpassengers) +
    labs(title='Arima (0,2,1)')

Write the model in terms of the backshift operator.

\(y_t = (1-B)^2 * (1+\theta B)\epsilon_t\)

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

fit_010 <- aus_airpassengers |>
    model(arima010 = ARIMA(Passengers ~ 1 + pdq(0,1,0)))
fit_010 |> forecast(h=10) |>
    autoplot(aus_airpassengers) +
    labs(title='Arima (0,1,0)')

The forecasts are fairly similar although part A has a larger forecast distribution which is harder to tell based on tick marks scaled. If the differencing is needed it might not be reasonable to use an ARIMA without stationarity despite its similarity.

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_212 <- aus_airpassengers |>
    model(arima010 = ARIMA(Passengers ~ 1 + pdq(2,1,2)))
fit_212 |> forecast(h=10) |>
    autoplot(aus_airpassengers) +
    labs(title='Arima (2,1,2)')

There appears to be more variability in the forecast that show a little uncertainty that’s reflected in the distribution shape as compared to the other models reviewed. The 80% interval seems to be slightly bigger although the full distribution looks to span a very similar area to part C

fit_212_noc <- aus_airpassengers |>
    model(arima010 = ARIMA(Passengers ~ 0 + pdq(2,1,2)))

## Warning: 1 error encountered for arima010
## [1] non-stationary AR part from CSS

This generates an invalid model likely stemming from the strong consistent trend that the constant is better able to account for even after differencing.

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

arima021c_fit <- aus_airpassengers |> 
    model(arima021c=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.

arima021c_fit |>
    forecast(h=10) |>
    autoplot(aus_airpassengers)

The fable package indicates that it is inadvisable to include 2 steps of differencing with a constant as it will induce higher polynomial trends. It is clear that the forecast generated appears to be growing with a steep slope.

For the United States GDP series (from global_economy):

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

global_economy |> filter(Country=='United States') |> autoplot(GDP)

us_gdp <- global_economy |> filter(Country=='United States') |> dplyr::select(GDP)

us_gdp_bc <- caret::preProcess(as.data.frame(us_gdp),method='BoxCox')

The best lambda is:

us_gdp_bc$bc$GDP$lambda

## [1] 0.2

us_gdp |> autoplot(GDP^0.2)

The transformation does seem to make the trend fairly linear as compared to the initial data although less interpretable.

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

us_gdp_fit <- us_gdp |> model(auto = ARIMA(GDP^0.2))

report(us_gdp_fit)

## Series: GDP 
## Model: ARIMA(1,1,0) w/ drift 
## Transformation: GDP^0.2 
## 
## Coefficients:
##          ar1  constant
##       0.4879    2.0618
## s.e.  0.1170    0.1724
## 
## sigma^2 estimated as 1.808:  log likelihood=-96.86
## AIC=199.72   AICc=200.17   BIC=205.84

The default model when using approximations and stepwise Hyndman algorithm indicates an ARIMA(1,1,0) with drift. The single level of differencing makes some sense given the clear and linear relationship of the trend.

us_gdp |> gg_tsdisplay(difference(GDP^0.2),plot_type='partial',lag=36)

## Warning: Removed 1 row containing missing values (`geom_line()`).

## Warning: Removed 1 rows containing missing values (`geom_point()`).

Plotting the residuals after applying one difference would indicate that an AR(1) or MA(1) might make sense as a model option which the automatic criteria selected AR(1).

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

other_usgdp_fit <- us_gdp |> 
    model(arima_011=ARIMA(GDP^0.2 ~ 0 + pdq(0,1,1)),
          arima_011_wc=ARIMA(GDP^0.2 ~ 1 + pdq(0,1,1)),
          arima_111=ARIMA(GDP^0.2 ~ 0 + pdq(1,1,1)))

report(other_usgdp_fit |> dplyr::select(arima_011))

## Series: GDP 
## Model: ARIMA(0,1,1) 
## Transformation: GDP^0.2 
## 
## Coefficients:
##          ma1
##       0.7757
## s.e.  0.0678
## 
## sigma^2 estimated as 7.872:  log likelihood=-139.64
## AIC=283.27   AICc=283.49   BIC=287.36

report(other_usgdp_fit |> dplyr::select(arima_011_wc))

## Series: GDP 
## Model: ARIMA(0,1,1) w/ drift 
## Transformation: GDP^0.2 
## 
## Coefficients:
##          ma1  constant
##       0.4124    4.0568
## s.e.  0.1083    0.2530
## 
## sigma^2 estimated as 1.915:  log likelihood=-98.46
## AIC=202.92   AICc=203.38   BIC=209.05

report(other_usgdp_fit |> dplyr::select(arima_111))

## Series: GDP 
## Model: ARIMA(1,1,1) 
## Transformation: GDP^0.2 
## 
## Coefficients:
##          ar1      ma1
##       0.9904  -0.6487
## s.e.  0.0122   0.1366
## 
## sigma^2 estimated as 1.936:  log likelihood=-99.9
## AIC=205.8   AICc=206.25   BIC=211.93

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

Considering that all 3 models have the same level of differencing it is permissible to use AICc as the criteria to evaluate the preferred model. Using that metric, the ARIMA(0,1,1) w/ drift appears to be the best model of the samples options chosen. It was slightly less accurate than the automated choice, but we will select based on the manual reviewed options.

gg_tsresiduals(other_usgdp_fit |> dplyr::select(arima_011_wc))

The residuals do not appear to violate relevant assumptions as they seem homoskedastic and do not show significant autocorrelation among lagged values. The distribution is close enough to approximately normal although there is a large outlier on the left which probably corresponds to the lowest residual from around 2010.

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

other_usgdp_fit |> dplyr::select(arima_011_wc) |> forecast(h=10) |>
    autoplot(us_gdp)

The forecasts seem reasonable enough and appear to track with the trend/cycle that exists in the data. The one caveat is that it reflects an expectation that the the same rate of growth with primarily continue unchanged and perhaps a damping trend will more accurately reflect the inconsistencies in economic cycles.

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

Let’s load two different ETS models using a damped trend based on the review of the last ARIMA model:

usgdp_ets_fit <- us_gdp |> model(aadn=ETS(GDP ~ error('A')+trend('Ad')+season('N')),
                                 amdn=ETS(GDP ~ error('A')+trend('Md')+season('N')))

usgdp_ets_fit |>
    forecast(h=10) |>
    autoplot(us_gdp)

The forecasts for each ETS model seem relatively similar to the values estimated by the ARIMA model. The primary difference will likely relate to the long term convergence of ARIMA versus the ETS methods. The exponential damping perhaps would be advantageous in this regard if forecasts needed to extend to a point farther out in the future. Both additive and multiplicative trend terms were included to compare if the appearance of a strong upward growth cycle would be better captured by a multiplicative trend component.

Let’s review the accuracy metrics for the two ETS methods:

knitr::kable(accuracy(usgdp_ets_fit) |> dplyr::select(c('.model','RMSE','MAE','MASE','RMSSE')))

.model	RMSE	MAE	MASE	RMSSE
aadn	167494825481	106487133842	0.3121186	0.4108994
amdn	154901893866	88142753333	0.2583504	0.3800063

Given the scale of the data is different than ARIMA the RMSE and likely most other error metrics using average error would indicate far more error values in these models even if the damping parameter does a better job for future unknown predictions farther out in the future. The range of the predicted estimates might prompt different model selection and would be an important consideration when evaluating the available models.

Data624_HW6

jfor: Josh Forster

2024-02-24