1 Introduction

The purpose of this project is to model and forecast U.S. Real Personal Consumption Expenditures (PCEC96) using Real Disposable Personal Income (DSPIC96) as the primary explanatory variable. The analysis will estimate the marginal propensity to consume (MPC), evaluate the stability of the income–consumption relationship over time, and assess the model’s forecasting performance using out-of-sample tests.

Specifically, I examine this consumption function from Keynesian consumption theory:

\[ C_t = \alpha + \beta Y_t + \varepsilon_t \] Where \(\beta\) is the Marginal Propensity to Consume (MPC).

This project combines benchmark time-series forecasting models, regression analysis, and other time-series methods. I began with the standard linear forecasting models, which include the mean, naive, drift, and a regression of consumption on disposable income. Additionally, I estimated multiple ARIMA specifications selected by AICc, and used an ETS exponential smoothing model as the other model. All models were evaluated by using both in-sample fit and out-of-sample forecast accuracy. The performance across these approaches determined the final model choice. Ultimately, I focused on the model that provides the most accurate and stable predictions over the short, medium, and long-term horizons.

2 Data

Previously, in the project proposal, I created a .csv file with both data sets. While creating that data set, I had to manually align the start and end dates of the data on the FRED website, then copy and paste it into a new .csv file, which became what I submitted.

A more efficient approach, which I have pursued here, is to connect the data to the R Markdown file using the FRED API. This way, I ensure the data aligns perfectly from a clear starting point. It also allows any change to be made efficiently in the future if needed.

## FRED API Key Recognition
fredr_set_key(Sys.getenv("FRED_API_KEY"))

## Retrieve consumption (PCEC96) and Income (DSPIC96) data from FRED
pcec96_raw <- fredr(
  series_id = "PCEC96",
  observation_start = as.Date("2007-01-01")
)

dspic96_raw <- fredr(
  series_id = "DSPIC96",
  observation_start = as.Date("2007-01-01")
)

## Merge data into one tsibble, arranged chronologically.
data <- pcec96_raw             |>
  select(date, consumption = value) |>
  full_join(
    dspic96_raw                |>
      select(date, income = value),
    by = "date"
  )                            |>
  arrange(date)                |>
  as_tsibble(index = date)

summary(data)
##       date             consumption        income     
##  Min.   :2007-01-01   Min.   :11068   Min.   :11957  
##  1st Qu.:2011-09-01   1st Qu.:11554   1st Qu.:12792  
##  Median :2016-05-01   Median :12870   Median :14135  
##  Mean   :2016-05-01   Mean   :13145   Mean   :14581  
##  3rd Qu.:2021-01-01   3rd Qu.:14205   3rd Qu.:16398  
##  Max.   :2025-09-01   Max.   :16575   Max.   :20520

2.1 Data Preparation

Using the FRED API, I was able to retrieve data on consumption and income at a specified starting point. I noticed while submitting the project proposal that the PCEC96 data did not start until 2007-01-01, which is far after the DSPIC96 data started.

When retrieving the data into my R Markdown file, I explicitly brought in data from both sets from that starting point. This way, I can compare both sets at an equal starting point. I then merged the income and consumption data into a single set, called ‘data’, for plotting.

2.1.1 Description of Variables

Real Personal Consumption Expenditures (PCEC96) measures the inflation-adjusted total amount of goods and services purchased by households in the United States.

  • Examples include spending on durable goods (e.g., cars, appliances), non-durable goods (e.g., food, clothing), and services (e.g., housing, healthcare, transportation).
  • Expressed in chained 2017 dollars - removes the effects of inflation and reflects changes in real consumption over time based on 2017 dollars.

Real Disposable Personal Income (DSPIC96) is the inflation-adjusted income available to households after taxes. It measures the real purchasing power that households have available after spending or saving.

  • Examples include wages, transfer payments, rental income, interest, and dividends, while subtracting personal taxes.
  • Expressed in chained 2017 dollars.

