Hierarchical ODE State Space Model for Unemployment Dynamics

A Bayesian Approach with Pooled Labor Market Parameters

WORK IN PROGRESS

This report is currently under development. Not all statements, interpretations, and conclusions have been reviewed or validated by a human researcher. Please treat all content as preliminary analysis subject to revision.

Overview

This report presents a hierarchical Bayesian state space model for unemployment dynamics across education levels. The model uses a time-varying equilibrium — constructed from a B-spline basis — that the ODE pulls unemployment toward. Shocks are cleanly identified as temporary deviations from this structural baseline.

Key Features

Economic Foundation - Explicit labor market dynamics: \(\frac{du}{dt} = \text{adj} \cdot (u_{\text{eq}}(t) - u)\) - Time-varying equilibrium \(u_{\text{eq}}(t)\) captures structural changes in education-specific labor markets - Adjustment speed (\(\text{adj}\)) controls convergence rate toward equilibrium - Shocks add to separation rate, creating transitory deviations

Hierarchical Structure - Population-level means with education-specific deviations (non-centered parameterization) - Pooling strengthens inference for small groups (e.g., PhD, professional degrees) - Hierarchical parameters: adjustment speed, shock magnitudes, recovery rates, seasonal effects

Shock Modeling - Asymmetric impulse functions: power-law rise (\(\rho\)) + stretched-exponential decay (\(\kappa\)) - 2008: convex rise (\(\rho=1.47\)) — gradual build then acceleration; stretched decay (\(\kappa=0.65\)) - 2020: near-linear rise (\(\rho=0.90\)); strongly stretched decay (\(\kappa=0.53\)) — disruption persisted - Education-specific magnitudes and recovery rates with hierarchical pooling

Spline-Based Equilibrium - B-spline basis (K=10) constructs smooth time-varying \(u_{\text{eq}}(t)\) - Random walk prior on spline coefficients for adaptive smoothness - Each education level has its own equilibrium trajectory with data-informed flexibility

Technical Approach - Non-centered parameterization for efficient MCMC sampling - Asymmetric shock shapes (power-law rise, stretched-exponential decay) for flexible impulse modeling - Education-level parallelization via Stan’s reduce_sum for multi-threaded sampling - Prior-centered initialization via make_init_at_prior_edu_parallel()

Recent Trends: December 2023 - April 2026

A close-up view of the most recent unemployment dynamics across education levels.

# Extract de-seasonalized trend for recent trends plot
trend_all <- extract_trend(result, summary = TRUE)

# Filter to December 2023 - April 2026
# Formula: year_frac = year + (month - 0.5) / 12
# December 2023 = 2023 + 11.5/12 ≈ 2023.958
# April 2026   = 2026 +  3.5/12 ≈ 2026.292
recent_trend <- data.frame(
  year_frac = trend_all$year_frac,
  education = trend_all$education,
  mean = trend_all$mean,
  lower = trend_all$q5,
  upper = trend_all$q95
)

# Add observed rates
for (i in seq_along(stan_data$education_levels)) {
  edu <- stan_data$education_levels[i]
  idx <- recent_trend$education == edu
  obs_rate <- stan_data$n_unemployed[, i] / stan_data$n_total[, i]
  recent_trend$observed[idx] <- obs_rate
}

# Filter to recent period
recent_trend <- recent_trend[recent_trend$year_frac >= 2023.958 &
                                recent_trend$year_frac <= 2026.292, ]

ggplot(recent_trend, aes(x = year_frac)) +
  geom_ribbon(aes(ymin = lower * 100, ymax = upper * 100, fill = education),
              alpha = 0.3) +
  geom_line(aes(y = mean * 100, color = education), linewidth = 1) +
  geom_point(aes(y = observed * 100), alpha = 0.5, size = 1.5) +
  facet_wrap(~education, scales = "free_y", ncol = 2) +
  scale_x_continuous(breaks = seq(2024, 2026.5, 0.5),
                    labels = function(x) {
                      year <- floor(x)
                      month <- round((x - year) * 12) + 1
                      sprintf("%d-%02d", year, month)
                    }) +
  scale_y_continuous(labels = scales::percent_format(scale = 1)) +
  labs(x = "Date",
       y = "Unemployment Rate",
       title = "Recent Unemployment Trends by Education Level",
       subtitle = "Points: observed data; Lines: model de-seasonalized trend; Bands: 90% credible intervals") +
  theme(legend.position = "none",
        strip.text = element_text(size = 11, face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1))

