Discussion_4

#libraries
library(fpp3)
library(tsibble)
library(feasts)
library(fable)
library(fredr)
library(ggplot2)
library(forecast)
library(tseries)
library(patchwork)
library(ggtime)

fredr_set_key ("2e2c37e9aafe0f62c76d6bc4ef83765f")

raw_GDP <- fredr(series_id = "GDP", 
      observation_start = as.Date("1947-01-01"),
      observation_end   = as.Date("2025-07-01")) |>
  transmute(Month = yearmonth(date), value) |>
  as_tsibble(index = Month)

raw_UNRATE <- fredr(series_id = "UNRATE",
  observation_start = as.Date("1948-01-01"),
  observation_end   = as.Date("2026-01-01")) |>
  transmute(Month = yearmonth(date), value) |>
  as_tsibble(index = Month)

raw_BEEF <- fredr(series_id = "APU0000703112",
  observation_start = as.Date("1984-01-01"),
  observation_end   = as.Date("2024-12-01")) |>
  transmute(Month = yearmonth(date), value) |>
  as_tsibble(index = Month)

autoplot(raw_GDP, value) +
  labs(
    title = "United States Gross Domestic Product", 
    subtitle = "1947-2025", 
    caption = "https://fred.stlouisfed.org/series/GDP",
    y = "Billions Of Dollars", 
    x = "Year"
)

autoplot(raw_UNRATE, value) +
  labs(
    title = "U.S. Unemployment Rate Throughout the Years",
    subtitle = "1948-2025",
    caption = "https://fred.stlouisfed.org/series/UNRATE",
    y = "Rate (%)",
    x = "Year"
)

autoplot(raw_BEEF, value) +
  labs(
    title = "Average Price: Ground Beef, 100% Beef in U.S. City Average",
    subtitle = "1984-2024",
    caption = "https://fred.stlouisfed.org/series/APU0000703112",
    y = "U.S. Dollars",
    x = "Year"
)

#ACF
feasts::ACF(raw_GDP, value) |> 
  autoplot() + 
  labs(subtitle = "GDP")

feasts::ACF(raw_UNRATE, value) |> 
  autoplot() + 
  labs(subtitle = "UNRATE")

feasts::ACF(raw_BEEF, value) |> 
  autoplot() + 
  labs(subtitle = "BEEF")

#Experimental Names
GDP <- raw_GDP
UNRATE <- raw_UNRATE
BEEF <- raw_BEEF

#GDP Transformations
fit_GDP <- GDP |> model(TSLM(value ~ trend()))
detrended_GDP <- augment(fit_GDP) |> select(Month, .resid)