2.2 Time-Series Plot of Consumption and Income

## Visualization of the merged data set
ggplot(data, aes(x = date)) +
  geom_line(aes(y = consumption, color = "Consumption")) +
  geom_line(aes(y = income, color = "Income")) +
  scale_color_manual(
    name = "",
    values = c( "Consumption" = "red", "Income" = "blue")
  ) +
  labs(
    title = "Trends in Real Consumption and Real Disposable Income",
    x = "",
    y = "Consumption & Income",
    caption = "Consumption: Billions of Chained 2017 Dollars \nIncome: Billions of Chained 2017 Dollars"
  ) + theme(
    plot.caption = element_text(hjust = 0, size = 10)
  )

Both Real Personal Consumption Expenditures (consumption) and Real Disposable Personal Income (income) display a clear upward trend over time. The trend reflects long-term economic growth driven by rising household income, population growth, and productivity improvements. Consumption remains consistently above income, as expected, since PCEC96 includes aggregate consumer spending rather than per-capita measures.

2.2.2 Abnormalities

  • The COVID-19 shock was a dramatic irregular movement in both series that is visible around 2020, which implies that income and consumption were heavily affected by the pandemic.

3 Seasonality

The above time-series plot shows a clear upward trend; however, it is unclear whether seasonality is present. I will implement STL decomposition, a seasonal subseries plot using ggseasonplot, ACF for seasonal lags, and Kruskal-Wallis across months to gather evidence on the existence of seasonality.

3.1 STL Decomposition

## Creating consumption and income monthly ts objects to use in STL decomposition
cons_ts <- data$consumption |>
  ts(start = c(2007, 1), frequency = 12)

inc_ts <- data$income |>
  ts(start = c(2007, 1), frequency = 12)

## Consumption STL decomposition
const_stl <- stl(cons_ts, s.window = "periodic")

autoplot(const_stl) +
  ggtitle("STL Decomposition of Real Consumption (PCEC96)")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the forecast package.
##   Please report the issue at <https://github.com/robjhyndman/forecast/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Income STL decomposition
inc_stl <- stl(inc_ts, s.window = "periodic")

autoplot(inc_stl) +
  ggtitle("STL Decomposition of Real Disposable Personal Income (DSPIC96)")

Both series show a strong, smooth upward trend from 2007 to 2025. This indicates long-run economic growth, rising household real income, and a modest increase in consumption capacity. The trend component is only interrupted by the COVID-19 shock around 2020. The STL decomposition of income and consumption confirms non-stationarity. Seasonality is clear, consistent, and stable over time. Both income and consumption exhibit nearly identical repeating patterns every 12 months with predictable peaks and troughs at the same time periods per year, and similar shapes across the sample. Consumption seasonality is likely driven by predictable holiday and post-holiday spending, while income seasonality is likely driven by pay-period cycles, bonuses, tax timing, and benefit disbursements. Ultimately, both series are dominated by trend; however, seasonality is still present. Interestingly, the remainder plot only spikes in activity during the COVID-19 shock.

While STL decomposition extracts seasonal components, I will implement ggseasonplot() to overlay each year’s monthly path. We do this as STL decomposition compresses all years into a single smooth curve. Moreover, this is valuable for my analysis, as it will aid in choosing the appropriate forecasting model to use.

3.2 ggseasonplot()

## Consumption seasonal plot
ggseasonplot(cons_ts,
             year.labels = TRUE,
             year.labels.left = TRUE) +
  labs(title = "Seasonal Plot of Real Consumption (PCEC96)",
    x = "Month",
    y = "PCEC96") +
  theme(legend.position = "none")
## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • na.rm = TRUE
## ℹ Did you misspell an argument name?