Recent unemployment trends by education level (Dec 2023 - Apr 2026)

Hierarchical Variance Analysis: Comparing Pooling Strength

The hierarchical structure pools information across education levels. The between-education standard deviation (\(\sigma\)) parameters quantify how much groups differ. Smaller \(\sigma\) = stronger pooling.

# Extract sigma parameters directly from v4 fit
N_edu <- result$stan_data$N_edu
edu_levels <- result$stan_data$education_levels

sigma_params <- c(
  "sigma_log_adj_speed",
  "sigma_log_shock_2008", "sigma_log_shock_2020",
  "sigma_decay_2008", "sigma_decay_2020",
  "sigma_seasonal"
)
sigma_summary <- result$fit$summary(variables = sigma_params)

cat("\n=== Between-Education Standard Deviations ===\n")

=== Between-Education Standard Deviations ===

cat("Smaller values = stronger pooling = more similar across education groups\n\n")

Smaller values = stronger pooling = more similar across education groups

sigma_display <- data.frame(
  Parameter = c("Adjustment speed", "2008 shock", "2020 shock",
                "2008 decay", "2020 decay", "Seasonal"),
  `Sigma median` = round(sigma_summary$median, 3),
  `90% CI` = sprintf("[%.3f, %.3f]", sigma_summary$q5, sigma_summary$q95),
  check.names = FALSE
)
knitr::kable(sigma_display, align = c("l", "r", "r"))

Parameter	Sigma median	90% CI
Adjustment speed	0.700	[0.470, 1.172]
2008 shock	0.102	[0.025, 0.236]
2020 shock	0.110	[0.060, 0.223]
2008 decay	2.844	[1.618, 4.942]
2020 decay	0.191	[0.015, 0.762]
Seasonal	0.048	[0.041, 0.057]

sigma_plot <- data.frame(
  parameter = factor(sigma_display$Parameter, levels = rev(sigma_display$Parameter)),
  median = sigma_summary$median,
  q5 = sigma_summary$q5,
  q95 = sigma_summary$q95
)

ggplot(sigma_plot, aes(x = median, y = parameter)) +
  geom_point(size = 3, color = "steelblue") +
  geom_errorbarh(aes(xmin = q5, xmax = q95), height = 0.2, color = "steelblue") +
  labs(x = "Between-Education SD", y = NULL,
       title = "Hierarchical Pooling Strength",
       subtitle = "Smaller values = education groups are more similar") +
  theme_minimal(base_size = 13)

Between-education standard deviations with 90% credible intervals

Education-Specific Parameter Values

Parameter estimates for key quantities across all seven education levels:

# Build per-edu parameter table
edu_params <- rbind(
  data.frame(result$fit$summary(variables = paste0("adj_speed[", 1:N_edu, "]")),
             param = "Adjustment Speed", edu = edu_levels),
  data.frame(result$fit$summary(variables = paste0("shock_2008_effect[", 1:N_edu, "]")),
             param = "2008 Shock Effect", edu = edu_levels),
  data.frame(result$fit$summary(variables = paste0("shock_2020_effect[", 1:N_edu, "]")),
             param = "2020 Shock Effect", edu = edu_levels),
  data.frame(result$fit$summary(variables = paste0("decay_2008[", 1:N_edu, "]")),
             param = "2008 Decay Rate", edu = edu_levels),
  data.frame(result$fit$summary(variables = paste0("decay_2020[", 1:N_edu, "]")),
             param = "2020 Decay Rate", edu = edu_levels),
  data.frame(result$fit$summary(variables = paste0("u_eq_mean[", 1:N_edu, "]")),
             param = "Mean Equilibrium", edu = edu_levels)
)

ggplot(edu_params, aes(x = median, y = edu)) +
  geom_point(size = 2.5, color = "steelblue") +
  geom_errorbarh(aes(xmin = q5, xmax = q95), height = 0.2, color = "steelblue") +
  facet_wrap(~ param, scales = "free_x", ncol = 2) +
  labs(y = NULL, x = "Posterior median with 90% CI",
       title = "Education-Specific Parameter Estimates") +
  theme_minimal(base_size = 12)

