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

## Loading required package: carData

library(eia)

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

library(fpp3)

## Warning: package 'fpp3' was built under R version 4.4.3

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

## ── Attaching packages ──────────────────────────────────────────── fpp3 1.0.1 ──

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

## ── 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()
## ✖ dplyr::recode()      masks car::recode()
## ✖ 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()
## ✖ dplyr::recode()          masks car::recode()
## ✖ purrr::set_names()       masks magrittr::set_names()
## ✖ purrr::some()            masks car::some()
## ℹ 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.

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

Time Series Features

So let’s look at some TS features. First, let’s include CPI in our analysis.

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

# Download CPI
cpi <- fredr(
  series_id = "CPIAUCSL",
  observation_start = as.Date("2020-04-01"),
  observation_end = as.Date("2025-03-31")
)

# Clean unemployment data
unemp$Month <- seq(as.Date("2020/4/1"), as.Date("2025/3/31"), by = "1 month")
unemp$Month <- yearmonth(unemp$Month)
unemp <- unemp %>% select(Month, Unemployment = value)

# Clean CPI data
cpi$Month <- seq(as.Date("2020/4/1"), as.Date("2025/3/31"), by = "1 month")
cpi$Month <- yearmonth(cpi$Month)
cpi <- cpi %>% select(Month, CPI = value)

# Merge into a single dataset
mydata <- left_join(unemp, cpi, by = "Month")

# Inspect structure
str(mydata)

## tibble [60 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Month       : mth [1:60] 2020 Apr, 2020 May, 2020 Jun, 2020 Jul, 2020 Aug, 2020 Sep,...
##  $ Unemployment: num [1:60] 14.8 13.2 11 10.2 8.4 7.8 6.9 6.7 6.7 6.4 ...
##  $ CPI         : num [1:60] 256 256 257 258 259 ...

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

# View first rows
print(myts)

## # A tsibble: 60 x 3 [1M]
##       Month Unemployment   CPI
##       <mth>        <dbl> <dbl>
##  1 2020 Apr         14.8  256.
##  2 2020 May         13.2  256.
##  3 2020 Jun         11    257.
##  4 2020 Jul         10.2  258.
##  5 2020 Aug          8.4  259.
##  6 2020 Sep          7.8  260.
##  7 2020 Oct          6.9  260.
##  8 2020 Nov          6.7  261.
##  9 2020 Dec          6.7  262.
## 10 2021 Jan          6.4  263.
## # ℹ 50 more rows

Pivot and Features

So now we pivot longer and generate features (i.e. quantiles).

myts_long <- myts %>%
  pivot_longer(
    cols = c(Unemployment, CPI),
    names_to = "Variable",
    values_to = "Value"
  )
myts_long <- myts_long %>%
  as_tsibble(
    key = Variable,
    index = Month
  )
myfeatures <- myts_long %>%
  features(Value, c(quantile,mean, feat_acf, feat_pacf, feat_stl, coef_hurst, 
                    feat_spectral,box_pierce,ljung_box,unitroot_kpss, unitroot_pp,
                    unitroot_ndiffs,unitroot_nsdiffs, var_tiled_mean, shift_level_max, 
                    shift_var_max, shift_kl_max, n_crossing_points,
                    stat_arch_lm, guerrero))

## New names:
## • `` -> `...6`

myfeatures <- myfeatures %>% ungroup()

temp <- myfeatures %>%
  pivot_longer(
    cols = -Variable,
    names_to = "Code",
    values_to = "Value"
  )


# 3. Pivot wider to get CPI and Unemployment columns
temp <- temp %>%
  distinct() %>%   # <= Important: remove any accidental duplicates
  pivot_wider(
    names_from = Variable,
    values_from = Value
  )


