R08：AAPL 的 ARCH、GARCH、非對稱波動與預測

AAPL 報酬的方向可能很難預測，波動幅度卻會連續幾天偏高。這種「報酬弱相關、平方報酬強相關」的資料特徵，是否能由 GARCH 描述？同樣幅度的負報酬是否對下一期波動帶來較大的條件關聯？本附錄對應第 11–12 章，以真實 AAPL 日對數報酬估計 GARCH(1,1) 與 GJR–GARCH(1,1)，並把答案放到最後的測試期比較。

資料源自原課程 S&P 500 價格檔的固定版本。每筆觀察值是一個交易日的 AAPL 對數報酬，有效樣本從 2019-01-03 到 2022-06-22，共 874 筆，原始單位為小數；圖中乘以 100 後才是百分比。價格調整與建檔方式見 data/DATA_SOURCES.md。

程式按日期排序後，把前 80% 定為訓練期，最後 20% 保留為一次性測試期；本例沒有另設驗證期。平均數、變異數參數與初始尺度都只由訓練期估計。進入測試期後，程式在目標日報酬出現以前先形成一步波動預測，再用實現報酬更新下一日狀態。手動版用 R 內建的 optim() 展開條件變異遞迴與目標函數；套件版沿用原課程的 fGarch::garchFit() 工作流程，並使用同一訓練期核對結果。

這些模型估計的是歷史條件變異動態與樣本外預測損失。正負報酬的不對稱參數可以描述條件關聯，若要把它命名為「壞消息造成波動」，還需要可信的消息度量與因果識別設計。

knitr::opts_chunk$set(
  echo = TRUE, message = FALSE, warning = FALSE,
  fig.width = 8, fig.height = 5,
  dev = "ragg_png", dpi = 144,
  dev.args = list(background = "white")
)

required_r <- "4.3.0"
stopifnot(getRversion() >= required_r)

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, ...)

stopifnot(
  requireNamespace("ragg", quietly = TRUE),
  requireNamespace("systemfonts", quietly = TRUE)
)
cwtex_file <- project_path("assets", "fonts", "cwTeXQKai-Medium.ttf")
stopifnot(file.exists(cwtex_file))
if (!"cwTeX Online" %in% systemfonts::registry_fonts()$family) {
  systemfonts::register_font("cwTeX Online", cwtex_file)
}
plot_family <- "cwTeX Online"

先確認資料與預測時間線#

波動預測最容易在切分時洩漏未來資訊。下列程式先排序日期、排除無效報酬，再一次固定訓練與測試分界。平均報酬只在訓練期計算；否則使用全樣本平均數中心化，已經讓測試期資訊進入模型。

aapl <- read.csv(project_path(
  "data", "processed", "aapl_adjusted_daily_2019_2022.csv"
))
aapl$date <- as.Date(aapl$date)
# 落後殘差與條件變異遞迴都假設列已按日期先後排列。
aapl <- aapl[order(aapl$date), ]
aapl <- aapl[is.finite(aapl$log_return), ]
row.names(aapl) <- NULL

stopifnot(
  !anyNA(aapl$date), !anyNA(aapl$log_return),
  all(diff(aapl$date) > 0)
)

n <- nrow(aapl)
train_end <- floor(0.80 * n)
train_dates <- aapl$date[seq_len(train_end)]
test_dates <- aapl$date[(train_end + 1L):n]

# 常數平均數只用訓練期估計，之後固定。
mu_train <- mean(aapl$log_return[seq_len(train_end)])
train <- aapl$log_return[seq_len(train_end)] - mu_train
test <- aapl$log_return[(train_end + 1L):n] - mu_train

split_table <- data.frame(
  區段 = c("訓練期", "測試期"),
  起日 = c(min(train_dates), min(test_dates)),
  迄日 = c(max(train_dates), max(test_dates)),
  觀察值 = c(length(train), length(test)),
  報酬單位 = "日對數報酬，小數",
  資料來源 = "原課程 S&P 500 價格檔的 AAPL 固定版本",
  check.names = FALSE
)
knitr::kable(split_table)

