---
title: "R11：S&P 500 平衡子樣本的 PCA 與跨股票落後值修正"
output:
  github_document:
    toc: true
    toc_depth: 3
---

在長格式股票資料中，若直接對整欄價格使用 `lag()`，會發生什麼事？又要如何在不讓未來資料進入中心化與負荷量估計的情況下，用 PCA 壓縮多檔股票報酬？本附錄先用一個小例子重現跨股票落後值的錯誤，再把同一套資料整理原則帶到真實股票報酬的時間切分、降維與跨期重建。

固定資料涵蓋 2013 年 1 月 3 日至 2022 年 6 月 22 日，共 2,384 個共同交易日與 89 檔股票。寬表中的一列是一個交易日，一欄是一檔股票；欄值是在各股票內按日期計算的日簡單報酬，採小數單位，例如 0.01 代表 1%。這 89 檔只是在共同日期都有觀察值的平衡子樣本，**不是完整的 S&P 500 指數成分股歷史**。平衡化可能帶來成分股選擇與存活者偏誤，因此以下結果用來說明資料整理與描述性降維，不衡量整體市場績效，也不構成交易策略或因果效果的證據。

原課程檔含價格、公司、權重與市場識別欄位，但原始供應商、成分股形成日與資料版本沒有完整留存。隨書提供 `sp500_returns_balanced_2013_2022.csv` 固定衍生面板，讓本 Rmd 可以直接重做；若要從原始市場資料重新建置，仍須補齊上游來源與形成日，再依「先按股票分組、組內排序、組內落後」的整理原則計算報酬。這項限制主要影響來源追溯、股票母體與經濟外推。

原課程的對應程式是
`slides/L09_Statistical_factor_models/W2L4_hands-on_R_factors/sp500/pca_sp500.R`：
第 38–48 行直接用 `stats::prcomp()` 做 PCA，第 50–66 行再用
`stats::factanal(..., factors = 3, rotation = "varimax", lower = 0.01)` 做三因子分析。
因此本附錄保留「先看清資料整理原則與時間切分，再呼叫成熟套件函數」的兩層教法；
`prcomp()` 代為處理中心化、尺度與矩陣分解，`factanal()` 代為處理最大概似估計與旋轉；學生不必重寫最佳化器，仍要自己決定時間切分、標準化方式、保留維度與解讀範圍。

```{r setup}
knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  fig.width = 7, fig.height = 4
)
stopifnot(getRversion() >= "4.3.0")
set.seed(1111)
```

## 先確認固定資料的期間、維度與版本

本檔不使用 `setwd()`、不安裝套件，也不在執行時下載資料。下列函數容許讀者從專案根目錄或 `online_appendix/` 執行，並在開始分析前確認日期順序、股票數、缺值與固定檔案版本。

```{r locate-files}
locate_project_file <- function(relative_path) {
  candidates <- c(
    relative_path,
    file.path("..", relative_path),
    file.path("../..", relative_path)
  )
  hit <- candidates[file.exists(candidates)]
  if (length(hit) == 0L) {
    stop("找不到專案檔案：", relative_path)
  }
  normalizePath(hit[1], mustWork = TRUE)
}

returns_file <- locate_project_file(
  "data/processed/sp500_returns_balanced_2013_2022.csv"
)
manifest_file <- locate_project_file("data/processed/manifest.csv")
```

```{r read-and-validate}
return_df <- read.csv(
  returns_file,
  check.names = FALSE,
  stringsAsFactors = FALSE
)
manifest <- read.csv(manifest_file, stringsAsFactors = FALSE)

dates <- as.Date(return_df$date)
R_all <- as.matrix(return_df[, -1, drop = FALSE])
storage.mode(R_all) <- "double"

manifest_key <- "data/processed/sp500_returns_balanced_2013_2022.csv"
manifest_row <- manifest[manifest$file == manifest_key, , drop = FALSE]
actual_md5 <- unname(tools::md5sum(returns_file))

stopifnot(
  nrow(manifest_row) == 1L,
  nrow(return_df) == manifest_row$rows,
  ncol(return_df) == manifest_row$columns,
  identical(actual_md5, manifest_row$md5),
  nrow(return_df) == 2384L,
  ncol(R_all) == 89L,
  !anyNA(dates),
  !anyDuplicated(dates),
  all(diff(dates) > 0),
  !anyNA(R_all),
  all(is.finite(R_all)),
  !anyDuplicated(colnames(R_all))
)

data.frame(
  first_date = min(dates),
  last_date = max(dates),
  trading_days = nrow(R_all),
  stocks = ncol(R_all),
  md5 = actual_md5
)
```