lookup <- tibble(
  Code = c("0%", "25%", "50%", "75%", "100%", "...6", "acf1", "acf10", "diff1_acf1", "diff1_acf10",
           "diff2_acf1", "diff2_acf10", "season_acf1", "pacf5", "diff1_pacf5", "diff2_pacf5",
           "season_pacf", "trend_strength", "seasonal_strength_year", "seasonal_peak_year",
           "seasonal_trough_year", "spikiness", "linearity", "curvature", "stl_e_acf1",
           "stl_e_acf10", "coef_hurst", "spectral_entropy", "bp_stat", "bp_pvalue", "lb_stat",
           "lb_pvalue", "kpss_stat", "kpss_pvalue", "pp_stat", "pp_pvalue", "ndiffs", "nsdiffs",
           "var_tiled_mean", "shift_level_max", "shift_level_index", "shift_var_max",
           "shift_var_index", "shift_kl_max", "shift_kl_index", "n_crossing_points",
           "stat_arch_lm", "lambda_guerrero"),
  Name = c(
    "Minimum", "25th Percentile", "Median", "75th Percentile", "Maximum", "Mean",
    "Lag-1 Autocorrelation", "Sum of First 10 ACFs",
    "ACF of 1st Difference", "Sum ACFs of 1st Diff",
    "ACF of 2nd Difference", "Sum ACFs of 2nd Diff",
    "Seasonal Lag-1 ACF", "PACF Lag-5", "PACF Lag-5 (1st Diff)",
    "PACF Lag-5 (2nd Diff)", "Seasonal PACF", "Trend Strength",
    "Seasonal Strength (Year)", "Seasonal Peak (Month)",
    "Seasonal Trough (Month)", "Spikiness", "Linearity",
    "Curvature", "Residual ACF Lag-1", "Residual ACF Sum (10)",
    "Hurst Exponent", "Spectral Entropy", "Box-Pierce Statistic",
    "Box-Pierce p-value", "Ljung-Box Statistic", "Ljung-Box p-value",
    "KPSS Statistic", "KPSS p-value", "Phillips-Perron Statistic",
    "Phillips-Perron p-value", "Required Differences",
    "Seasonal Differences", "Tiled Mean Variance", "Max Level Shift",
    "Level Shift Index", "Max Variance Shift", "Variance Shift Index",
    "Max KL Divergence", "KL Divergence Index", "Zero Crossing Count",
    "ARCH LM Statistic", "Guerrero Lambda"
  ),
  Description = c(
    "Smallest observed value.",
    "Value below which 25% of observations fall.",
    "Middle value of the distribution.",
    "Value below which 75% of observations fall.",
    "Largest observed value.",
    "Average value of the series.",
    "Correlation between series and itself lagged by 1 period.",
    "Cumulative autocorrelation over first 10 lags.",
    "Lag-1 autocorrelation of first differenced series.",
    "Sum of autocorrelations over first 10 lags (1st differenced).",
    "Lag-1 autocorrelation of second differenced series.",
    "Sum of autocorrelations over first 10 lags (2nd differenced).",
    "Lag-1 autocorrelation at seasonal period.",
    "Partial autocorrelation at lag 5.",
    "Partial autocorrelation at lag 5 (1st differenced).",
    "Partial autocorrelation at lag 5 (2nd differenced).",
    "Partial autocorrelation at seasonal lag.",
    "Proportion of variance explained by trend.",
    "Proportion of variance explained by yearly seasonality.",
    "Month of peak seasonal component.",
    "Month of lowest seasonal component.",
    "Measure of sudden large deviations in the series.",
    "Strength and direction of linear trend.",
    "Degree of nonlinearity in trend.",
    "Lag-1 autocorrelation of STL residuals.",
    "Sum of STL residual autocorrelations over 10 lags.",
    "Degree of long-term memory (0.5=random walk, >0.5 persistence).",
    "Complexity and unpredictability of frequency spectrum.",
    "Chi-squared test statistic for autocorrelation.",
    "p-value for Box-Pierce test.",
    "Adjusted chi-squared test statistic (Ljung-Box).",
    "p-value for Ljung-Box test.",
    "KPSS test statistic for stationarity.",
    "KPSS p-value.",
    "Phillips-Perron test statistic for unit root.",
    "Phillips-Perron p-value.",
    "Minimum number of differences needed for stationarity.",
    "Minimum number of seasonal differences needed.",
    "Variability of the mean across time tiles.",
    "Largest shift in level between adjacent windows.",
    "Index of maximum level shift.",
    "Largest variance change between windows.",
    "Index of maximum variance change.",
    "Largest KL divergence between windows.",
    "Index of maximum KL divergence.",
    "Number of times the series crosses its mean.",
    "ARCH LM test statistic for conditional heteroskedasticity.",
    "Box-Cox lambda estimated by Guerrero method."
  )
)
result_table <- lookup %>%
  left_join(temp, by = "Code") %>%
  select(Name, Description, everything(), -Code)