data.frame(
  訓練期平均日對數報酬 = mu_train,
  訓練期報酬標準差 = sd(train),
  check.names = FALSE
)

##   訓練期平均日對數報酬 訓練期報酬標準差
## 1          0.001879715       0.02195292

切分表給出了實際的訓練起迄日、測試起迄日與各段樣本數。以下定義 $a_t=r_t-\widehat{\mu}_{\mathrm{train}}$，並在後續模型中固定這個訓練期平均數。這種做法可以清楚交代形成波動預測時可用到哪些資料；它和同時估計平均數與變異數方程的聯合最大概似規格是兩個不同的比較基準。

報酬方向與波動幅度有不同的時間結構嗎？#

第一張圖把訓練與測試分界當作時間軸上的參考點，用報酬與絕對中心化報酬觀察波動群聚。ACF 只在訓練期計算，因為我們不希望根據測試期圖形回頭改選模型。

old_par <- par(
  mfrow = c(2, 1), mar = c(4.5, 4, 3, 1),
  family = plot_family
)
plot(
  aapl$date, 100 * aapl$log_return,
  type = "l", col = "#173B57",
  xlab = "日期", ylab = "日對數報酬（%）"
)
abline(v = max(train_dates), lty = 2, col = "#A34045")
plot(
  aapl$date, 100 * abs(aapl$log_return - mu_train),
  type = "l", col = "#A34045",
  xlab = "日期", ylab = "絕對中心化報酬（%）"
)
abline(v = max(train_dates), lty = 2, col = "#173B57")

par(old_par)

old_par <- par(
  mfrow = c(1, 2), mar = c(4.2, 4, 4, 1),
  family = plot_family, cex.main = 0.9
)
acf(train, lag.max = 30, main = "中心化報酬 ACF")
acf(train^2, lag.max = 30, main = "平方報酬 ACF")

par(old_par)

中心化報酬的 ACF 主要檢查方向性的線性關聯，平方報酬 ACF 則檢查幅度是否有持續性。若第一張 ACF 大多接近零，第二張卻有多個落後期明顯偏離零，接下來的合理決定是保持簡單平均數規格，並另建條件變異數模型。

ARCH–LM：殘差平方中還有可預測結構嗎？#

ARCH–LM 以訓練期殘差平方的落後值作為解釋變數。虛無假設是這些落後項的係數同時為零，也就是這項檢定沒有發現 ARCH 結構。

arch_lm <- function(residual, lags = 10L) {
  stopifnot(lags >= 1L, length(residual) > 5L * lags)
  x2 <- residual^2
  n <- length(x2)
  # 前 lags 期沒有完整的落後向量，因此從 lags + 1 期開始共同對齊。
  response <- x2[(lags + 1L):n]
  X <- sapply(seq_len(lags), function(j) {
    x2[(lags + 1L - j):(n - j)]
  })
  colnames(X) <- paste0("lag", seq_len(lags))
  fit <- lm(response ~ X)
  statistic <- nobs(fit) * summary(fit)$r.squared
  data.frame(
    落後階數 = lags,
    LM統計量 = statistic,
    卡方近似p值 = pchisq(
      statistic, df = lags, lower.tail = FALSE
    ),
    check.names = FALSE
  )
}

knitr::kable(arch_lm(train, lags = 10), digits = 6)

訓練期 ARCH–LM 統計量約為 166.93，卡方近似 $p$ 值在目前數值精度下顯示為 0。這是強烈的樣本證據，顯示常數條件變異沒有捕捉到殘差平方的落後關聯。因此下一步估計 GARCH 類模型是有資料依據的；但階數、不對稱項與創新分配仍需要以理論、殘差診斷和樣本外損失一起選擇。

手動作法：建立 GARCH 與 GJR–GARCH 條件變異遞迴#

對稱 GARCH 先回答「大幅度衝擊是否會讓波動持續」；GJR–GARCH 再加上衝擊符號，查看同樣幅度的負殘差是否對下一期變異數有額外關聯。以下將參數轉換、變異數濾波與概似目標分開寫，讓每個限制可以直接檢查。

