inst/replications/ipca_grunfeld.R

In sdim: An R Package for Supervised Dimension Reduction

## ============================================================================
## Validation of IPCA implementation against Python ipca package
## Data: Grunfeld (1958) investment dataset (statsmodels.datasets.grunfeld)
## Reference: Kelly, Pruitt, Su (2019), JFE
##
## Python reference (panel-based ALS, n_factors=1, no intercept):
##   Gamma = [0.99166014, 0.12888046]'
##   Converged in 15 iterations
## ============================================================================

devtools::load_all(".", quiet = TRUE)

# ---------- 1. Grunfeld data -------------------------------------------------
# 11 firms, 20 years (1935-1954), 3 variables: invest, value, capital

ret_data <- matrix(c(
  2.938,39.68,40.29,2.54,33.1,317.6,26.63,20.36,209.9,24.43,12.93,
  5.643,50.73,72.76,2,45,391.8,23.39,25.98,355.3,23.21,25.9,
  10.233,74.24,66.26,2.19,77.2,410.6,30.65,25.94,469.9,32.78,35.05,
  4.046,53.51,51.6,1.99,44.6,257.7,20.89,27.53,262.3,32.54,22.89,
  3.326,42.65,52.41,2.03,48.1,330.8,28.78,24.6,230.4,26.65,18.84,
  4.68,46.48,69.41,1.81,74.4,461.2,26.93,28.54,361.6,33.71,28.57,
  5.732,61.4,68.35,2.14,113,512,32.08,43.41,472.8,43.5,48.51,
  12.117,39.67,46.8,1.86,91.9,448,32.21,42.81,445.6,34.46,43.34,
  15.276,62.24,47.4,0.93,61.3,499.6,35.69,27.84,361.6,44.28,37.02,
  9.275,52.32,59.57,1.18,56.8,547.5,62.47,32.6,288.2,70.8,37.81,
  9.577,63.21,88.78,1.36,93.6,561.2,52.32,39.03,258.7,44.12,39.27,
  3.956,59.37,74.12,2.24,159.9,688.1,56.95,50.17,420.3,48.98,53.46,
  3.834,58.02,62.68,3.81,147.2,568.9,54.32,51.85,420.5,48.51,55.56,
  5.97,70.34,89.36,5.66,146.3,529.2,40.53,64.03,494.5,50,49.56,
  6.433,67.42,78.98,4.21,98.3,555.1,32.54,68.16,405.1,50.59,32.04,
  4.77,55.74,100.66,3.42,93.5,642.9,43.48,77.34,418.8,42.53,32.24,
  6.532,80.3,160.62,4.67,135.2,755.9,56.49,95.3,588.2,64.77,54.38,
  7.329,85.4,145,6,157.3,891.2,65.98,99.49,645.5,72.68,71.78,
  9.02,91.9,174.93,6.53,179.5,1304.4,66.11,127.52,641,73.86,90.08,
  6.281,81.43,172.49,5.12,189.6,1486.7,49.34,135.72,459.3,89.51,68.6
), nrow = 20, ncol = 11, byrow = TRUE)

