---
title: "Replicating Huang et al. (2022)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Replicating Huang et al. (2022)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

This vignette reproduces the industrial-production (IP) growth column
of Table 4 in Huang, Jiang, Li, Tong, and Zhou (2022). The exercise
forecasts one-month-ahead IP growth using factors extracted from 123
transformed FRED-MD macro variables, and compares two
dimension-reduction routes:

- **PCA**, the standard principal-component factors of the predictor
  matrix, which ignore the forecast target.
- **sPCA**, the *scaled* variant introduced by Huang et al. (2022).
  Before extracting components, each predictor is scaled by the OLS
  slope from regressing the target on that predictor. Predictors with
  stronger predictive content for the target receive larger weight; the
  irrelevant ones get shrunk toward zero. PCA on the scaled predictors
  therefore concentrates forecasting-relevant information into the
  first few factors.

The reported metric is the out-of-sample $R^2$, $R^2_{OS}$ (%),
relative to an AR benchmark with SIC-selected lag order.

## Data

The package ships the authors' original data:

- `huang2022_macro`: a $720 \times 123$ matrix of transformed FRED-MD
  predictors (January 1960 to December 2019).
- `huang2022_ip`: a 720-vector of monthly IP growth (log-differences of
  the IP index).

```{r}
library(sdim)

data(huang2022_macro)
data(huang2022_ip)

dim(huang2022_macro)
length(huang2022_ip)
```

## Methodology

The forecasting exercise is an expanding-window pseudo out-of-sample
experiment that mirrors the paper:

- **Initial estimation window**: January 1960 to December 1984 (300
  observations).
- **Out-of-sample period**: January 1985 to December 2019 (420
  forecasts).
- **Benchmark**: AR with lag order chosen by SIC (maximum lag 1).
- **Factor models**: ARDL with the AR lags plus one lag of the PCA or
  sPCA factors.

For sPCA, the scaling regression uses the *predictive* alignment
$y_{t+1} \sim X_{i,t}$ — that is, today's predictor is regressed on
next month's target. To dampen extreme slopes the absolute scaling
coefficients are winsorised at the 90th percentile, matching the
authors' MATLAB code. `spca_est()` supports this directly: when
`length(target) < nrow(X)`, the first `length(target)` rows are used
for the scaling regression while *all* rows of `X` are standardised and
used for factor extraction. This is what makes the predictive
alignment possible without losing observations from the factor panel.

## Out-of-sample loop

The function below runs the recursive forecast. At each step it:

1. selects an AR lag order via `select_ar_lag_sic()` and produces the
   AR-only forecast;
2. extracts up to `nfac_max` PCA factors from the standardised
   predictors and up to `nfac_max` sPCA factors with the predictive
   alignment and winsorisation just described;
3. builds an ARDL forecast for each number of factors using both factor
   sets, with the AR coefficients estimated jointly via
   `estimate_ardl_multi()`.

The loop runs ~420 iterations and takes a few minutes, so it is set to
`eval = FALSE` here:

