suppressPackageStartupMessages(library(forecast))
library(readxl)
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(dynlm))
suppressPackageStartupMessages(library(ggfortify))
suppressPackageStartupMessages(library(ggplot2))
library(stats)
suppressPackageStartupMessages(library(qcc))
suppressPackageStartupMessages(library(Metrics))
suppressPackageStartupMessages(library(vars))
suppressPackageStartupMessages(library(prophet))
library(tseries)
library(strucchange)
library(lmtest)
library(urca)

Introduction:

In this project, I will attempt to summarize the main concepts I learned in Econ 144, Economic Forecasting. This is mainly a reflection of the three projects we were assigned during the quarter, which included concepts like test of stationarity, ACF and PACF, trend and seasonal decomposition, utilizing forecasting methods, combining methods, and multivariate forecasting. The goal of this project is not necessarily to show the coding aspect of how I thought of all my models (as most of the code was generated through my previous projects and ChatGPT), but rather to show and analyze what the code was doing and how it is relevant to my goals of time series forecasting.

Part 1) Importing and Cleaning the Data

We start off by importing the two of the datasheets found from the Bureau of Labor Statistics. These data include unemployment, which is the main variable I will focus on, as well as inflation, which will only be used during multivariate forecasting in VAR models.

data1 <- read_excel("~/Downloads/UNRATE (2).xlsx", 
    sheet = "Monthly")
head(data1)
## # A tibble: 6 × 2
##   observation_date    UNRATE
##   <dttm>               <dbl>
## 1 1968-01-01 00:00:00    3.7
## 2 1968-02-01 00:00:00    3.8
## 3 1968-03-01 00:00:00    3.7
## 4 1968-04-01 00:00:00    3.5
## 5 1968-05-01 00:00:00    3.5
## 6 1968-06-01 00:00:00    3.7
tail(data1)
## # A tibble: 6 × 2
##   observation_date    UNRATE
##   <dttm>               <dbl>
## 1 2025-02-01 00:00:00    4.1
## 2 2025-03-01 00:00:00    4.2
## 3 2025-04-01 00:00:00    4.2
## 4 2025-05-01 00:00:00    4.2
## 5 2025-06-01 00:00:00    4.1
## 6 2025-07-01 00:00:00    4.2
data2 <- read_excel("~/Downloads/CORESTICKM159SFRBATL.xlsx", 
    sheet = "Monthly")
head(data2)
## # A tibble: 6 × 2
##   observation_date    CORESTICKM159SFRBATL
##   <dttm>                             <dbl>
## 1 1968-01-01 00:00:00                 3.65
## 2 1968-02-01 00:00:00                 3.67
## 3 1968-03-01 00:00:00                 4.14
## 4 1968-04-01 00:00:00                 4.16
## 5 1968-05-01 00:00:00                 4.09
## 6 1968-06-01 00:00:00                 4.55
tail(data2)
## # A tibble: 6 × 2
##   observation_date    CORESTICKM159SFRBATL
##   <dttm>                             <dbl>
## 1 2025-02-01 00:00:00                 3.52
## 2 2025-03-01 00:00:00                 3.26
## 3 2025-04-01 00:00:00                 3.18
## 4 2025-05-01 00:00:00                 3.16
## 5 2025-06-01 00:00:00                 3.29
## 6 2025-07-01 00:00:00                 3.42
data3 <- merge(data1, data2, by = "observation_date")

As we learned, the default class of the data we import is not in the time-series class, so we must create one ourselves before we do any forecasting.

is.ts(data3)
## [1] FALSE
ts_data <- ts(data3[, -1], start = c(1968, 1), end = c(2025, 7),frequency = 12)
unemployment <- ts_data[, "UNRATE"]
inflation <- ts_data[, "CORESTICKM159SFRBATL"]

This is a general overview of unemployment and inflation from the late 1960s to now. The movements seem very cyclical, with the two variables moving in a similar pattern. It is also important to notice the sudden increase in both variables during 2020. We will see how this can potentially affect our analysis and forecasting.

ts.plot(ts_data[,1], ts_data[,2],
        col = c("blue", "red"),
        lty = 1:1,
        xlab = "Time",
        ylab = "Numbers in Percentages",
        main = "Unemployment vs Inflation")