對稱 GARCH(1,1) 設定 $\gamma=0$；GJR 模型使用

\[ h_t=\omega+\alpha a_{t-1}^2 +\gamma\mathbf 1(a_{t-1}<0)a_{t-1}^2 +\beta h_{t-1}. \]

在對稱創新下，GJR 的二階持續性以 $\alpha+\beta+\gamma/2$ 衡量。下列參數轉換強制 $\omega>0$、$\alpha\geq0$、$\beta\geq0$、$\gamma\geq0$，並把持續性限制在 0.999 以下。

map_garch <- function(eta) {
  raw <- exp(eta[2:3])
  denominator <- 1 + sum(raw)
  c(
    omega = exp(eta[1]),
    alpha = 0.999 * raw[1] / denominator,
    beta = 0.999 * raw[2] / denominator,
    gamma = 0
  )
}

map_gjr <- function(eta) {
  raw <- exp(eta[2:4])
  denominator <- 1 + sum(raw)
  c(
    omega = exp(eta[1]),
    alpha = 0.999 * raw[1] / denominator,
    beta = 0.999 * raw[2] / denominator,
    gamma = 2 * 0.999 * raw[3] / denominator
  )
}

filter_variance <- function(a, par) {
  n <- length(a)
  h <- numeric(n)
  # 初始變異數只由傳入的訓練資料計算，避免用到未來尺度。
  h[1] <- var(a)
  for (t in 2:n) {
    h[t] <- par["omega"] +
      par["alpha"] * a[t - 1]^2 +
      par["gamma"] * as.numeric(a[t - 1] < 0) * a[t - 1]^2 +
      par["beta"] * h[t - 1]
  }
  h
}

gaussian_nll <- function(eta, a, model = c("garch", "gjr")) {
  model <- match.arg(model)
  par <- if (model == "garch") map_garch(eta) else map_gjr(eta)
  h <- filter_variance(a, par)
  if (any(!is.finite(h)) || any(h <= 0)) return(1e100)
  0.5 * sum(log(2 * pi) + log(h[-1]) + a[-1]^2 / h[-1])
}

這裡使用常態準最大概似（quasi-maximum likelihood, QML）來估計條件變異遞迴。常態密度在這裡是目標函數，不是對 AAPL 無條件報酬分配的描述。由於程式沒有計算穩健的三明治型標準誤（sandwich standard errors），下面將重點放在參數方向、條件變異遞迴、診斷與樣本外損失，不以這份輸出做係數顯著性推論。

為什麼要用多個起點做數值最佳化？#

條件變異模型的目標函數可能對初始值敏感。若只跑一次 optim()，「成功收斂」只表示演算法在那個起點附近停下，未必是可比較的最佳解。下面為每個模型使用三組具經濟意義的起點，檢查收斂碼與目標函數，再保留可接受解中的最小值。

garch_start <- function(v, alpha, beta) {
  persistence <- alpha + beta
  slack <- 0.999 - persistence
  stopifnot(slack > 0)
  c(
    log(v * (1 - persistence)),
    log(alpha / slack),
    log(beta / slack)
  )
}

gjr_start <- function(v, alpha, beta, gamma) {
  persistence <- alpha + beta + gamma / 2
  slack <- 0.999 - persistence
  stopifnot(slack > 0)
  c(
    log(v * (1 - persistence)),
    log(alpha / slack),
    log(beta / slack),
    log((gamma / 2) / slack)
  )
}

fit_multistart <- function(starts, a, model) {
  # 同一模型的所有起點共用同一訓練樣本與目標函數，才能直接比較。
  fits <- lapply(starts, function(start) {
    optim(
      start, gaussian_nll, a = a, model = model,
      method = "BFGS",
      control = list(maxit = 3000, reltol = 1e-10)
    )
  })
  acceptable <- vapply(fits, function(z) {
    z$convergence == 0 && is.finite(z$value)
  }, logical(1))
  stopifnot(any(acceptable))
  candidates <- which(acceptable)
  fits[[candidates[which.min(vapply(
    fits[candidates], function(z) z$value, numeric(1)
  ))]]]
}

