Naive, Basic, Dynamic, and Advanced Forecasting Models

Libraries

I often like to centralize my libraries, particularly if I am going to cite them later. (I don’t like to forget any citations of the packages.) Use the cite() function for full details.

library(eia)

## EIA_KEY not found. See `vignette("api", "eia")` for key storage options.

library(fpp3)

## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr

## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.0 ──

## ✔ tibble      3.2.1     ✔ tsibble     1.1.5
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.3.2
## ✔ lubridate   1.9.3     ✔ fable       0.3.4
## ✔ ggplot2     3.5.1     ✔ fabletools  0.4.2

## ── 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(fredr)
library(ggplot2)
library(httr)
library(jsonlite)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

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

library(latex2exp)
library(magrittr)

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:tidyr':
## 
##     extract

library(readxl)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0     ✔ readr   2.1.5
## ✔ purrr   1.0.2     ✔ stringr 1.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ magrittr::extract()      masks tidyr::extract()
## ✖ dplyr::filter()          masks stats::filter()
## ✖ purrr::flatten()         masks jsonlite::flatten()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ tsibble::interval()      masks lubridate::interval()
## ✖ dplyr::lag()             masks stats::lag()
## ✖ purrr::set_names()       masks magrittr::set_names()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Keys

I have API keys for foundational TS data (Fred and EIA). Obviously, I keep them hidden. :)

Data

So the easiest way to get FRED data is to invoke the fredr package.

unemp=fredr(
  series_id = "UNRATE",
  observation_start = as.Date("2020-04-01"),
  observation_end = as.Date("2025-03-31")
)

Look at the Data

A quick peak…

rbind(head(unemp),tail(unemp))

## # A tibble: 12 × 5
##    date       series_id value realtime_start realtime_end
##    <date>     <chr>     <dbl> <date>         <date>      
##  1 2020-04-01 UNRATE     14.8 2025-07-06     2025-07-06  
##  2 2020-05-01 UNRATE     13.2 2025-07-06     2025-07-06  
##  3 2020-06-01 UNRATE     11   2025-07-06     2025-07-06  
##  4 2020-07-01 UNRATE     10.2 2025-07-06     2025-07-06  
##  5 2020-08-01 UNRATE      8.4 2025-07-06     2025-07-06  
##  6 2020-09-01 UNRATE      7.8 2025-07-06     2025-07-06  
##  7 2024-10-01 UNRATE      4.1 2025-07-06     2025-07-06  
##  8 2024-11-01 UNRATE      4.2 2025-07-06     2025-07-06  
##  9 2024-12-01 UNRATE      4.1 2025-07-06     2025-07-06  
## 10 2025-01-01 UNRATE      4   2025-07-06     2025-07-06  
## 11 2025-02-01 UNRATE      4.1 2025-07-06     2025-07-06  
## 12 2025-03-01 UNRATE      4.2 2025-07-06     2025-07-06

Recoding / Variable Drops

I keep it easy and drop unnecessary variables. Tsibbles are also particularly about dates, so I force the months to be compliant.

unemp$Month=seq(as.Date("2020/4/1"), as.Date("2025/3/31"), by="1 month")
unemp$Month=yearmonth(unemp$Month)
unemp$series_id=NULL
mydata=unemp
str(mydata)

## tibble [60 × 5] (S3: tbl_df/tbl/data.frame)
##  $ date          : Date[1:60], format: "2020-04-01" "2020-05-01" ...
##  $ value         : num [1:60] 14.8 13.2 11 10.2 8.4 7.8 6.9 6.7 6.7 6.4 ...
##  $ realtime_start: Date[1:60], format: "2025-07-06" "2025-07-06" ...
##  $ realtime_end  : Date[1:60], format: "2025-07-06" "2025-07-06" ...
##  $ Month         : mth [1:60] 2020 Apr, 2020 May, 2020 Jun, 2020 Jul, 2020 Aug, 2020 Sep,...