Key parameter estimates by education level with 90% credible intervals

Key Insights

1. Adjustment speeds increase with education: less_than_hs has the slowest adjustment, PhD and professional degrees the fastest — consistent with more fluid labor markets at higher education levels.

2. Shock effects are tightly pooled across education groups — the Exponential(20) prior on \(\sigma\) reflects the economic prior that macro shocks affect all workers similarly. The model estimates similar shock magnitudes across groups.

3. Mean equilibrium unemployment ranges from ~1.3% (professional degrees) to ~10.2% (less than high school) — matching the raw data well while the time-varying equilibrium captures structural trends within each group.

Economic Parameter Estimates

Adjustment Speed and Equilibrium

The model estimates education-specific adjustment speeds and time-averaged equilibrium unemployment rates.

# Extract parameter summaries directly
N_edu <- result$stan_data$N_edu
edu_levels <- result$stan_data$education_levels

# Adjustment speeds
cat("\n=== Adjustment Speeds (adj = s + f, monthly) ===\n")

=== Adjustment Speeds (adj = s + f, monthly) ===

adj_summary <- result$fit$summary(variables = paste0("adj_speed[", 1:N_edu, "]"))
adj_summary$education <- edu_levels
print(adj_summary[, c("education", "median", "q5", "q95")], digits = 3)

# A tibble: 7 × 4
  education    median    q5   q95
  <chr>         <dbl> <dbl> <dbl>
1 bachelors     14.8  13.4  16.5 
2 high_school    8.25  7.60  8.99
3 less_than_hs   6.17  5.66  6.74
4 masters       22.8  20.3  25.4 
5 phd           37.3  30.0  45.2 
6 professional  30.9  26.5  35.9 
7 some_college  10.4   9.69 11.3 

# Time-averaged equilibrium
cat("\n=== Mean Equilibrium Unemployment Rates ===\n")

=== Mean Equilibrium Unemployment Rates ===

cat("(Time-average of u_eq[t] over the full period)\n\n")

(Time-average of u_eq[t] over the full period)

eq_summary <- result$fit$summary(variables = paste0("u_eq_mean[", 1:N_edu, "]"))
eq_summary$education <- edu_levels
print(eq_summary[, c("education", "median", "q5", "q95")], digits = 4)

# A tibble: 7 × 4
  education    median     q5    q95
  <chr>         <dbl>  <dbl>  <dbl>
1 bachelors    0.0247 0.0212 0.0269
2 high_school  0.0439 0.0408 0.0519
3 less_than_hs 0.103  0.0946 0.106 
4 masters      0.0223 0.0207 0.0229
5 phd          0.0154 0.0150 0.0158
6 professional 0.0147 0.0139 0.0151
7 some_college 0.0376 0.0327 0.0413

eq_df <- data.frame(
  education = edu_levels,
  median = eq_summary$median,
  lower = eq_summary$q5,
  upper = eq_summary$q95
)

eq_df$education <- factor(eq_df$education,
                          levels = eq_df$education[order(eq_df$median)])

ggplot(eq_df, aes(x = education, y = median * 100)) +
  geom_point(size = 4, color = "steelblue") +
  geom_errorbar(aes(ymin = lower * 100, ymax = upper * 100),
                width = 0.2, linewidth = 1, color = "steelblue") +
  coord_flip() +
  labs(x = "Education Level",
       y = "Mean Equilibrium Unemployment Rate (%)",
       title = "Structural Unemployment Rates by Education",
       subtitle = "Time-averaged u_eq(t) — the baseline toward which the ODE pulls") +
  theme(axis.text.y = element_text(size = 11))

Time-averaged equilibrium unemployment rates by education level

Shock Effects

The key advantage of the state space model: shock effects are positive and identifiable. With hierarchical priors, we estimate both population-level means and education-specific effects.

# Population-level shock parameters (hierarchical)
cat("\n=== Hierarchical Shock Parameters (Population Level) ===\n")

=== Hierarchical Shock Parameters (Population Level) ===

hier_params <- c("mu_log_shock_2008", "sigma_log_shock_2008",
                 "mu_log_shock_2020", "sigma_log_shock_2020")