## Income seasonal plot
ggseasonplot(inc_ts,
             year.labels = TRUE,
             year.labels.left = TRUE) +
  labs(title = "Seasonal Plot of Real Disposable Personal Income (DSPIC96)",
       x = "Month",
       y = "DSPIC96") +
  theme(legend.position = "none")
## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • na.rm = TRUE
## ℹ Did you misspell an argument name?

The purpose of using ggseasonplot() is to show evidence of strong, stable, and persistent monthly seasonality visually. The annual lines in both series have almost identical shapes, indicating stable seasonality, except for the pandemic years. Ultimately, the seasonal pattern remains consistent over time, with the timing of peaks and troughs remaining consistent across all years, and seasonality is an integral feature of both series.

3.3 Autocorrelation

## Income and consumption autocorrelation tests
par(mfrow = c(1, 2))

acf(cons_ts,
    lag.max = 48,
    main = "ACF of Real Consumption (PCEC96)")

acf(inc_ts,
    lag.max = 48,
    main = "ACF of Real Disposable Personal Income (DSPIC96)")

par(mfrow = c(1, 1))

The ACF plots for both series show an extremely strong and persistent autocorrelation. This is as expected when working with trending macroeconomic variables. Although at short lags (lags 1 to 6), both series exhibit auto-correlations near 1.0, indicating strong persistence, high memory, non-stationary behavior in levels, and the existence of a long-run trend. Lastly, the ACF confirms that both series exhibit stable and significant seasonal patterns that recur annually.

3.4 Kruskal-Wallis Test

Lastly, we will use a Kruskal-Wallis test across months to formally test seasonality. Since a strong trend is present, we must use the STL remainder to get an accurate p-value.

## Use STL to extract seasonal remainder components
cons_stl <- stl(cons_ts, s.window = "periodic")
inc_stl  <- stl(inc_ts,  s.window = "periodic")

## Extract remainder (trend + seasonal removed)
cons_rem <- cons_stl$time.series[, "remainder"]
inc_rem  <- inc_stl$time.series[, "remainder"]

## Extract months
cons_month <- cycle(cons_ts)
inc_month  <- cycle(inc_ts)

## Run KW tests on detrended, seasonally adjusted data
kw_cons_det <- kruskal.test(cons_rem ~ as.factor(cons_month))
kw_inc_det  <- kruskal.test(inc_rem  ~ as.factor(inc_month))

kw_cons_det
## 
##  Kruskal-Wallis rank sum test
## 
## data:  cons_rem by as.factor(cons_month)
## Kruskal-Wallis chi-squared = 53.204, df = 11, p-value = 1.648e-07
kw_inc_det
## 
##  Kruskal-Wallis rank sum test
## 
## data:  inc_rem by as.factor(inc_month)
## Kruskal-Wallis chi-squared = 41.776, df = 11, p-value = 1.771e-05

The Kruskal-Wallis results on the detrended STL residuals show that monthly medians differ significantly, with both p-values being less than 0.01. Once the long-run trend is removed, strong statistical evidence indicates stable and meaningful seasonality in both series.

4 Linear Models

Now we begin to use simple forecasts. We will use the last 24 months for testing.

4.1 Train/Test Split

## Number of observations
N <- length(cons_ts)

## Use last 24 months for test set
h <- 24

train_cons <- window(cons_ts, end = c(2023, 12))
test_cons  <- window(cons_ts, start = c(2024, 1))

train_inc  <- window(inc_ts,  end = c(2023, 12))
test_inc   <- window(inc_ts,  start = c(2024, 1))

4.2 Fit Benchmark Forecasting Models

## Mean Model
fit_mean <- meanf(train_cons, h = h)

## Naive Model
fit_naive <- naive(train_cons, h = h)

## Drift Model
fit_drift <- rwf(train_cons, drift = TRUE, h = h)

4.3 Fit Linear Regression Model (Consumption ~ Income)

## Linear regression model (consumption ~ income)

train_df <- tibble(
  cons = as.numeric(train_cons),
  inc  = as.numeric(train_inc)
)