value <- matrix(c(
  30.284,157.7,417.5,70.91,1170.6,3078.5,290.6,197,1362.4,138,191.5,
  43.909,167.9,837.8,87.94,2015.8,4661.7,291.1,210.3,1807.1,200.1,516,
  107.02,192.9,883.9,82.2,2803.3,5387.1,335,223.1,2676.3,210.1,729,
  68.306,156.7,437.9,58.72,2039.7,2792.2,246,216.7,1801.9,161.2,560.4,
  84.164,191.4,679.7,80.54,2256.2,4313.2,356.2,286.4,1957.3,161.7,519.9,
  69.157,185.5,727.8,86.47,2132.2,4643.9,289.8,298,2202.9,145.1,628.5,
  60.148,199.6,643.6,77.68,1834.1,4551.2,268.2,276.9,2380.5,110.6,537.1,
  49.332,189.5,410.9,62.16,1588,3244.1,213.3,272.6,2168.6,98.1,561.2,
  75.18,151.2,588.4,62.24,1749.4,4053.7,348.2,287.4,1985.1,108.8,617.2,
  62.05,187.7,698.4,61.82,1687.2,4379.3,374.2,330.3,1813.9,118.2,626.7,
  59.152,214.7,846.4,65.85,2007.7,4840.9,387.2,324.4,1850.2,126.5,737.2,
  68.424,232.9,893.8,69.54,2208.3,4900.9,347.4,401.9,2067.7,156.7,760.5,
  48.505,249,579,64.97,1656.7,3526.5,291.9,407.4,1796.7,119.4,581.4,
  40.507,224.5,694.6,68,1604.4,3254.7,297.2,409.2,1625.8,129.1,662.3,
  39.961,237.3,590.3,71.24,1431.8,3700.2,276.9,482.2,1667,134.8,583.8,
  36.494,240.1,693.5,69.05,1610.5,3755.6,274.6,673.8,1677.4,140.8,635.2,
  46.082,327.3,809,83.04,1819.4,4833,339.9,676.9,2289.5,179,723.8,
  57.616,359.4,727,74.42,2079.7,4924.9,474.8,702,2159.4,178.1,864.1,
  57.441,398.4,1001.5,63.51,2371.6,6241.7,496,793.5,2031.3,186.8,1193.5,
  47.165,365.7,703.2,58.12,2759.9,5593.6,474.5,927.3,2115.5,192.7,1188.9
), nrow = 20, ncol = 11, byrow = TRUE)

capital <- matrix(c(
  52.011,183.2,10.5,4.5,97.8,2.8,162,6.5,53.8,100.2,1.8,
  52.903,204,10.2,4.71,104.4,52.6,174,15.8,50.5,125,0.8,
  54.499,236,34.7,4.57,118,156.9,183,27.7,118.1,142.4,7.4,
  59.722,291.7,51.8,4.56,156.2,209.2,198,39.2,260.2,165.1,18.1,
  61.659,323.1,64.3,4.38,172.6,203.4,208,48.6,312.7,194.8,23.5,
  62.243,344,67.1,4.21,186.6,207.2,223,52.5,254.2,222.9,26.5,
  63.361,367.7,75.2,4.12,220.9,255.2,234,61.5,261.4,252.1,36.2,
  64.861,407.2,71.4,3.83,287.8,303.7,248,80.5,298.7,276.3,60.8,
  67.953,426.6,67.1,3.58,319.9,264.1,274,94.4,301.8,300.3,84.4,
  69.59,470,60.5,3.41,321.3,201.6,282,92.6,279.1,318.2,91.2,
  69.144,499.2,54.6,3.31,319.6,265,316,92.3,213.8,336.2,92.4,
  70.269,534.6,84.8,3.23,346,402.2,302,94.2,132.6,351.2,86,
  71.051,566.6,96.8,3.9,456.4,761.5,333,111.4,264.8,373.6,111.1,
  71.508,595.3,110.2,5.38,543.4,922.4,359,127.4,306.9,389.4,130.6,
  73.827,631.4,147.4,7.39,618.3,1020.1,370,149.3,351.1,406.7,141.8,
  75.847,662.3,163.2,8.74,647.4,1099,376,164.4,357.8,429.5,136.7,
  77.367,683.9,203.5,9.07,671.3,1207.7,391,177.2,342.1,450.6,129.7,
  78.631,729.3,290.6,9.93,726.1,1430.5,414,200,444.2,466.9,145.5,
  80.215,774.3,346.1,11.68,800.3,1777.3,443,211.5,623.6,486.2,174.8,
  83.788,804.9,414.9,14.33,888.9,2226.3,468,238.7,669.7,511.3,213.5
), nrow = 20, ncol = 11, byrow = TRUE)