hier_summary <- result$fit$summary(hier_params)
print(hier_summary[, c("variable", "mean", "sd", "q5", "q95", "rhat", "ess_bulk")], digits = 3)

# A tibble: 4 × 7
  variable               mean     sd      q5    q95  rhat ess_bulk
  <chr>                 <dbl>  <dbl>   <dbl>  <dbl> <dbl>    <dbl>
1 mu_log_shock_2008    -0.748 0.0832 -0.894  -0.624  1.00    1760.
2 sigma_log_shock_2008  0.113 0.0664  0.0253  0.236  1.00    1400.
3 mu_log_shock_2020     0.550 0.0689  0.436   0.650  1.00    2651.
4 sigma_log_shock_2020  0.121 0.0559  0.0596  0.223  1.00    3459.

cat("\n\nInterpretation on natural scale:\n")

Interpretation on natural scale:

cat(sprintf("  2008 mean effect: exp(%.2f) = %.1f%% increase in separation rate\n",
            hier_summary$mean[hier_summary$variable == "mu_log_shock_2008"],
            100 * exp(hier_summary$mean[hier_summary$variable == "mu_log_shock_2008"])))

  2008 mean effect: exp(-0.75) = 47.3% increase in separation rate

cat(sprintf("  2020 mean effect: exp(%.2f) = %.1f%% increase in separation rate\n",
            hier_summary$mean[hier_summary$variable == "mu_log_shock_2020"],
            100 * exp(hier_summary$mean[hier_summary$variable == "mu_log_shock_2020"])))

  2020 mean effect: exp(0.55) = 173.2% increase in separation rate

cat("\n=== 2008 Financial Crisis Shock Effects (Education-Specific) ===\n")

=== 2008 Financial Crisis Shock Effects (Education-Specific) ===

shock08 <- result$fit$summary(variables = paste0("shock_2008_effect[", 1:N_edu, "]"))
shock08$education <- edu_levels
print(shock08[, c("education", "median", "q5", "q95")], digits = 4)

# A tibble: 7 × 4
  education    median    q5   q95
  <chr>         <dbl> <dbl> <dbl>
1 bachelors     0.457 0.392 0.524
2 high_school   0.547 0.484 0.612
3 less_than_hs  0.494 0.430 0.558
4 masters       0.470 0.398 0.542
5 phd           0.456 0.344 0.544
6 professional  0.459 0.362 0.541
7 some_college  0.476 0.423 0.531

cat("\n=== 2020 COVID-19 Shock Effects (Education-Specific) ===\n")

=== 2020 COVID-19 Shock Effects (Education-Specific) ===

shock20 <- result$fit$summary(variables = paste0("shock_2020_effect[", 1:N_edu, "]"))
shock20$education <- edu_levels
print(shock20[, c("education", "median", "q5", "q95")], digits = 4)

# A tibble: 7 × 4
  education    median    q5   q95
  <chr>         <dbl> <dbl> <dbl>
1 bachelors      1.82  1.68  1.98
2 high_school    1.75  1.63  1.87
3 less_than_hs   1.62  1.49  1.75
4 masters        1.67  1.50  1.84
5 phd            1.58  1.28  1.84
6 professional   1.79  1.56  2.03
7 some_college   2.05  1.92  2.18

cat("\n=== 2008 Shock Half-Lives (Education-Specific) ===\n")

=== 2008 Shock Half-Lives (Education-Specific) ===

cat("(Time in years for shock effect to decay by 50%)\n\n")

(Time in years for shock effect to decay by 50%)

hl08 <- result$fit$summary(variables = paste0("halflife_2008[", 1:N_edu, "]"))
hl08$education <- edu_levels
print(hl08[, c("education", "median", "q5", "q95")], digits = 2)

# A tibble: 7 × 4
  education    median    q5   q95
  <chr>         <dbl> <dbl> <dbl>
1 bachelors     1.81  0.881 4.34 
2 high_school   5.14  2.23  6.88 
3 less_than_hs  1.33  0.746 3.07 
4 masters       0.676 0.368 1.64 
5 phd           0.151 0.139 0.286
6 professional  0.180 0.139 0.860
7 some_college  2.40  1.19  5.29 

cat("\n=== 2020 COVID Half-Lives (Education-Specific) ===\n")

=== 2020 COVID Half-Lives (Education-Specific) ===