Train / Test Split

Even when I am working on exploration, I like a train / test split. This is a foundational principal for ML. I define the TS and then set the split.

myts <- mydata %>%as_tsibble(index = Month)
train <- myts[1:48, ]
test <- myts[49:nrow(myts), ]

Viewing and LOESS

I am going to look at the data and fit a LOESS to the train set. I will then (for simplicity) use a linear forecast of the last 6 months of LOESS data and forecast the test set. Finally, I calculate metrics for this super simple forecasting method.

library(tsibble)
library(dplyr)
library(ggplot2)

# Create the tsibble
myts <- mydata %>%
  as_tsibble(index = Month)

# Split train/test
train <- myts[1:48, ]
test <- myts[49:nrow(myts), ]

# Convert Month to numeric for modeling
myts2 <- myts %>%
  mutate(x_numeric = as.numeric(Month))

train2 <- myts2[1:48, ]
test2 <- myts2[49:nrow(myts2), ]

# Fit LOESS on training data
loess_model <- loess(
  value ~ x_numeric,
  data = train2,
  span = 0.75
)

# Get LOESS fitted values for all train data
train2 <- train2 %>%
  mutate(
    loess_fit = predict(loess_model)
  )

# Select last 6 months of training for linear extension
last6 <- tail(train2, 6)

# Fit linear model on last 6 LOESS points
tail_model <- lm(
  loess_fit ~ x_numeric,
  data = last6
)

# Predict over test set with confidence intervals
predictions <- predict(
  tail_model,
  newdata = test2,
  interval = "confidence",
  level = 0.95
)

# Bind predictions into test2
test2 <- test2 %>%
  mutate(
    linear_extension = predictions[,"fit"],
    lower = predictions[,"lwr"],
    upper = predictions[,"upr"]
  )

# Compute forecast error metrics
errors <- test2$value - test2$linear_extension

ME <- mean(errors)
MPE <- mean(errors / test2$value) * 100
MAE <- mean(abs(errors))
RMSE <- sqrt(mean(errors^2))
MAPE <- mean(abs(errors) / test2$value) * 100

# For MASE: compute naive forecast MAE on training data
train_naive_errors <- train2$value[-1] - train2$value[-nrow(train2)]
naive_mae <- mean(abs(train_naive_errors))

MASE <- MAE / naive_mae

# Create a formatted metrics label
metrics_label <- sprintf(
  "Test Set Forecast Metrics:\nME = %.3f\nMPE = %.2f%%\nMAE = %.3f\nRMSE = %.3f\nMAPE = %.2f%%\nMASE = %.3f",
  ME, MPE, MAE, RMSE, MAPE, MASE
)

# Plot everything with metrics annotation
ggplot() +
  geom_line(data = train2, aes(x = Month, y = value, color = "Train Actual"), size = 1) +
  geom_line(data = test2, aes(x = Month, y = value, color = "Test Actual"), size = 1) +
  geom_line(data = train2, aes(x = Month, y = loess_fit, color = "LOESS Fit"), size = 1) +
  geom_line(data = test2, aes(x = Month, y = linear_extension, color = "Linear Extension"), size = 1) +
  geom_ribbon(
    data = test2,
    aes(x = Month, ymin = lower, ymax = upper, fill = "95% CI"),
    alpha = 0.3,
    color = "gray"
  ) +
  annotate(
    "text",
    x = min(test2$Month)-10 ,  # adjust x position if needed
    y = max(train2$value, na.rm=TRUE),  # adjust y position if needed
    label = metrics_label,
    hjust = 0,
    vjust = 1,
    size = 4,
    color = "black",
    fontface = "italic"
  ) +
  scale_color_manual(
    name = "Legend",
    values = c(
      "Train Actual" = "black",
      "Test Actual" = "red",
      "LOESS Fit" = "green",
      "Linear Extension" = "blue"
    )
  ) +
  scale_fill_manual(
    name = "",
    values = c("95% CI" = "gray")
  ) +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    title = "Monthly Unemployment Rate with LOESS, Linear Extension Forecast, and Error Metrics",
    x = "Month",
    y = "Unemployment Rate"
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

