1 Preprocessing: Converting Weekly to Daily then Monthly

When working with time-series data, such as weekly flu rates, directly aggregating to monthly values can introduce bias. Why? Weeks often straddle months, and blindly assigning an entire week to a single month over-or under-weights certain periods.

To avoid this, we use a two-step transformation.
1.Disaggregate weekly data into daily values Each weekly rate is spread evenly across its 7 days. This assumes the weekly rate applies uniformly across the week—reasonable when daily resolution is unavailable.

2.Aggregate daily values to monthly averages (or totals) By averaging the daily rates within each calendar month, we ensure each day contributes proportionally. This avoids artifacts from week boundaries and reflects the true monthly behavior.

This approach maintains temporal fidelity, enables accurate seasonal modeling, and aligns with best practices in the fpp3 framework.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tsibble)
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## 
## The following object is masked from 'package:lubridate':
## 
##     interval
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(lubridate)
library(readxl)
library(fable)
## Loading required package: fabletools
library(feasts)
library(fabletools)
library(purrr)

# Step 1: Load & prep weekly data
#C:\Users\lfult\OneDrive - bc.edu\_Courses\Boston College\Predictive analytics
flu_data <- read_excel("c:/users/lfult/Onedrive - bc.edu/_Courses/Boston College/Predictive Analytics/example for class.xlsx", sheet = "Sheet1") %>%
  filter(Measure == "Influenza Hospital Admissions Rate Per 100k") %>%
  mutate(
    week_start = ISOweek::ISOweek2date(paste0(Year, "-W", str_pad(Week, 2, pad = "0"), "-1"))
  ) %>%
  select(week_start, flu_rate = Measure_Value)

# Step 2: Expand to daily, convert to tsibble
flu_daily <- flu_data %>%
  rowwise() %>%
  mutate(
    date = list(seq(week_start, by = "1 day", length.out = 7)),
    daily_rate = list(rep(flu_rate, 7))
  ) %>%
  unnest(c(date, daily_rate)) %>%
  as_tsibble(index = date)

# Step 3: Aggregate daily to monthly average
flu_monthly <- flu_daily %>%
  index_by(month = yearmonth(date)) %>%
  summarise(monthly_avg_rate = mean(daily_rate, na.rm = TRUE)) %>%
  filter(!is.na(monthly_avg_rate))

2 Train / Test Split and Plot

A train/test split is essential for evaluating how well a time-series model can generalize to future data. We train models on historical data and then test their forecasts on unseen observationd This simulates real-world forecasting and helps prevent overfitting.

Plotting both the train/test data alongside model forecasts visually answers the following questions.

1. Does he model captures trends and seasonality?

2. Does the model perform well on unseen data?

3. Where do the model predictions fail?

This visual validation is critical before trusting the model for decisions or further analysis.

# Step 4: Train-test split
train_data <- flu_monthly %>%
  filter(month <= yearmonth("2024 May"))

test_data <- flu_monthly %>%
  filter(month > yearmonth("2024 May"))

# Step 5: Exploratory plots
# --- Seasonality
gg_season(train_data, monthly_avg_rate) +
  labs(title = "Seasonal Plot: Flu Rates by Month")

# --- Subseries (month-wise trend)
gg_subseries(train_data, monthly_avg_rate) +
  labs(title = "Subseries Plot: Flu Rates by Month")

# --- Lag Plot
gg_lag(train_data, monthly_avg_rate, geom = "point") +
  labs(title = "Lag Plot of Monthly Flu Rates")

# --- ACF & PACF
train_data %>%
  ACF(monthly_avg_rate) %>%
  autoplot() +
  labs(title = "ACF of Monthly Flu Rates")

train_data %>%
  PACF(monthly_avg_rate) %>%
  autoplot() +
  labs(title = "PACF of Monthly Flu Rates")

# --- Full display of time series (trend, ACF, seasonality)
train_data %>%
  gg_tsdisplay(monthly_avg_rate, plot_type = "partial") +
  labs(title = "TS Display: Trend, Seasonality, ACF")