test_df <- tibble(
  inc = as.numeric(test_inc)
)

fit_lm <- lm(cons ~ inc, data = train_df)

# Forecast using test income values
lm_forecast <- forecast(fit_lm, newdata = test_df)

4.4 Accuracy Tables

4.4.1 Out-of-Sample Accuracy Table

# Compare Forecast Accuracy using MAE and RMSE
acc_mean  <- accuracy(fit_mean,  test_cons)
acc_naive <- accuracy(fit_naive, test_cons)
acc_drift <- accuracy(fit_drift, test_cons)
acc_lm    <- accuracy(lm_forecast, test_cons)

# Combine into one table - extract only TEST rows (second row of each model)
accuracy_clean <- rbind(
  Mean       = acc_mean[2, ],
  Naive      = acc_naive[2, ],
  Drift      = acc_drift[2, ],
  Regression = acc_lm[2, ]
)

accuracy_clean <- as.data.frame(accuracy_clean)

knitr::kable(
  accuracy_clean,
  digits = 3,
  caption = "Out-of-Sample (Test-Set) Forecast Accuracy by Model"
)
Out-of-Sample (Test-Set) Forecast Accuracy by Model
ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Mean 3416.720 3424.673 3416.720 21.018 21.018 9.640 0.829 55.112
Naive 391.895 456.060 397.381 2.392 2.427 1.121 0.829 7.408
Drift 138.825 171.662 150.143 0.846 0.918 0.424 0.741 2.778
Regression 695.495 704.176 695.495 4.274 4.274 0.559 0.742 11.404

4.4.2 In-Sample Accuracy

accuracy_insample <- rbind(
  Mean       = acc_mean[1, ],
  Naive      = acc_naive[1, ],
  Drift      = acc_drift[1, ],
  Regression = acc_lm[1, ]
)

accuracy_insample <- accuracy_insample |>
  as.data.frame() |>
  tibble::rownames_to_column("Model")

knitr::kable(
  accuracy_insample,
  digits = 3,
  caption = "In-Sample (Train-Set) Forecast Accuracy by Model"
)
In-Sample (Train-Set) Forecast Accuracy by Model
Model ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
Mean 0.000 1439.974 1244.302 -1.215 9.649 3.511 0.980 NA
Naive 23.006 161.729 59.191 0.164 0.456 0.167 0.148 NA
Drift 0.000 160.084 53.491 -0.017 0.414 0.151 0.148 NA
Regression 0.000 522.909 294.728 -0.152 2.202 0.237 NA NA

4.4.3 Comparing Accuracy Tables and Assessing For Different Horizons

The in-sample and out-of-sample results show that different models perform best over different forecast horizons. In-sample, the naive and drift models fit the data extremely well, with RMSE values of 162 and 160, respectively. This is because real consumption is highly persistent and trending, making recent observations strong predictors of the recent past. The regression model fits reasonably well in-sample but exhibits larger training errors because it smooths out noise and explains consumption only through income, rather than perfectly tracking monthly fluctuations.

Out-of-sample performance tells a different story. Over the 24-month forecast horizon, the regression model produces substantially lower RMSE and MAE errors (699 and 690) than the mean and naive benchmarks and avoids the overly aggressive upward extrapolation seen in the drift model. The drift and naive models produce acceptable short-term forecasts but become increasingly optimistic over time and diverge from observed values. The mean model is the worst performer at all levels.

Ultimately, to forecast very short-term forecasts, use the naive model; for medium-term forecasts, use the regression model; and for forecasts far into the future, use the drift model.

4.5 Show Regression Model Parameters

## Show Model Parameters (regression model)
summary(fit_lm)
## 
## Call:
## lm(formula = cons ~ inc, data = train_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3865.6  -175.8    52.7   170.3  1028.3 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.110e+03  2.962e+02   7.125  1.8e-11 ***
## inc         7.524e-01  2.063e-02  36.467  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 525.5 on 202 degrees of freedom
## Multiple R-squared:  0.8681, Adjusted R-squared:  0.8675 
## F-statistic:  1330 on 1 and 202 DF,  p-value: < 2.2e-16