When we look at the data, there is clearly a leveling that occurs from a high of near 15% during the COVID spike to under 5%. Pushing a linear forecast forward from a LOESS fit of the previous 6 months is better than a simple naive forecast (MASE=0.8 = MAE Model Forecast / MAE Naive Model Forecast). The forecast is biased negatively (ME and MPE), and the MAE, RMSE, and MAPE suggest reasonably low variance. And we haven’t even decomposed the data yet.

Time Series Plots

Here, we take a look at some typical TS plots

myts%>%gg_season(value)

myts%>%gg_subseries(value)

myts%>%gg_lag(value, geom="point")

myts%>%ACF(value)%>%autoplot()+ylim(-1,1)

myts%>%PACF(value)%>%autoplot()+ylim(-1,1)

The seasonality ggseason plot reveals two things, a trend and a seasonal component. We can see that 2020 exhibited the highest unemployment followed by 2021, 2025, 2024, 2022, 2023. So we conclude that 2020 and 2021 exhibit true nonstationarity and that the trend effect is nonlinear. There is no obvious seasonal component from observing the data, but we will look further.

The ggsub_series plot appears to show a declining unemployment rate associated with the months April through December, but a word of caution…The outlier years have affected that appearance dramatically.

The lag plot shows obvious autocorrelation at lag 1 and some at 2 and 3 with it decreasing rapidly beyond that. The 45 degree reference line is the line if \(y_t=y_{t-k}\) where \(k\) is the lag. You can see how there is increasing spread and deviation from the 45 degree line as we increase lags. We will test to see what lags are appropriate at some point.

The ACF plot shows significant autocorrelation through 9 periods, finally dropping off at period 10. Recall, ACF is defined as \(r_{k} = \frac{\sum\limits_{t=k+1}^T (y_{t}-\bar{y})(y_{t-k}-\bar{y})}{\sum\limits_{t=1}^T (y_{t}-\bar{y})^2},\) This is just the correlation between \(y_t\) and \(y_{t-k}\)=$r_k = \(r_k=\frac{Cov(y_t, y_{t-k})}{Var(y_t)}\).

The PACF plot shows significant correlation only for 1 period. This plot removes the preceeding autocorrelation to see if there is remaining autocorrelation. If we want the PACF for \(y_{t-k}\), then we regress \(y_t\) on lags \(1, 2,...k-1\) and regress \(y_{t-k}\) on lags \(1, 2, k-1\). We correlate the residuals of both = PACF! What we don’t have is white noise,

Transformations

Our next step is transformation of the data prior to decomposition. In our case, we can use a Box-Cox transformation since the data are strictly non-negative. That transformation is defined as follows.

\[ w_t = \begin{cases} \log(y_t), & \text{if }\lambda = 0 \\[8pt] \displaystyle\frac{\operatorname{sign}(y_t)\,\lvert y_t\rvert^{\,\lambda}\,-\,1}{\lambda}, & \text{otherwise} \end{cases} \]

# Compute lambda using Guerrero method
lambda <- myts %>%
  features(value, features = guerrero) %>%
  pull(lambda_guerrero)

# Plot Box-Cox transformed series with LaTeX lambda label
myts %>%
  autoplot(box_cox(value, lambda)) +
  labs(
    y = "",
    title = latex2exp::TeX(paste0(
      "Transformed series with $\\lambda = ",
      round(lambda, 2),
      "$"
    ))
  )

We can see here that the lambda selected is almost a recipricol. But that lambda helped none.

Decomposition

