《財務時間序列分析》線上附錄
R01：R、tidyverse 與手動 OLS
本附錄對應第 1 章，要回答一個適合初學 R、也很適合檢查公式的問題：把普通最小平方法（ordinary least squares, OLS）的矩陣公式直接寫成程式，會不會得到與 lm() 相同的係數？我們以 AAPL 當日報酬為應變數，將 MSFT 與 NVDA 的等權平均當日報酬作為解釋變數，依序完成資料整理、手動 OLS、套件結果比較與殘差診斷。
資料共有 2,384 個共同交易日，期間為 2013-01-03 至 2022-06-22。每一列是一個交易日，每一欄是一檔股票的日簡單報酬，報酬以小數表示；例如 0.01 代表 1%。解釋變數使用的是與 AAPL 同一天收盤後才完整可得的報酬，因此這裡研究的是同日線性關係，不是盤前或即時預測。全樣本都用於這項代數練習，沒有另分訓練期、驗證期與測試期，也沒有足以支持因果解讀的識別設計。
knitr::opts_chunk$set(
echo = TRUE, message = FALSE, warning = FALSE,
fig.width = 7, fig.height = 4.5
)
stopifnot(requireNamespace("dplyr", quietly = TRUE))
stopifnot(requireNamespace("tibble", quietly = TRUE))
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tibble)
root_candidates <- c(".", "..")
is_root <- vapply(root_candidates, function(x) {
file.exists(file.path(x, "main.tex"))
}, logical(1))
stopifnot(any(is_root))
project_root <- root_candidates[which(is_root)[1]]
project_path <- function(...) file.path(project_root, ...)
先把觀察單位與排序確認清楚#
讀檔後先只留下本例需要的日期與三檔股票。日期轉成 R 的 Date 類別後再排序，才能確定第 $t-1$ 列真的早於第 $t$ 列。固定面板已保留三檔股票都有報酬的共同交易日，所以這一頁不做填補或刪除缺值。
data_path <- project_path(
"data", "processed", "sp500_returns_balanced_2013_2022.csv"
)
manifest_path <- project_path("data", "processed", "manifest.csv")
stopifnot(file.exists(data_path), file.exists(manifest_path))
raw <- read.csv(data_path, check.names = FALSE)
returns <- raw |>
transmute(
date = as.Date(date),
AAPL = as.numeric(AAPL),
MSFT = as.numeric(MSFT),
NVDA = as.numeric(NVDA)
) |>
arrange(date)
# 若日期重複、順序錯誤或任一欄缺值，先回到資料入口處理。
stopifnot(!anyNA(returns), !anyDuplicated(returns$date))
stopifnot(all(diff(returns$date) > 0))
glimpse(returns)
## Rows: 2,384
## Columns: 4
## $ date <date> 2013-01-03, 2013-01-04, 2013-01-07, 2013-01-08, 2013-01-09, 2013…
## $ AAPL <dbl> -0.012622072, -0.027854471, -0.005882822, 0.002691439, -0.0156289…
## $ MSFT <dbl> -0.0133961372, -0.0187154612, -0.0018699229, -0.0052453126, 0.005…
## $ NVDA <dbl> 0.0007860068, 0.0329928166, -0.0288975086, -0.0219260948, -0.0224…
glimpse() 應顯示 2,384 列與四個欄位；三個報酬欄都是數值，日期則為 Date。資料由原課程價格長表建立，建檔時先在每一檔股票內依日期計算報酬，再保留 89 檔都有觀察值的共同交易日。來源、建檔方式與 MD5 記錄見 data/DATA_SOURCES.md。
向量、資料框與管線#
接下來先熟悉兩個最常用的 R 物件。first_five 是只含五個數字的向量，returns 則是每列代表一個交易日的資料框；先分清楚物件的形狀，後面建立設計矩陣時才不會把列與欄顛倒。
first_five <- returns$AAPL[1:5]
object_summary <- tibble(
object = c("first_five", "returns"),
class = c(class(first_five)[1], class(returns)[1]),
length_or_rows = c(length(first_five), nrow(returns))
)
object_summary
## # A tibble: 2 × 3
## object class length_or_rows
## <chr> <chr> <int>
## 1 first_five numeric 5
## 2 returns data.frame 2384
現在用 mutate() 建立 MSFT 與 NVDA 的等權平均報酬 tech_equal，再用 summarise() 按年度整理 AAPL 的觀察值數、平均數與標準差。年度表的每一列是一個曆年；因原檔已是簡單報酬，mean_aapl 與 sd_aapl 仍以小數表示，而不是百分比。
analysis_df <- returns |>
mutate(
tech_equal = (MSFT + NVDA) / 2,
year = as.integer(format(date, "%Y"))
)
annual_summary <- analysis_df |>
group_by(year) |>
summarise(
observations = n(),
mean_aapl = mean(AAPL),
sd_aapl = sd(AAPL),
.groups = "drop"
)
annual_summary
## # A tibble: 10 × 4
## year observations mean_aapl sd_aapl
## <int> <int> <dbl> <dbl>
## 1 2013 251 0.000347 0.0179
## 2 2014 252 0.00145 0.0136
## 3 2015 252 0.0000199 0.0168
## 4 2016 252 0.000575 0.0147
## 5 2017 251 0.00164 0.0111
## 6 2018 251 -0.0000573 0.0181
## 7 2019 252 0.00266 0.0165
## 8 2020 253 0.00281 0.0294
## 9 2021 252 0.00131 0.0158
## 10 2022 118 -0.00202 0.0227
這張表先讓我們看出各年交易日數與波動是否大致合理。若某一年觀察值突然很少，應先檢查資料涵蓋與共同交易日規則，再進入迴歸；年度平均數的正負也只是樣本描述，不宜當作穩定的未來報酬率。
手動 OLS#
令 $y$ 為 AAPL 當日報酬，$x$ 為 MSFT、NVDA 的等權平均當日報酬。設計矩陣 $X$ 的每一列對應一個共同交易日，第一欄全為 1，用來估計截距，第二欄才是 tech_equal。在計算 $(X^\top X)^{-1}X^\top y$ 以前，先確認 $X$ 的列數與 $y$ 的長度相同，而且兩欄線性獨立。
y <- analysis_df$AAPL
X <- cbind(Intercept = 1, tech_equal = analysis_df$tech_equal)
# 維度不合或設計矩陣不滿秩時，OLS 公式不能照寫下去。
stopifnot(nrow(X) == length(y), qr(X)$rank == ncol(X))
beta_manual <- solve(crossprod(X), crossprod(X, y))
fitted_manual <- as.numeric(X %*% beta_manual)
residual_manual <- y - fitted_manual
tibble(
term = rownames(beta_manual),
estimate = as.numeric(beta_manual)
)
## # A tibble: 2 × 2
## term estimate
## <chr> <dbl>
## 1 Intercept 0.000124
## 2 tech_equal 0.569
crossprod(X) 等於 $X^\top X$，crossprod(X, y) 等於 $X^\top y$。輸出中的截距約為 0.000124，斜率約為 0.5693。就這份樣本而言，MSFT 與 NVDA 的等權平均當日報酬每增加 1 個百分點，AAPL 的線性投影值約增加 0.569 個百分點；這是同日共同波動的條件關聯，不表示前兩檔股票造成 AAPL 上漲。正式計算通常偏好 QR 分解，因為直接解正規方程式在近共線時較不穩定；此處保留矩陣公式，是為了讓課文中的推導與 R 運算一一對上。
用 lm() 與 QR 分解比較同一個答案#
手動結果算出後，再讓 qr.solve() 與 lm() 解同一個問題。三欄係數的名稱與次序必須一致，否則即使數字相近，也可能是在比較不同參數。
fit_lm <- lm(AAPL ~ tech_equal, data = analysis_df)
beta_qr <- qr.solve(X, y)
comparison <- tibble(
term = names(coef(fit_lm)),
manual = as.numeric(beta_manual),
qr = as.numeric(beta_qr),
lm = as.numeric(coef(fit_lm))
)
comparison
## # A tibble: 2 × 4
## term manual qr lm
## <chr> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.000124 0.000124 0.000124
## 2 tech_equal 0.569 0.569 0.569
# 三種演算法應在浮點計算容許範圍內得到相同係數。
stopifnot(
isTRUE(all.equal(as.numeric(beta_manual), as.numeric(beta_qr))),
isTRUE(all.equal(as.numeric(beta_manual), as.numeric(coef(fit_lm))))
)
comparison 的三個係數欄應重合。這項結果說明我們正確地把 OLS 公式翻成 R；它沒有比較模型好壞，也沒有替殘差的時間序列性質背書。
正規方程式與配適分解#
含截距的 OLS 殘差應與設計矩陣的每一欄正交。這裡直接計算 $X^\top \hat e$，也把總平方和分成模型解釋部分與殘差平方和，確認係數、配適值與殘差彼此一致。
orthogonality <- as.numeric(crossprod(X, residual_manual))
names(orthogonality) <- colnames(X)
orthogonality
## Intercept tech_equal
## 5.342948e-16 2.107689e-16
stopifnot(max(abs(orthogonality)) < 1e-10)
sst <- sum((y - mean(y))^2)
sse <- sum(residual_manual^2)
ssr <- sum((fitted_manual - mean(y))^2)
c(SST = sst, explained = ssr, SSE = sse, difference = sst - ssr - sse)
## SST explained SSE difference
## 7.771318e-01 2.989755e-01 4.781562e-01 -4.996004e-16
orthogonality 的兩個數應非常接近零，平方和表中的 difference 也應只剩浮點誤差。若差距明顯，應先檢查是否漏了截距、資料列是否錯位，或配適值與殘差是否來自同一個模型。
在計算以前先決定如何處理缺值#
固定面板目前沒有缺值，但換成別份資料時，刪除觀察值可能改變樣本期間與研究對象。以下函數的一般設定是遇到缺值就停止；只有呼叫者把 allow_drop 明確設為 TRUE 時，才保留指定欄位都完整的列。把決定集中在資料準備階段，比在每一個統計量中分別加入 na.rm = TRUE 更容易看出實際用了哪些觀察值。
prepare_complete <- function(data, columns, allow_drop = FALSE) {
stopifnot(all(columns %in% names(data)))
ok <- complete.cases(data[, columns, drop = FALSE])
# 先要求研究者選擇處理方式，避免各函數各自悄悄丟掉不同列。
if (!all(ok) && !allow_drop) {
stop("指定欄位含缺值；請決定刪除、填補或修正來源。")
}
data[ok, , drop = FALSE]
}
toy <- tibble(y = c(1, 2, NA), x = c(0, 1, 2))
expected_error <- try(
prepare_complete(toy, c("y", "x")),
silent = TRUE
)
inherits(expected_error, "try-error")
## [1] TRUE
prepare_complete(toy, c("y", "x"), allow_drop = TRUE)
## # A tibble: 2 × 2
## y x
## <dbl> <dbl>
## 1 1 0
## 2 2 1
第一次呼叫應得到預期錯誤，第二次則只留下前兩列。實際研究若允許刪除，還要報告刪除了多少列、落在哪些日期，以及缺值是否可能與市場狀況有關。
診斷圖#
係數相同只代表三種計算法一致，還不能回答線性模型是否充分。殘差對配適值圖用來找非線性或變異程度改變，常態 Q–Q 圖則讓我們觀察尾端是否偏離常態；兩張圖的每一點仍對應一個共同交易日。
par(mfrow = c(1, 2))
plot(
fitted_manual, residual_manual,
xlab = "配適值", ylab = "殘差",
pch = 16, cex = 0.45, col = "#173B57"
)
abline(h = 0, lty = 2, col = "#A34045")
qqnorm(residual_manual, pch = 16, cex = 0.45, col = "#173B57")
qqline(residual_manual, col = "#A34045", lwd = 2)
par(mfrow = c(1, 1))
圖形若呈現厚尾、極端點或波動群聚，OLS 係數的代數結果仍成立，但傳統獨立同分配標準誤便未必合適。下一步應檢查殘差的自相關與條件異質變異，並視研究目的採用時間序列穩健推論或重新指定模型。
這份練習回答了什麼#
這份固定資料顯示 AAPL 與另外兩檔大型科技股的同日報酬具有正向線性關聯，而且手動公式、QR 分解與 lm() 給出相同係數。它還沒有告訴我們這項關聯在未來是否穩定，也沒有因果含意。若要把問題改成預測，解釋變數必須在預測起點已經可得，例如用第 $t$ 日資訊預測第 $t+1$ 日報酬，並依第 8 章以時間排序的訓練期、驗證期與測試期評估。