This regression output estimates this function:

\[\widehat{C}_t = 2110 + 0.752\, Y_t\]

The coefficient on income (0.752) is highly significant, indicating that a one-unit increase in real disposable income is associated with a 0.75-unit increase in real consumption, holding trend and seasonality effects implicit in the series constant. We have now estimated the marginal propensity to consume (MPC), which is equal to 0.752. The R-squared value is significantly high, meaning over 86% of the variation in consumption is explained by the income data alone. The F-stat is 1330, which is statistically significant, meaning the model explains more variation in consumption than a model with no predictors.

4.6 Create TS Object for LM Forecast and Plot All Forecasts

## Convert LM forecast to a time-series object
lm_forecast_ts <- ts(lm_forecast$mean, 
                     start = start(test_cons), 
                     frequency = 12)

## Forecast Plot
autoplot(cons_ts) +
  autolayer(fit_mean,  series = "Mean") +
  autolayer(fit_naive, series = "Naive") +
  autolayer(fit_drift, series = "Drift") +
  autolayer(lm_forecast_ts, series = "Regression") +
  labs(
    title = "Forecast Comparison for Real Consumption (PCEC96)",
    y = "PCEC96", x = "Year"
  ) +
  guides(color = guide_legend(title = "Models"))

The forecast comparison plot shows clear differences in how each model projects future consumption. - Naive: tracks the last observed value and produces stable short-run forecasts. - Drift: grows linearly and generates the most optimistic forecasts, widening as the horizon extends due to extrapolation of the long-run trend. - Mean: poor performs by flattening consumption at a level far below recent observations as it failed to account for trend or persistence. - Regression: incorporates disposable income, produces a smoother and more realistic forecast than the other models.

5 Forecasting Using ARIMA Models

We’ll estimate three models:

  • Random walk with drift
  • Non-seasonal ARIMA chosen by AICc
  • Seasonal ARIMA chosen by AICc

5.1 Fitting ARIMA Models

## Random walk with drift
fit_rw <- Arima(train_cons, order = c(0,1,0), include.drift = TRUE)

## Best non-seasonal ARIMA by AICc
fit_arima_ns <- auto.arima(
  train_cons,
  seasonal = FALSE,
  stepwise = FALSE,
  approximation = FALSE
)

## Best seasonal ARIMA by AICc
fit_arima_s <- auto.arima(
  train_cons,
  seasonal = TRUE,
  stepwise = FALSE,
  approximation = FALSE
)

fit_rw
## Series: train_cons 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##         drift
##       23.0064
## s.e.  11.2357
## 
## sigma^2 = 25754:  log likelihood = -1318.41
## AIC=2640.82   AICc=2640.88   BIC=2647.45
fit_arima_ns
## Series: train_cons 
## ARIMA(0,1,3) with drift 
## 
## Coefficients:
##          ma1      ma2      ma3    drift
##       0.1558  -0.3458  -0.1668  22.9199
## s.e.  0.0686   0.0686   0.0658   6.7451
## 
## sigma^2 = 22412:  log likelihood = -1302.95
## AIC=2615.9   AICc=2616.2   BIC=2632.46
fit_arima_s
## Series: train_cons 
## ARIMA(0,1,3) with drift 
## 
## Coefficients:
##          ma1      ma2      ma3    drift
##       0.1558  -0.3458  -0.1668  22.9199
## s.e.  0.0686   0.0686   0.0658   6.7451
## 
## sigma^2 = 22412:  log likelihood = -1302.95
## AIC=2615.9   AICc=2616.2   BIC=2632.46

5.2 Forecast for the test period