We can decompose our data into different components: Seasonal, Trend, Residual. Additive decomposition assumes \(y_{t} = S_{t} + T_{t} + R_t\) whereas multiplicative decomposition is additive in the logs: \(y_{t} = S_{t} \times T_{t} \times R_t\). One of the common techniques is STL: Seasonal decomposition of time series by LOESS. The trend and seasonal components are measured using LOESS smoothing. Think of fitting a curve through all Januaries, all Februaries, etc. Think of fitting a nonlinear function for trend. Recall, LOESS fits a locally weighted polynomial regression:

\[ \hat{y}(x_0) = \sum_{i=1}^{n} w_i(x_0)\, y_i \]

\[ \sum_{i=1}^{n} K\left( \frac{x_i - x_0}{h} \right) \Bigl( y_i \;-\; p(x_i) \Bigr)^2 \] where

\[ K(u) = (1 - |u|^3)^3 \quad \text{for } |u| \le 1 \] - \(h\) is the bandwidth parameter controlling the span of the weights - \(p(x)\) is a local polynomial (usually linear or quadratic) - \(w_i(x_0)\) are the values of the kernel function applied to the scaled distance between \(x_i\) and \(x_0\): \(K \frac{x_i-x_0}{h}\) - \(x_0\) is the point at which the fitted value is estimated

comp=myts%>%model(stl=STL(value))%>%components()
comp2=myts%>%mutate(logvalue=log(value))%>%
  model(stl=STL(logvalue))%>%components()

comp%>%as_tsibble() %>%autoplot(value)+geom_line(aes(y=trend), colour = "red")+
  geom_line(aes(y=season_adjust), colour = "blue")

comp2%>%as_tsibble() %>%autoplot(exp(logvalue))+geom_line(aes(y=trend), colour = "red")+
  geom_line(aes(y=season_adjust), colour = "blue")

comp%>%autoplot()

comp2%>%autoplot()+labs(
  title = "log(value) = trend + seasonality + remainder",
  subtitle = "All components are in log scale"
)

So we have run both additive and multiplicative decomposition. Which is more appropriate for the model?

Decomposition Comparison

m1 <- myts %>% model(add = STL(value))
m2 <- myts %>% model(mult = STL(log(value)))
rbind(accuracy(m1),accuracy(m2))

## # A tibble: 2 × 10
##   .model .type         ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>  <chr>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 add    Training -0.0856 0.491 0.311 -2.03  5.77 0.217 0.209 0.613
## 2 mult   Training -0.0252 0.403 0.208 -1.09  3.47 0.146 0.172 0.543

Looking at the full data,the multiplicative has the lower ME and MPE (bias measures). It has the lower RMSE, MAE, and MAPE (variance measures). It has the smaller comparative measures (MASE and RMSSE). RMSSE is the ratio of mean squared error of the model versus a naive. Multiplicative wins. We use the log going forward.

Now, we cannnot forecast with a traditional STL. We have to use other decomposition methods. Why? The LOESS component…We would need more traditional decompositions to do so.

Forecasting Models (Basic)

So next we move to same basic forecasting models. We will use log(value) to help us since that appears to be a better solution by multiplicative decomposition.

#Log
train$logvalue=log(train$value)

#Naive Models
m1=train |>model(Mean_Model=MEAN(logvalue))
m2=train |>model(Drift_Model=RW(logvalue~drift()))
m3=train |>model(Naive_Model=NAIVE(logvalue))
m4=train |>model(Seasonal_Naive_Model=SNAIVE(logvalue))

#Basic Models
m5=train |>model(TS_Model = TSLM(logvalue))
m6=train |>model(ETS_Model=ETS(logvalue))
m7=train |>model(ARIMA_Model=ARIMA(logvalue))

#Dynamic Ensemble Models
m8=train|>model(Ensemble_Model=(ETS(logvalue)+ARIMA(logvalue))/2)
m9=train|>model(Ensemble_Model_2=(ETS(logvalue)+TSLM(logvalue~fourier(K=2))+ARIMA(logvalue))/3)

The first model is nothing more than a mean model. We acquire the mean of the TS and use that for foreasting. Simple.

