《財務時間序列分析》線上附錄
R09：美國十產業投資組合、CAPM 與 Fama–French 三因子
本附錄用美國十個產業的月報酬回答兩個實際問題。第一，估計共變異數矩陣時加入收縮，能否讓全域最小變異數（global minimum variance, GMV）投資組合在保留期更穩定？第二，加入 SMB 與 HML 後，Fama–French 三因子模型能否比 CAPM 更貼近同月已實現的產業超額報酬？前一題是投資組合風險比較，後一題是同期條件解釋，兩者的資訊時間與結論不能混在一起。
資料涵蓋 1967 年 1 月至 2021 年 11 月，共 659 個月、10 個產業。長資料中的一列是一個「產業—月份」觀察，所以共有 6,590 列；轉成矩陣後，一列是一個月份，一欄是一個產業。ret 是產業投資組合的月超額報酬，factor_ff_rf 是月無風險報酬，市場、SMB 與 HML 都是月因子報酬。所有報酬都以小數表示，例如 0.01 代表 1%。原課程程式透過 Kenneth French Data Library 取得因子與十產業投資組合，再合併 global-q 與 Welch–Goyal 總體變數；本附錄只使用其中的十產業、無風險利率與 FF3 欄位。
隨書提供的 data/processed CSV 是固定資料快照，讓程式在沒有網路時也能重做。若讀者改從上游來源下載，請記錄資料庫版本、下載日、FIZ/CIZ 方法版本、百分比轉小數，以及產業報酬減去無風險利率的步驟。CRSP 從 FIZ 改為 CIZ 時，同時調整了資料格式與月報酬累積方式，因此資料庫更新後，數字不必與這份固定快照完全相同。
knitr::opts_chunk$set(
echo = TRUE, message = FALSE, warning = FALSE,
fig.width = 7, fig.height = 4.3
)
stopifnot(getRversion() >= "4.3.0")
先確認一列資料代表什麼#
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)
}
panel_file <- locate_project_file(
"data/processed/ff_qf_macro_industries_1967_2021.csv"
)
manifest_file <- locate_project_file("data/processed/manifest.csv")
d <- read.csv(panel_file, stringsAsFactors = FALSE, check.names = FALSE)
manifest <- read.csv(manifest_file, stringsAsFactors = FALSE)
d$month <- as.Date(d$month)
# 排序後，同一月份的十個產業會相鄰，後面的矩陣轉換才有明確順序。
d <- d[order(d$month, d$industry), ]
manifest_key <- "data/processed/ff_qf_macro_industries_1967_2021.csv"
manifest_row <- manifest[manifest$file == manifest_key, , drop = FALSE]
actual_md5 <- unname(tools::md5sum(panel_file))
stopifnot(
nrow(manifest_row) == 1L,
nrow(d) == 6590L,
ncol(d) == 24L,
identical(actual_md5, manifest_row$md5),
!anyDuplicated(d[c("month", "industry")]),
length(unique(d$month)) == 659L,
length(unique(d$industry)) == 10L
)
data.frame(
first_month = min(d$month),
last_month = max(d$month),
months = length(unique(d$month)),
industries = length(unique(d$industry)),
return_unit = "decimal per month",
md5 = actual_md5
)
## first_month last_month months industries return_unit
## 1 1967-01-01 2021-11-01 659 10 decimal per month
## md5
## 1 69563611584d8a2dfd984ec6a53822a4
輸出中的日期與筆數先回答最基本的資料問題：每個「月份—產業」鍵只出現一次，而且 659 個月都各有 10 個產業。MD5 值只用來辨認本附錄使用哪一份固定快照；它不代表資料內容正確，也不取代來源與轉換紀錄。
下面把資料來源與本附錄實際使用的欄位並列，方便讀者分清楚原始合併檔與本題分析範圍。
data.frame(
component = c("十產業與 FF3", "固定合併面板", "本附錄輸出"),
source = c(
"Kenneth French Data Library；原課程以 frenchdata 下載",
"原課程 fffqmacro.R；另含 global-q 與 Welch–Goyal 欄位",
"隨書提供的 processed 固定快照"
),
use_here = c(
"ret、RF、MKT-RF、SMB、HML",
"只讀取 FF3 相關欄位",
"方法教學，不作投資建議或因果解讀"
)
)
## component source
## 1 十產業與 FF3 Kenneth French Data Library；原課程以 frenchdata 下載
## 2 固定合併面板 原課程 fffqmacro.R；另含 global-q 與 Welch–Goyal 欄位
## 3 本附錄輸出 隨書提供的 processed 固定快照
## use_here
## 1 ret、RF、MKT-RF、SMB、HML
## 2 只讀取 FF3 相關欄位
## 3 方法教學，不作投資建議或因果解讀
如何安排訓練、驗證與測試月份？#
先確認同一月份的因子與無風險利率在十個產業列完全一致，再把長資料轉成「月份乘產業」矩陣。這一步也把資訊單位改清楚：共變異數與因子迴歸都以月份為觀察，而不是把同月的十個產業誤當成十個獨立時間點。
common_names <- c(
"factor_ff_rf", "factor_ff_mkt_excess",
"factor_ff_smb", "factor_ff_hml"
)
common_check <- aggregate(
d[, common_names],
by = list(month = d$month),
FUN = function(x) max(x) - min(x)
)
stopifnot(max(as.matrix(common_check[, -1])) < 1e-12)
industries <- sort(unique(d$industry))
months <- sort(unique(d$month))
R_excess <- xtabs(ret ~ month + industry, data = d)
R_excess <- R_excess[, industries, drop = FALSE]
factor_by_month <- d[!duplicated(d$month), c("month", common_names)]
factor_by_month <- factor_by_month[match(months, factor_by_month$month), ]
rf <- factor_by_month$factor_ff_rf
R_total <- sweep(R_excess, 1, rf, "+")
stopifnot(
nrow(R_excess) == 659L,
ncol(R_excess) == 10L,
!anyNA(R_excess),
all(is.finite(R_excess)),
max(abs(rowMeans(R_total - R_excess) - rf)) < 1e-12
)
依時間順序把前 60% 設為訓練期、接續 20% 設為驗證期、最後 20% 留作測試期。訓練期估計初始共變異數；驗證期只選收縮強度；選定後再用訓練與驗證期重估最終權重。測試期只用一次，不參與權重、收縮強度或因子模型的選擇。
n_month <- length(months)
train_end <- floor(0.60 * n_month)
validation_end <- floor(0.80 * n_month)
train_id <- seq_len(train_end)
validation_id <- (train_end + 1L):validation_end
test_id <- (validation_end + 1L):n_month
stopifnot(max(train_id) < min(validation_id), max(validation_id) < min(test_id))
data.frame(
sample = c("訓練", "驗證", "測試"),
first_month = months[c(min(train_id), min(validation_id), min(test_id))],
last_month = months[c(max(train_id), max(validation_id), max(test_id))],
observations = c(length(train_id), length(validation_id), length(test_id))
)
## sample first_month last_month observations
## 1 訓練 1967-01-01 1999-11-01 395
## 2 驗證 1999-12-01 2010-11-01 132
## 3 測試 2010-12-01 2021-11-01 132
請把表中的三段日期視為研究設計的一部分。若先看過測試期再改候選網格或模型規格，最後一段就已參與調校，不能再當作未見的保留期。
收縮能否改善 GMV 的保留期表現？#
GMV 權重為
\[ \widehat w_{GMV}= \frac{\widehat\Sigma^{-1}\mathbf 1} {\mathbf 1^\top\widehat\Sigma^{-1}\mathbf 1}. \]
GMV 會反轉共變異數矩陣，因此很小、又估得不準的特徵值可能把權重放大。為降低這種不穩定，候選共變異數為
\[ \widehat\Sigma_\lambda=(1-\lambda)\widehat\Sigma+ \lambda\operatorname{diag}(\widehat\Sigma), \]
且只用驗證期實現波動選擇 \(\lambda\)。當 \(\lambda=0\) 時保留完整樣本共變異數；當 \(\lambda=1\) 時只保留各產業自己的變異數。這裡比較的是一小組事先列出的候選值，不是在測試期反覆調整。
gmv_weights <- function(Sigma, lambda = 0) {
stopifnot(lambda >= 0, lambda <= 1)
# 收縮只改變共變異數估計；所有候選權重仍滿足加總為一。
Sigma_use <- (1 - lambda) * Sigma + lambda * diag(diag(Sigma))
one <- rep(1, nrow(Sigma_use))
raw <- solve(Sigma_use, one)
drop(raw / sum(raw))
}
Sigma_train <- cov(R_total[train_id, , drop = FALSE])
lambda_grid <- c(0, 0.10, 0.25, 0.50, 0.75, 1)
validation_sd <- vapply(lambda_grid, function(lambda) {
w <- gmv_weights(Sigma_train, lambda)
sd(drop(R_total[validation_id, , drop = FALSE] %*% w))
}, numeric(1))
shrinkage_table <- data.frame(
lambda = lambda_grid,
validation_monthly_sd = validation_sd
)
selected_lambda <- lambda_grid[which.min(validation_sd)]
shrinkage_table
## lambda validation_monthly_sd
## 1 0.00 0.04387532
## 2 0.10 0.04185025
## 3 0.25 0.04111311
## 4 0.50 0.04110520
## 5 0.75 0.04189224
## 6 1.00 0.04337890
selected_lambda
## [1] 0.5
驗證期月標準差在 \(\lambda=0.50\) 時最低，因此後續固定使用 0.50。這不是「真實最佳參數」，而是這個候選網格、這段驗證期與這組產業下的選擇。選定 \(\lambda\) 後，合併訓練與驗證期重新估計一次；測試期仍未參與。
development_id <- c(train_id, validation_id)
Sigma_development <- cov(R_total[development_id, , drop = FALSE])
n_asset <- ncol(R_total)
weights <- rbind(
Equal = rep(1 / n_asset, n_asset),
GMV = gmv_weights(Sigma_development, 0),
Validation_Shrunk_GMV = gmv_weights(Sigma_development, selected_lambda)
)
colnames(weights) <- industries
stopifnot(max(abs(rowSums(weights) - 1)) < 1e-10)
round(weights, 3)
## Durbl Enrgy HiTec Hlth Manuf NoDur Other Shops Telcm
## Equal 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100
## GMV 0.019 0.128 0.000 0.165 0.084 0.276 -0.465 0.077 0.268
## Validation_Shrunk_GMV 0.013 0.115 0.016 0.126 0.056 0.157 -0.002 0.051 0.169
## Utils
## Equal 0.100
## GMV 0.447
## Validation_Shrunk_GMV 0.300
未收縮 GMV 出現較大的正負權重，反映矩陣反轉會放大產業間細微的共變異數差異；收縮後的權重較接近零，也較少依賴大幅放空。這是穩定化的直觀效果，但是否真的降低保留期風險，仍要看下一張表。
max_drawdown <- function(r) {
wealth <- c(1, cumprod(1 + r))
min(wealth / cummax(wealth) - 1)
}
test_performance <- t(vapply(seq_len(nrow(weights)), function(j) {
total <- drop(R_total[test_id, , drop = FALSE] %*% weights[j, ])
excess <- total - rf[test_id]
c(
annualized_mean_total = 12 * mean(total),
annualized_sd = sqrt(12) * sd(total),
annualized_sharpe = sqrt(12) * mean(excess) / sd(excess),
maximum_drawdown = max_drawdown(total)
)
}, numeric(4)))
rownames(test_performance) <- rownames(weights)
round(test_performance, 4)
## annualized_mean_total annualized_sd annualized_sharpe
## Equal 0.1433 0.1412 0.9768
## GMV 0.1045 0.1166 0.8512
## Validation_Shrunk_GMV 0.1178 0.1198 0.9389
## maximum_drawdown
## Equal -0.2291
## GMV -0.1861
## Validation_Shrunk_GMV -0.2224
在這個測試期，未收縮 GMV 的年化標準差約為 11.66%，收縮 GMV 約為 11.98%，都低於等權重的 14.12%；最低風險並沒有自動帶來最高平均報酬或 Sharpe 比率。收縮 GMV 的表現介於等權重與未收縮 GMV 之間，說明驗證期選出的穩定化程度不保證在每個保留期都勝過未收縮估計。
wealth <- sapply(seq_len(nrow(weights)), function(j) {
cumprod(1 + drop(R_total[test_id, , drop = FALSE] %*% weights[j, ]))
})
matplot(
months[test_id], wealth, type = "l", lty = 1, lwd = 1.5,
col = c("#173B57", "#A34045", "#1D6D73"),
xlab = "月份", ylab = "累積財富",
main = "十產業投資組合：固定權重的測試期表現"
)
legend(
"topleft", c("等權重", "未收縮 GMV", "驗證期收縮 GMV"),
col = c("#173B57", "#A34045", "#1D6D73"),
lty = 1, lwd = 1.5, bty = "n"
)
這是歷史保留期的描述性比較。程式假設在測試期持有固定權重，沒有納入交易成本、做空限制、估計誤差、多重嘗試，也沒有處理產業投資組合在實務上如何取得，因此不能直接視為可交易策略的績效。
FF3 是否比 CAPM 更能重建同月產業報酬？#
以下先以製造業（Manuf）為例，所有係數只用訓練期估計。測試月份的因子與產業報酬都是當月結束後才完整觀察，因此這裡問的是：給定同月已實現的市場、SMB 與 HML，模型能重建多少產業超額報酬？這是同期條件重建，不是月初可以執行的事前預測。
focus_industry <- "Manuf"
reg_data <- data.frame(
month = months,
y = as.numeric(R_excess[, focus_industry]),
MKT = factor_by_month$factor_ff_mkt_excess,
SMB = factor_by_month$factor_ff_smb,
HML = factor_by_month$factor_ff_hml
)
fit_capm <- lm(y ~ MKT, data = reg_data, subset = train_id)
fit_ff3 <- lm(y ~ MKT + SMB + HML, data = reg_data, subset = train_id)
newey_west_vcov <- function(model, lag = 6L) {
# 月資料可能同時有異質變異與短期序列相關，故以 Bartlett 權重計算 HAC。
X <- model.matrix(model)
e <- residuals(model)
n <- nrow(X)
Xe <- X * as.numeric(e)
meat <- crossprod(Xe)
if (lag > 0L) {
for (ell in seq_len(lag)) {
weight <- 1 - ell / (lag + 1)
Gamma <- crossprod(
Xe[(ell + 1L):n, , drop = FALSE],
Xe[seq_len(n - ell), , drop = FALSE]
)
meat <- meat + weight * (Gamma + t(Gamma))
}
}
bread <- solve(crossprod(X))
bread %*% meat %*% bread
}
V_ff3 <- newey_west_vcov(fit_ff3, lag = 6L)
data.frame(
term = names(coef(fit_ff3)),
estimate = unname(coef(fit_ff3)),
HAC_se_L6 = sqrt(diag(V_ff3)),
row.names = NULL
)
## term estimate HAC_se_L6
## 1 (Intercept) -0.001245651 0.0008102358
## 2 MKT 1.055043308 0.0188703733
## 3 SMB 0.031724972 0.0251161499
## 4 HML 0.049466265 0.0430464166
製造業的市場係數約為 1.06，表示在這段訓練樣本中，製造業超額報酬與市場超額報酬大致同向且幅度略高。SMB 與 HML 的點估計較小；HAC 標準誤反映序列相依與異質變異後的不確定性。這些係數描述條件關聯，不等於因子的因果效果。
score_reconstruction <- function(model, rows) {
prediction <- predict(model, newdata = reg_data[rows, ])
actual <- reg_data$y[rows]
c(
RMSE = sqrt(mean((actual - prediction)^2)),
MAE = mean(abs(actual - prediction)),
conditional_R2 = 1 - sum((actual - prediction)^2) /
sum((actual - mean(reg_data$y[train_id]))^2)
)
}
conditional_test <- rbind(
CAPM = score_reconstruction(fit_capm, test_id),
FF3 = score_reconstruction(fit_ff3, test_id)
)
round(conditional_test, 4)
## RMSE MAE conditional_R2
## CAPM 0.0157 0.0120 0.8826
## FF3 0.0152 0.0117 0.8905
製造業測試期的 CAPM RMSE 約為 0.0157，FF3 約為 0.0152；FF3 的同期重建誤差略低，條件 \(R^2\) 也由約 0.883 提高到 0.891。差距雖然方向一致，幅度不大，不能只憑這一個產業宣稱三因子模型全面勝出。
all_industry_rmse <- do.call(rbind, lapply(industries, function(industry) {
one <- transform(reg_data, y = as.numeric(R_excess[, industry]))
capm <- lm(y ~ MKT, data = one, subset = train_id)
ff3 <- lm(y ~ MKT + SMB + HML, data = one, subset = train_id)
data.frame(
industry = industry,
CAPM_test_RMSE = sqrt(mean((one$y[test_id] - predict(capm, one[test_id, ]))^2)),
FF3_test_RMSE = sqrt(mean((one$y[test_id] - predict(ff3, one[test_id, ]))^2))
)
}))
all_industry_rmse$FF3_minus_CAPM <-
all_industry_rmse$FF3_test_RMSE - all_industry_rmse$CAPM_test_RMSE
all_industry_rmse
## industry CAPM_test_RMSE FF3_test_RMSE FF3_minus_CAPM
## 1 Durbl 0.05813085 0.05715474 -0.0009761165
## 2 Enrgy 0.05894024 0.05692949 -0.0020107531
## 3 HiTec 0.01990691 0.01798862 -0.0019182973
## 4 Hlth 0.02467100 0.02520409 0.0005330911
## 5 Manuf 0.01574446 0.01519984 -0.0005446143
## 6 NoDur 0.02414882 0.02488379 0.0007349714
## 7 Other 0.01724476 0.01375786 -0.0034869024
## 8 Shops 0.02093788 0.02194264 0.0010047576
## 9 Telcm 0.02558006 0.02507557 -0.0005044857
## 10 Utils 0.03220611 0.03544833 0.0032422134
十個產業的結果並不一致：FF3 在耐久財、能源、高科技、製造業、其他產業與電信的 RMSE 較低，但在醫療、非耐久財、商店與公用事業較高。加入因子提高模型彈性，是否改善條件重建仍取決於產業與保留期間。
capm_test <- predict(fit_capm, newdata = reg_data[test_id, ])
ff3_test <- predict(fit_ff3, newdata = reg_data[test_id, ])
matplot(
months[test_id],
cbind(reg_data$y[test_id], capm_test, ff3_test),
type = "l", lty = c(1, 2, 3), lwd = c(1.6, 1.2, 1.2),
col = c("black", "#A34045", "#1D6D73"),
xlab = "月份", ylab = "月超額報酬（小數）",
main = "製造業：測試期的同期因子重建"
)
legend(
"topleft", c("實際值", "CAPM", "FF3"),
col = c("black", "#A34045", "#1D6D73"),
lty = c(1, 2, 3), lwd = c(1.6, 1.2, 1.2), bty = "n"
)
因此，FF3 測試誤差低於 CAPM 時，合宜的結論是：在這個固定切分下，加入同期 SMB 與 HML 改善了部分產業的條件解釋。這個結果沒有證明未來報酬可預測，也沒有識別因子的因果效果，更不是策略獲利的證據。
從結果回到研究問題#
投資組合部分顯示，GMV 的確降低了這段測試期的波動，但驗證期選出的收縮 GMV 並未在測試期進一步勝過未收縮 GMV。這不是收縮方法失效的普遍證據，而是提醒我們：估計穩定、權重不極端與單一保留期風險最低，是三個相關但不相同的目標。若要做更接近實務的比較，下一步可加入非負權重、交易成本與滾動重估。
因子模型部分顯示，FF3 對製造業以及部分產業的同期重建略有改善，但各產業方向不一。由於測試期因子是同月已實現值，這個練習回答的是條件解釋，而不是事前預測或因果問題。若研究目標改為報酬預測，就必須把解釋變數延遲到預測形成時已知的資訊。
重做本附錄時，請保留資料來源、版本、單位轉換與時間切分。若改用最新資料，應把新結果視為新的資料版本，而不是期待逐位數複製這份固定快照。
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.5.1
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
##
## time zone: Asia/Tokyo
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] compiler_4.5.2 cli_3.6.5 tools_4.5.2 otel_0.2.0 knitr_1.51
## [6] xfun_0.57 rlang_1.1.7 evaluate_1.0.5