## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(factorselect)

## ----simulate-----------------------------------------------------------------
set.seed(42)
X <- simulate_factor_model(N = 100, TT = 200, k = 3, sd = 0.5)
dim(X)

## ----ahn_horenstein-----------------------------------------------------------
result <- select_factors(X, method = "ahn_horenstein", kmax = 8)
print(result)

## ----all_methods--------------------------------------------------------------
result_all <- select_factors(
  X,
  method = c("ahn_horenstein", "bai_ng", "abc",
             "lam_yao", "onatski_2009", "onatski_2010"),
  kmax   = 8
)
print(result_all)

## ----scree, fig.width = 6, fig.height = 4-------------------------------------
result_ah <- select_factors(X, method = "ahn_horenstein", kmax = 8)
plot(result_ah, main = "Scree Plot — Ahn & Horenstein (2013)")

## ----simulation, cache = TRUE-------------------------------------------------
set.seed(123)
n_reps  <- 100
k_true  <- 3
configs <- list(
  large_both  = list(N = 100, TT = 200),
  small_N     = list(N = 25,  TT = 200),
  small_T     = list(N = 200, TT = 25)
)

results <- lapply(configs, function(cfg) {
  estimates <- replicate(n_reps, {
    X <- simulate_factor_model(N = cfg$N, TT = cfg$TT,
                               k = k_true, sd = 0.5)
    res <- select_factors(X,
                          method = c("ahn_horenstein", "bai_ng",
                                     "onatski_2010"),
                          kmax   = 8)
    res$k
  })
  rowMeans(estimates == k_true)
})

# Percentage correct for each configuration
do.call(rbind, lapply(names(results), function(nm) {
  data.frame(
    config         = nm,
    ahn_horenstein = round(results[[nm]]["ahn_horenstein"] * 100),
    bai_ng         = round(results[[nm]]["bai_ng"] * 100),
    onatski_2010   = round(results[[nm]]["onatski_2010"] * 100)
  )
}))

## ----lam_yao------------------------------------------------------------------
result_ly <- select_factors(X, method = "lam_yao", kmax = 8, h = 1)
print(result_ly)

## ----onatski_2009-------------------------------------------------------------
result_o09 <- select_factors(X, method = "onatski_2009",
                              kmax = 8, alpha = 0.05)
print(result_o09)

## ----onatski_2010-------------------------------------------------------------
result_o10 <- select_factors(X, method = "onatski_2010", kmax = 8)
print(result_o10)