輸出應顯示 2,384 個遞增且不重複的交易日、89 檔股票，而且沒有缺值。檢查碼只確認本次分析讀到的是隨書固定版本；它不能替代上游供應者、成分股形成日、資料版本或經濟意義的查核。

## 為什麼不能在整張長表直接做 `lag()`？

若長格式價格表依日期與股票代碼排序後，直接對整欄價格取落後值，某一列的「前一期」往往會是另一檔股票。下列微型資料刻意重現此錯誤。

```{r cross-ticker-bug}
toy_price <- data.frame(
  date = rep(as.Date("2024-01-01") + 0:2, each = 2),
  symbol = rep(c("AAA", "BBB"), times = 3),
  price = c(100, 50, 110, 45, 121, 49.5)
)
toy_price <- toy_price[order(toy_price$date, toy_price$symbol), ]

# 錯誤：完全沒有按股票分組。
toy_price$return_wrong <- c(
  NA_real_,
  diff(toy_price$price) / head(toy_price$price, -1)
)
toy_price
```

輸出會讓錯誤一目了然：`BBB` 第一天拿 `AAA` 第一天的價格當分母，`AAA` 第二天又拿 `BBB` 第一天當分母。這些數字看似都是有限值，卻沒有任何金融意義。正確做法是先按股票分組、在每一組內按日期排序，然後才計算報酬。

```{r grouped-return-function}
within_symbol_return <- function(data) {
  required <- c("date", "symbol", "price")
  stopifnot(all(required %in% names(data)))

  groups <- split(data, data$symbol, drop = TRUE)
  corrected <- lapply(groups, function(one_stock) {
    one_stock <- one_stock[order(one_stock$date), , drop = FALSE]
    one_stock$return <- c(
      NA_real_,
      diff(one_stock$price) / head(one_stock$price, -1)
    )
    one_stock
  })
  corrected <- do.call(rbind, corrected)
  rownames(corrected) <- NULL
  corrected[order(corrected$date, corrected$symbol), ]
}

toy_correct <- within_symbol_return(toy_price[, c("date", "symbol", "price")])
toy_correct

# 可執行的方向與分組單元測試。
first_in_group <- !duplicated(toy_correct$symbol)
stopifnot(
  all(is.na(toy_correct$return[first_in_group])),
  isTRUE(all.equal(
    toy_correct$return[toy_correct$symbol == "AAA"][-1],
    c(0.10, 0.10), tolerance = 1e-12
  )),
  isTRUE(all.equal(
    toy_correct$return[toy_correct$symbol == "BBB"][-1],
    c(-0.10, 0.10), tolerance = 1e-12
  ))
)
```

單元測試確認 AAA 的兩期報酬都是 10%，BBB 則先跌 10%、再漲 10%；每檔股票第一筆價格沒有前一期，因此報酬保留為 `NA`，不硬填成 0。隨書固定面板也是依「先按股票代碼分組，再計算報酬」的程序建立。

## 如何安排估計、驗證與測試交易日？

先查看每檔股票日報酬的平均數、標準差與極端值。這些摘要不等同於投資績效比較，尤其資料是經平衡化後的子樣本。

```{r descriptive-statistics}
stock_summary <- data.frame(
  symbol = colnames(R_all),
  mean = colMeans(R_all),
  sd = apply(R_all, 2, sd),
  minimum = apply(R_all, 2, min),
  maximum = apply(R_all, 2, max)
)

stock_summary[order(stock_summary$sd, decreasing = TRUE)[1:10], ]
```

為避免未來資料影響中心化、尺度與主成分負荷量，本例固定切成連續三段：前 65% 為估計期、接著 15% 為驗證期、最後 20% 為測試期。估計期建立 PCA 與 80% 維度規則；驗證期只觀察固定負荷量的重建；選定維度後，估計期與驗證期合稱發展期，用來重估最後一組負荷量；測試期直到最後才使用。日資料不能先隨機打散再切分。

```{r chronological-split}
T_total <- nrow(R_all)
train_end <- floor(0.65 * T_total)
validation_end <- floor(0.80 * T_total)

train_id <- seq_len(train_end)
validation_id <- seq.int(train_end + 1L, validation_end)
test_id <- seq.int(validation_end + 1L, T_total)

split_table <- data.frame(
  sample = c("估計期", "驗證期", "測試期"),
  first_date = dates[c(min(train_id), min(validation_id), min(test_id))],
  last_date = dates[c(max(train_id), max(validation_id), max(test_id))],
  observations = c(length(train_id), length(validation_id), length(test_id))
)
split_table
```