cat("(Time in years for shock effect to decay by 50%)\n\n")

(Time in years for shock effect to decay by 50%)

hl20 <- result$fit$summary(variables = paste0("halflife_2020[", 1:N_edu, "]"))
hl20$education <- edu_levels
print(hl20[, c("education", "median", "q5", "q95")], digits = 2)

# A tibble: 7 × 4
  education    median    q5   q95
  <chr>         <dbl> <dbl> <dbl>
1 bachelors     0.170 0.154 0.197
2 high_school   0.168 0.153 0.190
3 less_than_hs  0.163 0.147 0.181
4 masters       0.167 0.152 0.190
5 phd           0.169 0.153 0.203
6 professional  0.165 0.150 0.187
7 some_college  0.165 0.151 0.185

shock_df <- rbind(
  data.frame(
    shock = "2008 Financial Crisis",
    education = edu_levels,
    median = shock08$median,
    lower = shock08$q5,
    upper = shock08$q95
  ),
  data.frame(
    shock = "2020 COVID-19",
    education = edu_levels,
    median = shock20$median,
    lower = shock20$q5,
    upper = shock20$q95
  )
)

ggplot(shock_df, aes(x = education, y = median * 100, color = shock)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = lower * 100, ymax = upper * 100),
                position = position_dodge(width = 0.5), width = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  coord_flip() +
  scale_color_manual(values = c("2008 Financial Crisis" = "#E74C3C",
                                "2020 COVID-19" = "#3498DB")) +
  labs(x = "Education Level",
       y = "Additional Separation Rate (%)",
       color = "Economic Shock",
       title = "Shock Effects on Job Separation by Education",
       subtitle = "Positive values indicate increased job loss during crisis") +
  theme(legend.position = "bottom")

Comparison of shock effects across education levels

Latent Unemployment Trajectories

# Extract latent unemployment rates
latent <- extract_latent_rates(result, summary = TRUE)

# Create data for plotting
plot_data <- data.frame(
  year_frac = latent$year_frac,
  education = latent$education,
  mean = latent$mean,
  lower = latent$q5,
  upper = latent$q95
)

# Add observed rates for comparison
for (i in seq_along(stan_data$education_levels)) {
  edu <- stan_data$education_levels[i]
  idx <- plot_data$education == edu
  obs_rate <- stan_data$n_unemployed[, i] / stan_data$n_total[, i]
  plot_data$observed[idx] <- obs_rate
}

ggplot(plot_data, aes(x = year_frac)) +
  geom_ribbon(aes(ymin = lower * 100, ymax = upper * 100, fill = education),
              alpha = 0.3) +
  geom_line(aes(y = mean * 100, color = education), linewidth = 0.8) +
  geom_point(aes(y = observed * 100), alpha = 0.3, size = 0.5) +
  geom_vline(xintercept = c(2008.75, 2020.25), linetype = "dashed",
             color = "red", alpha = 0.5) +
  facet_wrap(~education, scales = "free_y", ncol = 2) +
  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  scale_y_continuous(labels = scales::percent_format(scale = 1)) +
  labs(x = "Year",
       y = "Unemployment Rate",
       title = "Latent Unemployment Trajectories by Education Level",
       subtitle = "Points: observed data; Lines: model estimates; Bands: 90% CI; Red lines: shock onsets") +
  theme(legend.position = "none",
        strip.text = element_text(size = 10))

Latent unemployment rates from state space model (with 90% credible bands)

Non-Seasonal Trend

The state space model allows us to extract the underlying trend by removing seasonal effects. This shows the unemployment dynamics driven only by:

• Baseline separation and finding rates
• Economic shock effects (2008, 2020)
• Stochastic innovations

# Extract the non-seasonal trend
trend <- extract_trend(result, summary = TRUE)

# Create data for plotting
trend_data <- data.frame(
  year_frac = trend$year_frac,
  education = trend$education,
  mean = trend$mean,
  lower = trend$q5,
  upper = trend$q95
)