legend("topright", legend = c("Unemployment", "Inflation"),
       col = c("blue", "red"), lty = 1)

Part Two) Testing for Stationarity

The ADF test is a very familar test we have used in previous classes when testing if a data is stationary. I have also included the KPSS test to ensure that our results are accurate. Correcting the data for stationarity is necessary because it makes our data’s structure more consistent and stable over time. Here the H0 is rejected for both the ADF and KPSS test. However, it is important to know that the H0 for both tests are actually the opposites, the null for ADF says that the series is non-stationary, and the null for KPSS says that the series is stationary. Therefore only the KPSS test says that our data for unemployment is not stationary, but ADF test barely rejects the null as well. We should probably adjust the model in some way:

adf.test(unemployment)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  unemployment
## Dickey-Fuller = -3.4304, Lag order = 8, p-value = 0.049
## alternative hypothesis: stationary
kpss.test(unemployment)
## 
##  KPSS Test for Level Stationarity
## 
## data:  unemployment
## KPSS Level = 0.73896, Truncation lag parameter = 6, p-value = 0.01
diff_unemployment <- diff(unemployment)
adf.test(diff_unemployment)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_unemployment
## Dickey-Fuller = -9.0868, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(diff_unemployment)
## 
##  KPSS Test for Level Stationarity
## 
## data:  diff_unemployment
## KPSS Level = 0.059506, Truncation lag parameter = 6, p-value = 0.1

This is the visual data after we differentiate the data of unemployment. I have made the data more stable by differentiating it, but I could have utilized other methods to make the movements more stable, such as logging (shrinking) the data to reduce the raw numbers. Anyway, the both the ADF and KPSS test has confirmed that the new data called diff_unemployment is now stationary.

ggtsdisplay(diff_unemployment, main = "Overview of Unemployment Data", )

Briefly, we also learned the concept of burst movements by analyzing what the CUSUM test does. The CUSUM test checks for sudden burst movements within the data that can potentially harm our ability to later analyze and forecast our data. We already made a remark of how there was a sudden spike of unemployment during covid. Although we can also see the spike here, we would only worry about it if the spike crosses the absolute value threshold that is displayed on the chart (over the absolute value of one). Based on the formalities of this test, we shouldn’t have to worry about it.

cusum <- efp(diff_unemployment ~ 1, type = "OLS-CUSUM")
plot(cusum)

Part Three) Fitting an Estimated Model to Analyze Residuals and Variance

We will now attempt to fit both a linear and non-linear model that best fits the past values of our unemployment data. In other words, we will do a regressive analysis of our residuals and variance to see our properties of the data, and to see if any of the data’s components needs a change. Although this may not directly impact our later work of forecasting the model, it is imperative that we are familiar with the properties of the data we are working with so that we understand our forecasting results better.

time_index <- 1:length(diff_unemployment)

df <- data.frame(
  t = time_index,
  y = as.numeric(diff_unemployment)
)


fit_lin <- lm(y ~ t, data = df)

fit_quad <- lm(y ~ poly(t, 2, raw = TRUE), data = df)

fit_cub <- lm(y ~ poly(t, 3, raw = TRUE), data = df)

fit_quart <- lm(y ~ poly(t, 4, raw = TRUE), data = df)

models <- list(Linear = fit_lin, Quadratic = fit_quad, Cubic = fit_cub, Quartic = fit_quart)

for (name in names(models)) {
  cat("\nModel:", name, "\n")
  print(summary(models[[name]])$adj.r.squared)   
  print(AIC(models[[name]]))                     
}
## 
## Model: Linear 
## [1] -0.0009457434
## [1] 872.2577
## 
## Model: Quadratic 
## [1] -0.002271492
## [1] 874.1674
## 
## Model: Cubic 
## [1] -0.003450984
## [1] 875.9738
## 
## Model: Quartic 
## [1] -0.004896096
## [1] 877.9602
par(mfrow = c(2,2))  
for (name in names(models)) {
  fit <- models[[name]]
  plot(df$t, df$y, main = paste(name, "Fit"), xlab = "Time", ylab = "diff_unemployment", pch = 20)
  lines(df$t, fitted(fit), col = "blue", lwd = 2)
}