表中的截止日決定每個統計量可以看見哪些資料。之後若改動分割比例，也應把它視為另一個研究設計，而不是在看過測試結果後微調。

## 套件作法：用 `prcomp()` 只在估計期建立 PCA

個股波動尺度不同，因此本例先用**估計期**平均數與標準差把各欄標準化，再做 PCA。驗證期與測試期不得各自重新標準化。

原課程直接使用 `prcomp()`。本附錄沿用這個作法，但把函數可見的資料明確限制在
估計期。函數代為計算中心、尺度、主成分方向與分數；研究者仍須決定是否標準化、
何時截止估計，以及後面如何選維度。其後用 `eigen()` 核對所有特徵值，而不是另造一個 PCA 估計器。

```{r fit-training-pca}
pca_train <- prcomp(
  R_all[train_id, , drop = FALSE],
  center = TRUE,
  scale. = TRUE
)

eigenvalues <- pca_train$sdev^2
pve <- eigenvalues / sum(eigenvalues)
explained <- data.frame(
  component = seq_along(pve),
  eigenvalue = eigenvalues,
  PVE = pve,
  cumulative_PVE = cumsum(pve)
)
head(explained, 15)
```

前 15 個特徵值讓我們看見共同變動並非集中在一、兩個方向。陡坡圖適合觀察下降速度，但本例不靠肉眼事後選維度，而是在下一節用事先寫好的 80% 累積解釋比例規則。

```{r eigenvalue-check}
Z_train <- scale(
  R_all[train_id, , drop = FALSE],
  center = pca_train$center,
  scale = pca_train$scale
)
direct_eigen <- eigen(cor(Z_train), symmetric = TRUE, only.values = TRUE)$values

stopifnot(isTRUE(all.equal(
  unname(eigenvalues), unname(direct_eigen), tolerance = 1e-9
)))
```

```{r scree-plot}
plot(
  explained$component[1:30],
  explained$eigenvalue[1:30],
  type = "b", pch = 19, col = "#173B57",
  xlab = "主成分", ylab = "特徵值",
  main = "估計期相關矩陣的前 30 個特徵值"
)
abline(h = 1, lty = 2, col = "#A34045")
```

### 負荷量與正負號不定性

特徵向量整欄乘以 \(-1\) 仍是同一個主成分。因此，不能把 PC1 負荷量的正負號當成可識別的經濟方向。下表依絕對值列出前兩個主成分的重要股票，並保留原符號以利重現。

```{r top-loadings}
largest_loadings <- function(fit, component, n = 12L) {
  loading <- fit$rotation[, component]
  keep <- order(abs(loading), decreasing = TRUE)[seq_len(n)]
  data.frame(
    component = paste0("PC", component),
    symbol = names(loading)[keep],
    loading = unname(loading[keep]),
    absolute_loading = abs(unname(loading[keep]))
  )
}

rbind(
  largest_loadings(pca_train, 1),
  largest_loadings(pca_train, 2)
)
```

將估計期 PC1 分數與 89 檔股票的等權平均報酬相比，可作為「共同移動」的描述性核對；這個等權組合不是 S&P 500 指數。

```{r pc1-equal-weight-check}
equal_weight_train <- rowMeans(R_all[train_id, , drop = FALSE])
pc1_train <- pca_train$x[, 1]

data.frame(
  correlation = cor(pc1_train, equal_weight_train),
  absolute_correlation = abs(cor(pc1_train, equal_weight_train))
)
```

PC1 與 89 檔股票等權平均報酬的絕對相關約為 0.995，支持「第一主成分捕捉廣泛共同移動」的描述。相關的正負號可隨主成分整欄反號而改變；而且這個等權平均不是市值加權 S&P 500 指數。

### 套件作法：用 `factanal()` 估計三個共同因子

原課程也提供統計因子模型的套件作法。以下完全沿用原程式的
`factanal()`、三因子、varimax 與 `lower = 0.01` 設定，但只把估計期送入函數。
它和 PCA 的目標不同：PCA 排序總變異的正交方向；最大概似因子分析則以
「共同因子加個別變異」近似相關矩陣。因此應比較共同性、個別變異與重建相關
矩陣的誤差，不要求兩者負荷量逐格相同。`factanal()` 代為進行最大概似估計與旋轉；
研究者仍須決定因子數、旋轉方式、個別變異下界，並檢查邊界解與整體配適。
`factanal()` 回報的 `PVAL` 是概似比檢定的 $p$ 值；只有在觀察值可視為獨立同分配的
多變量常態樣本、最大概似解位於參數空間內部等參考條件下，才可使用通常的
卡方近似。以下因此把它標成「名目概似比 $p$ 值」，並另外檢查是否有個別變異
停在 `lower = 0.01` 的人為下界。