The second model is a drift model (not a linear fit). We acquire the first and last observations, and use these to generate a drift component. To easy.

The third model is a naive model. Next period forecast is last period observation. Obviously, it runs out of juice for longer forecasts.

The fourth model is a seasonal naive. Next January forecast is this January’s observation. This type of forecast can work well for data that only has seasonality (no trend).

The fifth model is a time series linear regression model, a model that works well for data that only has a trend.

The sixth model is an Error, Trend, and Seasonality Model. We allow the software to auto-choose the terms. More on this later in the course.

The seventh model is an ARIMA, autoregressive integrated moving average. Again, later in the course.

Model 8 and 9 are simple ensembles. We will discuss later.

Reports

Of course, we want complete reports for the estimates.

for (i in 1:9){report(eval(parse(text=paste0('m',i))))}

## Series: logvalue 
## Model: MEAN 
## 
## Mean: 1.5606 
## sigma^2: 0.1472 
## Series: logvalue 
## Model: RW w/ drift 
## 
## Drift: -0.0284 (se: 0.0079)
## sigma^2: 0.003 
## Series: logvalue 
## Model: NAIVE 
## 
## sigma^2: 0.003 
## Series: logvalue 
## Model: SNAIVE 
## 
## sigma^2: 0.081 
## Series: logvalue 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3369 -0.2797 -0.1997  0.2477  1.1340 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.56064    0.05537   28.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3836 on 47 degrees of freedom
## Series: logvalue 
## Model: ETS(M,Ad,N) 
##   Smoothing parameters:
##     alpha = 0.5870324 
##     beta  = 0.3298261 
##     phi   = 0.8427385 
## 
##   Initial states:
##      l[0]       b[0]
##  2.922793 -0.2324905
## 
##   sigma^2:  6e-04
## 
##       AIC      AICc       BIC 
## -121.7800 -119.7312 -110.5528 
## Series: logvalue 
## Model: ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1     ma2
##       -0.8881  0.4431
## s.e.   0.1357  0.1359
## 
## sigma^2 estimated as 0.001706:  log likelihood=81.89
## AIC=-157.78   AICc=-157.21   BIC=-152.3
## Series: logvalue 
## Model: COMBINATION 
## Combination: logvalue * 0.5
## 
## ===========================
## 
## Series: logvalue 
## Model: COMBINATION 
## Combination: logvalue + logvalue
## 
## ================================
## 
## Series: logvalue 
## Model: ETS(M,Ad,N) 
##   Smoothing parameters:
##     alpha = 0.5870324 
##     beta  = 0.3298261 
##     phi   = 0.8427385 
## 
##   Initial states:
##      l[0]       b[0]
##  2.922793 -0.2324905
## 
##   sigma^2:  6e-04
## 
##       AIC      AICc       BIC 
## -121.7800 -119.7312 -110.5528 
## 
## Series: logvalue 
## Model: ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1     ma2
##       -0.8881  0.4431
## s.e.   0.1357  0.1359
## 
## sigma^2 estimated as 0.001706:  log likelihood=81.89
## AIC=-157.78   AICc=-157.21   BIC=-152.3
## 
## 
## Series: logvalue 
## Model: COMBINATION 
## Combination: logvalue * 0.333333333333333
## 
## =========================================
## 
## Series: logvalue 
## Model: COMBINATION 
## Combination: logvalue + logvalue
## 
## ================================
## 
## Series: logvalue 
## Model: COMBINATION 
## Combination: logvalue + logvalue
## 
## ================================
## 
## Series: logvalue 
## Model: ETS(M,Ad,N) 
##   Smoothing parameters:
##     alpha = 0.5870324 
##     beta  = 0.3298261 
##     phi   = 0.8427385 
## 
##   Initial states:
##      l[0]       b[0]
##  2.922793 -0.2324905
## 
##   sigma^2:  6e-04
## 
##       AIC      AICc       BIC 
## -121.7800 -119.7312 -110.5528 
## 
## Series: logvalue 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4506 -0.2499 -0.1057  0.1795  1.0351 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.56064    0.05583  27.951   <2e-16 ***
## fourier(K = 2)C1_12 -0.11102    0.07896  -1.406    0.167    
## fourier(K = 2)S1_12  0.06551    0.07896   0.830    0.411    
## fourier(K = 2)C2_12 -0.03336    0.07896  -0.423    0.675    
## fourier(K = 2)S2_12 -0.04850    0.07896  -0.614    0.542    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3868 on 43 degrees of freedom
## Multiple R-squared: 0.06969, Adjusted R-squared: -0.01685
## F-statistic: 0.8053 on 4 and 43 DF, p-value: 0.52862
## 
## 
## Series: logvalue 
## Model: ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1     ma2
##       -0.8881  0.4431
## s.e.   0.1357  0.1359
## 
## sigma^2 estimated as 0.001706:  log likelihood=81.89
## AIC=-157.78   AICc=-157.21   BIC=-152.3