v_train <- var(train)
garch_starts <- list(
  garch_start(v_train, 0.08, 0.87),
  garch_start(v_train, 0.10, 0.80),
  garch_start(v_train, 0.05, 0.92)
)
gjr_starts <- list(
  gjr_start(v_train, 0.05, 0.88, 0.08),
  gjr_start(v_train, 0.05, 0.84, 0.10),
  gjr_start(v_train, 0.03, 0.90, 0.10)
)

fit_garch <- fit_multistart(garch_starts, train, "garch")
fit_gjr <- fit_multistart(gjr_starts, train, "gjr")
par_garch <- map_garch(fit_garch$par)
par_gjr <- map_gjr(fit_gjr$par)

parameter_table <- rbind(
  GARCH = c(
    par_garch,
    persistence = unname(
      par_garch["alpha"] + par_garch["beta"]
    ),
    Gaussian_QML_objective = fit_garch$value
  ),
  GJR = c(
    par_gjr,
    persistence = unname(
      par_gjr["alpha"] + par_gjr["beta"] +
        par_gjr["gamma"] / 2
    ),
    Gaussian_QML_objective = fit_gjr$value
  )
)
knitr::kable(parameter_table, digits = 7)

stopifnot(
  par_garch["alpha"] + par_garch["beta"] < 0.999,
  par_gjr["alpha"] + par_gjr["beta"] + par_gjr["gamma"] / 2 < 0.999
)

GJR 的不對稱係數估計值約為 0.1614，訓練期高斯 QML 目標函數也較低；這是條件關聯與配適結果，沒有穩健標準誤時不作顯著性宣稱，更不能作因果解讀。

兩個模型使用同一訓練資料與常態 QML 目標，因此目標函數可直接比較。GJR 多了一個不對稱參數，訓練期彈性也較高；是否對預測有實際幫助，要留到固定測試期，用共同損失函數回答。

套件作法：用 garchFit() 完成平均數與波動模型估計#

原課程在 slides/L07_ARCH_GARCH/W2L2_R_template_GARCH.R 使用這套工作流程。程式以 garchFit(~ arma(1,0) + garch(1,1), cond.dist = "QMLE") 聯合估計 AR(1) 平均數與 GARCH(1,1)，並以 garchFit(~ arma(1,0) + aparch(1,1), delta = 2) 示範非對稱規格。 以下保留兩個公式，但把即時下載的 AAPL 換成同一份固定 CSV，且一律只用訓練期。garchFit() 代為處理條件變異遞迴、參數限制與數值最佳化；研究者仍須決定平均數階數、波動規格、創新分配、訓練截止日，以及哪些診斷結果足以支持下一步。

為了與前面的手動 GARCH 作公平核對，另外配適一個「對齊規格」：先扣除已固定的 $\widehat{\mu}_{\mathrm{train}}$，再設定 include.mean = FALSE。這個版本與手動 GARCH 使用相同的固定平均數與變異數方程；兩者仍可能因初始變異數與概似實作細節而有 小幅差異。APARCH 的 gamma1 則屬於套件的符號與冪次參數化，不能直接當成前面 GJR 指標函數中的 $\gamma$ 比較。

stopifnot(requireNamespace("fGarch", quietly = TRUE))

raw_train <- aapl$log_return[seq_len(train_end)]

# 先重現原課程的對稱規格：平均數與 GARCH 方程一起估計。
fgarch_lecture_garch <- fGarch::garchFit(
  ~ arma(1, 0) + garch(1, 1),
  data = raw_train,
  cond.dist = "QMLE",
  trace = FALSE
)

# 再重現原課程的非對稱 APARCH 規格。這裡的 gamma1 採 APARCH
# 參數化，不能直接和前面 GJR 指標函數中的 gamma 比大小。
fgarch_lecture_aparch <- fGarch::garchFit(
  ~ arma(1, 0) + aparch(1, 1),
  data = raw_train,
  delta = 2,
  include.delta = FALSE,
  cond.dist = "norm",
  trace = FALSE
)