#GDP Detrending
p1_GDP <- autoplot(GDP, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Quarter")
p2_GDP <- autoplot(detrended_GDP, .resid) + 
  labs(title = "Detrended", 
       y = "Residuals",
       x="Quarter")
#GDP First Difference
differenced_GDP <- GDP |> mutate(diff_value = difference(value))

p3_GDP <- autoplot(GDP, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Quarter")
p4_GDP <- autoplot(differenced_GDP, diff_value) + 
  labs(title = "First Difference", 
       y = "Δ Value",
       x="Quarter")
#GDP Log Transformed
log_transformed_GDP <- GDP |> mutate(log_value = log(value))

p5_GDP <- autoplot(GDP, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Quarter")
p6_GDP <- autoplot(log_transformed_GDP, log_value) + 
  labs(title = "Log Transformed", 
       y = "log(Value)",
       x="Quarter")
#GDP Seasonal Difference
seasonal_diff_GDP <- GDP |> mutate(seas_diff = difference(value, lag = 4))

p7_GDP <- autoplot(GDP, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Quarter")
p8_GDP <- autoplot(seasonal_diff_GDP, seas_diff) + 
  labs(title = "Seasonal Difference (lag=4)", 
       y = "Change in y_{4} Value",
       x="Quarter")

(p1_GDP / p2_GDP) | (p3_GDP / p4_GDP)

(p5_GDP / p6_GDP) | (p7_GDP / p8_GDP)

#GDP Log Difference
logdiff_GDP <- GDP |> 
  mutate(dlog_value = difference(log(value))) |>
  filter(!is.na(dlog_value))

p9_GDP <- autoplot(GDP, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Quarter")
p10_GDP <- autoplot(logdiff_GDP, dlog_value) +
  labs(title = "Log First Difference", 
       y = "Δ log(Value)",
       x="Quarter")

(p9_GDP / p10_GDP)

#UNRATE Transformations

#Detrending
fit_UNRATE <- UNRATE |> model(TSLM(value ~ trend()))
detrended_UNRATE <- augment(fit_UNRATE) |> select(Month, .resid)

p1_UNRATE <- autoplot(UNRATE, value) + 
  labs(title = "Original Series", 
       y = "Rate(%)",
       x="Month")
p2_UNRATE <- autoplot(detrended_UNRATE, .resid) + 
  labs(title = "Detrended", 
       y = "Residuals",
       x="Month")

#Differncing
differenced_UNRATE <- UNRATE |> mutate(diff_value = difference(value))

p3_UNRATE <- autoplot(UNRATE, value) + 
  labs(title="Original Series", 
       y="Rate(%)",
       x="Month")
p4_UNRATE <- autoplot(differenced_UNRATE, diff_value) + 
  labs(title="First Difference", 
       y="Δ Rate (pp)",
       x="Month")

(p1_UNRATE / p2_UNRATE) | (p3_UNRATE / p4_UNRATE)

#BEEF Transformations
fit_BEEF <- BEEF |> model(TSLM(value ~ trend()))
detrended_BEEF <- augment(fit_BEEF) |> select(Month, .resid)

#BEEF Detrending
p1_BEEF <- autoplot(BEEF, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Month")
p2_BEEF <- autoplot(detrended_BEEF, .resid) + 
  labs(title = "Detrended", 
       y = "Residuals",
       x="Month")

#BEEF First Difference
differenced_BEEF <- BEEF |> mutate(diff_value = difference(value))

p3_BEEF <- autoplot(BEEF, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Month")
p4_BEEF <- autoplot(differenced_BEEF, diff_value) + 
  labs(title = "First Difference", 
       y = "Δ Value",
       x="Month")

#GDP Log Transformed
log_transformed_BEEF <- BEEF |> mutate(log_value = log(value))

p5_BEEF <- autoplot(BEEF, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Month")
p6_BEEF <- autoplot(log_transformed_BEEF, log_value) + 
  labs(title = "Log Transformed", 
       y = "log(Value)",
       x="Month")

#Beef Seasonal Difference
seasonal_diff_BEEF <- BEEF |> mutate(seas_diff = difference(value, lag = 12))

p7_BEEF <- autoplot(BEEF, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Month")
p8_BEEF <- autoplot(seasonal_diff_BEEF, seas_diff) + 
  labs(title = "Seasonal Difference (lag=12)", 
       y = "Change in y_{12} Value",
       x="Month")

(p1_BEEF / p2_BEEF) | (p3_BEEF / p4_BEEF)

(p5_BEEF / p6_BEEF) | (p7_BEEF / p8_BEEF)

#BEEF Log Difference
logdiff_BEEF <- BEEF |>
  tsibble::fill_gaps() |>
  mutate(dlog_value = difference(log(value))) |>
  slice(-1)

p9_BEEF <- autoplot(BEEF, value) + 
  labs(title = "Original Series", 
       y = "Value",
       x="Month")
p10_BEEF <- autoplot(logdiff_BEEF, dlog_value) +
  labs(title = "Log First Difference", 
       y = "Δ log(Value)",
       x="Month")

(p9_BEEF / p10_BEEF)

#Fix Knit
UNRATE_clean <- UNRATE |>
  tsibble::fill_gaps() |>
  tidyr::fill(value, .direction = "downup")

stats_UNRATE <- UNRATE |>
  mutate(diff_value = difference(value)) |>
  slice(-1)

BEEF_clean <- BEEF |>
  tsibble::fill_gaps() |>
  tidyr::fill(value, .direction = "downup")

#GDP KPSS
GDP |> features(value, unitroot_kpss)

# A tibble: 1 × 2
  kpss_stat kpss_pvalue
      <dbl>       <dbl>
1      4.72        0.01

logdiff_GDP |> features(dlog_value, unitroot_kpss)

# A tibble: 1 × 2
  kpss_stat kpss_pvalue
      <dbl>       <dbl>
1     0.865        0.01

#UNRATE ADF & KPSS
UNRATE_clean |> features(value, unitroot_kpss)

# A tibble: 1 × 2
  kpss_stat kpss_pvalue
      <dbl>       <dbl>
1     0.960        0.01

stats_UNRATE |> features(diff_value, unitroot_kpss)

# A tibble: 1 × 2
  kpss_stat kpss_pvalue
      <dbl>       <dbl>
1    0.0400         0.1

#BEEF KPSS
BEEF_clean |> features(value, unitroot_kpss)

# A tibble: 1 × 2
  kpss_stat kpss_pvalue
      <dbl>       <dbl>
1      7.32        0.01

logdiff_BEEF |> features(dlog_value, unitroot_kpss)

# A tibble: 1 × 2
  kpss_stat kpss_pvalue
      <dbl>       <dbl>
1     0.238         0.1

alpha <- 0.05

kpss_decision <- function(p) ifelse(p < alpha, "Reject H0 (NOT stationary)", "Fail to reject H0 (stationary)")

#KPSS GDP untransformed
kpss_GDP <- GDP |> features(value, unitroot_kpss)
kpss_decision(kpss_GDP$kpss_pvalue)

[1] "Reject H0 (NOT stationary)"

#KPSS UNRATE untransformed
kpss_UNRATE <- UNRATE |> features(value, unitroot_kpss)
kpss_decision(kpss_UNRATE$kpss_pvalue)

[1] "Reject H0 (NOT stationary)"

#KPSS BEEF untransformed
kpss_BEEF <- BEEF |> features(value, unitroot_kpss)
kpss_decision(kpss_BEEF$kpss_pvalue)

[1] "Reject H0 (NOT stationary)"

#KPSS GDP transformed
kpss_logdiff_GDP <- logdiff_GDP |> features(dlog_value, unitroot_kpss)
kpss_decision(kpss_logdiff_GDP$kpss_pvalue)

[1] "Reject H0 (NOT stationary)"

#KPSS UNRATE transformed
kpss_diff_UNRATE <- stats_UNRATE |> features(diff_value, unitroot_kpss)
kpss_decision(kpss_diff_UNRATE$kpss_pvalue)

[1] "Fail to reject H0 (stationary)"

#KPSS BEEF transformed
kpss_logdiff_BEEF <- logdiff_BEEF |> features(dlog_value, unitroot_kpss)
kpss_decision(kpss_logdiff_BEEF$kpss_pvalue)

[1] "Fail to reject H0 (stationary)"

#ADF
adf_decision <- function(p) ifelse(
  p < alpha,
  "Reject H0 (NO unit root -> stationary)",
  "Fail to reject H0 (unit root -> NOT stationary)"
)

#ADF GDP
adf.test(na.omit(GDP$value), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(GDP$value)
Dickey-Fuller = 2.9862, Lag order = 6, p-value = 0.99
alternative hypothesis: stationary

adf.test(na.omit(logdiff_GDP$dlog_value),  alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(logdiff_GDP$dlog_value)
Dickey-Fuller = -5.9411, Lag order = 6, p-value = 0.01
alternative hypothesis: stationary

#ADF UNRATE
adf.test(na.omit(UNRATE_clean$value), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(UNRATE_clean$value)
Dickey-Fuller = -3.9726, Lag order = 9, p-value = 0.01037
alternative hypothesis: stationary

adf.test(na.omit(stats_UNRATE$diff_value), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(stats_UNRATE$diff_value)
Dickey-Fuller = -9.7886, Lag order = 9, p-value = 0.01
alternative hypothesis: stationary

#ADF BEEF
adf.test(na.omit(BEEF_clean$value), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(BEEF_clean$value)
Dickey-Fuller = -0.63105, Lag order = 7, p-value = 0.9758
alternative hypothesis: stationary

adf.test(na.omit(logdiff_BEEF$dlog_value), alternative = "stationary")

    Augmented Dickey-Fuller Test

data:  na.omit(logdiff_BEEF$dlog_value)
Dickey-Fuller = -8.033, Lag order = 7, p-value = 0.01
alternative hypothesis: stationary

#GDP
adf_GDP <- adf.test(na.omit(GDP$value), alternative="stationary")
adf_decision(adf_GDP$p.value)

[1] "Fail to reject H0 (unit root -> NOT stationary)"

adf_logdiff_GDP <- adf.test(na.omit(logdiff_GDP$dlog_value), alternative="stationary")
adf_decision(adf_logdiff_GDP$p.value)

[1] "Reject H0 (NO unit root -> stationary)"

#UNRATE
adf_UNRATE <- adf.test(na.omit(UNRATE_clean$value), alternative = "stationary")
adf_decision(adf_UNRATE$p.value)

[1] "Reject H0 (NO unit root -> stationary)"

adf_diff_UNRATE <- adf.test(na.omit(stats_UNRATE$diff_value), alternative = "stationary")
adf_decision(adf_diff_UNRATE$p.value)

[1] "Reject H0 (NO unit root -> stationary)"

#BEEF
adf_BEEF <- adf.test(na.omit(BEEF_clean$value), alternative = "stationary")
adf_decision(adf_BEEF$p.value)

[1] "Fail to reject H0 (unit root -> NOT stationary)"

adf_logdiff_BEEF <- adf.test(na.omit(logdiff_BEEF$dlog_value), alternative = "stationary")
adf_decision(adf_logdiff_BEEF$p.value)

[1] "Reject H0 (NO unit root -> stationary)"

#ACF & PCF For GDP
logdiff_GDP |> ggtime::gg_tsdisplay(dlog_value, plot_type = "partial")

#ACF & PCF For UNRATE
stats_UNRATE <- UNRATE |>
  mutate(diff_value = difference(value)) |>
  slice(-1)

stats_UNRATE |> ggtime::gg_tsdisplay(diff_value, plot_type = "partial")

#ACF & PCF For BEEF
logdiff_BEEF |> ggtime::gg_tsdisplay(dlog_value, plot_type = "partial")

#GDP Decomposition Plots
dcmp_GDP <- GDP |>
  model(stl = STL(value ~ season(window = "periodic"), robust = TRUE)) |>
  components()

autoplot(dcmp_GDP) + labs(
  title = "GDP: STL decomposition",
  x= "Quarterly")

#UNRATE Decomposition Plots
UNRATE_clean <- UNRATE |>
  tsibble::fill_gaps() |>
  tidyr::fill(value, .direction = "downup")

dcmp_UNRATE <- UNRATE_clean |>
  model(stl = STL(value ~ season(window = "periodic"), robust = TRUE)) |>
  components()

autoplot(dcmp_UNRATE) + labs(title = "UNRATE: STL decomposition")

#BEEF Decomposition Plots
BEEF_clean <- BEEF |>
  tsibble::fill_gaps() |>
  tidyr::fill(value, .direction = "downup")

dcmp_BEEF <- BEEF_clean |>
  model(stl = STL(value ~ season(window = "periodic"), robust = TRUE)) |>
  components()

autoplot(dcmp_BEEF) + labs(title = "BEEF: STL decomposition")

#Optional Decomposition Plots
dcmp1 <- dcmp_GDP |> 
  gg_season(season_year) + 
  labs(title = "GDP seasonal component (gg_season)",
       x= "Quarter")
dcmp2 <- dcmp_GDP |> 
  gg_subseries(season_year) + 
  labs(title = "GDP seasonal subseries",
       x= "Quarter")
dcmp3 <- dcmp_BEEF |> gg_season(season_year) + labs(title = "BEEF seasonal component (gg_season)")
dcmp4 <- dcmp_BEEF |> gg_subseries(season_year) + labs(title = "BEEF seasonal subseries")

(dcmp1 / dcmp2) 

(dcmp3 / dcmp4)

Random Citation (Hyndman, n.d.)

References

Hyndman, Rob. n.d. “Fpp3: Data for "Forecasting: Principles and Practice" (3rd Edition).” https://doi.org/10.32614/CRAN.package.fpp3.