fc_rw        <- forecast(fit_rw,        h = h)
fc_arima_ns  <- forecast(fit_arima_ns,  h = h)
fc_arima_s   <- forecast(fit_arima_s,   h = h)

5.3 Accuracy: in-sample and out-of-sample

## In-sample (train)
acc_rw_train       <- accuracy(fc_rw)       # 1st row = train
acc_ns_train       <- accuracy(fc_arima_ns)
acc_s_train        <- accuracy(fc_arima_s)

acc_arima_train <- rbind(
  RW_Drift   = acc_rw_train[1, ],
  ARIMA_NS   = acc_ns_train[1, ],
  ARIMA_Seasonal = acc_s_train[1, ]
) |>
  as.data.frame() |>
  tibble::rownames_to_column("Model")

knitr::kable(
  acc_arima_train,
  digits = 3,
  caption = "In-Sample (Train-Set) Forecast Accuracy for ARIMA Models"
)
In-Sample (Train-Set) Forecast Accuracy for ARIMA Models
Model ME RMSE MAE MPE MAPE MASE ACF1
RW_Drift 0.055 159.693 53.283 -0.017 0.412 0.150 0.148
ARIMA_NS -0.068 147.859 56.420 -0.024 0.436 0.159 0.007
ARIMA_Seasonal -0.068 147.859 56.420 -0.024 0.436 0.159 0.007 
## Out-of-sample (test)
acc_rw_test      <- accuracy(fc_rw,       test_cons)
acc_ns_test      <- accuracy(fc_arima_ns, test_cons)
acc_s_test       <- accuracy(fc_arima_s,  test_cons)

acc_arima_test <- rbind(
  RW_Drift   = acc_rw_test[2, ],
  ARIMA_NS   = acc_ns_test[2, ],
  ARIMA_Seasonal = acc_s_test[2, ]
) |>
  as.data.frame() |>
  tibble::rownames_to_column("Model")

knitr::kable(
  acc_arima_test,
  digits = 3,
  caption = "Out-of-Sample (Test-Set) Forecast Accuracy for ARIMA Models"
)
Out-of-Sample (Test-Set) Forecast Accuracy for ARIMA Models
Model ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
RW_Drift 138.825 171.662 150.143 0.846 0.918 0.424 0.741 2.778
ARIMA_NS 175.681 204.787 182.947 1.073 1.119 0.516 0.725 3.324
ARIMA_Seasonal 175.681 204.787 182.947 1.073 1.119 0.516 0.725 3.324

The model with the lowest test RMSE / MAE is the best ARIMA forecaster. The in-sample results indicate that the ARIMA(0,1,3) models provide a better fit to the training data than the random walk with drift. Out-of-sample performance tells a different story. The random walk with drift clearly produces the most accurate forecasts for the 24-month test period, with the lowest RMSE (169) and MAE (147). The ARIMA(0,1,3) models perform worse out of sample, despite having a stronger in-sample fit. The results imply that the best ARIMA forecast model is the random walk with drift, especially for short-term and medium-term forecasts.

5.4 Plotting ARIMA Forecasts

I’m displaying the full forecast plot and a zoomed-in test period plot to see the differences more clearly.

## Full plot
autoplot(cons_ts) +
  autolayer(fc_rw,        series = "RW + Drift") +
  autolayer(fc_arima_ns,  series = "ARIMA(0,1,3)") +
  autolayer(fc_arima_s,   series = "ARIMA Seasonal") +
  labs(
    title = "ARIMA Forecast Comparison for Real Consumption (PCEC96)",
    x     = "Year",
    y     = "PCEC96"
  ) +
  guides(colour = guide_legend(title = "Models"))