numeric_cols <- c("CPI", "Unemployment")

# Round and format
result_table <- result_table %>%
  mutate(across(all_of(numeric_cols), ~format(round(.x, 3), nsmall = 3, scientific = FALSE)))
knitr::kable(result_table, align = "l", caption = "Summary of Time Series Features")

Summary of Time Series Features

Name	Description	CPI	Unemployment
Minimum	Smallest observed value.	255.802	3.400
25th Percentile	Value below which 25% of observations fall.	271.651	3.675
Median	Middle value of the distribution.	297.200	4.000
75th Percentile	Value below which 75% of observations fall.	309.000	5.500
Maximum	Largest observed value.	319.775	14.800
Mean	Average value of the series.	291.204	4.967
Lag-1 Autocorrelation	Correlation between series and itself lagged by 1 period.	0.959	0.828
Sum of First 10 ACFs	Cumulative autocorrelation over first 10 lags.	5.811	2.276
ACF of 1st Difference	Lag-1 autocorrelation of first differenced series.	0.368	0.603
Sum ACFs of 1st Diff	Sum of autocorrelations over first 10 lags (1st differenced).	0.528	0.981
ACF of 2nd Difference	Lag-1 autocorrelation of second differenced series.	-0.231	-0.712
Sum ACFs of 2nd Diff	Sum of autocorrelations over first 10 lags (2nd differenced).	0.516	0.951
Seasonal Lag-1 ACF	Lag-1 autocorrelation at seasonal period.	0.426	0.171
PACF Lag-5	Partial autocorrelation at lag 5.	0.925	0.697
PACF Lag-5 (1st Diff)	Partial autocorrelation at lag 5 (1st differenced).	0.403	0.540
PACF Lag-5 (2nd Diff)	Partial autocorrelation at lag 5 (2nd differenced).	0.534	0.639
Seasonal PACF	Partial autocorrelation at seasonal lag.	-0.044	-0.023
Trend Strength	Proportion of variance explained by trend.	0.998	0.953
Seasonal Strength (Year)	Proportion of variance explained by yearly seasonality.	0.000	0.324
Seasonal Peak (Month)	Month of peak seasonal component.	0.000	1.000
Seasonal Trough (Month)	Month of lowest seasonal component.	2.000	8.000
Spikiness	Measure of sudden large deviations in the series.	0.001	0.000
Linearity	Strength and direction of linear trend.	157.158	-12.132
Curvature	Degree of nonlinearity in trend.	-21.860	10.849
Residual ACF Lag-1	Lag-1 autocorrelation of STL residuals.	0.675	0.613
Residual ACF Sum (10)	Sum of STL residual autocorrelations over 10 lags.	0.687	0.771
Hurst Exponent	Degree of long-term memory (0.5=random walk, >0.5 persistence).	0.997	0.995
Spectral Entropy	Complexity and unpredictability of frequency spectrum.	0.226	0.501
Box-Pierce Statistic	Chi-squared test statistic for autocorrelation.	55.129	41.158
Box-Pierce p-value	p-value for Box-Pierce test.	0.000	0.000
Ljung-Box Statistic	Adjusted chi-squared test statistic (Ljung-Box).	57.932	43.250
Ljung-Box p-value	p-value for Ljung-Box test.	0.000	0.000
KPSS Statistic	KPSS test statistic for stationarity.	1.566	0.931
KPSS p-value	KPSS p-value.	0.010	0.010
Phillips-Perron Statistic	Phillips-Perron test statistic for unit root.	-1.564	-12.743
Phillips-Perron p-value	Phillips-Perron p-value.	0.100	0.010
Required Differences	Minimum number of differences needed for stationarity.	1.000	2.000
Seasonal Differences	Minimum number of seasonal differences needed.	0.000	0.000
Tiled Mean Variance	Variability of the mean across time tiles.	1.185	0.833
Max Level Shift	Largest shift in level between adjacent windows.	21.721	6.825
Level Shift Index	Index of maximum level shift.	21.000	2.000
Max Variance Shift	Largest variance change between windows.	49.840	8.224
Variance Shift Index	Index of maximum variance change.	29.000	9.000
Max KL Divergence	Largest KL divergence between windows.	0.231	11.536
KL Divergence Index	Index of maximum KL divergence.	18.000	4.000
Zero Crossing Count	Number of times the series crosses its mean.	1.000	4.000
ARCH LM Statistic	ARCH LM test statistic for conditional heteroskedasticity.	0.996	0.865
Guerrero Lambda	Box-Cox lambda estimated by Guerrero method.	2.000	-0.900

