Midterm

1. Select a Monthly Time Series

Name of series: Unemployment Rate
FRED Series Identifier: UNRATENSA
BLS Source Code: LNU04000000
Frequency: Monthly
Units of measurement: Percent, Not Seasonally Adjusted
Period of Data: 01/1997 - 12/2005

# Relative path: the CSV must sit in the same folder as this .Rmd file
df <- read.csv("UNRATENSA.csv")

Data-access statement/ Reason I selected this data

I downloaded the CSV data from the Federal Reserve Economic Data database for series UNRATENSA. https://fred.stlouisfed.org/series/unratensa

library(fpp3)
library(patchwork)

2. Protect the Test Period

# Convert data.frame into a tsibble
df <- df |>
  mutate(Date = as.Date(observation_date)) |>  
  mutate(Month = yearmonth(Date)) |>
  as_tsibble(index = Month) |>
  filter_index("1990 Jan" ~ "2005 Dec")

# Training Data 
train <- df |> filter_index("1990 Jan" ~ "2002 Dec")

# Test Set 
test <- df |> filter_index("2003 Jan" ~ "2005 Dec")

nrow(train)  # should be 156

## [1] 156

nrow(test)   # should be 36

## [1] 36

The split is 156/36 (81%/19%) rather than exactly 154/38 so that the training and test periods both cover complete calendar years.

3. Diagnose the Training Series

# plot data
train |>
  autoplot(UNRATENSA) +
  labs(title = "U.S. Unemployment Rate",x = "Year",y = "Percent")

When examining this plot, the trend is gradual rises to above 8% in 1992 and then it falls to below 4% by 2000 and then continues upward again in 2001. There is roughly constant seasonality as well.

train |>
  gg_season(UNRATENSA, labels = "both") +
  labs(title = "Seasonal Plot: U.S. Unemployment Rate (NSA)",y = "Percent")

There is some seasonality here as unemployment tends to dip in early spring, rise into summer, and decline again during autumn. The highest years of average unemployment were 1992-1993.

train |> gg_tsdisplay(UNRATENSA, plot_type = "partial", lag_max = 36)

# Box-Cox lambda to consider a transformation
train |>
  features(UNRATENSA, features = guerrero)

A transformation like log or Box-Cox is not needed for this data. Transformations are used when the size of the seasonal swings gets bigger or smaller as the overall level changes. In our data, the seasonal ups and downs stay roughly the same size (approx 1 to 1.5 percentage points) whether unemployment is high or low. Keeping the data in its original percentage format also makes the model’s results easy to interpret without needing a complex back-transformations.

4. Establish a Benchmark

# Seasonal Naive (SNAIVE)
fit_benchmark <- train |>
  model(SNaive = SNAIVE(UNRATENSA))

# Residual diagnostics
fit_benchmark |>
  gg_tsresiduals() +
  labs(title = "ETS(A,Ad,A) Residuals")

# Ljung-Box test
fit_benchmark |>
  augment() |>
  features(.innov, ljung_box, lag = 24)

The Seasonal Naive (SNAIVE) model was chosen as the benchmark because the U.S. Unemployment Rate exhibits strong annual seasonality (ex. recurring surge in June/July due to academic graduates and seasonal workers entering the workforce), a simple Naive or Mean forecast would perform poorly. The SNAIVE method sets each future month’s forecast equal to the observed value from the exact same month of the previous year.

5. Fit an ETS Model

# Automatically selected ETS model 
fit_auto <- train |>
  model(ETS(UNRATENSA))
report(fit_auto)

## Series: UNRATENSA 
## Model: ETS(A,Ad,A) 
##   Smoothing parameters:
##     alpha = 0.6941198 
##     beta  = 0.1781651 
##     gamma = 0.0001681018 
##     phi   = 0.8377267 
## 
##   Initial states:
##      l[0]       b[0]       s[0]      s[-1]      s[-2]      s[-3]       s[-4]
##  5.369058 0.06988356 -0.3514921 -0.3193674 -0.4029559 -0.2343927 -0.06383893
##      s[-5]     s[-6]      s[-7]      s[-8]     s[-9]    s[-10]    s[-11]
##  0.1642968 0.1846738 -0.2197351 -0.1813709 0.2933087 0.4929547 0.6379189
## 
##   sigma^2:  0.0251
## 
##      AIC     AICc      BIC 
## 230.6481 235.6408 285.5455