Plots

And we want plots.

myts$logvalue=log(myts$value)
test$logvalue=log(test$value)

myplot=function(model){
  model|>forecast(test)|>autoplot(myts)+geom_line(aes(y=.fitted,col=.model),
                                                   data=augment(model))+ggtitle(names(model))}

for (i in 1:9){print(myplot(eval(parse(text=paste0('m',i)))))} #NNETAR takes long to plot

## 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_line()`).

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

Accuracy Metrics

And of course, accuracy metrics.

# Initialize an empty tibble to collect all results
all_metrics <- tibble()

# Loop over your 9 models
for (i in 1:9) {
  # Get the model object
  mod <- eval(parse(text = paste0("m", i)))
  
  # Forecast over the test set
  fc <- forecast(mod, new_data = test)
  
  # Compute accuracy using the *full dataset* so MASE/RMSSE are available
  acc <- accuracy(fc, data = myts)
  
  # Add a model name column
  acc <- acc %>% mutate(Model = names(mod))
  
  # Append to the collection tibble
  all_metrics <- bind_rows(all_metrics, acc)
}

# Clean up the output to only keep the metrics you care about
metrics_table <- all_metrics %>%
  select(Model, ME, MPE, RMSE, MAE, MAPE, MASE, RMSSE)

# Arrange by MAPE for example
metrics_table %>% arrange(MAPE)

## # A tibble: 9 × 8
##   Model                     ME    MPE   RMSE    MAE  MAPE   MASE  RMSSE
##   <chr>                  <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>
## 1 ETS_Model             0.0172   1.20 0.0258 0.0213  1.50 0.0712 0.0664
## 2 Ensemble_Model       -0.0411  -2.91 0.0598 0.0454  3.22 0.152  0.154 
## 3 Naive_Model           0.0498   3.50 0.0546 0.0498  3.50 0.166  0.140 
## 4 Ensemble_Model_2     -0.0773  -5.50 0.0804 0.0773  5.50 0.258  0.207 
## 5 ARIMA_Model          -0.0994  -7.03 0.125  0.0994  7.03 0.332  0.322 
## 6 Seasonal_Naive_Model  0.101    7.16 0.108  0.101   7.16 0.338  0.278 
## 7 Mean_Model           -0.150  -10.7  0.152  0.150  10.7  0.501  0.390 
## 8 TS_Model             -0.150  -10.7  0.152  0.150  10.7  0.501  0.390 
## 9 Drift_Model           0.234   16.5  0.258  0.234  16.5  0.783  0.664

The metrics show a clear winner in terms of bias. The lowest bias (ME and MPE measures) belong to the ETS model. In terms of variance, the ETS wins RMSE, MPE, and MAPE. Comparatively, the ETS wins MASE and RMSSE as well. ETS it is. Now, looking at the components of ETS, we have multiplicate error, additive trend, and no seasonality. More on this later.