We readily generate the 5-number summary here along with the mean. We also look at a variety of TS features.

Time Series Diagnostic Statistics: Definitions, Formulae, and Interpretations

This appendix provides the formulae and interpretation for selected metrics.

---

Autocorrelation and Partial Autocorrelation

Lag-1 Autocorrelation

Formula:

\[ \rho_1 = \frac{\sum_{t=2}^{n} (x_t - \bar{x})(x_{t-1} - \bar{x})}{\sum_{t=1}^{n} (x_t - \bar{x})^2} \]

Interpretation: Measures the correlation of the series with itself lagged by 1 period. High values indicate strong persistence.

---

Sum of First 10 ACFs

Formula:

\[ \sum_{k=1}^{10} \rho_k \]

Interpretation: Cumulative autocorrelation over the first 10 lags. Higher sums imply strong persistence.

---

ACF of 1st Difference

Formula: Same as Lag-1 autocorrelation but applied to the differenced series:

\[ \rho_1(\Delta x_t) \]

Interpretation: Autocorrelation after first differencing. Lower values indicate stationarity achieved by differencing.

---

Sum ACFs of 1st Diff

Formula:

\[ \sum_{k=1}^{10} \rho_k(\Delta x_t) \]

Interpretation: Aggregate autocorrelation of the first differenced series.

---

ACF of 2nd Difference

Formula:

\[ \rho_1(\Delta^2 x_t) \]

Interpretation: Lag-1 autocorrelation after second differencing.

---

Sum ACFs of 2nd Diff

Formula:

\[ \sum_{k=1}^{10} \rho_k(\Delta^2 x_t) \]

Interpretation: Cumulative autocorrelation of the twice differenced series.

---

Seasonal Lag-1 ACF

Formula:

\[ \rho_s \]

where \(s\) = seasonal lag (e.g., 12).

Interpretation: Measures seasonality correlation.

---

PACF Lag-5

Formula: Partial autocorrelation at lag 5, defined recursively via the Yule-Walker equations.

Interpretation: Correlation at lag 5 controlling for lags 1–4.

---

Trend and Seasonality Strength

Trend Strength

Formula: Proportion of variance explained by the trend component:

\[ 1 - \frac{\text{Var}(R)}{\text{Var}(X)} \]

where \(R\) = remainder (residual).

Interpretation: Near 1 means almost all variance comes from the trend.

---

Seasonal Strength

Formula: Proportion of variance explained by seasonality:

\[ 1 - \frac{\text{Var}(R + T)}{\text{Var}(X)} \]

where \(T\) = trend.

Interpretation: Higher values mean stronger seasonality.

---

Spikiness

Definition:
Spikiness measures how much the local variance of residuals fluctuates over time. It helps detect sudden bursts of volatility (e.g., outliers, volatility clustering).

Formula:

Let residuals be: \[ e_1, e_2, \dots, e_n. \]

