# \donttest{
#########################################################
######### Forecasting quarterly U.S. inflation ##########
#### Please see Koop & Korobilis (2023) for further ####
#### details regarding the data & external forecasts ####
#########################################################
# Load Package
library("hdflex")
library("ggplot2")
library("cowplot")
########## Get Data ##########
# Load Package Data
inflation_data <- inflation_data
# Set Target Variable
y <- inflation_data[, 1]
# Set 'P-Signals'
X <- inflation_data[, 2:442]
# Set 'F-Signals'
Ext_F <- inflation_data[, 443:462]
# Get Dates and Number of Observations
tdates <- rownames(inflation_data)
tlength <- length(tdates)
# First complete observation (no missing values)
first_complete <- which(complete.cases(inflation_data))[1]
########## Rolling AR2-Benchmark ##########
# Set up matrix for predictions
benchmark <- matrix(NA, nrow = tlength,
ncol = 1, dimnames = list(tdates, "AR2"))
# Set Window-Size (15 years of quarterly data)
window_size <- 15 * 4
# Time Sequence
t_seq <- seq(window_size, tlength - 1)
# Loop with rolling window
for (t in t_seq) {
# Split Data for Training Train Data
x_train <- cbind(int = 1, X[(t - window_size + 1):t, 1:2])
y_train <- y[(t - window_size + 1):t]
# Split Data for Prediction
x_pred <- cbind(int = 1, X[t + 1, 1:2, drop = FALSE])
# Fit AR-Model
model_ar <- .lm.fit(x_train, y_train)
# Predict and store in benchmark matrix
benchmark[t + 1, ] <- x_pred %*% model_ar$coefficients
}
########## STSC ##########
# Set TV-C-Parameter
init <- 5 * 4
lambda_grid <- c(0.90, 0.95, 1.00)
kappa_grid <- c(0.94, 0.96, 0.98)
bias <- TRUE
# Set DSC-Parameter
gamma_grid <- c(0.40, 0.50, 0.60, 0.70, 0.80, 0.90,
0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00)
n_tvc <- (ncol(X) + ncol(Ext_F)) * length(lambda_grid) * length(kappa_grid)
psi_grid <- c(1:100, sapply(1:4, function(i) floor(i * n_tvc / 4)))
delta <- 0.95
burn_in <- first_complete + init / 2
burn_in_dsc <- 1
metric <- 5
equal_weight <- TRUE
incl <- NULL
parallel <- FALSE
n_threads <- NULL
# Apply STSC-Function
results <- hdflex::stsc(y,
X,
Ext_F,
init,
lambda_grid,
kappa_grid,
bias,
gamma_grid,
psi_grid,
delta,
burn_in,
burn_in_dsc,
metric,
equal_weight,
incl,
parallel,
n_threads,
NULL)
########## Evaluation ##########
# Define Evaluation Period (OOS-Period)
eval_period <- which(tdates >= "1991-04-01" & tdates <= "2021-12-01")
# Get Evaluation Summary for STSC
eval_results <- summary(obj = results, eval_period = eval_period)
# Calculate (Mean-)Squared-Errors for AR2-Benchmark
se_ar2 <- (y[eval_period] - benchmark[eval_period, 1])^2
mse_ar2 <- mean(se_ar2)
# Create CSSED-Plot
cssed <- cumsum(se_ar2 - eval_results$MSE[[2]])
plot_cssed <- ggplot(
data = data.frame(eval_period, cssed),
aes(x = eval_period, y = cssed)
) +
geom_line() +
ylim(-0.0008, 0.0008) +
ggtitle("Cumulative Squared Error Differences") +
xlab("Time Index") +
ylab("CSSED") +
geom_hline(yintercept = 0, linetype = "dashed", color = "darkgray") +
theme_minimal(base_size = 15) +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
axis.ticks = element_line(colour = "black"),
plot.title = element_text(hjust = 0.5)
)
# Show Plots
options(repr.plot.width = 15, repr.plot.height = 15)
plots_list <- eval_results$Plots
plots_list <- c(list(plot_cssed), plots_list)
cowplot::plot_grid(plotlist = plots_list,
ncol = 2,
nrow = 3,
align = "hv")
# Relative MSE
print(paste("Relative MSE:", round(eval_results$MSE[[1]] / mse_ar2, 4)))
# }
Run the code above in your browser using DataLab