--- title: "R01：R、tidyverse 與手動 OLS" output: github_document: toc: true toc_depth: 3 --- 本附錄對應第 1 章，要回答一個適合初學 R、也很適合檢查公式的問題：把普通最小平方法（ordinary least squares, OLS）的矩陣公式直接寫成程式，會不會得到與 `lm()` 相同的係數？我們以 AAPL 當日報酬為應變數，將 MSFT 與 NVDA 的等權平均當日報酬作為解釋變數，依序完成資料整理、手動 OLS、套件結果比較與殘差診斷。資料共有 2,384 個共同交易日，期間為 2013-01-03 至 2022-06-22。每一列是一個交易日，每一欄是一檔股票的日簡單報酬，報酬以小數表示；例如 0.01 代表 1%。解釋變數使用的是與 AAPL 同一天收盤後才完整可得的報酬，因此這裡研究的是同日線性關係，不是盤前或即時預測。全樣本都用於這項代數練習，沒有另分訓練期、驗證期與測試期，也沒有足以支持因果解讀的識別設計。 ``` r 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 ``` ``` r 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$ 列。固定面板已保留三檔股票都有報酬的共同交易日，所以這一頁不做填補或刪除缺值。 ``` r 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 2013-01-03, 2013-01-04, 2013-01-07, 2013-01-08, 2013-01-09, 2013… ## $ AAPL -0.012622072, -0.027854471, -0.005882822, 0.002691439, -0.0156289… ## $ MSFT -0.0133961372, -0.0187154612, -0.0018699229, -0.0052453126, 0.005… ## $ NVDA 0.0007860068, 0.0329928166, -0.0288975086, -0.0219260948, -0.0224… ``` `glimpse()` 應顯示 2,384 列與四個欄位；三個報酬欄都是數值，日期則為 `Date`。資料由原課程價格長表建立，建檔時先在每一檔股票內依日期計算報酬，再保留 89 檔都有觀察值的共同交易日。來源、建檔方式與 MD5 記錄見 `data/DATA_SOURCES.md`。 ## 向量、資料框與管線接下來先熟悉兩個最常用的 R 物件。`first_five` 是只含五個數字的向量，`returns` 則是每列代表一個交易日的資料框；先分清楚物件的形狀，後面建立設計矩陣時才不會把列與欄顛倒。 ``` r 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 ## ## 1 first_five numeric 5 ## 2 returns data.frame 2384 ``` 現在用 `mutate()` 建立 MSFT 與 NVDA 的等權平均報酬 `tech_equal`，再用 `summarise()` 按年度整理 AAPL 的觀察值數、平均數與標準差。年度表的每一列是一個曆年；因原檔已是簡單報酬，`mean_aapl` 與 `sd_aapl` 仍以小數表示，而不是百分比。 ``` r 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 ## ## 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$ 的長度相同，而且兩欄線性獨立。 ``` r 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 ## ## 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()` 解同一個問題。三欄係數的名稱與次序必須一致，否則即使數字相近，也可能是在比較不同參數。 ``` r 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 ## ## 1 (Intercept) 0.000124 0.000124 0.000124 ## 2 tech_equal 0.569 0.569 0.569 ``` ``` r # 三種演算法應在浮點計算容許範圍內得到相同係數。 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$，也把總平方和分成模型解釋部分與殘差平方和，確認係數、配適值與殘差彼此一致。 ``` r orthogonality <- as.numeric(crossprod(X, residual_manual)) names(orthogonality) <- colnames(X) orthogonality ``` ``` ## Intercept tech_equal ## 5.342948e-16 2.107689e-16 ``` ``` r 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` 更容易看出實際用了哪些觀察值。 ``` r 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 ``` ``` r prepare_complete(toy, c("y", "x"), allow_drop = TRUE) ``` ``` ## # A tibble: 2 × 2 ## y x ## ## 1 1 0 ## 2 2 1 ``` 第一次呼叫應得到預期錯誤，第二次則只留下前兩列。實際研究若允許刪除，還要報告刪除了多少列、落在哪些日期，以及缺值是否可能與市場狀況有關。 ## 診斷圖係數相同只代表三種計算法一致，還不能回答線性模型是否充分。殘差對配適值圖用來找非線性或變異程度改變，常態 Q–Q 圖則讓我們觀察尾端是否偏離常態；兩張圖的每一點仍對應一個共同交易日。 ``` r 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) ``` ![OLS 殘差與配適值；此圖只診斷線性投影，不代表預測或因果。](../R01_r_tidyverse_manual_ols_files/figure-gfm/diagnostics-1.png) ``` r par(mfrow = c(1, 1)) ``` 圖形若呈現厚尾、極端點或波動群聚，OLS 係數的代數結果仍成立，但傳統獨立同分配標準誤便未必合適。下一步應檢查殘差的自相關與條件異質變異，並視研究目的採用時間序列穩健推論或重新指定模型。 ## 這份練習回答了什麼這份固定資料顯示 AAPL 與另外兩檔大型科技股的同日報酬具有正向線性關聯，而且手動公式、QR 分解與 `lm()` 給出相同係數。它還沒有告訴我們這項關聯在未來是否穩定，也沒有因果含意。若要把問題改成預測，解釋變數必須在預測起點已經可得，例如用第 $t$ 日資訊預測第 $t+1$ 日報酬，並依第 8 章以時間排序的訓練期、驗證期與測試期評估。