Compute the leave-one-out variance for each observation \(i\): \[ v_i = \frac{1}{n - 2}\sum_{j \neq i}\left(e_j - \bar{e}_{(-i)}\right)^2, \] where \(\bar{e}_{(-i)}\) is the mean of all residuals excluding observation \(i\).

The spikiness is: \[ \text{Spikiness} = \text{Var}\left\{v_1, v_2, \dots, v_n\right\}. \]

Interpretation:

• Low spikiness means residuals have stable variance.
• High spikiness indicates periods of volatility shifts or sudden outliers.

---

Linearity and Curvature

Linearity

Formula: Estimated slope coefficient in a linear model:

\[ y_t = \beta_0 + \beta_1 t + \varepsilon_t \]

Interpretation: Positive indicates upward trend; negative indicates downward.

---

Curvature

Formula: Quadratic coefficient from:

\[ y_t = \beta_0 + \beta_1 t + \beta_2 t^2 + \varepsilon_t \]

Interpretation: Curvature in trend (concave or convex).

---

Residual ACF

Residual ACF Lag-1

Formula: Lag-1 autocorrelation of residuals:

\[ \rho_1(\hat{\varepsilon}) \]

Interpretation: Should be near zero if model is adequate.

---

Residual ACF Sum (10)

Formula: Cumulative ACF of residuals over 10 lags:

\[ \sum_{k=1}^{10} \rho_k(\hat{\varepsilon}) \]

---

Hurst Exponent

Hurst Exponent

Interpretation: Measures long-term memory.

• \(H=0.5\): random walk.
• \(H >0.5\): persistence.
• \(H<0.5\): mean reversion.

---

Spectral Entropy

Spectral Entropy

Formula: Entropy of normalized power spectrum:

\[ H = -\sum_{i} p_i \log(p_i) \]

Interpretation: Lower entropy = more predictable series.

---

Box-Pierce and Ljung-Box

Box-Pierce Statistic

Formula:

\[ Q = n \sum_{k=1}^{h} \hat{\rho}_k^2 \]

Interpretation: Large \(Q\) indicates autocorrelation.

---

Ljung-Box Statistic

Formula:

\[ Q = n(n+2)\sum_{k=1}^{h} \frac{\hat{\rho}_k^2}{n-k} \]

Interpretation: Adjusted for small samples.

---

KPSS Test

KPSS Statistic

Formula: Sum of squared partial sums of residuals divided by variance estimate.

Interpretation: Higher values indicate nonstationarity.

---

Phillips-Perron

Phillips-Perron Test

Interpretation: Tests for unit root. Negative statistics with low p-values reject unit root hypothesis.

---

Required Differences

Required Differences

Interpretation: Minimum differencing needed to achieve stationarity.

---

Guerrero Lambda

Guerrero Lambda

Interpretation: Box-Cox transformation parameter estimate.

# Time Series Diagnostic Statistics: Definitions, Formulae, and Interpretations

This appendix describes each diagnostic statistic, shows formulas where applicable, and explains how to interpret the results.

Autocorrelation and Partial Autocorrelation

Lag-1 Autocorrelation

Formula: \(\rho_1 = \frac{\sum_{t=2}^{n} (x_t - \bar{x})(x_{t-1} - \bar{x})}{\sum_{t=1}^{n} (x_t - \bar{x})^2}\)

Interpretation: - Measures persistence from one observation to the next. - CPI: 0.959 (very strong persistence) - Unemployment: 0.828 (strong persistence)

---

Sum of First 10 ACFs

Formula: \(\sum_{k=1}^{10} \rho_k\)

Interpretation: - Captures cumulative autocorrelation. - CPI: 5.811 (high persistence) - Unemployment: 2.276 (moderate)

---

ACF of 1st Difference

Interpretation: - Autocorrelation after differencing once. - Lower is better if goal is stationarity. - CPI: 0.368 (residual correlation) - Unemployment: 0.603 (notable persistence)

---

Sum ACFs of 1st Diff

Interpretation: - Aggregated persistence post-differencing. - CPI: 0.528 - Unemployment: 0.981

---

ACF of 2nd Difference