# 最後另估一個固定平均數的對齊規格，目的在和手動 GARCH
# 比較同一條變異數方程，而不是拿不同平均數模型硬作比較。
fgarch_aligned <- fGarch::garchFit(
  ~ garch(1, 1),
  data = train,
  include.mean = FALSE,
  cond.dist = "QMLE",
  trace = FALSE
)

coefficient_or_na <- function(coefficient, term) {
  if (term %in% names(coefficient)) {
    unname(coefficient[term])
  } else {
    NA_real_
  }
}

coef_lecture_garch <- fgarch_lecture_garch@fit$coef
coef_lecture_aparch <- fgarch_lecture_aparch@fit$coef
coef_aligned <- fgarch_aligned@fit$coef

# 先把三種規格的共同參數排在同一張表；不適用的欄位保留 NA，
# 免得學生誤以為三種模型使用完全相同的參數化。
fgarch_specification_table <- data.frame(
  規格 = c(
    "原課程 AR(1)-GARCH；QMLE",
    "原課程 AR(1)-APARCH；delta=2",
    "對齊手動版：固定平均 GARCH；QMLE"
  ),
  平均數處理 = c(
    "聯合估計 AR(1)",
    "聯合估計 AR(1)",
    "固定訓練期平均數"
  ),
  mu = c(
    coefficient_or_na(coef_lecture_garch, "mu"),
    coefficient_or_na(coef_lecture_aparch, "mu"),
    mu_train
  ),
  ar1 = c(
    coefficient_or_na(coef_lecture_garch, "ar1"),
    coefficient_or_na(coef_lecture_aparch, "ar1"),
    NA_real_
  ),
  omega = c(
    coefficient_or_na(coef_lecture_garch, "omega"),
    coefficient_or_na(coef_lecture_aparch, "omega"),
    coefficient_or_na(coef_aligned, "omega")
  ),
  alpha1 = c(
    coefficient_or_na(coef_lecture_garch, "alpha1"),
    coefficient_or_na(coef_lecture_aparch, "alpha1"),
    coefficient_or_na(coef_aligned, "alpha1")
  ),
  beta1 = c(
    coefficient_or_na(coef_lecture_garch, "beta1"),
    coefficient_or_na(coef_lecture_aparch, "beta1"),
    coefficient_or_na(coef_aligned, "beta1")
  ),
  gamma1_APARCH = c(
    NA_real_,
    coefficient_or_na(coef_lecture_aparch, "gamma1"),
    NA_real_
  ),
  check.names = FALSE
)
knitr::kable(fgarch_specification_table, digits = 7)

# 套件對齊版必須滿足正變異數與定態限制，才值得進入下一步比較。
par_fgarch_aligned <- c(
  omega = coefficient_or_na(coef_aligned, "omega"),
  alpha = coefficient_or_na(coef_aligned, "alpha1"),
  beta = coefficient_or_na(coef_aligned, "beta1"),
  gamma = 0
)
stopifnot(
  all(is.finite(par_fgarch_aligned)),
  par_fgarch_aligned["omega"] > 0,
  par_fgarch_aligned["alpha"] >= 0,
  par_fgarch_aligned["beta"] >= 0,
  par_fgarch_aligned["alpha"] + par_fgarch_aligned["beta"] < 1
)

# fGarch 與手動 optim 的內部目標函數定義未必完全相同；因此兩組
# 估計值都代回同一個高斯目標函數，才有可直接比較的數值基準。
common_gaussian_objective <- function(a, par) {
  h <- filter_variance(a, par)
  0.5 * sum(
    log(2 * pi) + log(h[-1]) + a[-1]^2 / h[-1]
  )
}