# Residual diagnostics
fit_auto |>
  gg_tsresiduals() +
  labs(title = "ETS(A,Ad,A) Residuals")

# Ljung-Box test
fit_auto |>
  augment() |>
  features(.innov, ljung_box, lag = 24)

The model selected ETS(A,Ad,A): this model has additive errors, an additive damped trend, and additive seasonality. Additive errors and seasonality fit this series because the seasonal swings stay about the same size at any level of unemployment. The damped trend makes sense as over time without one the unemployment rate would cross zero and become negative which is not possible.

The Ljung-Box test gives p = 0.416, so no significant autocorrelation is left in the residuals. The ETS model captured the trend and seasonal structure of the training data.

6. Fit an ARIMA Model

fit_arima <- train |> model(ARIMA(UNRATENSA))
report(fit_arima)

## Series: UNRATENSA 
## Model: ARIMA(1,1,3)(2,1,0)[12] 
## 
## Coefficients:
##           ar1      ma1     ma2     ma3     sar1     sar2
##       -0.0944  -0.0360  0.0871  0.2217  -0.4733  -0.2370
## s.e.   0.3434   0.3378  0.0926  0.0871   0.0902   0.0964
## 
## sigma^2 estimated as 0.03141:  log likelihood=45.85
## AIC=-77.7   AICc=-76.87   BIC=-56.96

# Residual diagnostics
fit_arima |>
  gg_tsresiduals() +
  labs(title = "ARIMA(1,1,3)(2,1,0)[12] Residuals")

# Ljung-Box test 
fit_arima |>
  augment() |>
  features(.innov, ljung_box, lag = 24, dof = 6)

The model selected ARIMA(1,1,3)(2,1,0)[12]. One regular difference (d = 1) removes the wandering level, and one seasonal difference (D = 1) removes the yearly pattern. The model has one AR term, three MA terms,and two seasonal AR terms that link each month to the same month one and two years earlier. Overall this is consistent with a series that has both a changing level and strong seasonality.

With 6 degrees of freedom removed for the estimated parameters, the Ljung-Box test gives Q = 22.84 with p = 0.197. The innovations look like white noise, so the model captured the temporal structure of the training data.

7. Fit a Times-Series Regression Model

fit_regression <- train |>
  model(
    TSLM_Trend_Season = TSLM(UNRATENSA ~ trend() + season())
  )

# Report estimated coefficients and model fit statistics
fit_regression |> report()

## Series: UNRATENSA 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5129 -0.6191 -0.1494  0.5633  1.6129 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     7.367579   0.249338  29.549  < 2e-16 ***
## trend()        -0.016205   0.001451 -11.167  < 2e-16 ***
## season()year2  -0.114564   0.319211  -0.359  0.72020    
## season()year3  -0.306052   0.319221  -0.959  0.33930    
## season()year4  -0.774462   0.319237  -2.426  0.01651 *  
## season()year5  -0.804411   0.319260  -2.520  0.01285 *  
## season()year6  -0.388206   0.319290  -1.216  0.22605    
## season()year7  -0.410462   0.319326  -1.285  0.20073    
## season()year8  -0.617334   0.319369  -1.933  0.05521 .  
## season()year9  -0.754976   0.319418  -2.364  0.01944 *  
## season()year10 -0.877232   0.319474  -2.746  0.00681 ** 
## season()year11 -0.768719   0.319537  -2.406  0.01742 *  
## season()year12 -0.752514   0.319606  -2.355  0.01991 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8138 on 143 degrees of freedom
## Multiple R-squared: 0.5115,  Adjusted R-squared: 0.4705
## F-statistic: 12.48 on 12 and 143 DF, p-value: < 2.22e-16