Interpretation: - Autocorrelation after differencing twice. - Negative values may indicate over-differencing. - CPI: -0.231 - Unemployment: -0.712

---

Sum ACFs of 2nd Diff

Interpretation: - CPI: 0.516 - Unemployment: 0.951

---

Seasonal Lag-1 ACF

Interpretation: - Autocorrelation at seasonal lag (e.g., lag 12). - CPI: 0.426 (moderate seasonality) - Unemployment: 0.171 (low)

---

PACF Lag-5

Interpretation: - Direct dependence at lag 5. - CPI: 0.925 - Unemployment: 0.697

---

PACF Lag-5 (1st Diff)

Interpretation: - CPI: 0.403 - Unemployment: 0.540

---

PACF Lag-5 (2nd Diff)

Interpretation: - CPI: 0.534 - Unemployment: 0.639

---

Seasonal PACF

Interpretation: - CPI: -0.044 - Unemployment: -0.023

---

Trend and Seasonality Strength

Trend Strength

Interpretation: - Proportion of variance explained by trend. - Close to 1 means very strong trend. - CPI: 0.998 - Unemployment: 0.953

---

Seasonal Strength (Year)

Interpretation: - Proportion explained by yearly seasonality. - CPI: 0.000 - Unemployment: 0.324

---

Seasonal Peak (Month)

Interpretation: - Month of peak seasonal component (0–11). - CPI: 0 (January) - Unemployment: 1 (February)

---

Seasonal Trough (Month)

Interpretation: - Month of lowest seasonal component. - CPI: 2 (March) - Unemployment: 8 (September)

---

Spikiness, Linearity, and Curvature

Spikiness

Interpretation: - Measures sudden large deviations. - Close to 0 means smooth. - CPI: 0.001 - Unemployment: 0.000

---

Linearity

Interpretation: - Positive indicates upward trend. - Negative indicates downward. - CPI: 157.158 (very strong upward trend) - Unemployment: -12.132 (weak negative trend)

---

Curvature

Interpretation: - Non-linear trends. - Positive: convex; Negative: concave. - CPI: -21.860 - Unemployment: 10.849

---

Residual Diagnostics

Residual ACF Lag-1

Interpretation: - Should be near 0 for good model fit. - CPI: 0.675 - Unemployment: 0.613

---

Residual ACF Sum (10)

Interpretation: - Sum over 10 lags. - CPI: 0.687 - Unemployment: 0.771

---

Memory and Complexity

Hurst Exponent

Interpretation: - \(H = 0.5\): random walk. - \(H > 0.5\): long memory. - CPI: 0.997 - Unemployment: 0.995

---

Spectral Entropy

Interpretation: - Lower entropy: more predictable. - CPI: 0.226 (predictable) - Unemployment: 0.501 (less predictable)

---

Box-Pierce and Ljung-Box Tests

Box-Pierce Statistic

Interpretation: - Large values = autocorrelation present. - CPI: 55.129 - Unemployment: 41.158

---

Box-Pierce p-value

Interpretation: - p < 0.05: reject null of no autocorrelation. - Both p ≈ 0.

---

Ljung-Box Statistic

Interpretation: - Adjusted test for autocorrelation. - CPI: 57.932 - Unemployment: 43.250

---

Ljung-Box p-value

Interpretation: - Both p ≈ 0: autocorrelation detected.

---

Stationarity Tests

KPSS Statistic

Interpretation: - High: non-stationarity. - CPI: 1.566 - Unemployment: 0.931

---

KPSS p-value

Interpretation: - p = 0.01: evidence of non-stationarity.

---

Phillips-Perron Statistic

Interpretation: - Negative: reject unit root. - CPI: -1.564 - Unemployment: -12.743

---

Phillips-Perron p-value

Interpretation: - CPI: 0.10 (weak evidence of stationarity) - Unemployment: 0.01 (stationary)

---

Required Differences

Interpretation: - Min differences to stationarity. - CPI: 1 - Unemployment: 2

---

Seasonal Differences

Interpretation: - Seasonal differencing needed. - Both: 0