# 這張表只比較「固定平均數」的兩個對齊版本；原課程的 AR(1)
# 平均數規格留在上表，不混入這項數值核對。
aligned_comparison <- rbind(
  手動_optim固定平均 = c(
    par_garch[c("omega", "alpha", "beta")],
    persistence = unname(par_garch["alpha"] + par_garch["beta"]),
    共同目標函數 = common_gaussian_objective(train, par_garch)
  ),
  fGarch固定平均 = c(
    par_fgarch_aligned[c("omega", "alpha", "beta")],
    persistence = unname(
      par_fgarch_aligned["alpha"] + par_fgarch_aligned["beta"]
    ),
    共同目標函數 = common_gaussian_objective(
      train, par_fgarch_aligned
    )
  )
)
knitr::kable(aligned_comparison, digits = 7)

前一張規格表的第一、二列回答原課程的聯合平均數—波動問題；第三列才是用來核對本附錄手動固定平均版本的橋梁。緊接著的兩列對照表只比較後者與手動 optim()。若兩個對齊版本的參數不完全相同，應先查看初始化與概似定義，而不是把差異誤認為資料或公式錯誤。

標準化殘差診斷#

估計完成後，先問兩件事：標準化殘差 $z_t=a_t/\sqrt{h_t}$ 是否仍有線性相依？其平方 $z_t^2$ 是否仍保留波動群聚？前者若顯著，應回頭檢討平均數模型；後者若顯著，表示條件變異方程還沒有吸收完可預測的波動結構。

diagnose_variance <- function(a, par, label) {
  h <- filter_variance(a, par)
  z <- a / sqrt(h)
  q_z <- Box.test(z, lag = 20, type = "Ljung-Box")
  q_z2 <- Box.test(z^2, lag = 20, type = "Ljung-Box")
  data.frame(
    模型 = label,
    Q20_標準化殘差 = unname(q_z$statistic),
    p_標準化殘差 = q_z$p.value,
    Q20_平方標準化殘差 = unname(q_z2$statistic),
    p_平方標準化殘差 = q_z2$p.value,
    check.names = FALSE
  )
}

diagnostic_table <- rbind(
  diagnose_variance(train, par_garch, "GARCH"),
  diagnose_variance(train, par_gjr, "GJR"),
  diagnose_variance(
    train, par_fgarch_aligned, "fGarch 固定平均 GARCH"
  )
)
knitr::kable(diagnostic_table, digits = 6)

三個規格的 $z_t$ 與 $z_t^2$ Ljung–Box $p$ 值都高於 0.20；在落後 20 期的這項檢查中，沒有明顯證據顯示線性相依或平方相依仍未被吸收。表中使用 fitdf = 0，沒有扣除平均數與波動參數的估計自由度，因此應把它讀成方便比較的近似診斷，而不是已完整校正的規格檢定。這不代表模型已經正確，仍要查看尾端形狀、符號不對稱、參數是否貼近邊界，以及結果在不同樣本期間是否穩定。若其他樣本出現 $z_t$ 顯著，應先檢討平均數模型；若 $z_t^2$ 顯著，則應修改波動規格。

無資料洩漏的一步波動預測#

現在把問題從樣本內配適改成真正的預測：在目標日 $t$ 的報酬尚未出現時，只用截至 $t-1$ 日的資訊形成 $h_{t\mid t-1}$。第一個測試日以訓練期最後一日為預測起點，後續各日才依序納入已實現的測試期報酬；歷史 60 期變異數則提供一個簡單、同樣遵守時間邊界的基準。

one_step_h <- function(previous_a, previous_h, par) {
  par["omega"] +
    par["alpha"] * previous_a^2 +
    par["gamma"] * as.numeric(previous_a < 0) * previous_a^2 +
    par["beta"] * previous_h
}

forecast_locked <- function(train, test, par) {
  h_train <- filter_variance(train, par)
  previous_a <- tail(train, 1)
  previous_h <- tail(h_train, 1)
  forecast <- numeric(length(test))

  for (j in seq_along(test)) {
    forecast[j] <- one_step_h(previous_a, previous_h, par)
    # test[j] 在 forecast[j] 形成後，才進入下一期狀態。
    previous_a <- test[j]
    previous_h <- forecast[j]
  }
  forecast
}