# Residual diagnostics
fit_regression |>
  gg_tsresiduals() +
  labs(title = "Time-Series Regression Residuals (Trend + Season)")

# Ljung-Box test
fit_regression |>
  augment() |>
  features(.innov, ljung_box, lag = 24)

The model uses two predictors, a time trend to show the long-term direction of unemployment, and month dummies to show the yearly seasonal pattern, with January as the base month. I picked these because the plots in Section 3 showed a trend and strong seasonality. The trend shows unemployment falling about 0.016 points per month, or about 0.19 points per year. Every month is lower than January. October (−0.88) and May (−0.80) are the lowest, and June is about 0.42 points higher than May because graduates and students start looking for work in the summer.

The Ljung-Box test has at p < 0.001, so the residuals are strongly autocorrelated. That means the model’s prediction intervals are unreliable and its trend is likely to predict poorly.

8. Construct a Simple Forecast Combination

fit_combo <- train |>
  model(
    combination = (ETS(UNRATENSA) +
                   ARIMA(UNRATENSA) +
                   TSLM(UNRATENSA ~ trend() + season())) / 3)

# 2. Forecast
fc_combo <- fit_combo |> forecast(h = 36)

# 3. Plot the combined forecast
fc_combo |>
  autoplot(df) +
  labs(title = "Equal-Weight Combined Forecast (ETS + ARIMA + Regression)",
       y = "UNRATENSA")

# 4. View the combined forecast values
fc_combo

9. Compare 5 Forecasting Approaches

# Fit all five approaches on the training data
fit_all <- train |>
  model(
    benchmark   = SNAIVE(UNRATENSA),
    ets         = ETS(UNRATENSA),
    arima       = ARIMA(UNRATENSA),
    regression  = TSLM(UNRATENSA ~ trend() + season())
  ) |>
  mutate(combination = (ets + arima + regression) / 3)

# Forecast the same test observations
fc_all <- fit_all |> forecast(new_data = test)

# Plot all five forecasts against the actual data
fc_all |>
  autoplot(df, level = NULL) +
  labs(title = "Comparison of Five Forecasting Approaches",
       y = "UNRATENSA")

fc_all |>
  autoplot(test, level = NULL) +
  labs(title = "Forecasts vs Actual: Test Period (2003–2005)",
       y = "UNRATENSA")

Evaluate Out-of-Sample Performance

# Compare test-set accuracy, round 3 decimal places
fc_all |>
  accuracy(df) |>
  select(.model,MAE, RMSE, MASE, ME) |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  arrange(RMSE)

# Residual plots (residuals over time, ACF, histogram)
fit_all |> select(benchmark) |> gg_tsresiduals() + labs(title = "Seasonal Naive residuals")

fit_all |> select(ets) |> gg_tsresiduals() + labs(title = "ETS residuals")

fit_all |> select(arima) |> gg_tsresiduals() + labs(title = "ARIMA residuals")

fit_all |> select(regression) |> gg_tsresiduals() + labs(title = "Regression residuals")

fit_all |> select(combination) |> gg_tsresiduals(type = "response") +
  labs(title = "Combination residuals")

# Ljung-Box tests (p < 0.05 = autocorrelation left in residuals)
fit_all |>
  select(benchmark, ets, regression) |>
  augment() |>
  features(.innov, ljung_box, lag = 24)

fit_all |>
  select(arima) |>
  augment() |>
  features(.innov, ljung_box, lag = 24, dof = 6)

fit_all |>
  select(combination) |>
  augment() |>
  features(.resid, ljung_box, lag = 24)

11. Interpret the Results

Which model performs best overall: The seasonal naive benchmark performed best on every metric (MAE 0.406, RMSE 0.485, MASE 0.670). None of the more complex models beat it on the test period.

