Naive, Basic, Dynamic, and Advanced Forecasting Models

basic_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)           # for EIA data (not used below, but kept for parity with the original)

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

library(fpp3)          # tsibble + feasts + fable + ggplot2

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

── Attaching packages ──────────────────────────────────────────── fpp3 1.0.2 ──

✔ tibble      3.3.0     ✔ tsibble     1.1.6
✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
✔ tidyr       1.3.1     ✔ feasts      0.4.2
✔ lubridate   1.9.4     ✔ fable       0.4.1
✔ ggplot2     3.5.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)         # FRED API
library(ggplot2)       # plotting
library(httr)
library(jsonlite)
library(kableExtra)

Attaching package: 'kableExtra'

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

    group_rows

library(latex2exp)     # LaTeX labels in ggplot
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.1.0     ✔ 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.

# If you have a FRED key, uncomment and set it here:

# fredr_set_key(Sys.getenv("FRED_API_KEY"))

fredr_set_key("8a9ec1330374c1696f05cc8e526233b5") # replace with your own key please

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 peek…

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

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-09-07" "2025-09-07" ...
 $ realtime_end  : Date[1:60], format: "2025-09-07" "2025-09-07" ...
 $ 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 principle 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"), linewidth = 0.6) +
  geom_line(data = test2,  aes(x = Month, y = value, color = "Test Actual"),  linewidth = 0.6) +
  geom_line(data = train2, aes(x = Month, y = loess_fit,  color = "LOESS Fit"), linewidth = 0.6) +
  geom_line(data = test2,  aes(x = Month, y = linear_extension, color = "Linear Extension"), linewidth = 0.8) +
  geom_ribbon(
    data = test2,
    aes(x = Month, ymin = lower, ymax = upper, fill = "95% CI"),
    alpha = 0.3
  ) +
  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
  ) +
  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"
  )

We can see the 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). 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)

Warning: `gg_season()` was deprecated in feasts 0.4.2.
ℹ Please use `ggtime::gg_season()` instead.

myts %>% gg_subseries(value)

Warning: `gg_subseries()` was deprecated in feasts 0.4.2.
ℹ Please use `ggtime::gg_subseries()` instead.

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

Warning: `gg_lag()` was deprecated in feasts 0.4.2.
ℹ Please use `ggtime::gg_lag()` instead.

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

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

The gg_season 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 gg_subseries 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 the correlation of (y_t) with (y_{t-k}) (after removing the mean). The PACF plot shows significant correlation only for 1 period. This plot removes the preceding 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:

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

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 reciprocal. But that lambda helped none.

Decomposition

We can decompose our data into different components: Seasonal, Trend, and Residual.

• Additive decomposition assumes: \[ y_t = S_t + T_t + R_t \]

• Multiplicative decomposition (which becomes additive in log-scale) assumes: \[ \log(y_t) = S_t + T_t + R_t \]

One of the most commonly used methods is STL (Seasonal-Trend decomposition using LOESS), which estimates the trend and seasonal components via LOESS smoothing.

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

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

Plot variable not specified, automatically selected `.vars = value`

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

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 cannot 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 some 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_1 = (ETS(logvalue) + ARIMA(logvalue)) / 2)

The first model is nothing more than a mean model. We acquire the mean of the TS and use that for forecasting. 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. 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.

The sixth model is an Error, Trend, and Seasonality Model. We allow the software to auto‑choose the terms.

The seventh model is an ARIMA, autoregressive integrated moving average.

Models 8 is a simple ensembles.

Reports

Of course, we want complete reports for the estimates.

for (i in 1:8) {
  fit <- get(paste0("m", i))  # cleaner than eval(parse(...))
  report(fit)
}

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

Plots

Next, we visualize each model’s forecast against the original data. The custom plotting function compares the fitted values and forecasts over the test period.

# Ensure log-transformed values are present
myts$logvalue  <- log(myts$value)
test$logvalue  <- log(test$value)

# Define plotting function
myplot <- function(model) {
  forecast(model, new_data = test) |>  # Fix: forecast should use model, not test
    autoplot(myts) +
    geom_line(aes(y = .fitted, color = .model), data = augment(model)) +
    ggtitle(names(model))
}

# Generate plots for models m1 to m8
for (i in 1:8) {
  print(myplot(get(paste0("m", i))))
}

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

We calculate common forecast accuracy metrics for each model using the test set. Results are stored and compared by MAPE.

# Initialize empty matrix for storing metrics
metrics_all <- NULL

for (i in 1:8) {
  model_forecast <- get(paste0("m", i)) |> forecast(new_data = test)
  model_accuracy <- model_forecast |> accuracy(test)
  metrics_all    <- rbind(metrics_all, model_accuracy)
}

# Clean and rank by MAPE
metrics_all <- metrics_all[metrics_all$RMSE > 0, ]
metrics_all$MASE <- metrics_all$RMSSE <- metrics_all$.type <- NULL

# Display ranked accuracy table
metrics_all |> arrange(MAPE)