par(mfrow = c(2,2))  
for (name in names(models)) {
  fit <- models[[name]]
  hist(residuals(fit), main = paste(name, "Residuals"), xlab = "Residuals", col = "lightgray", border = "white")
}

I have modeled four different fits of our unemployment data, including a linear one, as well as a qudratic, cubic and quartic fit. Visually, it is difficult to see which fit is the most accurate, mainly because the data is differentiated and so the changes are too stablized. We are going to need to conduct a formal test of error to see which one is the most preferred:

check_residuals <- function(models, y) {
  res <- residuals(models)
  rmse_val <- rmse(y, fitted(models))
  mae_val  <- mae(y, fitted(models))
  return(c(RMSE = rmse_val, MAE = mae_val))
}

results <- sapply(models, check_residuals, y = df$y)
results
##         Linear Quadratic     Cubic   Quartic
## RMSE 0.4532935 0.4532638 0.4532002 0.4531957
## MAE  0.1531012 0.1533018 0.1531933 0.1530419

Unsurprisingly the Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE), which are the two most common error metrics, show marginal differences between all four of the fits. For simplicity, I will go with the linear fit.

The Ljung box test is used to see if the data’s values are autocorrelated. In other words, we must see if the numbers are independently distributed, to ensure that the data is random/has white noise. In the two results shown above, it is clear that with a lag value of 4, the H0 of no autocorrelation is rejected, which means there is autocorrelation. However, after we add a lag, the H0 is accepted. This means that in order to forecast our data, we must use a lag value of 5 or more to avoid autocorrelation. If we use a lag value of 4 or lower, our new values are too dependent on our past values (in this case 4 or less) and would not be accurate.

Box.test(residuals(fit_lin), lag = 4, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(fit_lin)
## X-squared = 9.7767, df = 4, p-value = 0.04436
Box.test(residuals(fit_lin), lag = 5, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  residuals(fit_lin)
## X-squared = 9.8154, df = 5, p-value = 0.08064

The Shapiro-Wilk test checks for the normality of our residuals/errors. Because it only checks our residuals, a non-normal pattern of residuals will not affect our point forecasts, but can affect the interval forecasts that we ensure for safety. If we have non-normal residuals, it means that there are extreme values/outliers within the data. Our H0 of normality is rejected, so we now make a new object called diff_unemployment 2 where we rescale the unemployment numbers so that the variance can be stabilized and in Gaussian form.

shapiro.test(residuals(fit_lin))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fit_lin)
## W = 0.29416, p-value < 2.2e-16
lambda <- BoxCox.lambda(diff_unemployment)
diff_unemployment2 <- BoxCox(diff_unemployment, lambda)

The final test we will conduct is for heteroskedasticity, which means a non-constant variance over time. Heteroskedasticity within the data means that there are errors and that our results may be biased, so we must change the error variance model by using GARCH/robust forecast intervals. Fortunately, the H0 for the BP test is not rejected, meaning that our data is homeskedastistic.

bptest(fit_lin)
## 
##  studentized Breusch-Pagan test
## 
## data:  fit_lin
## BP = 2.3767, df = 1, p-value = 0.1232

Another concept we learned is being able to decompose our data into different three different components: the trend, seasonality, and the remainder/noise. The trend is the most generalized movement of the data. Because we differentiated the data, the trend is mostly flat. The seasonality is very consistent, almost looking like a radio sound signal wave. Because the swings are mostly the same, we conclude that an additive decomposition is preferred. Finally, there is very little remainder left after dissecting the trend and seasonality, meaning that the trend and seasonality explains most of the movements within the data.

autoplot(decompose(diff_unemployment2), type = c("additive"))

autoplot(decompose(diff_unemployment2), type = c("multiplicative"))

Part Four) Forecasting

We finally arrive to the main part of our project, which is forecasting. Here are the main forecast methods we learned, and how they operate:

ARIMA: Uses the auto regressive (AR), moving average (MA), as well as the integration aspect to forecast. The forecast relies on its past values and errors to make predictions.

Holt-Winters: Utilizes the trend, seasonality, and remainder (what we did in the decomposition). This method may be less flexible than ARIM due to relying on consistent seasonality.