Does the ranking change across metrics: The top three are the same on MAE, RMSE and MASE: benchmark, combination, ETS. The bottom two switch places as regression has the lower RMSE (1.56 vs 1.72), but ARIMA has the lower MAE (1.42 vs 1.54) and MASE (2.34 vs 2.54).

Does either model show systemic bias: The regression has a large positive bias (+1.54), so it under forecast by about 1.5 points on average. Its linear trend (−0.19 points per year) was estimated mostly from the 1992–2000 decline, so it kept pushing unemployment down after the trend had reversed. ARIMA (−1.41) and ETS (−0.73) over-forecast. Both picked up the 2001–2002 recession rise and carried it forward, but unemployment peaked in mid-2003 and then fell. The benchmark (−0.24) and combination (−0.20) have little bias.

How much does the best model improve on the benchmark? None. The best non benchmark method, the combination, had a 21% higher RMSE than seasonal naive (0.589 vs 0.485) and a 24% higher MAE.

Does the forecast combination improve on its individual components: The combination beat all three of its components. Its RMSE was lower than the best single component, ETS (0.589 vs 0.854). This happened because the components’ errors went in opposite directions: the regression was too low and ARIMA/ETS were too high, so averaging them cancelled each other out. This is the main reason to combine forecasts.

Did residual diagnostics correspond with out of sample performance: Residual diagnostics did not line up well with out-of-sample performance. ETS (p = 0.416) and ARIMA (p = 0.197) were the only models whose residuals passed the Ljung-Box test, yet both forecast worse than the benchmark. ARIMA had the worst RMSE of all five methods. The seasonal naive benchmark failed the test badly (Q = 918.94, p < 0.001) but had the best test accuracy. The regression was the one case where they agreed: it had heavily autocorrelated residuals (Q = 979.39) and large, biased forecast errors. The combination also failed (Q = 649.97) because it averages in the regression’s autocorrelated residuals, but it still forecast better than any of its components. Passing a residual test means a model fits the training data well. It does not mean the model will forecast well when the series changes direction after the forecast origin. In this case, ETS and ARIMA fit the 2001–2002 rise closely and extrapolated it, while unemployment peaked in 2003 and fell.

What characteristics of the series appear to explain the winning model’s performance? Unemployment has a stable seasonal pattern but no lasting trend. It follows the business cycle, which none of these models can predict. Over 2003–2005 the level changed slowly, so “same month last year” was a good guess. The models that extrapolated a trend (regression, ARIMA, ETS) were hurt when the cycle turned

12. Select and Defend the Final Forecasting Approach

I recommend the seasonal naive method. It had the lowest MAE, RMSE and MASE, almost no bias, and no parameters to estimate. It needs no future predictor values and is easy to explain to a business audience. The more complex models did not add enough accuracy to justify their extra structure. Each of them was worse than the benchmark on the test set. The benchmark’s failed residual test means its prediction intervals are too narrow, so its point forecasts can be trusted more than its intervals.

The equal-weight combination is a reasonable second choice. It had low bias and was much more accurate than any single model. If the organization expects a turning point in the business cycle, the combination may be more robust than relying on one specification. One limit of seasonal naive is that it can’t anticipate a change in level, so its prediction intervals should be treated as a rough guide during substantial changes in employment.

Visuals for poster

# Forecasts vs actual
fc_all |>
  autoplot(df |> filter_index("2001 Jan" ~ .), level = NULL) +
  labs(title = "Forecasts vs Actual, 2003–2005", x = NULL, y = "Percent")

# Data with test period marked
df |> autoplot(UNRATENSA) +
  geom_vline(xintercept = as.Date("2003-01-01"), linetype = "dashed") +
  labs(title = "U.S. Unemployment Rate (NSA), 1990–2005", x = NULL, y = "Percent")

# Benchmark with prediction intervals
fc_all |> filter(.model == "benchmark") |>
  autoplot(df |> filter_index("2001 Jan" ~ .)) +
  labs(title = "Seasonal Naive Forecast with 80% and 95% Intervals", x = NULL, y = "Percent")