```{r original-course-factanal-shortcut}
fa_lower_bound <- 0.01
fa_train <- stats::factanal(
  R_all[train_id, , drop = FALSE],
  factors = 3,
  method = "mle",
  rotation = "varimax",
  lower = fa_lower_bound
)
fa_loadings <- unclass(fa_train$loadings)
fa_at_lower_bound <-
  abs(fa_train$uniquenesses - fa_lower_bound) < 1e-7
fa_model_correlation <- tcrossprod(fa_loadings)
diag(fa_model_correlation) <-
  diag(fa_model_correlation) + fa_train$uniquenesses
fa_empirical_correlation <- cor(R_all[train_id, , drop = FALSE])

fa_top_loadings <- do.call(rbind, lapply(seq_len(ncol(fa_loadings)), function(j) {
  keep <- order(abs(fa_loadings[, j]), decreasing = TRUE)[1:10]
  data.frame(
    factor = paste0("Factor", j),
    symbol = rownames(fa_loadings)[keep],
    loading = fa_loadings[keep, j],
    absolute_loading = abs(fa_loadings[keep, j]),
    row.names = NULL
  )
}))

data.frame(
  estimator = c("stats::prcomp", "stats::factanal"),
  fitted_on = "training period only",
  main_output = c(
    "ordered principal-component directions and scores",
    "rotated common-factor loadings and uniquenesses"
  )
)
data.frame(
  factors = fa_train$factors,
  correlation_reconstruction_RMSE = sqrt(mean(
    (fa_empirical_correlation - fa_model_correlation)^2
  )),
  minimum_uniqueness = min(fa_train$uniquenesses),
  maximum_uniqueness = max(fa_train$uniquenesses),
  minimum_uniqueness_at_0_01_boundary = any(fa_at_lower_bound),
  boundary_symbols = paste(
    names(fa_train$uniquenesses)[fa_at_lower_bound],
    collapse = ", "
  ),
  nominal_likelihood_ratio_p_value = unname(fa_train$PVAL)
)
fa_top_loadings
```

三因子重建相關矩陣的 RMSE 約為 0.045，名目概似比 $p$ 值接近 0。在獨立同分配的多變量常態觀察且最大概似解位於參數空間內部的參考條件下，這會拒絕「三因子足以精確描述母體相關矩陣」的虛無假設。不過，最小個別變異正好碰到 0.01 下界，已是邊界解；日報酬的厚尾與時間相依也使獨立同分配多變量常態的參考條件可疑。因此這個數字只能稱為名目 $p$ 值，不能當成一般條件下仍有效的正式模型檢定。此時負荷量仍可用來做探索性描述，卻不宜把三因子模型說成完整或已通過配適檢查。

旋轉後的因子可交換次序或整欄反號；前十大負荷量是描述性標籤線索，不是產業
分類真值，也不能把統計共同因子自動命名為可交易風險因子。

## 要保留多少個主成分？

若以同一組變數的重建誤差選維度，保留全部 89 個主成分必然達到零誤差。為避免看過測試期才改規則，本例在分析前指定：選擇估計期累積解釋比例首次達 80% 的最小主成分個數。

```{r locked-dimension-rule}
target_pve <- 0.80
r_selected <- which(explained$cumulative_PVE >= target_pve)[1]

data.frame(
  target_PVE = target_pve,
  selected_components = r_selected,
  achieved_training_PVE = explained$cumulative_PVE[r_selected]
)
```

這份估計樣本需要 39 個主成分才首次達到 80%，累積比例約為 80.06%。89 檔股票並沒有被壓縮成極少數方向，這本身就是值得報告的資料特徵。80% 不是自然法則，而是一條事先寫明的教學規則；研究者也可改用經濟損失函數或驗證期準則，但必須在查看最終測試期以前決定。

## 固定的 39 維空間能否延伸到較晚期間？

下列函數固定使用估計樣本所得到的中心、尺度與負荷量，把新資料投影到前 \(r\) 個主成分後再重建。