ggplot(trend_data, aes(x = year_frac)) +
  geom_ribbon(aes(ymin = lower * 100, ymax = upper * 100, fill = education),
              alpha = 0.3) +
  geom_line(aes(y = mean * 100, color = education), linewidth = 1) +
  geom_vline(xintercept = c(2008.75, 2020.25), linetype = "dashed",
             color = "red", alpha = 0.7) +
  annotate("text", x = 2009.5, y = Inf, label = "2008 Crisis",
           vjust = 2, hjust = 0, size = 3, color = "red") +
  annotate("text", x = 2021, y = Inf, label = "COVID-19",
           vjust = 2, hjust = 0, size = 3, color = "red") +
  facet_wrap(~education, scales = "free_y", ncol = 2) +
  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  scale_y_continuous(labels = scales::percent_format(scale = 1)) +
  labs(x = "Year",
       y = "Unemployment Rate (Trend)",
       title = "Non-Seasonal Unemployment Trend by Education Level",
       subtitle = "Seasonal effects removed; shows baseline dynamics + shock impacts") +
  theme(legend.position = "none",
        strip.text = element_text(size = 10))

Non-seasonal unemployment trend (shock + baseline dynamics only)

Seasonal Oscillation: Direct Visualization

The seasonal effect is the difference between the full model (with seasonality) and the trend (without). This directly shows how much unemployment oscillates due to seasonal hiring patterns.

# Extract seasonal effects directly computed by the model
seasonal_effects <- extract_seasonal_effects(result, summary = TRUE)

# Focus on key education levels
sample_edu <- c("phd", "bachelors", "less_than_hs")
seas_subset <- seasonal_effects[seasonal_effects$education %in% sample_edu, ]

# Plot the seasonal oscillation over time
ggplot(seas_subset, aes(x = year_frac)) +
  geom_ribbon(aes(ymin = q5 * 100, ymax = q95 * 100, fill = education),
              alpha = 0.3) +
  geom_line(aes(y = mean * 100, color = education), linewidth = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
  facet_wrap(~education, scales = "free_y", ncol = 1) +
  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  labs(x = "Year",
       y = "Seasonal Effect (percentage points)",
       title = "Seasonal Oscillation in Unemployment by Education",
       subtitle = "Positive = unemployment above trend; Negative = below trend; Band = 90% CI") +
  theme(legend.position = "none",
        strip.text = element_text(size = 12, face = "bold"))

Seasonal oscillation in unemployment rates (Full model minus Trend)

Monthly Seasonal Pattern

The model includes a direct seasonal effect on unemployment (on the logit scale). This captures observed seasonal patterns like academic hiring calendars, summer employment fluctuations, and annual hiring cycles.

# Extract direct seasonal unemployment parameters
seasonal_u_summary <- result$fit$summary(variables = "seasonal_u")

# Parse indices
seasonal_u_summary$month <- as.integer(
  gsub("seasonal_u\\[(\\d+),\\d+\\]", "\\1", seasonal_u_summary$variable)
)
seasonal_u_summary$edu_index <- as.integer(
  gsub("seasonal_u\\[\\d+,(\\d+)\\]", "\\1", seasonal_u_summary$variable)
)
seasonal_u_summary$education <- stan_data$education_levels[seasonal_u_summary$edu_index]

# Create month labels
month_labels <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun",
                  "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
seasonal_u_summary$month_label <- factor(month_labels[seasonal_u_summary$month],
                                          levels = month_labels)

# Plot for selected education levels
seasonal_u_plot <- seasonal_u_summary[seasonal_u_summary$education %in% sample_edu, ]

ggplot(seasonal_u_plot, aes(x = month_label, y = mean, group = education)) +
  geom_ribbon(aes(ymin = q5, ymax = q95, fill = education),
              alpha = 0.3) +
  geom_line(aes(color = education), linewidth = 1.2) +
  geom_point(aes(color = education), size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
  facet_wrap(~education, ncol = 1, scales = "free_y") +
  labs(x = "Month",
       y = "Seasonal Effect (logit scale)",
       title = "Direct Seasonal Effects on Unemployment",
       subtitle = "Positive = higher unemployment; Negative = lower unemployment; Band = 90% CI") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1),
        strip.text = element_text(size = 12, face = "bold"))

Direct monthly seasonal effects on unemployment by education level

Model Decomposition: Equilibrium, Trend, and Full Trajectory

This visualization compares three components:

1. Equilibrium (green): The time-varying \(u_{\text{eq}}[t]\) — where the ODE pulls unemployment. No shocks, no seasonal.
2. Trend (orange): Equilibrium + shock effects (2008 and 2020). The de-seasonalized dynamics.
3. Full model (blue): Trend + seasonal effects. The complete model prediction.