# --- Plot full series with test set as dashed
flu_monthly %>%
  mutate(Source = if_else(month <= yearmonth("2024 May"), "Train", "Test")) %>%
  ggplot(aes(x = month, y = monthly_avg_rate, color = Source)) +
  geom_line() +
  labs(title = "Monthly Flu Rate with Train/Test Split",
       x = "Month", y = "Rate per 100k") +
  scale_color_manual(values = c("Train" = "blue", "Test" = "red")) +
  theme_minimal()

These six graphs each represent flu forecast models with varying predictive intervals and smoothing characteristics. The forecasts are shown with ribbons indicating uncertainty bounds, and they generally track the seasonal surge of flu rates during the winter months. In several of the plots, we observe a well-calibrated alignment between the forecast mean and observed values early on, but widening uncertainty further into the forecast horizon, a typical signal of increasing variance over time. Some models (e.g., in the second and fifth images) show sharper prediction intervals, suggesting tighter confidence around the central tendency, though this may come at the cost of underestimating variance in peak flu seasons. Others (such as the fourth and sixth plots) exhibit broader intervals, capturing more of the volatility but risking lower specificity. Overall, the model forecasts are reasonable, with varying degrees of sharpness and calibration, reflecting the trade-off between confidence and coverage.

3 Models

ETS (Error, Trend, Seasonal) models describe time series using combinations of:

Error type: Additive (A) or Multiplicative (M) or None (‘N’) Trend: None (N), Additive (A), or Multiplicative (M) Seasonality: None (N), Additive (A), or Multiplicative (M)

This yields a 3 × 3 grid of 9 common ETS models.