hhat_garch <- forecast_locked(train, test, par_garch)
hhat_gjr <- forecast_locked(train, test, par_gjr)
hhat_fgarch <- forecast_locked(train, test, par_fgarch_aligned)

# 60 期歷史變異數基準：每一期只用目標日以前的觀察。
all_centered <- c(train, test)
hhat_hist60 <- numeric(length(test))
for (j in seq_along(test)) {
  target_index <- train_end + j
  past_index <- (target_index - 60L):(target_index - 1L)
  hhat_hist60[j] <- var(all_centered[past_index])
}

stopifnot(
  all(hhat_garch > 0),
  all(hhat_gjr > 0),
  all(hhat_fgarch > 0),
  all(hhat_hist60 > 0)
)

QLIKE 與平方損失#

qlike <- function(realized_square, variance_forecast) {
  log(variance_forecast) + realized_square / variance_forecast
}

loss_row <- function(h, label) {
  data.frame(
    模型 = label,
    QLIKE = mean(qlike(test^2, h)),
    平方損失_x1e8 = 1e8 * mean((test^2 - h)^2),
    平均預測日波動_pct = 100 * mean(sqrt(h)),
    check.names = FALSE
  )
}

loss_table <- rbind(
  loss_row(hhat_hist60, "歷史 60 期"),
  loss_row(hhat_garch, "GARCH"),
  loss_row(hhat_gjr, "GJR"),
  loss_row(hhat_fgarch, "fGarch 固定平均 GARCH")
)
knitr::kable(loss_table, digits = 8)

兩種損失都是越小越好。QLIKE 由 GJR 最低（約 $-6.807$），平方損失則由手動 GARCH 最低（放大後約 42.29）；兩個指標沒有選出同一模型。套件對齊版與手動 GARCH 的結果很接近，因為兩者使用相同固定平均數、逐期資訊集合與損失函數；這一列主要反映估計實作與初始化差異。

test^2 是不可觀察條件變異數的高雜訊代理，因此 QLIKE 與平方損失給出不同排序並不意外。合宜的結論是：GJR 的不對稱項改善了 QLIKE，卻沒有同時改善平方損失；不能事後只挑其中一個指標宣告全面勝出。

old_par <- par(family = plot_family)
plot(
  test_dates, 100 * abs(test),
  type = "l", col = "gray70",
  xlab = "目標日期", ylab = "幅度或預測標準差（%）"
)
lines(test_dates, 100 * sqrt(hhat_garch), col = "#173B57", lwd = 1.5)
lines(test_dates, 100 * sqrt(hhat_gjr), col = "#A34045", lwd = 1.5)
lines(test_dates, 100 * sqrt(hhat_fgarch), col = "#7A5C99", lwd = 1.3)
lines(test_dates, 100 * sqrt(hhat_hist60), col = "#3F7158", lwd = 1.3)
legend(
  "topright",
  c(
    "|中心化報酬|", "GARCH", "GJR",
    "fGarch 固定平均 GARCH", "歷史 60 期"
  ),
  col = c(
    "gray70", "#173B57", "#A34045", "#7A5C99", "#3F7158"
  ),
  lty = 1, lwd = c(1, 1.5, 1.5, 1.3, 1.3),
  bty = "n", cex = 0.8
)

par(old_par)

用日期確認每個預測真的只看過去#

origin_dates <- c(max(train_dates), head(test_dates, -1))
timing_table <- data.frame(
  目標日期 = test_dates,
  可用到的最後日期 = origin_dates,
  時序正確 = origin_dates < test_dates,
  check.names = FALSE
)
knitr::kable(head(timing_table, 8))

stopifnot(all(timing_table$時序正確))

模擬只作單元檢查#

下列固定種子模擬不產生任何 AAPL 實證數字，只檢查 GJR 遞迴是否保持正變異數，以及同幅度負面衝擊是否在 $\gamma>0$ 時給出較高的下一期變異數。