## Zoomed-in
autoplot(window(cons_ts, start = c(2018, 1))) +
  autolayer(fc_rw,        series = "RW + Drift") +
  autolayer(fc_arima_ns,  series = "ARIMA(0,1,3)") +
  autolayer(fc_arima_s,   series = "ARIMA Seasonal") +
  labs(
    title = "ARIMA Forecasts for Real Consumption (Zoomed on Recent Years)",
    x     = "Year",
    y     = "PCEC96"
  ) +
  guides(colour = guide_legend(title = "Models"))

6 Error-Trend-Seasonal (ETS) Exponential Smoothing Model

I chose this model in addition to the simple and ARIMA forecasts because the data already shows a strong trend and stable monthly seasonality which the ETS Exponential Smoothing Model is designed for.

6.1 ETS Model Equation

The ETS model decomposes the series into level, trend, and seasonal components, more formally: ETS(A, A, A) – Additive Error, Additive Trend, Additive Seasonality.

Mathematically the equation looks like this:

\[\hat{y}_t = \ell_{t-1} + b_{t-1} + s_{t-m}\] ## ETS Model Execution

## Fit the model
fit_ets <- ets(train_cons)

fit_ets
## ETS(M,N,N) 
## 
## Call:
## ets(y = train_cons)
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 11179.2587 
## 
##   sigma:  0.0124
## 
##      AIC     AICc      BIC 
## 3153.383 3153.503 3163.337
## Forecast consumption using ETS
fc_ets <- forecast(fit_ets, h = h)

## Add ETS to ARIMA in-sample accuracy tables
acc_ets_train <- accuracy(fc_ets) # train is first row

acc_arima_train <- rbind(
  RW_Drift       = acc_rw_train[1, ],
  ARIMA_NS       = acc_ns_train[1, ],
  ARIMA_Seasonal = acc_s_train[1, ],
  ETS            = acc_ets_train[1, ]
) |>
  as.data.frame() |>
  tibble::rownames_to_column("Model")

knitr::kable(
  acc_arima_train,
  digits  = 3,
  caption = "In-Sample (Train-Set) Forecast Accuracy for ARIMA and ETS Models"
)
In-Sample (Train-Set) Forecast Accuracy for ARIMA and ETS Models
Model ME RMSE MAE MPE MAPE MASE ACF1
RW_Drift 0.055 159.693 53.283 -0.017 0.412 0.150 0.148
ARIMA_NS -0.068 147.859 56.420 -0.024 0.436 0.159 0.007
ARIMA_Seasonal -0.068 147.859 56.420 -0.024 0.436 0.159 0.007
ETS 22.904 161.335 58.911 0.163 0.454 0.166 0.148 
## Add ETS to ARIMA out-of-sample accuracy tables
acc_ets_test <- accuracy(fc_ets, test_cons)

acc_arima_test <- rbind(
  RW_Drift       = acc_rw_test[2, ],
  ARIMA_NS       = acc_ns_test[2, ],
  ARIMA_Seasonal = acc_s_test[2, ],
  ETS            = acc_ets_test[2, ]
) |>
  as.data.frame() |>
  tibble::rownames_to_column("Model")

knitr::kable(
  acc_arima_test,
  digits  = 3,
  caption = "Out-of-Sample (Test-Set) Forecast Accuracy for ARIMA and ETS Models"
)
Out-of-Sample (Test-Set) Forecast Accuracy for ARIMA and ETS Models
Model ME RMSE MAE MPE MAPE MASE ACF1 Theil’s U
RW_Drift 138.825 171.662 150.143 0.846 0.918 0.424 0.741 2.778
ARIMA_NS 175.681 204.787 182.947 1.073 1.119 0.516 0.725 3.324
ARIMA_Seasonal 175.681 204.787 182.947 1.073 1.119 0.516 0.725 3.324
ETS 391.902 456.065 397.387 2.392 2.427 1.121 0.829 7.408 
## Plot ETS vs. best ARIMA
autoplot(cons_ts) +
  autolayer(fc_arima_ns, series = "ARIMA(0,1,3)") +
  autolayer(fc_ets,      series = "ETS") +
  labs(
    title = "ARIMA(0,1,3) and ETS Forecasts for Real Consumption (PCEC96)",
    x     = "Year",
    y     = "PCEC96"
  ) +
  guides(colour = guide_legend(title = "Models"))