---

Change Point Statistics

Tiled Mean Variance

Interpretation: - Variability over time tiles. - CPI: 1.185 - Unemployment: 0.833

---

Max Level Shift

Interpretation: - Largest shift in level. - CPI: 21.721 - Unemployment: 6.825

---

Level Shift Index

Interpretation: - Index of max shift. - CPI: 21 - Unemployment: 2

---

Max Variance Shift

Interpretation: - Largest change in variance. - CPI: 49.840 - Unemployment: 8.224

---

Variance Shift Index

Interpretation: - Index of change. - CPI: 29 - Unemployment: 9

---

Max KL Divergence

Interpretation: - Largest divergence between windows. - CPI: 0.231 - Unemployment: 11.536

---

KL Divergence Index

Interpretation: - Index of divergence. - CPI: 18 - Unemployment: 4

---

Zero Crossing Count

Interpretation: - Number of mean crossings. - CPI: 1 - Unemployment: 4

---

Volatility and Transformation

ARCH LM Statistic

Interpretation: - Test for conditional heteroskedasticity. - CPI: 0.996 - Unemployment: 0.865

---

Guerrero Lambda

Interpretation: - Box-Cox transformation parameter. - CPI: 2.000 - Unemployment: -0.900

---

---

Note: Metrics like Max Level Shift, Variance Shift, and KL Divergence are calculated numerically, so no standard formula is shown.

Forecaster’s Toolbox

We can use a simple algorithm to prepare data for forecasting. First, we tidy the data. Then we visualize, estimate, and generate forecasts which are evaluated. The visualization, estimation, and forecast generation process is repeated. Hyndman illustrates.

Process

Let’s do that with our Unemployment and CPI data. We will create a proxy for the misery index.

Tidy the Data

myts <- myts %>%
  arrange(Month) %>%
  mutate(
    CPI_inflation_mom = 100 * (CPI / lag(CPI, 1) - 1)
  )

Plot the Data

myts$Misery=myts$Unemployment+myts$CPI_inflation_mom
myts |>
  autoplot(Misery) +
  labs(y = "Misery Index", title = "Misery Index")

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

myts%>%gg_season(Misery)

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

myts%>%gg_subseries(Misery)

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

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

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

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

## Guerrero

# Compute lambda using Guerrero method
lambda <- myts %>%filter(!is.na(Misery), Misery > 0) %>%
  features(Misery, features = guerrero) %>%
  pull(lambda_guerrero)

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

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

Nearly inverse! But we will ignore for now (and for interpretability).

Estimate the Models

Train and Test Set

train2 <- myts %>%
  slice(2:48)

test2 <- myts %>%
  slice(49:n())

Models

#Log
train2$logMisery=log(train2$Misery)

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

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

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

Evaluate

Plot

myts$logMisery=log(myts$Misery)
test2$logMisery=log(test2$Misery)

myplot=function(model){
  model|>forecast(test2)|>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()`).
## 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()`).
## 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()`).

## 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 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 1 row containing missing values or values outside the scale range
## (`geom_line()`).

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 = test2)
  
  # 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 Ensemble_Model       -0.0160  -1.13  0.0334 0.0246  1.71 0.0880 0.101
## 2 Naive_Model           0.0110   0.706 0.0350 0.0317  2.17 0.114  0.106
## 3 ETS_Model             0.0265   1.77  0.0433 0.0378  2.57 0.135  0.131
## 4 ARIMA_Model          -0.0584  -4.04  0.0687 0.0584  4.04 0.209  0.208
## 5 Ensemble_Model_2     -0.0671  -4.67  0.0821 0.0671  4.67 0.240  0.248
## 6 Seasonal_Naive_Model  0.0744   5.07  0.0920 0.0784  5.35 0.281  0.278
## 7 Drift_Model           0.170   11.6   0.200  0.171  11.6  0.611  0.604
## 8 TS_Model             -0.169  -11.6   0.172  0.169  11.6  0.604  0.520
## 9 Mean_Model           -0.169  -11.6   0.172  0.169  11.6  0.604  0.520