simulate_gjr <- function(
    n, omega, alpha, beta, gamma,
    burn = 300L, seed = 808L) {
  persistence <- alpha + beta + gamma / 2
  stopifnot(
    omega > 0, alpha >= 0, beta >= 0,
    gamma >= 0, persistence < 1
  )
  set.seed(seed)
  total <- n + burn
  z <- rnorm(total)
  a <- numeric(total)
  h <- numeric(total)
  h[1] <- omega / (1 - persistence)
  a[1] <- sqrt(h[1]) * z[1]
  for (t in 2:total) {
    h[t] <- omega + alpha * a[t - 1]^2 +
      gamma * as.numeric(a[t - 1] < 0) * a[t - 1]^2 +
      beta * h[t - 1]
    a[t] <- sqrt(h[t]) * z[t]
  }
  keep <- (burn + 1L):total
  data.frame(a = a[keep], h = h[keep])
}

unit_truth <- c(
  omega = 2e-6, alpha = 0.05, beta = 0.88, gamma = 0.08
)
unit_sim <- simulate_gjr(
  n = 500,
  omega = unit_truth["omega"],
  alpha = unit_truth["alpha"],
  beta = unit_truth["beta"],
  gamma = unit_truth["gamma"]
)

h_reference <- mean(unit_sim$h)
positive_next <- one_step_h(0.02, h_reference, unit_truth)
negative_next <- one_step_h(-0.02, h_reference, unit_truth)

unit_table <- data.frame(
  檢查 = c("模擬最小變異數", "+2% 衝擊的下一期變異數", "-2% 衝擊的下一期變異數"),
  數值 = c(min(unit_sim$h), positive_next, negative_next),
  check.names = FALSE
)
knitr::kable(unit_table, digits = 8)

stopifnot(
  all(is.finite(unit_sim$h)),
  all(unit_sim$h > 0),
  negative_next > positive_next
)

從結果回到波動預測決定#

AAPL 的估計、殘差診斷與測試期損失都來自同一份固定資料；最後的模擬只用來確認遞迴方向與正變異數，不參與實證模型排名。手動 optim() 與原課程採用的 garchFit() 並列後，可以分辨結果差異究竟來自模型規格，還是初始化與概似實作。

實際選模時，先依標準化殘差診斷排除明顯遺漏，再以固定測試期的 QLIKE 與平方損失比較預測。若損失函數給出不同排序，應如實報告，而不是事後挑選有利指標。這個練習仍只涵蓋單一股票與一次時間切割；要判斷結論能否推廣，還需要滾動視窗、多個資產與不同市場期間。

若要在報告中重做同一結果，請一併記錄 R 與 fGarch 版本、訓練截止日、多起點設定、所有收斂碼、損失函數及逐期預測值。這些資訊讓讀者知道模型在何時、看過哪些資料後形成每一個預測。

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] Matrix_1.7-4        gtable_0.3.6        dplyr_1.2.1        
##  [4] compiler_4.5.2      gbutils_0.5.1       fBasics_4052.98    
##  [7] tidyselect_1.2.1    Rcpp_1.1.0          cvar_0.6           
## [10] parallel_4.5.2      systemfonts_1.3.2   scales_1.4.0       
## [13] timeSeries_4052.112 textshaping_1.0.5   lattice_0.22-7     
## [16] ggplot2_4.0.3       R6_2.6.1            generics_0.1.4     
## [19] fGarch_4052.93      knitr_1.51          rbibutils_2.4.1    
## [22] tibble_3.3.0        spatial_7.3-18      forecast_9.0.2     
## [25] timeDate_4052.112   pillar_1.11.1       RColorBrewer_1.1-3 
## [28] rlang_1.1.7         urca_1.3-4          xfun_0.57          
## [31] S7_0.2.2            otel_0.2.0          cli_3.6.5          
## [34] magrittr_2.0.4      Rdpack_2.6.6        grid_4.5.2         
## [37] lifecycle_1.0.5     nlme_3.1-168        fracdiff_1.5-4     
## [40] vctrs_0.7.2         evaluate_1.0.5      glue_1.8.0         
## [43] farver_2.1.2        ragg_1.5.2          zoo_1.8-15         
## [46] colorspace_2.1-3    tools_4.5.2         pkgconfig_2.0.3