ETS(A,N,N) – Simple exponential smoothing ’ETS(M,A,A)* – Multiplicative Holt-Winters with additive trend *ETS(A,A,A)`() – Fully additive model with seasonal pattern

Each structure fits different types of time series:

Additive models assume constant variation over time. Multiplicative models capture series where variation grows with the level (e.g., flu rates that spike higher in bad seasons). Trend and seasonal components model structural behaviors, such as yearly flu cycles or persistent increases/decreases.

*How fpp3 Chooses Automatically

In the fpp3 framework, you can let ETS() select the best model using ETS(monthly_avg_rate ~ error(“ZZZ”)) The “ZZZ” tells fpp3 to automatically search across all valid combinations of error, trend, and seasonality. Additionally, ffp3 will fit them using maximum likelihood and rank them using AICc (Akaike Information Criterion with correction). The result is an automatically selected model that balances fit and complexity—giving you a strong, data-driven starting point for forecasting.

# Step 1: Fit all possible ETS models on the training data
ets_models <- train_data %>%
  model(
    ETS_ANN = ETS(monthly_avg_rate ~ error("A") + trend("N") + season("N")),
    ETS_AAN = ETS(monthly_avg_rate ~ error("A") + trend("A") + season("N")),
    ETS_AAA = ETS(monthly_avg_rate ~ error("A") + trend("A") + season("A")),
    ETS_AMN = ETS(monthly_avg_rate ~ error("A") + trend("M") + season("N")),
    ETS_AMA = ETS(monthly_avg_rate ~ error("A") + trend("M") + season("A")),
    ETS_MNN = ETS(monthly_avg_rate ~ error("M") + trend("N") + season("N")),
    ETS_MAN = ETS(monthly_avg_rate ~ error("M") + trend("A") + season("N")),
    ETS_MAA = ETS(monthly_avg_rate ~ error("M") + trend("A") + season("A")),
    ETS_MAM = ETS(monthly_avg_rate ~ error("M") + trend("A") + season("M"))
  )

# Step 2: Generate forecasts over the test period
library(fabletools)

ets_fc <- forecast(ets_models, h = nrow(test_data)) %>%
  hilo() %>%
  unpack_hilo(`95%`)

ets_fc2=forecast(ets_models, h = nrow(test_data), level = c(80, 95))

4 Model Reports

Of course, running the models doesn’t help us unless we run the reports!

for (name in names(ets_models)) {
  cat("\n\n### Report for", name, "\n\n")
  print(report(ets_models[[name]][[1]]))
}
## 
## 
## ### Report for ETS_ANN 
## 
## Series: monthly_avg_rate 
## Model: ETS(A,N,N) 
##   Smoothing parameters:
##     alpha = 0.9998999 
## 
##   Initial states:
##          l[0]
##  0.0007126544
## 
##   sigma^2:  1.1826
## 
##      AIC     AICc      BIC 
## 182.8007 183.3861 188.2207 
##      AIC     AICc      BIC 
## 182.8007 183.3861 188.2207 
## 
## 
## ### Report for ETS_AAN 
## 
## Series: monthly_avg_rate 
## Model: ETS(A,A,N) 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.0001000259 
## 
##   Initial states:
##       l[0]        b[0]
##  -1.397439 -0.02347322
## 
##   sigma^2:  1.2905
## 
##      AIC     AICc      BIC 
## 188.5871 190.1255 197.6204 
##      AIC     AICc      BIC 
## 188.5871 190.1255 197.6204 
## 
## 
## ### Report for ETS_AAA 
## 
## Series: monthly_avg_rate 
## Model: ETS(A,A,A) 
##   Smoothing parameters:
##     alpha = 0.0001079829 
##     beta  = 0.000103126 
##     gamma = 0.0001103105 
## 
##   Initial states:
##        l[0]       b[0]       s[0]      s[-1]      s[-2]      s[-3]      s[-4]
##  -0.1778538 0.04352008 -0.7462467 -0.5874411 -0.1470452 -0.3342339 -0.5820749
##        s[-5]     s[-6]     s[-7]    s[-8]     s[-9]    s[-10]     s[-11]
##  -0.03707391 0.1346033 0.5485954 2.459856 0.3576924 -0.466304 -0.6003276
## 
##   sigma^2:  1.4891
## 
##      AIC     AICc      BIC 
## 203.4469 226.1136 234.1602 
##      AIC     AICc      BIC 
## 203.4469 226.1136 234.1602 
## 
## 
## ### Report for ETS_AMN 
## 
## Series: monthly_avg_rate 
## Model: ETS(A,M,N) 
##   Smoothing parameters:
##     alpha = 0.9996546 
##     beta  = 0.0001000001 
## 
##   Initial states:
##      l[0]      b[0]
##  0.217303 0.7082392
## 
##   sigma^2:  1.1041
## 
##      AIC     AICc      BIC 
## 181.5678 183.1063 190.6012 
##      AIC     AICc      BIC 
## 181.5678 183.1063 190.6012 
## 
## 
## ### Report for ETS_AMA 
## 
## Series: monthly_avg_rate 
## Model: ETS(A,M,A) 
##   Smoothing parameters:
##     alpha = 0.08077987 
##     beta  = 0.0007414703 
##     gamma = 0.0001187467 
## 
##   Initial states:
##       l[0]      b[0]       s[0]      s[-1]      s[-2]      s[-3]      s[-4]
##  0.7817316 0.7989138 -0.6295529 -0.7474424 -0.1808402 -0.1300691 -0.3344528
##       s[-5]      s[-6]     s[-7]    s[-8]     s[-9]     s[-10]     s[-11]
##  -0.1980208 -0.1906247 0.6406311 2.639687 0.6849273 -0.5544393 -0.9998035
## 
##   sigma^2:  1.8162
## 
##      AIC     AICc      BIC 
## 212.3825 235.0492 243.0958 
##      AIC     AICc      BIC 
## 212.3825 235.0492 243.0958 
## 
## 
## ### Report for ETS_MNN 
## 
## Series: monthly_avg_rate 
## Model: ETS(M,N,N) 
##   Smoothing parameters:
##     alpha = 0.9716986 
## 
##   Initial states:
##       l[0]
##  0.5927189
## 
##   sigma^2:  1.4569
## 
##      AIC     AICc      BIC 
## 93.73594 94.32131 99.15593 
##      AIC     AICc      BIC 
## 93.73594 94.32131 99.15593 
## 
## 
## ### Report for ETS_MAN 
## 
## Series: monthly_avg_rate 
## Model: ETS(M,A,N) 
##   Smoothing parameters:
##     alpha = 0.03630488 
##     beta  = 0.03630486 
## 
##   Initial states:
##       l[0]        b[0]
##  0.1127431 -0.01274659
## 
##   sigma^2:  3.2893
## 
##      AIC     AICc      BIC 
## 137.4494 138.9879 146.4827 
##      AIC     AICc      BIC 
## 137.4494 138.9879 146.4827 
## 
## 
## ### Report for ETS_MAA 
## 
## Series: monthly_avg_rate 
## Model: ETS(M,A,A) 
##   Smoothing parameters:
##     alpha = 0.01608937 
##     beta  = 0.002170562 
##     gamma = 0.05280681 
## 
##   Initial states:
##        l[0]       b[0]       s[0]      s[-1]      s[-2]      s[-3]      s[-4]
##  -0.1683132 0.04674163 -0.7081568 -0.5958517 -0.3672397 -0.1744883 -0.1111187
##       s[-5]      s[-6]     s[-7]    s[-8]     s[-9]     s[-10]     s[-11]
##  -0.1873498 -0.3270449 0.4971496 2.662934 0.6446474 -0.5749135 -0.7585675
## 
##   sigma^2:  1.4521
## 
##      AIC     AICc      BIC 
## 158.2073 180.8740 188.9206 
##      AIC     AICc      BIC 
## 158.2073 180.8740 188.9206 
## 
## 
## ### Report for ETS_MAM 
## 
## Series: monthly_avg_rate 
## Model: ETS(M,A,M) 
##   Smoothing parameters:
##     alpha = 0.2172262 
##     beta  = 0.001824539 
##     gamma = 0.03205124 
## 
##   Initial states:
##         l[0]       b[0]      s[0]     s[-1]     s[-2]     s[-3]    s[-4]
##  -0.03395065 0.03395065 0.2677047 0.3288935 0.4772897 0.8835812 1.246323
##     s[-5]     s[-6]    s[-7]    s[-8]    s[-9]    s[-10]    s[-11]
##  1.128106 0.7775542 1.509588 3.353253 1.432243 0.4190449 0.1764187
## 
##   sigma^2:  0.7427
## 
##      AIC     AICc      BIC 
## 48.72180 71.38847 79.43506 
##      AIC     AICc      BIC 
## 48.72180 71.38847 79.43506

Among the ETS family of models fitted to the monthly_avg_rate series, several models offered a good in-sample fit, with variations in structural complexity and implied component behavior (trend, seasonality, and error form).

The ETS(M,A,M) model achieved the lowest AIC (48.7), AICc (71.4), and BIC (79.4), with a relatively low estimated variance (σ² = 0.7427), indicating a strong in-sample fit. It includes multiplicative error, additive trend, and multiplicative seasonality, suggesting the data may exhibit percentage-based fluctuations and a trend that grows or declines at a constant rate. However, the inclusion of both trend and seasonality increases model complexity.

By contrast, simpler models like ETS(ANN) and ETS(AMN) also performed well (AIC ≈ 182), indicating good fits with fewer components. Their comparable training performance suggests that gains from including trend or seasonality may be modest, especially if overfitting is a concern.

Some other models, such as ETS(AAA) and ETS(AMA), had much higher AIC and variance, suggesting poorer fits on the training data, possibly due to overparameterization or poor alignment with the underlying signal.

While ETS(M,A,M) appears to be the most promising model based on training fit, this conclusion should be hedged until evaluated on validation data. Overfitting is a risk—especially with more complex seasonal and multiplicative models—so cross-validation or out-of-sample testing is needed to assess robustness.

5 Accuracy Comparison

The accuracy() function in fpp3 computes forecast accuracy metrics by comparing model predictions to actual values in the test set. This allows us to rank and evaluate models using multiple metrics.

The key metrics reported include:

RMSE (Root Mean Squared Error):
Penalizes larger errors more heavily. Useful for general error magnitude.

MAE (Mean Absolute Error):
Simple average of absolute forecast errors. Less sensitive to outliers than RMSE.

MAPE (Mean Absolute Percentage Error):
Expresses forecast error as a percentage. Can be unstable if actual values are close to zero.

MASE (Mean Absolute Scaled Error):
Scales MAE relative to a naïve seasonal forecast. Values > 1 indicate worse-than-naïve performance.
Note: Can return NA if the training data has no variation (e.g., all values are flat), making the denominator zero.

RMSSE (Root Mean Squared Scaled Error):
Similar to MASE, but based on RMSE instead of MAE.
Also requires a non-zero benchmark error from the training set, or it returns NA.

5.1 Why MASE and RMSSE Might Be NA

MASE and RMSSE compare model errors to a naïve seasonal benchmark, typically using a seasonal naïve forecast on the training set. If the training data has very low variation, contains many flat or repeated values, or lacks sufficient seasonal structure, then the benchmark’s error (denominator) may be zero or undefined, leading to NA values in MASE and RMSSE.

This is especially common in public health time series like flu rates, where entire seasons (e.g., summers) may have consistently zero or near-zero values.

5.2 Recommendation

When MASE and RMSSE are NA, rely instead on RMSE and MAE for error magnitude, visual diagnostics (residual plots, forecast overlays), and comparisons across multiple metrics, not just one.

ets_accuracy <- ets_fc2 %>%
  accuracy(test_data) %>%
  arrange(RMSE)

ets_accuracy
## # A tibble: 9 × 10
##   .model  .type       ME     RMSE      MAE         MPE    MAPE  MASE RMSSE  ACF1
##   <chr>   <chr>    <dbl>    <dbl>    <dbl>       <dbl>   <dbl> <dbl> <dbl> <dbl>
## 1 ETS_AMA Test  -3.83e-5     3.90     2.89    -6.62e+2 6.89 e2   NaN   NaN 0.599
## 2 ETS_AAA Test   6.65e-1     3.91     2.52    -5.15e+2 5.46 e2   NaN   NaN 0.615
## 3 ETS_MAM Test   1.08e+0     3.96     2.20    -1.72e+2 2.04 e2   NaN   NaN 0.522
## 4 ETS_MAA Test   6.47e-1     3.97     2.57    -5.07e+2 5.38 e2   NaN   NaN 0.607
## 5 ETS_MAN Test   1.58e+0     4.47     2.68    -3.93e+2 4.43 e2   NaN   NaN 0.691
## 6 ETS_MNN Test   2.45e+0     4.83     2.56    -1.03e+1 1.000e2   NaN   NaN 0.689
## 7 ETS_ANN Test   2.47e+0     4.84     2.56    -8.71e-1 9.39 e1   NaN   NaN 0.689
## 8 ETS_AAN Test   2.65e+0     4.95     2.70     4.29e+1 9.57 e1   NaN   NaN 0.692
## 9 ETS_AMN Test  -1.27e+4 41138.   12672.      -4.07e+6 4.08 e6   NaN   NaN 0.133

Several models struggled with high percentage errors—MAPE exceeding 500% in most cases—despite relatively low RMSE and MAE. This suggests volatility in the scale of the response variable, likely due to small or near-zero values in parts of the series. A log or Box-Cox transformation could help stabilize variance and reduce the distortion of percentage-based errors. However, applying such transformations requires care, especially if zeros are present, as log-transforms are undefined for non-positive values. The particularly poor performance of ETS(AMN) may also be due to the multiplicative error structure interacting poorly with untransformed, noisy data. These results indicate that without a transformation to regularize scale and variance, model performance may remain unstable, particularly in out-of-sample forecasts.

6 Forecast Plots

Forecast plots are a vital part of time series analysis because they reveal how a model behaves in time, beyond what numerical accuracy metrics alone can show.

6.0.0.1 1. Visualizing Fit and Forecast Behavior

Plots let us see how well a model tracks the historical data and how it projects into the future. We can visually assess overfitting or underfitting, Whether trends or seasonality are captured, and sudden shifts or poor adaptation in the forecast

6.0.0.2 2. Understanding Forecast Uncertainty

By showing prediction intervals (e.g., 80% or 95%), forecast plots help us understand how much confidence to place in the forecast, whether the uncertainty grows over time (as it often should), and if forecasts are overly narrow or unrealistically wide

6.0.0.3 3. Communicating Results Effectively

Plots are intuitive. Stakeholders and students may not interpret RMSE or MASE easily, but a well-labeled plot makes it immediately clear whether a model’s predictions make sense.

6.0.0.4 4. Detecting Model Pathology

Outliers, bias, delayed response, or unrealistic extrapolations can often be spotted at a glance in a forecast plot—things metrics might miss.

In short, forecast plots bring transparency, intuition, and accountability to your modeling. They’re not just for aesthetics—they’re essential for evaluation.

plot_data <- flu_monthly %>%
  mutate(Source = case_when(
    month <= yearmonth("2024 May") ~ "Train",
    TRUE ~ "Test"
  ))

ggplot() +
  geom_ribbon(data = ets_fc,
              aes(x = month, ymin = `95%_lower`, ymax = `95%_upper`, fill = "95% PI"),
              alpha = 0.2, inherit.aes = FALSE) +

  geom_line(data = plot_data, aes(x = month, y = monthly_avg_rate, color = Source), size = 1) +

  geom_line(data = plot_data,
            aes(x = month, y = monthly_avg_rate, color = Source), size = 1) +
  
  geom_line(data = ets_fc2, aes(x = month, y = .mean, color = "Forecast"), size = 1)+

  scale_color_manual(values = c("Train" = "blue", "Test" = "black", "Forecast" = "red")) +
  scale_fill_manual(values = c("95% PI" = "pink")) +

  facet_wrap(~ .model, ncol = 3, scales = "free_y") +

  labs(title = "ETS Forecasts vs Actual Flu Rates",
       subtitle = "Training (blue), Test (black), Forecast (red), 95% PI (pink)",
       x = "Month", y = "Flu Rate per 100k",
       color = "Series", fill = "Prediction Interval") +

  theme_minimal(base_size = 13) +
  theme(
  legend.position = "bottom",
  axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
)+
  theme(legend.position = "bottom")
## 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.

What’s clear across all models is that the forecasts consistently fail to anticipate the sharp increase in flu rates during the test period, particularly around early 2025. This surge is evident in the black line (test set), which climbs steeply and exceeds not just the point forecasts but also the upper bound of the 95% prediction intervals, highlighting a systematic underprediction of the spike.

Models like ETS(MAM) and ETS(MAA) offer somewhat better alignment with the upward trend just before the peak, while models like ETS(AMN) and ETS(MNN) appear unstable or overreactive, with exaggerated prediction intervals and poor directional accuracy. Notably, ETS(AMN) shows a forecast range that is wildly off-scale—this matches its extreme error metrics and renders it unusable for practical forecasting.

Taken together, the chart suggests that these exponential smoothing models, while structurally diverse, lack the flexibility to accommodate sudden regime changes or outbreak dynamics. The consistent underperformance across specifications, despite relatively stable prediction intervals, hints at potential value in either transforming the data (e.g., log or Box-Cox to stabilize spikes), integrating leading indicators, or pivoting toward models that can more explicitly handle shock or threshold behavior.