Error-Trend Seasonal (ETS): Similar to Holt-Winter, relies on exponential smoothing built upon the components of trend, seasonality, and remainder.

Neural Networks (NNETAR): Creates a series of feed-forward version of networks of lagged values. This method is more complex and can result in more accurate results but often requires more data and may have more variability.

Trigonometric, Box-Cox, ARMA, Trend, Seasonal Components (TBATS): We learned about Fourier terms and how it uses trigonometry to forecast. This method relies on long seasonal cycles (which we have) but is also more complex and requires a lot of computational data to use, just like Neural Networks.

fit_arima <- auto.arima(diff_unemployment2)
fit_hw    <- HoltWinters(diff_unemployment2)
fit_ets   <- ets(diff_unemployment2)
fit_nnet  <- nnetar(diff_unemployment2)
fit_tbats <- tbats(diff_unemployment2)


start_year  <- start(diff_unemployment2)[1]
start_month <- start(diff_unemployment2)[2]

dates <- seq.Date(
  from = as.Date(sprintf("%d-%02d-01", start_year, start_month)),
  by = "month",
  length.out = length(diff_unemployment2)
)

df <- data.frame(ds = dates, y = as.numeric(diff_unemployment2))
fit_prophet <- prophet(df)
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
future <- make_future_dataframe(fit_prophet, periods = 12, freq = "month")
fcast_prophet <- predict(fit_prophet, future)


fcast_arima <- forecast(fit_arima, h = 12)
fcast_hw    <- forecast(fit_hw, h = 12)
fcast_ets   <- forecast(fit_ets, h = 12)
fcast_nnet  <- forecast(fit_nnet, h = 12)
fcast_tbats <- forecast(fit_tbats, h = 12)

future <- make_future_dataframe(fit_prophet, periods = 12, freq = "month")
fcast_prophet <- predict(fit_prophet, future)


n <- length(diff_unemployment2)

train <- ts(diff_unemployment2[1:(n - 12)], 
            start = start(diff_unemployment2), 
            frequency = frequency(diff_unemployment2))

test <- ts(diff_unemployment2[(n - 12 + 1):n], 
           start = time(diff_unemployment2)[n - 12 + 1], 
           frequency = frequency(diff_unemployment2))


models <- list(
  ARIMA = auto.arima(train),
  HW    = HoltWinters(train),
  ETS   = ets(train),
  NNET  = nnetar(train),
  TBATS = tbats(train)
)


errors <- lapply(models, function(m) {
  fc <- forecast(m, h = 12)
  c(RMSE = rmse(test, fc$mean),
    MAE  = mae(test, fc$mean))
})

results <- do.call(rbind, errors)
print(results)
##             RMSE        MAE
## ARIMA 0.08060121 0.06637752
## HW    0.10569756 0.09621599
## ETS   0.08165215 0.06688084
## NNET  0.13708017 0.10603820
## TBATS 0.08413264 0.07342961

As the results show above, we have calculated the RMSE and MAE for each of the forecasting methods. Since there are error values, it would make sense that the forecast methods with the lowest values would be more ideal to go with. However, the differences here are very marginal between all five of the methods, perhaps because we have made the data stationary, thus the variance is too stabilized. We will choose the ARIMA and ETS forecasts for now since they have the lowest errors, and these two methods were the ones most emphasized during class.

Visual Overview of the Five Forecasts:

plot(fcast_arima, main="ARIMA Forecast")

plot(fcast_hw, main="Holt-Winters Forecast")

plot(fcast_ets, main="ETS Forecast")

plot(fcast_nnet, main="NNETAR Forecast")

plot(fcast_tbats, main="TBATS Forecast")

plot(fit_prophet, fcast_prophet) + ggplot2::ggtitle("Prophet Forecast")

Part Five) Multivariate Forecasting

Around the second half of the class, we have also learned to consider other variables that may potentially affect the movement of the main variable we are examining. We assume ceteris paribus when we forecast a variable by just examining its own past values and natures, but in the real world, that is not very realistic. As we said in the beginning, inflation was another variable that was part of our initial data. Before we dissect the variable of inflation, it must go through the stationarity test, just like unemployment did.