# Build Z array (T x N x L)
Z <- array(NA_real_, dim = c(20, 11, 2))
Z[, , 1] <- value
Z[, , 2] <- capital

cat("ret:", nrow(ret_data), "x", ncol(ret_data), "\n")
cat("Z:  ", paste(dim(Z), collapse = " x "), "\n\n")

# ---------- 2. Fit IPCA with sdim --------------------------------------------

fit <- ipca_est(ret_data, Z, nfac = 1)

cat("=== sdim Gamma (lambda) ===\n")
print(fit$lambda)

cat("\n=== sdim Factors ===\n")
print(fit$factors)

cat("\n=== sdim Eigvals ===\n")
print(fit$eigvals)

# ---------- 3. Python reference values ----------------------------------------

py_gamma <- c(0.99166014, 0.12888046)

py_factors <- c(
  0.1031968381, 0.0884489515, 0.0838496628, 0.0845069923, 0.0722523449,
  0.0995068155, 0.1228840058, 0.1422623752, 0.1197532025, 0.1179724004,
  0.1087561863, 0.1357521189, 0.1579348267, 0.1660545375, 0.1484923276,
  0.1586634303, 0.1596007400, 0.1759379247, 0.1921695585, 0.2111065868
)

# ---------- 4. Compare (rotation-invariant) -----------------------------------
# Gamma and F are identified only up to sign flip; compare absolute values
# and check that fitted values match.

cat("\n=== Comparison ===\n")

# Sign-align: if correlation is negative, flip R result
r_gamma <- as.numeric(fit$lambda)
if (cor(r_gamma, py_gamma) < 0) {
  r_gamma  <- -r_gamma
  r_factors <- -as.numeric(fit$factors)
} else {
  r_factors <- as.numeric(fit$factors)
}

cat("Gamma (Python):  ", sprintf("%.8f", py_gamma), "\n")
cat("Gamma (R sdim):  ", sprintf("%.8f", r_gamma), "\n")
cat("Gamma max |diff|:", sprintf("%.2e", max(abs(r_gamma - py_gamma))), "\n\n")

cat("Factors max |diff|:", sprintf("%.2e", max(abs(r_factors - py_factors))), "\n")
cat("Factors correlation:", sprintf("%.10f", cor(r_factors, py_factors)), "\n\n")

# Fitted values comparison: y_hat = Z_{t} %*% Gamma %*% f_t
py_fitted <- numeric(20 * 11)
r_fitted  <- numeric(20 * 11)
idx <- 1
for (t in seq_len(20)) {
  for (i in seq_len(11)) {
    z_ti <- Z[t, i, ]
    py_fitted[idx] <- sum(z_ti * py_gamma) * py_factors[t]
    r_fitted[idx]  <- sum(z_ti * r_gamma) * r_factors[t]
    idx <- idx + 1
  }
}

cat("Fitted values max |diff|:", sprintf("%.2e", max(abs(r_fitted - py_fitted))), "\n")
cat("Fitted values R-squared: ", sprintf("%.10f", cor(r_fitted, py_fitted)^2), "\n")

# ---------- 5. Pass/fail summary ---------------------------------------------

tol_gamma   <- 1e-4
tol_factors <- 1e-4
tol_fitted  <- 1e-3

pass_gamma   <- max(abs(r_gamma - py_gamma)) < tol_gamma
pass_factors <- max(abs(r_factors - py_factors)) < tol_factors
pass_fitted  <- max(abs(r_fitted - py_fitted)) < tol_fitted

cat("\n=== VALIDATION RESULTS ===\n")
cat("Gamma match:   ", ifelse(pass_gamma, "PASS", "FAIL"),
    sprintf("(tol = %.0e)\n", tol_gamma))
cat("Factors match: ", ifelse(pass_factors, "PASS", "FAIL"),
    sprintf("(tol = %.0e)\n", tol_factors))
cat("Fitted match:  ", ifelse(pass_fitted, "PASS", "FAIL"),
    sprintf("(tol = %.0e)\n", tol_fitted))