## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  cache = TRUE,
  fig.width = 7,
  fig.height = 5,
  comment = "#>"
)

opt <- options(max.print = 70)

## ----setup, message=FALSE-----------------------------------------------------
library(dfms)
library(xts)

## -----------------------------------------------------------------------------
# Using the monthly series from BM14
dim(BM14_M)
range(index(BM14_M))
head(colnames(BM14_M))
plot(scale(BM14_M), lwd = 1)

## -----------------------------------------------------------------------------
head(BM14_Models, 3)

# Using only monthly data
BM14_Models_M <- subset(BM14_Models, freq == "M")

## -----------------------------------------------------------------------------
library(magrittr)
# log-transforming and first-differencing the data
BM14_M[, BM14_Models_M$log_trans] %<>% log()
BM14_M_diff <- diff(BM14_M)
plot(scale(BM14_M_diff), lwd = 1)

## -----------------------------------------------------------------------------
ic <- ICr(BM14_M_diff)
print(ic)
plot(ic)

## -----------------------------------------------------------------------------
screeplot(ic)

## -----------------------------------------------------------------------------
# Using vars::VARselect() with 4 principal components to estimate the VAR lag order
vars::VARselect(ic$F_pca[, 1:4])

## -----------------------------------------------------------------------------
# Estimating the model with 4 factors and 3 lags using BM14's EM algorithm
model_m <- DFM(BM14_M_diff, r = 4, p = 3)
print(model_m)
plot(model_m)

## -----------------------------------------------------------------------------
dfm_summary <- summary(model_m)
print(dfm_summary) # Large model with > 40 series: defaults to compact = 2

# Can request more detailed printouts
# print(dfm_summary, compact = 1)
# print(dfm_summary, compact = 0) 

## -----------------------------------------------------------------------------
plot(resid(model_m, orig.format = TRUE), lwd = 1)
plot(fitted(model_m, orig.format = TRUE), lwd = 1)

## -----------------------------------------------------------------------------
plot(model_m, method = "all", type = "individual")

## -----------------------------------------------------------------------------
# Default: all estimates in long format
head(as.data.frame(model_m, time = index(BM14_M_diff)))

## -----------------------------------------------------------------------------
# 12-period ahead DFM forecast
fc <- predict(model_m, h = 12)
print(fc)

## -----------------------------------------------------------------------------
# Setting an appropriate plot range to see the forecast
plot(fc, xlim = c(320, 370))

## -----------------------------------------------------------------------------
# Factor forecasts in wide format
head(as.data.frame(fc, pivot = "wide"))

## ----MQ-----------------------------------------------------------------------
# Quarterly series from BM14
head(BM14_Q, 3)
# Pre-processing the data
BM14_Q[, BM14_Models$log_trans[BM14_Models$freq == "Q"]] %<>% log()
BM14_Q_diff <- diff(BM14_Q)
# Merging to monthly data
BM14_diff <- merge(BM14_M_diff, BM14_Q_diff)

# Estimating the model with 5 factors and 3 lags using BM14's EM algorithm
model_mq <- DFM(BM14_diff, r = 5, p = 3, quarterly.vars = colnames(BM14_Q))
print(model_mq)
plot(model_mq)

## ----MQ_AR1-------------------------------------------------------------------
BM14_med_diff <- BM14_diff[, BM14_Models$medium]
Q_vars_med <- intersect(colnames(BM14_med_diff), colnames(BM14_Q))

# Mixed-frequency model with AR(1) idiosyncratic errors
model_mq_ar1 <- DFM(BM14_med_diff, r = 5, p = 3, idio.ar1 = TRUE,
                    quarterly.vars = Q_vars_med)
print(model_mq_ar1)

## -----------------------------------------------------------------------------
# AR(1) coefficients for the first 10 series
head(model_mq_ar1$rho, 10)

# Summary of AR(1) coefficients
summary(model_mq_ar1$rho)

## -----------------------------------------------------------------------------
# Dimension of estimated errors
dim(model_mq_ar1$e)

# Plot estimated errors for first 5 monthly series
matplot(model_mq_ar1$e[-(1:120), 1:5], type = "l", lty = 1,
        main = "Estimated Idiosyncratic Errors (First 5 Series)",
        xlab = "Time", ylab = "Error")

## ----news---------------------------------------------------------------------
# Create old vintage with ragged edge
model_mq_ar1$rm.rows # Previous model removed first row
X_old <- BM14_med_diff[-model_mq_ar1$rm.rows, ]
tail(X_old, 3) # See current ragged edge
(obscurr <- colnames(X_old)[is.finite(tail(X_old, 1))]) # Observe current
# Set some observed values to NA
set.seed(1)
X_old[nrow(X_old), sample(obscurr, 10)] <- NA
# Also check series missing the previous observation and set some to NA
(obsprev <- colnames(X_old)[is.finite(tail(X_old, 2)[1, ])] |> setdiff(obscurr)) 
set.seed(1)
X_old[nrow(X_old), sample(obsprev, 4)] <- NA

# Estimating same DFM on old vintage
model_mq_ar1_old <- model_mq_ar1$call
model_mq_ar1_old$X <- quote(X_old)
model_mq_ar1_old$max.missing = 1 # No further row removals
print(model_mq_ar1_old)
model_mq_ar1_old <- eval(model_mq_ar1_old)

# Compute news for "gdp" variable at the last period (t.fcst = nrow(X_old))
news_gdp <- news(model_mq_ar1_old, model_mq_ar1, target.vars = "gdp")
print(news_gdp)
# Verifying the news decomposition
na.omit(news_gdp$news_df) |> 
  transform(test1 = news - (actual - forecast), test2 = impact - news * gain)

# Now computing news for all quarterly variables
news_qvars <- news(model_mq_ar1_old, model_mq_ar1, target.vars = Q_vars_med)
print(news_qvars)
# Data frame with details
as.data.frame(news_qvars) |> na.omit() |> head()

## ----include=FALSE------------------------------------------------------------
options(opt)