Forecast Accuracy & Decomposition: HOUST

Setup

library(fpp3)
library(fredr)
library(tidyverse)
library(patchwork)
library(knitr)
library(kableExtra)
library(writexl)

fredr_set_key("bdfddd4beeaf18c36abe8754e4f3929b")

# Import data
myts <- fredr(series_id = "HOUST",
              observation_start = as.Date("1990-01-01"),
              observation_end   = as.Date("2024-12-01")) |>
  transmute(Month = yearmonth(date), value) |>
  as_tsibble(index = Month)

# Visualize raw data
autoplot(myts, value) +
  labs(title = "US Housing Starts (Monthly)", y = "Units")

# Calculate split point
split_index <- floor(0.8 * nrow(myts))

# Create sets
train <- myts[1:split_index, ]
test  <- myts[(split_index + 1):nrow(myts), ]

print(paste("Train ends:", max(train$Month)))

[1] "Train ends: 2017 12月"

print(paste("Test starts:", min(test$Month)))

[1] "Test starts: 2018 1月"

# Fit models
models <- train |>
  model(
    MEAN   = MEAN(value),
    NAIVE  = NAIVE(value),
    SNAIVE = SNAIVE(value),
    DRIFT  = RW(value ~ drift()), 
    WAVG   = TSLM(value ~ trend())
  )

# Forecast for the length of the test set
fc <- models |> forecast(h = nrow(test))

# Visualize forecasts
fc |> 
  autoplot(train, level = 95) + 
  labs(title = "Forecasts vs Actuals", y = "Units") +
  facet_wrap(~ .model, ncol = 2) +
  theme_minimal()

# Compute accuracy on test set
acc <- accuracy(fc, test) |>
  select(.model, ME, MPE, RMSE, MAE, MAPE)

# Display formatted table
kable(acc, caption = "Forecast Accuracy Metrics", digits = 2) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Forecast Accuracy Metrics

.model	ME	MPE	RMSE	MAE	MAPE
DRIFT	281.66	18.65	335.00	291.77	19.66
MEAN	107.99	6.25	201.92	156.42	10.59
NAIVE	234.21	15.33	289.77	246.52	16.55
SNAIVE	206.55	13.36	270.69	224.43	15.04
WAVG	423.42	28.82	462.41	425.17	29.01

# Prepare data: Join Forecasts with Actual values
df_export <- fc |> 
  as_tibble() |> 
  select(Month, .model, .mean) |> 
  pivot_wider(names_from = .model, values_from = .mean) |> 
  left_join(test |> rename(Actual = value), by = "Month")

# Export to Excel
write_xlsx(df_export, "accuracy_check.xlsx")
print("Data exported to 'accuracy_check.xlsx' for manual verification.")

[1] "Data exported to 'accuracy_check.xlsx' for manual verification."

# Additive Decomposition
p1 <- myts |>
  model(STL(value ~ season(window="periodic"))) |>
  components() |>
  autoplot() + ggtitle("Additive Decomposition")

# Multiplicative Decomposition (via Log transformation)
p2 <- myts |>
  model(STL(log(value) ~ season(window="periodic"))) |>
  components() |>
  autoplot() + ggtitle("Multiplicative Decomposition (Log)")

# Display plots
p1

p2