```{r, eval = FALSE}
run_oos <- function(y, Z, h = 1, p_max = 1, nfac_max = 5) {

  TT <- length(y)
  M  <- (1984 - 1959) * 12
  NN <- TT - M

  FC_AR    <- rep(NA, NN - (h - 1))
  FC_PCA   <- matrix(NA, NN - (h - 1), nfac_max)
  FC_sPCA  <- matrix(NA, NN - (h - 1), nfac_max)
  actual_y <- rep(NA, NN - (h - 1))

  for (n in seq_len(NN - (h - 1))) {

    actual_y[n] <- mean(y[(M + n):(M + n + h - 1)])

    y_n  <- y[1:(M + n - 1)]
    Z_n  <- Z[1:(M + n - 1), ]
    Zs_n <- oos_standardize(Z_n)
    T_n  <- length(y_n)

    y_n_h <- vapply(
      seq_len(T_n - (h - 1)),
      function(t) mean(y_n[t:(t + h - 1)]),
      numeric(1)
    )

    # --- AR benchmark with SIC lag selection ---
    p_ar <- select_ar_lag_sic(y_n, h, p_max)

    if (p_ar > 0L) {

      ar_out   <- estimate_ar_res(y_n, h, p_ar)
      y_n_last <- rev(y_n[(T_n - p_ar + 1):T_n])
      FC_AR[n] <- sum(c(1, y_n_last) * ar_out$a_hat)

    } else {

      FC_AR[n] <- mean(y_n)

    }

    # --- PCA factors ---
    pca_fit <- pca_est(X = Zs_n, nfac = nfac_max)
    z_pc_n  <- predict(pca_fit, Zs_n)

    # --- sPCA factors (predictive alignment + winsorization) ---
    spca_fit <- spca_est(
      target       = y_n_h[2:length(y_n_h)],
      X            = Z_n,
      nfac         = nfac_max,
      winsorize    = TRUE,
      winsor_probs = c(0, 90)
    )

    z_trans_n <- predict(spca_fit, Z_n)

    # --- ARDL forecast for each number of factors ---
    for (cc in seq_len(nfac_max)) {

      for (jj in 1:2) {

        z_f <- if (jj == 1) {

          z_pc_n[, 1:cc, drop = FALSE]

        } else {

          z_trans_n[, 1:cc, drop = FALSE]

        }

        p_ardl <- c(p_ar, 1)

        if (p_ar > 0L) {

          c_hat    <- estimate_ardl_multi(y_n, z_f, h, p_ardl)
          y_n_last <- rev(y_n[(T_n - p_ar + 1):T_n])
          fc       <- sum(c(1, y_n_last, z_f[T_n, ]) * c_hat)

        } else {

          dep   <- y_n_h[2:length(y_n_h)]
          reg   <- cbind(1, z_f[1:(length(y_n_h) - 1 - (h - 1)), 1:cc])
          c_hat <- lm.fit(x = reg, y = dep)$coefficients
          fc    <- sum(c(1, z_f[T_n, 1:cc]) * c_hat)

        }

        if (jj == 1) FC_PCA[n, cc]  <- fc
        if (jj == 2) FC_sPCA[n, cc] <- fc

      }

    }

  }

  # R²_OS for each number of factors
  r2_pca <- r2_spca <- numeric(nfac_max)
  sse_ar <- sum((actual_y - FC_AR)^2)

  for (cc in seq_len(nfac_max)) {

    r2_pca[cc]  <- 100 * (1 - sum((actual_y - FC_PCA[, cc])^2)  / sse_ar)
    r2_spca[cc] <- 100 * (1 - sum((actual_y - FC_sPCA[, cc])^2) / sse_ar)

  }

  data.frame(K = seq_len(nfac_max), PCA = round(r2_pca, 2), sPCA = round(r2_spca, 2))

}

# Run
res <- run_oos(huang2022_ip, huang2022_macro, h = 1, p_max = 1, nfac_max = 5)
print(res)
```

## Results

Running the code above produces:

```
  K   PCA  sPCA
  1  8.97  9.65
  2  8.06 10.68
  3  8.22 11.09
  4  7.99 11.97
  5  7.88 13.17
```

With five factors, PCA reaches $R^2_{OS}$ = **7.88%** and sPCA reaches
$R^2_{OS}$ = **13.17%** — both matching the paper to two decimals. sPCA
dominates PCA at every factor count, which is precisely the result the
paper highlights: when many candidate predictors are weak or
irrelevant, weighting each one by its target-predictive slope lets PCA
focus on the components that actually carry forecasting power.

## Key `spca_est()` features used

1. **Predictive alignment**. Passing a `target` that is one observation
   shorter than `X` (i.e. $T-1$ versus $T$) makes the scaling regression
   pair $X_{i,t}$ with $y_{t+1}$, while factors are still extracted
   from the full $T$-row predictor matrix.

2. **Winsorisation**. `winsorize = TRUE` with `winsor_probs = c(0, 90)`
   caps the absolute scaling slopes at their 90th percentile,
   reproducing the trimming used in the authors' MATLAB code.

3. **`predict()`**. Projects the training `X` onto the estimated sPCA
   loadings using the training-window standardisation and scaling, so
   in-sample and out-of-sample factor draws are constructed
   consistently.

## References

Huang, D., Jiang, F., Li, K., Tong, G., and Zhou, G. (2022). Scaled
PCA: A New Approach to Dimension Reduction. *Management Science*,
68(3), 1678--1695.