# Extract time-varying equilibrium
eq_trajectory <- result$fit$summary(variables = "u_eq")
eq_trajectory$time_index <- as.integer(gsub("u_eq\\[(\\d+),\\d+\\]", "\\1", eq_trajectory$variable))
eq_trajectory$edu_index <- as.integer(gsub("u_eq\\[\\d+,(\\d+)\\]", "\\1", eq_trajectory$variable))
eq_trajectory$year_frac <- result$stan_data$year_frac[eq_trajectory$time_index]
eq_trajectory$education <- result$stan_data$education_levels[eq_trajectory$edu_index]

# Create plot data
eq_df <- data.frame(
  year_frac = eq_trajectory$year_frac,
  education = eq_trajectory$education,
  rate = eq_trajectory$median,
  type = "Equilibrium (u_eq[t])"
)

trend_df <- data.frame(
  year_frac = trend_data$year_frac,
  education = trend_data$education,
  rate = trend_data$mean,
  type = "Trend (equilibrium + shocks)"
)

full_df <- data.frame(
  year_frac = plot_data$year_frac,
  education = plot_data$education,
  rate = plot_data$mean,
  type = "Full Model (+ seasonal)"
)

decomp_data <- rbind(eq_df, trend_df, full_df)
decomp_data$type <- factor(decomp_data$type,
                           levels = c("Equilibrium (u_eq[t])", "Trend (equilibrium + shocks)", "Full Model (+ seasonal)"))

decomp_subset <- decomp_data[decomp_data$education %in% sample_edu, ]

ggplot(decomp_subset, aes(x = year_frac, y = rate * 100, color = type)) +
  geom_line(linewidth = 0.8) +
  geom_vline(xintercept = c(2008.75, 2020.25), linetype = "dashed",
             color = "gray40", alpha = 0.5) +
  facet_wrap(~education, scales = "free_y", ncol = 1) +
  scale_color_manual(values = c("Equilibrium (u_eq[t])" = "#27AE60",
                                "Trend (equilibrium + shocks)" = "#F39C12",
                                "Full Model (+ seasonal)" = "#3498DB")) +
  scale_x_continuous(breaks = seq(2000, 2025, 5)) +
  labs(x = "Year",
       y = "Unemployment Rate (%)",
       color = "Model Component",
       title = "Decomposition of Unemployment Dynamics",
       subtitle = "Green = time-varying equilibrium; Orange = equilibrium + shocks; Blue = full with seasonality") +
  theme(legend.position = "bottom",
        strip.text = element_text(size = 12, face = "bold"))

Model decomposition: equilibrium, trend (ODE + shocks), and full trajectory

Zoom: Seasonal Oscillation Detail (2022-2026)

A closer look at recent years to see the seasonal oscillation clearly:

# Zoom into 2022-2026 to see seasonal pattern clearly
zoom_years <- decomp_subset[decomp_subset$year_frac >= 2022 &
                             decomp_subset$year_frac <= 2026, ]

cat("Zoom rows:", nrow(zoom_years), "\n")

Zoom rows: 423 

cat("Types in zoom:", unique(as.character(zoom_years$type)), "\n")

Types in zoom: Equilibrium (u_eq[t]) Trend (equilibrium + shocks) Full Model (+ seasonal) 

ggplot(zoom_years, aes(x = year_frac, y = rate * 100, color = type)) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~education, scales = "free_y", ncol = 1) +
  scale_color_manual(values = c("Equilibrium (u_eq[t])" = "#27AE60",
                                "Trend (equilibrium + shocks)" = "#F39C12",
                                "Full Model (+ seasonal)" = "#3498DB")) +
  scale_x_continuous(breaks = seq(2022, 2026, 1),
                     minor_breaks = seq(2022, 2026, 0.25)) +
  labs(x = "Year",
       y = "Unemployment Rate (%)",
       color = "Component",
       title = "Seasonal Pattern Detail (2022-2026)",
       subtitle = "Blue oscillation around orange trend shows seasonal hiring patterns") +
  theme(legend.position = "bottom",
        strip.text = element_text(size = 12, face = "bold"),
        panel.grid.minor.x = element_line(color = "gray90"))