```{r reconstruction-functions}
standardized_reconstruction <- function(fit, newdata, r) {
  stopifnot(r >= 0L, r <= ncol(fit$rotation))
  # 新交易日沿用舊樣本的中心、尺度與負荷量，不在保留期重新標準化。
  Z <- scale(newdata, center = fit$center, scale = fit$scale)
  if (r == 0L) {
    Z_hat <- matrix(0, nrow = nrow(Z), ncol = ncol(Z))
  } else {
    V <- fit$rotation[, seq_len(r), drop = FALSE]
    Z_hat <- (Z %*% V) %*% t(V)
  }
  list(actual = Z, reconstructed = Z_hat)
}

reconstruction_diagnostics <- function(object) {
  sse <- sum((object$actual - object$reconstructed)^2)
  sst <- sum(object$actual^2)
  c(
    standardized_MSE = mean((object$actual - object$reconstructed)^2),
    out_of_sample_fraction_reconstructed = 1 - sse / sst
  )
}
```

先用驗證期確認所選維度的表現，但不再改動 80% 規則。

```{r validation-diagnostics}
r_grid <- sort(unique(c(1L, 2L, 5L, 10L, 20L, r_selected, 40L, 60L, 89L)))
r_grid <- r_grid[r_grid <= ncol(R_all)]

validation_table <- do.call(rbind, lapply(r_grid, function(r) {
  diagnostic <- reconstruction_diagnostics(
    standardized_reconstruction(
      pca_train,
      R_all[validation_id, , drop = FALSE],
      r
    )
  )
  data.frame(r = r, t(diagnostic), row.names = NULL)
}))
validation_table
```

驗證表顯示，39 個成分重建約 84.9% 的標準化變動；保留全部 89 個成分當然會得到近乎零誤差，因此不能用同一組變數的重建誤差無限制地增加維度。這張表用來理解壓縮與重建的取捨，不再改動先前選定的 39。

選定 \(r\) 後，把估計期與驗證期合併重新估計中心、尺度與負荷量；最終測試期仍完全未參與估計或選擇。

```{r locked-final-test}
development_id <- c(train_id, validation_id)
pca_development <- prcomp(
  R_all[development_id, , drop = FALSE],
  center = TRUE,
  scale. = TRUE
)

test_object <- standardized_reconstruction(
  pca_development,
  R_all[test_id, , drop = FALSE],
  r_selected
)
test_diagnostic <- reconstruction_diagnostics(test_object)

data.frame(
  selected_components = r_selected,
  test_standardized_MSE = unname(test_diagnostic[1]),
  test_fraction_reconstructed = unname(test_diagnostic[2])
)
```

測試期標準化 MSE 約為 0.254，重建比例約為 78.8%，接近但略低於原先的 80% 目標。這是「同一期觀察值的降維重建」：程式先觀察交易日 $t$ 的 89 維報酬，才計算當日分數，不是用 $t-1$ 期預測 $t$ 期報酬。若研究問題是預測，還須另外建立只使用當時可取得資料的因子動態，並以真正的樣本外預測損失評量。

## 測試期分數的時間圖

測試期分數用發展期的固定負荷量計算，不會重新估計 PCA。圖形只呈現共同波動的相對尺度，PC 的正負號仍可整欄翻轉。

```{r test-score-plot}
Z_test <- scale(
  R_all[test_id, , drop = FALSE],
  center = pca_development$center,
  scale = pca_development$scale
)
test_scores <- Z_test %*% pca_development$rotation[, 1:2, drop = FALSE]

matplot(
  dates[test_id], test_scores,
  type = "l", lty = 1, col = c("#173B57", "#A34045"),
  xlab = "日期", ylab = "主成分分數",
  main = "以發展期負荷量計算的測試期分數"
)
legend(
  "topright", legend = c("第一主成分", "第二主成分"),
  col = c("#173B57", "#A34045"), lty = 1, bty = "n"
)
```

## 從降維結果回到資料與研究問題

這個案例最重要的第一步發生在 PCA 以前：報酬必須在股票內計算，否則後面再精密的矩陣分解也只是在整理錯誤資料。固定面板讓程式可以重做，但 89 檔平衡子樣本的形成方式仍限制了結果能否代表 S&P 500 全體成分股；若要作市場層面的推論，應補齊供應者、歷史成分股形成日與資料版本。

在降維部分，39 個主成分於測試期保留約 78.8% 的同日標準化變動。這表示資料可壓縮，但不是只靠一、兩個成分就足夠，也不表示未來報酬可預測。PC 的正負號可翻轉，旋轉因子的次序也可交換；經濟命名需要額外資料與理論，而不是只看負荷量最大的股票。

若下一步要做預測，應在每個預測起點重新確認可用股票母體，只用過去資料估計中心、尺度、負荷量與因子動態，再用未來實現值評量。若目標只是同日共變動描述，則可以把重點放在低秩近似的穩定性、股票宇宙敏感度與不同維度規則的比較。

```{r session-info}
sessionInfo()
```