adf.test(inflation)
## Warning in adf.test(inflation): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  inflation
## Dickey-Fuller = -4.2567, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(inflation)
## Warning in kpss.test(inflation): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  inflation
## KPSS Level = 4.6296, Truncation lag parameter = 6, p-value = 0.01
diff_inflation <- diff(inflation)

Inflation did not pass for stationarity, but we have noticed and fixed that now.

We will now combine the two variables into one data. When doing a VAR analysis, it is necessary for us to know how many lags is ideal before we actually forecast the data. Previously for unemmployment, we did a similar process where conducted an autocorrelation test and found out that any lag below 5 will result in inaccuracy. This is a similar process, except we are now doing it for two variables:

new_data <- cbind(diff_unemployment2, diff_inflation)
colnames(new_data) <- c("Unemployment", "Inflation")
VARselect(new_data, lag.max = 30, type = "const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     25     13      1     25 
## 
## $criteria
##                  1           2           3           4           5           6
## AIC(n) -4.47280182 -4.49867687 -4.49016309 -4.48864809 -4.51249939 -4.51036601
## HQ(n)  -4.45697259 -4.47229482 -4.45322822 -4.44116041 -4.45445888 -4.44177269
## SC(n)  -4.43196327 -4.43061263 -4.39487316 -4.36613246 -4.36275806 -4.33339899
## FPE(n)  0.01141529  0.01112371  0.01121883  0.01123586  0.01097107  0.01099455
##                  7           8           9          10          11           12
## AIC(n) -4.50448477 -4.50908238 -4.53704711 -4.53348453 -4.53569484 -4.714204000
## HQ(n)  -4.42533862 -4.41938341 -4.43679533 -4.42267993 -4.41433741 -4.582293758
## SC(n)  -4.30029205 -4.27766396 -4.27840300 -4.24761472 -4.22259933 -4.373882800
## FPE(n)  0.01105946  0.01100881  0.01070531  0.01074364  0.01072006  0.008967651
##                  13           14           15           16           17
## AIC(n) -4.781707970 -4.771291607 -4.760420293 -4.752315719 -4.745600530
## HQ(n)  -4.639244909 -4.618275727 -4.596851594 -4.578194200 -4.560926193
## SC(n)  -4.414161074 -4.376519015 -4.338422006 -4.303091735 -4.269150851
## FPE(n)  0.008382437  0.008470393  0.008563195  0.008633124  0.008691573
##                 18           19           20           21           22
## AIC(n) -4.75256305 -4.749009006 -4.737599239 -4.734640405 -4.745971073
## HQ(n)  -4.55733590 -4.543229030 -4.521266443 -4.507754790 -4.508532639
## SC(n)  -4.24888768 -4.218107935 -4.179472471 -4.149287941 -4.133392914
## FPE(n)  0.00863158  0.008662661  0.008762458  0.008788856  0.008690305
##                  23          24         25           26           27
## AIC(n) -4.734820656 -4.80122530 -4.8373288 -4.827557366 -4.819266525
## HQ(n)  -4.486829402 -4.54268123 -4.5682319 -4.547907655 -4.529063994
## SC(n)  -4.095016800 -4.13419575 -4.1430736 -4.106076423 -4.070559886
## FPE(n)  0.008788268  0.00822417  0.0079331  0.008011604  0.008078964
##                  28           29           30
## AIC(n) -4.811577346 -4.805316097 -4.800173247
## HQ(n)  -4.510821996 -4.494007928 -4.478312258
## SC(n)  -4.035645011 -4.002158066 -3.969789520
## FPE(n)  0.008142039  0.008193952  0.008237032

Based on the selection options (AIC,HQ,SC,FPE), I will conclude that 13 lags is appropriate. Although the AIC and FPE says 25, the SC is much lower, so HQ’s lag (13) is a good compromise.

From there on we could do some fun tests that relate to multivariable regression analysis, such as Granger-casuality. The Granger-casuality tests which variable is essentally causing the other, hence the name, “casuality”. For example, we would expect the amount of hours someone studies for a test to cause what their grade outcome is. However, the cases are not always so simple and straight-forward, like our case. We should do a formal test to see which variable is the causer.

var_model <- VAR(new_data, p = 13, type = "const")

causality(var_model, cause = "Inflation")
## $Granger
## 
##  Granger causality H0: Inflation do not Granger-cause Unemployment
## 
## data:  VAR object var_model
## F-Test = 0.67578, df1 = 13, df2 = 1300, p-value = 0.7884
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Inflation and Unemployment
## 
## data:  VAR object var_model
## Chi-squared = 9.3495, df = 1, p-value = 0.00223
causality(var_model, cause = "Unemployment")
## $Granger
## 
##  Granger causality H0: Unemployment do not Granger-cause Inflation
## 
## data:  VAR object var_model
## F-Test = 1.8831, df1 = 13, df2 = 1300, p-value = 0.02802
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Unemployment and Inflation
## 
## data:  VAR object var_model
## Chi-squared = 9.3495, df = 1, p-value = 0.00223

It is clear after looking at the results that it is unemployment that causes inflation. Our H0 of “inflation does not cause unemployment” has a p-value of 0.7884, meaning that our H0 is accepted. In contrast, the H0 of “Unemployment does not cause inflation” has a p-value of 0.02802, meaning that the H0 is rejected. Realistically, both variables may affect each other, as our economy is not always so straight-forward. But for the sake of this project, we will just conclude one of the variables as the causer. This result is a bit ironic to what we’re trying to do, since we ideally want inflation to cause unemployment, as that is the variable we are forecasting.

The Impulse Response Function is another way to see how one variable is affected by another. When the x of the variable (the impulse) is derived (increased by a unit of 1), we can see how the other variable responds to that increase over time. It is important to note that we are mostly concenred with the short-term changes here, hence the name “impulse”.

plot(irf(var_model, impulse = "Inflation", response = "Unemployment", n.ahead = 13, boot = TRUE))

plot(irf(var_model, impulse = "Unemployment", response = "Inflation", n.ahead = 13, boot = TRUE))

The first graph shows how unemployment changes after a burst increase in inflation. For four periods (four months), there isn’t much of a response at all, until an increase that persists until the end of the year (12 months). The confidence intervals are labeled as the dotted lines, meaning that there is possibility that there may be a decrease in unemployment when inflation increases as well.

The second graph where unemployment acts as the impulse looks a lot different. Even initially, the value of inflation is under 0 when unemployment increases. In terms of the derivative, this is not that important, as that only represents the intercept. Near the end of the 12-13 month period, inflation seems to be back around 0, which is a relatively long time to recover after an impulse change.

The first graph is more relevant to our study since we are concerned with how unemployment moves. We can now infer that if there is a sudden shock of increase in inflation, the increase in unemployment will also persist for about a year, though the increase won’t start for about four months. This may or may not be helpful to our forecast, since we don’t know if there actually will be a sudden change in inflation when we forecast unemployment.

Another test that is commonly applied with two or more variables is the test for co-integration. We use the Johansan test to see if there is co-integration between unemployment and inflation. If two variables are co-integrated, it means they are moving together, indicating some relationship. Fortunately, we have already made both of the data stationary, so we don’t have to worry about this. I have tested for it below nevertheless:

summary(ca.jo(new_data, type = "trace", K = 13, ecdet = "const", spec = "transitory"))
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 8.590445e-02 6.198138e-02 1.387779e-17
## 
## Values of teststatistic and critical values of test:
## 
##            test 10pct  5pct  1pct
## r <= 1 |  43.32  7.52  9.24 12.97
## r = 0  | 104.13 17.85 19.96 24.60
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##                 Unemployment.l1 Inflation.l1  constant
## Unemployment.l1       1.0000000    1.0000000  1.000000
## Inflation.l1          0.4837856   -1.5641378 -1.621987
## constant              1.0006540    0.9945516 20.030259
## 
## Weights W:
## (This is the loading matrix)
## 
##                Unemployment.l1 Inflation.l1      constant
## Unemployment.d      -0.6688915   -0.4070893 -1.854847e-18
## Inflation.d         -0.3901764    0.1306444 -1.962716e-18

Finally, the VAR forecast that applies inflation onto the unemployment forecast is finally made. Visually, it is hard to tell the differences between this and the five previous forecasts we made, so we go back to the errors:

var_forecast <- predict(var_model, n.ahead = 13)
plot(var_forecast, names = "Unemployment")

n <- nrow(new_data)

train <- window(new_data, end = time(new_data)[n - 13])   
test  <- window(new_data, start = time(new_data)[n - 13 + 1])  


var_model_train <- VAR(train, p = 13, type = "const")


fc <- predict(var_model_train, n.ahead = 13)

preds <- fc$fcst$Unemployment[, "fcst"]

actual <- test[, "Unemployment"]

rmse_val <- rmse(actual, preds)
mae_val  <- mae(actual, preds)

c(RMSE = rmse_val, MAE = mae_val)
##       RMSE        MAE 
## 0.08914236 0.07780328

The errors are also similar. However, the VAR forecast is conceptually different enough where it incorporates a different variable to predict the movement of the main one, so presumably it would cover up a lot of the weaknesses that the single-variable forecasting won’t be able to see. For our combination forecasts, we will go with ARIMA, ETS, and VAR.

Combination forecasts are a good way to incorporate multiple aspects of each of the forecast methods we are using. Each forecast method has its own weakness, so by combining them their weaknesses are more likely covered up and balanced through there being multiple forecasts. Of course, this isn’t always ideal, since a single forecast may explain the model far better than any other method. But in this case it is necessary, since all of our forecasts showed marginal differences in errors.

h <- 12
ts_data <- diff_unemployment   


fit_arima <- auto.arima(ts_data)
fc_arima  <- forecast(fit_arima, h = h)

fit_ets <- ets(ts_data)
fc_ets  <- forecast(fit_ets, h = h)

data_mat <- cbind(
  Unemployment = ts_data,
  Inflation    = diff_inflation
)
lag_sel <- VARselect(data_mat, lag.max = 25, type = "const")
p_chosen <- lag_sel$selection[["AIC(n)"]]   
var_model <- VAR(data_mat, p = p_chosen, type = "const")
fc_var <- predict(var_model, n.ahead = h)

var_mean  <- fc_var$fcst$Unemployment[, "fcst"]
var_lower <- fc_var$fcst$Unemployment[, "lower"]
var_upper <- fc_var$fcst$Unemployment[, "upper"]


point_mat <- cbind(
  ARIMA = as.numeric(fc_arima$mean),
  ETS   = as.numeric(fc_ets$mean),
  VAR   = as.numeric(var_mean)
)
avg_forecast <- rowMeans(point_mat)

lower_mat <- cbind(fc_arima$lower[,2], fc_ets$lower[,2], var_lower)
upper_mat <- cbind(fc_arima$upper[,2], fc_ets$upper[,2], var_upper)
avg_lower <- rowMeans(lower_mat)
avg_upper <- rowMeans(upper_mat)

fc_mean_ts  <- ts(avg_forecast, start = time(ts_data)[length(ts_data)] + 1/12, frequency = frequency(ts_data))
fc_lower_ts <- ts(avg_lower,  start = time(ts_data)[length(ts_data)] + 1/12, frequency = frequency(ts_data))
fc_upper_ts <- ts(avg_upper,  start = time(ts_data)[length(ts_data)] + 1/12, frequency = frequency(ts_data))


ts.plot(ts_data,
        col = "black", lwd = 1.5,
        main = "Combined Forecast (ARIMA + ETS + VAR)",
        ylab = "Differenced Unemployment (Δ %)")
lines(fc_mean_ts,  col = "blue", lwd = 2)
lines(fc_lower_ts, col = "red",  lty = 2)
lines(fc_upper_ts, col = "red",  lty = 2)
legend("topleft",
       legend = c("History (Diff)", "Combined Forecast", "95% Interval"),
       col = c("black", "blue", "red"),
       lty = c(1,1,2), bty = "n")

fc_mean_ts
##                Jan           Feb           Mar           Apr           May
## 2025                                                                      
## 2026  9.270338e-04  6.387096e-03 -2.818944e-05  9.137832e-03  1.899044e-04
##                Jun           Jul           Aug           Sep           Oct
## 2025                             -1.647829e-02 -3.181440e-02 -3.030097e-03
## 2026  2.343254e-03  8.759024e-04                                          
##                Nov           Dec
## 2025 -7.597461e-03 -1.781612e-03
## 2026

Visually, there is almost no difference between the combination forecasts and previous three. This is because there was almost no difference between the three already before the forecasts were combined. In real life scenarios, this is necessary, since forecasting trends and changes is difficult, and we want to ensure that our forecasts are within a band width that we are confident with before we report the data. The raw numbers of the forecasts do not mean much either since we adjusted the data for normality, so they don’t help much either. Instead, what we can do is try forecasting before the raw data was adjusted:

Here is the combination forecast of the original data of employment before it was differentiated and adjusted for normality:

h <- 12   
ts_data <- unemployment  



fit_arima <- auto.arima(ts_data)
fc_arima  <- forecast(fit_arima, h = h)

fit_ets <- ets(ts_data)
fc_ets  <- forecast(fit_ets, h = h)


data_mat <- cbind(Unemployment = ts_data, Inflation = inflation)
var_model <- VAR(data_mat, p = 2, type = "const")   
fc_var <- predict(var_model, n.ahead = h)


fc_var_mean  <- fc_var$fcst$Unemployment[, "fcst"]
fc_var_lower <- fc_var$fcst$Unemployment[, "lower"]
fc_var_upper <- fc_var$fcst$Unemployment[, "upper"]


point_mat <- cbind(
  ARIMA = fc_arima$mean,
  ETS   = fc_ets$mean,
  VAR   = fc_var_mean
)
avg_forecast <- rowMeans(point_mat)


lower_mat <- cbind(fc_arima$lower[,2], fc_ets$lower[,2], fc_var_lower)
upper_mat <- cbind(fc_arima$upper[,2], fc_ets$upper[,2], fc_var_upper)
avg_lower <- rowMeans(lower_mat)
avg_upper <- rowMeans(upper_mat)

combined_fc <- list(
  mean  = ts(avg_forecast, start = time(ts_data)[length(ts_data)] + 1/12, frequency = 12),
  lower = ts(avg_lower, start = time(ts_data)[length(ts_data)] + 1/12, frequency = 12),
  upper = ts(avg_upper, start = time(ts_data)[length(ts_data)] + 1/12, frequency = 12)
)


ts.plot(ts_data, col = "black", lwd = 1.5,
        main = "Combined Forecast (ARIMA + ETS + VAR)",
        ylab = "Unemployment Rate (%)")
lines(combined_fc$mean, col = "blue", lwd = 2)
lines(combined_fc$lower, col = "red", lty = 2)
lines(combined_fc$upper, col = "red", lty = 2)

legend("topleft", legend = c("History", "Combined Forecast", "95% Interval"),
       col = c("black", "blue", "red"), lty = c(1,1,2), bty = "n")

combined_fc$mean
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2025                                                                4.218808
## 2026 4.308168 4.325063 4.341589 4.357744 4.373530 4.388951 4.404007         
##           Sep      Oct      Nov      Dec
## 2025 4.237148 4.255349 4.273291 4.290906
## 2026

Here, both the visuals and numbers of the forecast are a lot more noticeable. It seems like the combination forecast predicts unemployment to slightly increase, but we can never be sure since the 95% interval is too large. If we were to look at the raw numbers, the forecast predicts the unemployment rate during September to be 4.2% (Currently in September it is at 4.3%, which turned out to be fairly accurate). However, we have to remember that this is a forecast of the original data before all the modifications, and in real life scenarios we would report the former data that hovers around the mean, since the original data did not pass the test for stationarity.

Part Six) Conclusion and References:

Our main objective in this project was to utilize a relevant time series data and analyze the necessary components of it before we ultimately performed a forecast. Taking the Unemployment data, we tested for stationarity, which did not pass. This adjustment of differentiation made the data a lot more mean recurring, thus the forecast also showed a lot more general results. When we forecasted the original data, the results were a lot more noticeable, but unreliable since the original data did not pass for stationarity. We also did a other computations that gave information to the qualities and characteristics of unemployment, such as its additive decomposition, its linear fit and residuals, and finding the variable’s relationship to another crucial time series data (inflation) that most likely had an effect on unemployment.

Sources: Bureau Labor of Statistics, ChatGPT, and previous Econ 144 projects