6.2 Accuracy Table Results

The ETS model performs worst overall, which could suggest that exponential smoothing is too restrictive for the stochastic trend present in consumption. Out of the options presented in the table, it is implied that the Random Walk with Drift provides the most reliable medium-term forecasts for real consumption.

7 Final Forecast

Based on the in-sample and out-of-sample forecast accuracy tables, the regression model is the best overall forecasting model, which confirms the original question of this analysis. It is the best because it has the lowest RMSE and MAE in the out-of-sample models, is economically meaningful by incorporating the theoretical driver of consumption, and is the most stable medium- and long-horizon forecaster. ARIMA relies heavily on past fluctuations, whereas regression generalizes the past, which is beneficial in this case.

7.1 Forecasts for 1, 2, 5, 10, 20, 30, 40, 50 Months Ahead

# Forecast for 50 months ahead using regression model
fc_reg <- forecast(fit_lm, newdata = data.frame(inc = tail(data$income, 50)), h = 50)

# Extract selected horizons
horizons <- c(1, 2, 5, 10, 20, 30, 40, 50)
reg_forecast_values <- fc_reg$mean[horizons]

# Put into a table
forecast_table <- tibble(
  Horizon = horizons,
  Forecast = as.numeric(reg_forecast_values)
)

knitr::kable(
  forecast_table,
  digits = 2,
  caption = "Regression Model Forecasts for Selected Horizons"
)
Regression Model Forecasts for Selected Horizons
Horizon Forecast
1 14837.97
2 14696.07
5 14607.52
10 14225.47
20 14974.68
30 15333.25
40 15540.91
50 15740.06

7.2 Plotting Final Forecast (Regression Model)

full_df <- tibble(
  cons = as.numeric(cons_ts),
  inc  = as.numeric(inc_ts)
)

fit_lm_full <- lm(cons ~ inc, data = full_df)

fit_inc_full <- auto.arima(inc_ts)
h_final      <- 50  # or however many months ahead you want

fc_inc_full <- forecast(fit_inc_full, h = h_final)
future_inc  <- as.numeric(fc_inc_full$mean)

fc_reg_full <- forecast(
  fit_lm_full,
  newdata = data.frame(inc = future_inc),
  h = h_final
)

last <- end(cons_ts)  # e.g. c(2025, 9)

next_period <- if (last[2] == 12) c(last[1] + 1, 1) else c(last[1], last[2] + 1)

fc_reg_full_ts <- ts(
  fc_reg_full$mean,
  start     = next_period,
  frequency = 12
)

autoplot(cons_ts) +
  autolayer(fc_reg_full_ts, series = "Regression Forecast") +
  labs(
    title = "Final Regression Model Forecast for Real Consumption (PCEC96)",
    x = "Year",
    y = "PCEC96"
  ) +
  guides(colour = guide_legend(title = "Model"))

8 Conclusion

The goal of the project was to model and forecast U.S. real personal consumption expenditures using real disposable personal income as the primary explanatory variable. The analysis evaluated the in-sample and out-of-sample predictive accuracy of several forecasting methods used to determine which approach best captures and forecasts consumption behavior over time.

The regression model provided the most accurate medium-run forecasts and produced a meaningful estimate of the marginal propensity to consume. The random walk with drift was the strongest time-series predictor, but was not chosen because it ignores the fundamental economic driver of consumption. The purpose of the project was to estimate the marginal propensity to consume (MPC) and the broader income-consumption relationship, in addition to forecasting consumption, which a random walk with drift cannot do. Despite strong ARIMA performance in shorter horizons, the regression model was chosen as it was the more appropriate and theoretically meaningful choice for the final long-horizon consumption forecast.