GARCH bất đối xứng và đuôi dày: EGARCH, GJR và Student-t

Trong Phần 1 của loạt bài này, chúng ta đã xây dựng GARCH(1,1) từ đầu: trực giác về hiện tượng phân cụm biến động (volatility clustering), phép truy hồi phương sai có điều kiện (conditional variance), ước lượng hợp lý cực đại (maximum likelihood), dự báo, và các chẩn đoán phần dư chuẩn với thư viện arch. Nếu bạn chưa đọc, hãy bắt đầu từ đó — bài viết này giả định rằng bạn đã có thể khớp và diễn giải một mô hình GARCH(1,1) thuần túy, và sẽ không suy dẫn lại các phần cơ bản.

GARCH(1,1) thuần túy là một mô hình nền tốt nhưng lại là một câu trả lời cuối cùng tồi. Nó có hai khiếm khuyết về mặt cấu trúc mà việc bỏ qua thì rẻ trong một backtest nhưng đắt khi làm với vốn thật. Thứ nhất, nó đối xứng: mô hình phản ứng với một ngày $+5\%$ y hệt như phản ứng với một ngày $-5\%$, bởi vì cú sốc đi vào phép truy hồi phương sai chỉ thông qua bình phương của nó, $\varepsilon_{t-1}^2$. Việc bình phương làm mất đi dấu. Thứ hai, nó giả định các cải tiến (innovations) tuân theo phân phối Gaussian: ngay cả sau khi GARCH đã hấp thụ hết hiện tượng phân cụm biến động, phần dư chuẩn hóa của BTC và ETH vẫn có đuôi dày rõ rệt, và hàm hợp lý Gaussian định giá đuôi thấp một cách có hệ thống. VaR 99% của GARCH(1,1)-Normal sẽ bị vi phạm nhiều hơn nhiều so với 1% số ngày.

Bài viết này khắc phục cả hai khiếm khuyết. Chúng ta bổ sung tính bất đối xứng bằng GJR-GARCH và EGARCH, và đuôi dày bằng các cải tiến Student-$t$ và skew-$t$ của Hansen. Sau đó chúng ta làm điều thực sự sinh lời: biến phân phối có điều kiện đã khớp thành một dự báo Value-at-Risk và Expected Shortfall một bước, rồi backtest dự báo đó một cách trung thực bằng các kiểm định Kupiec và Christoffersen. Một mô hình biến động mà bạn không bao giờ kiểm định rủi ro chỉ là vật trang trí.

Hiệu ứng đòn bẩy, và vì sao crypto lộn xộn hơn

Trong thị trường cổ phiếu, tính bất đối xứng có một cái tên và một câu chuyện. Hiệu ứng đòn bẩy (leverage effect) (Black, 1976): khi cổ phiếu của một công ty giảm, tỷ lệ nợ trên vốn chủ sở hữu tăng, vốn cổ phần trở nên rủi ro hơn về mặt cơ học, và biến động tăng lên. Tin xấu làm tăng biến động tương lai nhiều hơn tin tốt có cùng độ lớn. Về mặt thực nghiệm, đây là một trong những sự kiện cách điệu (stylized fact) vững chắc nhất trong tài liệu về biến động cổ phiếu.

Crypto không có vốn cổ phần và không có đòn bẩy bảng cân đối kế toán theo nghĩa doanh nghiệp, thế nhưng một tính bất đối xứng giống hiệu ứng đòn bẩy vẫn xuất hiện trong phần lớn thời gian — được thúc đẩy bởi việc giảm đòn bẩy bắt buộc (forced deleveraging) chứ không phải bởi kế toán. Khi BTC giảm mạnh, các khoản vay được thế chấp vượt mức bị thanh lý, các vị thế mua (long) hợp đồng tương lai vĩnh cửu bị đóng cưỡng bức, funding đảo chiều, và chuỗi phản ứng dây chuyền nuôi thêm biến động. Vậy nên cơ chế thì khác nhưng dấu thường trùng với thị trường cổ phiếu: các đợt giảm giá làm biến động vọt lên nhiều hơn.

Lưu ý quan trọng: crypto lộn xộn hơn, và bạn nên xem tính bất đối xứng như một câu hỏi thực nghiệm chứ không phải một quy luật. Các đợt tăng giá dữ dội — short squeeze, một đợt bùng nổ tăng giá nhờ đòn bẩy, một cú gap khi ETF được phê duyệt — cũng có thể làm biến động thực hiện (realized volatility) vọt lên. Tùy vào tài sản và cửa sổ mẫu, tính bất đối xứng ước lượng được có thể mạnh, yếu, hoặc đôi khi mang dấu "sai". Kỷ luật mà bài viết này nhấn mạnh: khớp mô hình bất đối xứng, xem tham số bất đối xứng có ý nghĩa thống kê và có đúng chiều kỳ vọng hay không, và chỉ giữ lại tham số bổ sung nếu nó xứng đáng với vị trí của mình. Đừng giả định rằng câu chuyện của thị trường cổ phiếu chuyển giao được sang crypto; hãy kiểm định nó.

Kiểm định tính bất đối xứng trước khi mô hình hóa nó

Phần trên nói "hãy xem tính bất đối xứng là thực nghiệm" — vậy nên trước khi khớp một mô hình bất đối xứng, hãy chạy một kiểm định hình thức chi phí thấp để xem tính bất đối xứng có thực sự hiện diện hay không. Các kiểm định sai lệch dấu (sign-bias) của Engle-Ng (1993) làm chính xác điều này. Trước tiên khớp một GARCH(1,1) đối xứng, lấy bình phương phần dư chuẩn hóa $z_t^2$ của nó, và hồi quy chúng theo các biến chỉ báo về dấu và độ lớn của cú sốc trước đó:

$z_t^2 = a_0 + a_1 S_{t-1}^- + a_2 S_{t-1}^- \varepsilon_{t-1} + a_3 S_{t-1}^+ \varepsilon_{t-1} + u_t$

trong đó $S_{t-1}^- = \mathbf{1}\{\varepsilon_{t-1} < 0\}$ và $S_{t-1}^+ = 1 - S_{t-1}^-$. Logic là: nếu mô hình đối xứng đã nắm bắt hết mọi thứ, thì dấu và độ lớn của cú sốc hôm qua không nên dự báo được bình phương phần dư hôm nay, nên $a_1 = a_2 = a_3 = 0$. Các kiểm định $t$ riêng lẻ là kiểm định sai lệch dấu ($a_1$), sai lệch độ lớn âm ($a_2$), và sai lệch độ lớn dương ($a_3$); một kiểm định $F$ đồng thời trên cả ba là kiểm định tổng quát. Một $a_1$ hoặc $a_2$ có ý nghĩa cho thấy các cú sốc âm bị mô hình đối xứng định giá sai một cách có hệ thống — đó là dấu hiệu cho thấy GJR hoặc EGARCH sẽ có ích.

import statsmodels.api as sm

def sign_bias_test(symmetric_res):
    """Engle-Ng sign-bias tests on a fitted symmetric GARCH result."""
    z = symmetric_res.std_resid.dropna()
    z2 = (z ** 2).values[1:]
    eps_lag = z.values[:-1]                      # standardized shock proxy
    neg = (eps_lag < 0).astype(float)
    X = np.column_stack([
        np.ones_like(eps_lag),                   # intercept
        neg,                                     # sign bias
        neg * eps_lag,                           # negative size bias
        (1 - neg) * eps_lag,                     # positive size bias
    ])
    ols = sm.OLS(z2, X).fit()
    names = ["const", "sign_bias", "neg_size_bias", "pos_size_bias"]
    for nm, coef, t, p in zip(names, ols.params, ols.tvalues, ols.pvalues):
        print(f"{nm:16s} coef={coef:+.4f}  t={t:+.2f}  p={p:.3f}")
    print(f"Joint F p-value: {ols.f_pvalue:.4f}")
    return ols

sign_bias_test(models["GARCH-N"])

Nếu kiểm định $F$ đồng thời không có ý nghĩa, bạn có căn cứ thực nghiệm để giữ mô hình đối xứng và tiết kiệm hai tham số. Nếu nó có ý nghĩa — kết quả thường gặp đối với BTC/ETH — hãy tiến tới GJR/EGARCH với lương tâm thanh thản, biết rằng bạn đang mô hình hóa một đặc trưng có thật chứ không phải đuổi theo nhiễu. Đây chính là kỷ luật thực nghiệm mà phần mở đầu đòi hỏi: đừng giả định câu chuyện đòn bẩy của thị trường cổ phiếu, hãy kiểm định nó.

GJR-GARCH: Tính bất đối xứng qua một số hạng ngưỡng

Mô hình Glosten-Jagannathan-Runkle (1993) — đôi khi được gọi là TGARCH hay threshold GARCH — là chỉnh sửa nhỏ nhất có thể đối với GARCH(1,1) để tin xấu và tin tốt có tác động khác nhau. Nhớ lại phép truy hồi phương sai có điều kiện đối xứng từ Phần 1:

$\sigma_t^2 = \omega + \alpha\,\varepsilon_{t-1}^2 + \beta\,\sigma_{t-1}^2$

GJR thêm một số hạng ngưỡng: một liều phương sai bổ sung chỉ được kích hoạt sau một cú sốc âm.

$\sigma_t^2 = \omega + \alpha\,\varepsilon_{t-1}^2 + \gamma\,I_{t-1}\,\varepsilon_{t-1}^2 + \beta\,\sigma_{t-1}^2$

trong đó $I_{t-1}$ là biến chỉ báo

$I_{t-1} = \begin{cases} 1 & \text{if } \varepsilon_{t-1} < 0 \\ 0 & \text{if } \varepsilon_{t-1} \geq 0 \end{cases}$

Hãy đọc phép truy hồi theo từng trường hợp. Sau một cú sốc dương ($\varepsilon_{t-1} \geq 0$), biến chỉ báo bằng không và tác động của cú sốc bình phương lên phương sai kỳ tiếp theo chỉ là $\alpha$. Sau một cú sốc âm, biến chỉ báo bằng một và tác động là $\alpha + \gamma$. Tham số $\gamma$ chính là toàn bộ câu chuyện bất đối xứng gói gọn trong một con số:

• $\gamma > 0$: các cú sốc âm làm tăng biến động nhiều hơn các cú sốc dương cùng độ lớn. Đây là hiệu ứng đòn bẩy, và là điều bạn kỳ vọng tìm thấy ở BTC/ETH trong phần lớn thời gian.
• $\gamma = 0$: mô hình thu về GARCH(1,1) đối xứng. Do đó, một kiểm định tỷ số hợp lý (likelihood-ratio) hoặc kiểm định $t$ trên $\gamma$ là một kiểm định trực tiếp về việc tính bất đối xứng có tồn tại hay không.
• $\gamma < 0$: các cú sốc dương làm tăng biến động nhiều hơn — chế độ bùng nổ tăng giá thi thoảng gặp trong crypto. Hiếm, nhưng đừng loại trừ nó một cách tiên nghiệm.

Tính dương và tính dừng

Bởi vì $\sigma_t^2$ vẫn được xây dựng theo kiểu cộng, ta cần mỗi số hạng đều không âm. Các điều kiện đủ để đảm bảo tính dương là

$\omega > 0, \quad \alpha \geq 0, \quad \alpha + \gamma \geq 0, \quad \beta \geq 0.$

Lưu ý rằng bản thân $\gamma$ được phép âm miễn là $\alpha + \gamma \geq 0$, nên tác động sau tin xấu không bao giờ trở nên âm.

Để có tính dừng theo hiệp phương sai (covariance stationarity), giả định các cải tiến $z_t = \varepsilon_t/\sigma_t$ được chuẩn hóa với một phân phối đối xứng quanh không, sao cho $P(\varepsilon_{t-1} < 0) = 1/2$ và trung bình biến chỉ báo đóng góp $\gamma/2$. Điều kiện dừng trở thành

$\alpha + \beta + \tfrac{1}{2}\gamma < 1.$

Khi đó phương sai vô điều kiện (dài hạn) là

$\bar{\sigma}^2 = \frac{\omega}{1 - \alpha - \beta - \tfrac{1}{2}\gamma}.$

Đây là phiên bản GJR tương ứng của kết quả $\bar{\sigma}^2 = \omega/(1-\alpha-\beta)$ ở Phần 1, với số hạng $\tfrac{1}{2}\gamma$ bổ sung tính đến đóng góp trung bình của phần bán biên đòn bẩy. Nếu phân phối cải tiến của bạn bị lệch (skew-$t$ của Hansen, bên dưới), thì $1/2$ được thay bằng xác suất thực tế mà $z_t < 0$, nhưng $1/2$ là mốc tham chiếu chuẩn được dùng cho độ dai dẳng (persistence) được báo cáo.

EGARCH: Mô hình hóa log-phương sai, không ràng buộc tính dương

GJR giữ bạn trong một cái áo bó buộc về tính dương của phương sai: mọi tổ hợp tham số đều phải được kiểm tra theo các ràng buộc bất đẳng thức, điều này gây phiền toái trong quá trình tối ưu và tệ hơn trong quá trình tái ước lượng lăn (rolling re-estimation) khi một cửa sổ đôi lúc lang thang vào vùng bất khả thi. Exponential GARCH của Nelson (1991) né tránh điều này hoàn toàn bằng cách mô hình hóa logarit của phương sai có điều kiện. Bởi vì $\log \sigma_t^2$ có thể là bất kỳ số thực nào, nên $\sigma_t^2 = \exp(\log \sigma_t^2)$ tự động dương bất kể tham số là gì. Không có ràng buộc nào phải áp đặt.

Viết phép truy hồi theo cải tiến chuẩn hóa $z_{t-1} = \varepsilon_{t-1}/\sigma_{t-1}$:

$\log \sigma_t^2 = \omega + \beta \log \sigma_{t-1}^2 + \alpha\bigl(|z_{t-1}| - \mathbb{E}|z_{t-1}|\bigr) + \gamma\, z_{t-1}$

Hai số hạng mang cú sốc, và việc tách chúng ra chính là toàn bộ ý tưởng:

• Số hạng độ lớn $\alpha(|z_{t-1}| - \mathbb{E}|z_{t-1}|)$ phản ứng với kích thước của cú sốc, đã bỏ dấu. Việc trừ đi $\mathbb{E}|z_{t-1}|$ căn giữa nó sao cho một cú sốc có độ lớn trung bình không đóng góp gì. Với phân phối chuẩn tắc, $\mathbb{E}|z| = \sqrt{2/\pi} \approx 0.7979$; với Student-$t$ chuẩn hóa thì giá trị tuyệt đối kỳ vọng nhỏ hơn và phụ thuộc vào $\nu$, nhưng arch xử lý điều này nội bộ.
• Số hạng dấu $\gamma\, z_{t-1}$ là tính bất đối xứng. Nó tuyến tính theo cải tiến có dấu, nên một $z_{t-1}$ âm đẩy $\log\sigma_t^2$ theo hướng ngược lại với một $z_{t-1}$ dương.

Quy ước dấu rất quan trọng và hay làm người ta vấp. Trong tham số hóa này, hiệu ứng đòn bẩy (tin xấu làm tăng biến động) tương ứng với $\gamma < 0$: khi đó một cú sốc âm $z_{t-1} < 0$ làm cho $\gamma z_{t-1} > 0$, làm tăng log-phương sai. Đây là dấu ngược với $\gamma > 0$ của GJR. Hãy luôn đọc tài liệu của chính mô hình để biết quy ước thay vì giả định; arch báo cáo EGARCH với quy ước dấu của riêng nó, và chúng ta kiểm tra điều đó bằng đường cong tác động của tin tức (news impact curve) bên dưới thay vì tin vào trí nhớ.

Bởi vì mọi thứ đều cộng theo dạng log, độ dai dẳng của một EGARCH(1,1) được điều khiển bởi hệ số tự hồi quy duy nhất $\beta$ trên $\log\sigma_{t-1}^2$; tính dừng chỉ yêu cầu $|\beta| < 1$. Đó là một điều kiện gọn hơn nhiều so với bất đẳng thức của GJR, và là một lợi thế thực tiễn có thật khi bạn khớp lại trên các cửa sổ lăn.

Một điểm tinh tế đáng nêu: phản ứng của EGARCH với các cú sốc là hàm mũ theo cải tiến (bạn lấy hàm mũ ở bước cuối), trong khi GJR là hàm bậc hai. Do đó EGARCH phản ứng dữ dội hơn với các cú sốc lớn — một đặc điểm hữu ích trong crypto, nơi các sự kiện đuôi mới là điều quan trọng, nhưng cũng là lý do EGARCH đôi khi có thể tạo ra các dự báo phương sai lớn một cách không hợp lý sau một ngày ngoại lai. Không cái nào tốt hơn một cách phổ quát; bạn chọn dựa trên độ khớp ngoài mẫu (out-of-sample) và các backtest rủi ro, vốn là mục đích của cả loạt bài này.

Đường cong tác động của tin tức

Cách sạch sẽ nhất để thấy sự khác biệt giữa GARCH đối xứng, GJR và EGARCH là đường cong tác động của tin tức (Engle và Ng, 1993): giữ $\sigma_{t-1}$ cố định ở mức dài hạn của nó và vẽ phương sai có điều kiện kỳ tiếp theo $\sigma_t^2$ như một hàm của cú sốc gần nhất $\varepsilon_{t-1}$. Nó trả lời câu hỏi "với một cú sốc có kích thước và dấu như thế này, mô hình làm tăng biến động ngày mai lên bao nhiêu?"

• GARCH đối xứng tạo ra một parabol đối xứng căn giữa tại không. Một cú sốc $-5\%$ và một cú sốc $+5\%$ rơi vào cùng độ cao. Đây chính xác là khiếm khuyết mà chúng ta đang khắc phục.
• GJR tạo ra một parabol có nếp gấp tại không — dốc hơn ở bên trái (các cú sốc âm) so với bên phải khi $\gamma > 0$. Hai nửa có độ cong lần lượt là $\alpha+\gamma$ và $\alpha$.
• EGARCH tạo ra một hình chữ V hàm mũ, bất đối xứng: hai nhánh có độ dốc khác nhau vì số hạng $\gamma z$, và toàn bộ đường cong đi lên nhanh hơn parabol vì phép lấy hàm mũ cuối cùng.

Chúng ta sẽ vẽ cả ba từ các tham số đã khớp ở phần sau, trong phần triển khai — đó là chẩn đoán hữu ích nhất để truyền đạt xem tính bất đối xứng mang lại điều gì.

Đuôi dày: Các cải tiến Student-t và skew-t

Tính bất đối xứng khắc phục phản ứng của mô hình với dấu của các cú sốc. Nó không làm gì với phân phối của bản thân các cú sốc. GARCH thuần túy giả định $z_t \sim \mathcal{N}(0,1)$, và giả định đó gần như luôn sai đối với crypto. Ngay cả sau khi GARCH loại bỏ hiện tượng phân cụm biến động, phần dư chuẩn hóa $z_t = \varepsilon_t / \sigma_t$ vẫn giữ độ nhọn vượt trội (excess kurtosis) — chúng có đuôi dày. Một hàm hợp lý Gaussian, khi khớp phần vai của phân phối, đánh giá thấp tần suất thực sự xảy ra của một ngày chuẩn hóa $4$, $5$ hay $6$ sigma.

Hệ quả đối với rủi ro là trực tiếp. VaR 99% Gaussian dùng phân vị $\Phi^{-1}(0.01) \approx -2.326$, nên nó dự báo $\text{VaR}_{0.99} \approx 2.326\,\sigma_t$. Nếu phân phối chuẩn hóa thực sự là Student-$t$ với, chẳng hạn, $\nu = 5$ bậc tự do, thì phân vị 1% thực sự nằm gần $-3.36$ — VaR Gaussian lạc quan quá mức khoảng $44\%$ tại mức tin cậy đó. Bạn sẽ vi phạm nó nhiều hơn 1% số ngày và liên tục bị bất ngờ bởi những ngày "bất khả". Đây không phải là một điều kỳ quặc riêng của crypto; Bollerslev (1987) đã giới thiệu $t$-GARCH chính vì phần dư của cổ phiếu và ngoại hối cũng cho thấy cùng loại đuôi dày. Crypto chỉ là phiên bản cực đoan hơn của cùng vấn đề đó.

Student-t chuẩn hóa

Mật độ Student-$t$ có một tham số bậc tự do $\nu > 2$ điều khiển độ dày của đuôi: $\nu$ nhỏ nghĩa là đuôi dày, và khi $\nu \to \infty$ thì $t$ hội tụ về Gaussian. Điểm mắc là phân phối $t_\nu$ thô có phương sai $\nu/(\nu-2) \neq 1$, nên ta phải chuẩn hóa nó về phương sai đơn vị trước khi dùng làm cải tiến — nếu không thì "$\sigma_t$" trong phép truy hồi GARCH sẽ không thực sự là độ lệch chuẩn có điều kiện.

Cải tiến Student-$t$ chuẩn hóa với phương sai đơn vị có mật độ

$f(z;\nu) = \frac{\Gamma\!\left(\frac{\nu+1}{2}\right)}{\sqrt{\pi(\nu-2)}\;\Gamma\!\left(\frac{\nu}{2}\right)}\left(1 + \frac{z^2}{\nu-2}\right)^{-\frac{\nu+1}{2}}, \qquad \nu > 2.$

Chú ý $(\nu-2)$ ở bên trong — đó là phần chuẩn hóa, tái tỷ lệ sao cho $\mathrm{Var}(z)=1$. Đóng góp log-hợp lý của một quan sát, cho trước phương sai có điều kiện GARCH $\sigma_t^2$ và $z_t = \varepsilon_t/\sigma_t$, là

$\ell_t = \log\Gamma\!\left(\tfrac{\nu+1}{2}\right) - \log\Gamma\!\left(\tfrac{\nu}{2}\right) - \tfrac{1}{2}\log\bigl(\pi(\nu-2)\bigr) - \tfrac{1}{2}\log\sigma_t^2 - \frac{\nu+1}{2}\log\!\left(1 + \frac{z_t^2}{\nu-2}\right).$

Số hạng $-\tfrac{1}{2}\log\sigma_t^2$ là Jacobian của phép biến đổi từ $\varepsilon_t$ sang $z_t$ — chính số hạng bạn đã thấy trong hàm hợp lý GARCH Gaussian ở Phần 1. Chỉ có dạng (shape) thay đổi. Cực đại hóa $\sum_t \ell_t$ đồng thời theo các tham số GARCH và $\nu$ chính là điều mà arch làm khi bạn truyền dist='t'.

Bản thân $\nu$ ước lượng được cũng mang thông tin. Với lợi suất hằng ngày của BTC/ETH, bạn thường rơi vào khoảng $\nu \approx 3\text{–}6$ — đuôi dày, nhưng với phương sai hữu hạn (cần $\nu > 2$) và thường là độ nhọn hữu hạn (cần $\nu > 4$). Nếu $\nu$ khớp được của bạn tụt xuống dưới 4, hãy lưu ý rằng độ nhọn mẫu về mặt kỹ thuật là vô hạn trong mô hình và một số ước lượng trở nên bất ổn; đó là tín hiệu để xem xét kỹ các điểm ngoại lai và chất lượng dữ liệu.

Skew-t của Hansen

Student-$t$ có đuôi dày nhưng vẫn đối xứng — đuôi trái và đuôi phải nặng như nhau. Phần dư lợi suất crypto thường còn bị lệch nữa: đuôi trái (các đợt sụp đổ) nặng hơn đuôi phải. Skew-$t$ của Hansen (1994) tổng quát hóa $t$ chuẩn hóa với một tham số độ lệch (skewness) $\lambda \in (-1, 1)$ bên cạnh $\nu$:

$f(z;\nu,\lambda) = \begin{cases} b\,c\left(1 + \dfrac{1}{\nu-2}\left(\dfrac{bz+a}{1-\lambda}\right)^2\right)^{-\frac{\nu+1}{2}} & z < -a/b \\[2mm] b\,c\left(1 + \dfrac{1}{\nu-2}\left(\dfrac{bz+a}{1+\lambda}\right)^2\right)^{-\frac{\nu+1}{2}} & z \geq -a/b \end{cases}$

trong đó các hằng số $a = 4\lambda c\,\frac{\nu-2}{\nu-1}$, $b^2 = 1 + 3\lambda^2 - a^2$, và $c = \frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\pi(\nu-2)}\,\Gamma(\frac{\nu}{2})}$ được chọn sao cho $z$ có trung bình bằng không và phương sai đơn vị với mọi $(\nu,\lambda)$ hợp lệ. Phân phối tách tại $z = -a/b$, dùng một hệ số tỷ lệ khác nhau trong mỗi mảnh để dồn thêm khối lượng vào một đuôi.

Diễn giải: $\lambda < 0$ cho một phân phối lệch trái (rủi ro giảm nặng hơn), vốn là phát hiện thường gặp đối với crypto và là điều bạn kỳ vọng đi kèm với một hiệu ứng đòn bẩy. $\lambda = 0$ khôi phục Student-$t$ đối xứng, nên một kiểm định $\lambda = 0$ cho bạn biết số hạng lệch có mang lại gì hay không. Trong arch đây là dist='skewt', ước lượng cả $\nu$ và $\lambda$. Phần thưởng là một VaR có phân vị đuôi trái nặng hơn một cách trung thực so với phân vị đuôi phải — chính xác là điều bạn muốn khi các khoản lỗ mà bạn đang cố sống sót qua có tính bất đối xứng. Điều này liên hệ trực tiếp tới tính bất đối xứng giữa lỗ và lãi trong kết quả vị thế: một khoản sụt giảm $x\%$ cần hơn $x\%$ để hồi phục, nên việc mô hình hóa sai đuôi trái tốn kém hơn mô hình hóa sai đuôi phải.

Triển khai bằng Python

Bây giờ chúng ta khớp tất cả những điều này bằng thư viện arch. Thiết lập phản chiếu Phần 1: lấy lợi suất hằng ngày, tỷ lệ hóa nhân 100 để cải thiện điều kiện số (các bộ tối ưu GARCH hoạt động kém khi lợi suất ở cỡ $O(0.01)$), và khớp với một trung bình hằng số. Nếu bạn muốn dữ liệu trong ngày (intraday) hay một mô hình trung bình khác, bộ máy vẫn y hệt.

Thiết lập và dữ liệu

import numpy as np
import pandas as pd
from arch import arch_model
from arch.univariate import GARCH, EGARCH, ConstantMean, StudentsT, SkewStudent
from scipy import stats

def fetch_returns(symbol="BTC-USD", start="2019-01-01", end="2025-12-31"):
    """Daily log returns, in percent (scaled x100 for the optimizer)."""
    import yfinance as yf
    px = yf.download(symbol, start=start, end=end)["Close"].dropna()
    ret = 100.0 * np.log(px / px.shift(1)).dropna()
    ret.name = symbol
    return ret

r = fetch_returns("BTC-USD")
print(f"{len(r)} daily observations, "
      f"annualized vol ~ {r.std() * np.sqrt(365):.0f}%")

Crypto giao dịch 24/7, nên chúng ta quy về hằng năm với 365 chứ không phải 252 — một nguồn gây nhầm lẫn nhỏ nhưng lặp lại khi bạn so sánh Sharpe hoặc biến động của crypto với các con số của một bàn giao dịch cổ phiếu.

Khớp bốn mô hình

Mẫu hình trong arch: vol='Garch' với p=1, q=1 là GARCH đối xứng; thêm o=1 bật số hạng ngưỡng GJR; vol='EGARCH' chuyển sang mô hình log-phương sai. Phân phối cải tiến được đặt bằng dist: 'normal', 't', 'skewt'.

def fit(returns, vol="Garch", p=1, o=0, q=1, dist="normal"):
    am = arch_model(returns, mean="Constant",
                    vol=vol, p=p, o=o, q=q, dist=dist)
    res = am.fit(disp="off")
    return res

models = {
    "GARCH-N":    fit(r, vol="Garch",  o=0, dist="normal"),  # Part 1 baseline
    "GJR-t":      fit(r, vol="Garch",  o=1, dist="t"),        # asymmetry + fat tails
    "EGARCH-t":   fit(r, vol="EGARCH", o=1, dist="t"),        # log-variance asymmetry
    "GJR-skewt":  fit(r, vol="Garch",  o=1, dist="skewt"),    # + skew
}

for name, res in models.items():
    print(f"\n===== {name} =====")
    print(res.summary().tables[1])   # parameter table

Với vol='EGARCH', đối số o điều khiển số hạng bất đối xứng ($\gamma z$) còn p/q điều khiển các số hạng độ lớn và độ trễ; o=1, p=1, q=1 là EGARCH(1,1) chuẩn. Một điểm dễ vấp: tên tham số của EGARCH trong arch cũng là các chữ cái đó nhưng quy ước dấu trên số hạng bất đối xứng là của Nelson, nên một ước lượng âm chính là hiệu ứng đòn bẩy. Chúng ta xác minh điều này từ đường cong tác động của tin tức chứ không từ trí nhớ.

Đọc kết quả khớp GJR

Một bảng tham số GJR-$t$ trông đại khái như thế này (giá trị minh họa, không phải một thí nghiệm được báo cáo — hãy khớp lại trên dữ liệu của chính bạn):

                  coef    std err       t      P>|t|
omega           0.0480     0.017     2.82     0.005
alpha[1]        0.0620     0.018     3.44     0.001
gamma[1]        0.0910     0.028     3.25     0.001
beta[1]         0.8850     0.021    42.1      0.000
nu              4.30       0.55      7.82     0.000

Cách đọc nó:

• gamma[1] = 0.091 với thống kê $t$ trên 3 là một hiệu ứng đòn bẩy có ý nghĩa thống kê. Sau một cú sốc âm, tác động của cú sốc bình phương là $\alpha + \gamma = 0.062 + 0.091 = 0.153$; sau một cú sốc dương thì chỉ là $\alpha = 0.062$. Tin xấu làm biến động của mô hình này thay đổi khoảng gấp $2.5\times$ so với tin tốt cùng độ lớn.
• nu = 4.3 xác nhận đuôi dày — cách xa Gaussian ($\nu \to \infty$), và đủ thấp để mô men bậc bốn chỉ vừa vặn hữu hạn. Một VaR Gaussian trên chuỗi này sẽ lạc quan quá mức một cách tệ hại.
• Độ dai dẳng là $\alpha + \beta + \tfrac{1}{2}\gamma = 0.062 + 0.885 + 0.0455 \approx 0.993$ — rất cao, như thường thấy đối với crypto hằng ngày: các cú sốc phân rã chậm và biến động phân cụm mạnh.

Dòng quan trọng nhất cần kiểm tra là dòng $\gamma$. Nếu $p$-value của nó lớn, số hạng bất đối xứng không xứng đáng với vị trí của mình trên tài sản và cửa sổ này, và bạn nên ưu tiên mô hình đối xứng đơn giản hơn. Đây là kỷ luật lựa chọn mô hình, không phải vật trang trí — sẽ nói thêm bên dưới.

So sánh các mô hình bằng tiêu chuẩn thông tin

Log-hợp lý luôn cải thiện khi bạn thêm tham số, nên bạn không thể chọn dựa trên log-hợp lý đơn thuần. Hãy dùng AIC/BIC, vốn phạt số lượng tham số (BIC phạt mạnh hơn):

def compare(models: dict) -> pd.DataFrame:
    rows = []
    for name, res in models.items():
        rows.append({
            "model": name,
            "n_params": len(res.params),
            "loglik": res.loglikelihood,
            "AIC": res.aic,
            "BIC": res.bic,
        })
    df = pd.DataFrame(rows).set_index("model")
    return df.sort_values("BIC")

print(compare(models))

Quy tắc kinh nghiệm khi diễn giải: một cải thiện BIC lớn hơn ~6 so với mô hình nền là bằng chứng mạnh rằng cấu trúc bổ sung là thực; một chênh lệch 1–2 là nhiễu. Nếu GJR-t thắng GARCH-N 30+ điểm BIC nhưng GJR-skewt chỉ thắng GJR-t 1 điểm, hãy giữ $t$ và bỏ phần lệch — tham số lệch không tự trang trải được trên dữ liệu này. Đừng đọc AIC/BIC như một sự thay thế cho kiểm định ngoài mẫu; chúng tưởng thưởng độ khớp trong mẫu đã điều chỉnh theo độ phức tạp, điều cần thiết nhưng chưa đủ. Bài kiểm định thực sự là backtest VaR và, rốt cuộc, đánh giá walk-forward.

Vẽ đường cong tác động của tin tức

Đây là biểu đồ phần thưởng — nó làm tính bất đối xứng trở nên hữu hình và xác minh quy ước dấu của EGARCH.

import matplotlib.pyplot as plt

def news_impact_curve(res, shock_grid):
    """
    Next-period conditional variance as a function of the last shock,
    holding sigma_{t-1} at the model's unconditional level.
    Works for symmetric GARCH, GJR (o=1), and EGARCH.
    """
    p = res.params
    vol_name = res.model.volatility.__class__.__name__
    sigma2_bar = np.mean(res.conditional_volatility ** 2)

    omega = p["omega"]
    alpha = p.get("alpha[1]", 0.0)
    gamma = p.get("gamma[1]", 0.0)
    beta  = p.get("beta[1]", 0.0)

    if vol_name == "EGARCH":
        sig_prev = np.sqrt(sigma2_bar)
        z = shock_grid / sig_prev
        E_abs = np.sqrt(2 / np.pi)   # approx; arch uses the fitted dist's E|z|
        log_s2 = (omega + beta * np.log(sigma2_bar)
                  + alpha * (np.abs(z) - E_abs) + gamma * z)
        return np.exp(log_s2)
    else:
        ind = (shock_grid < 0).astype(float)
        return omega + (alpha + gamma * ind) * shock_grid**2 + beta * sigma2_bar

shocks = np.linspace(-8, 8, 401)   # daily % shocks
plt.figure(figsize=(8, 5))
for name in ["GARCH-N", "GJR-t", "EGARCH-t"]:
    nic = news_impact_curve(models[name], shocks)
    plt.plot(shocks, nic, label=name, lw=2)
plt.axvline(0, color="gray", lw=0.8, ls="--")
plt.xlabel("Last shock  $\\varepsilon_{t-1}$  (daily %)")
plt.ylabel("Next-period conditional variance  $\\sigma_t^2$")
plt.title("News impact curve: symmetric vs GJR vs EGARCH (BTC)")
plt.legend()
plt.tight_layout()

Khi bạn chạy đoạn này, đường cong GARCH-N đối xứng là một parabol sạch sẽ căn giữa tại không — một cú sốc $-6\%$ và $+6\%$ cho phương sai giống hệt nhau. GJR-t là một parabol có nếp gấp tại gốc, cao hơn ở nhánh trái. EGARCH-t là hình chữ V hàm mũ, và nếu nhánh trái của nó nằm trên nhánh phải thì bạn đã xác nhận hiệu ứng đòn bẩy và quy ước dấu chỉ trong một cái liếc. Nếu nhánh trái của EGARCH nằm dưới nhánh phải, thì hoặc $\gamma$ ước lượng dương (một chế độ tăng biến động khi giá lên) hoặc bạn đã hiểu ngược dấu — biểu đồ cho bạn biết là cái nào mà không cần đoán.

So sánh trực tiếp bốn mô hình

Trước khi chuyển sang rủi ro, sẽ có ích khi đặt bốn mô hình cạnh nhau. Mỗi hàng là một quyết định thiết kế, và các cột cho thấy quyết định đó tốn gì và mua được gì.

Mẫu hình cần thấm nhuần: bước nhảy từ cột 1 sang cột 2 — thêm đồng thời cả tính bất đối xứng lẫn đuôi dày — là nơi chứa gần như toàn bộ phần cải thiện trong hiệu chỉnh rủi ro. Các tinh chỉnh tiếp theo (dạng hàm của EGARCH, số hạng lệch) là có thật nhưng bậc hai, và trên nhiều chuỗi crypto chúng nằm trong nhiễu. Hãy dồn ngân sách mô hình hóa của bạn vào bước nhảy đầu tiên và hoài nghi phần còn lại.

Ứng dụng rủi ro: VaR và Expected Shortfall

Khớp một mô hình biến động phức tạp hơn chỉ đáng công nếu nó cải thiện một quyết định. Quyết định sạch sẽ nhất để cải thiện là dự báo rủi ro đuôi một bước: ngày mai có thể tệ đến đâu? Chúng ta tạo ra một dự báo Value-at-Risk và Expected Shortfall một ngày tới (còn gọi là Conditional VaR, thứ mà quy trình danh mục HRP/CVaR dùng làm hàm mục tiêu) trực tiếp từ dự báo GARCH-$t$/skew-$t$ đã khớp.

Từ phân phối có điều kiện tới VaR

Bộ máy GARCH cho một dự báo một bước của trung bình có điều kiện $\mu_{t+1}$ và độ lệch chuẩn có điều kiện $\sigma_{t+1}$. Lợi suất được mô hình hóa là $r_{t+1} = \mu_{t+1} + \sigma_{t+1} z_{t+1}$ với $z_{t+1}$ được lấy từ phân phối chuẩn hóa đã khớp (Gaussian, $t$, hoặc skew-$t$). Vậy nên phân vị $\alpha$ của lợi suất chỉ là một phép biến đổi affine của phân vị $\alpha$ của phân phối chuẩn hóa:

$\text{VaR}_{\alpha}(t+1) = -\Bigl(\mu_{t+1} + \sigma_{t+1}\, F_z^{-1}(1-\alpha)\Bigr)$

trong đó $F_z^{-1}$ là phân vị (CDF nghịch đảo) của cải tiến chuẩn hóa và dấu trừ đứng đầu tuân theo quy ước rằng VaR là một con số lỗ dương. Với VaR 99%, $\alpha = 0.99$ và bạn thế vào $F_z^{-1}(0.01)$. Toàn bộ lợi ích của $t$/skew-$t$ hiện ra ở đây: $F_z^{-1}(0.01)$ âm hơn so với $-2.326$ của Gaussian, nên VaR lớn hơn một cách trung thực.

Expected Shortfall

VaR cho bạn biết ngưỡng; nó không nói gì về việc vi phạm tệ đến đâu khi nó xảy ra. Expected Shortfall — khoản lỗ trung bình có điều kiện vượt ngưỡng VaR — thì có, và nó là một độ đo nhất quán (coherent, có tính cộng dưới - subadditive), đó là lý do nó là độ đo rủi ro đứng sau tối ưu hóa CVaR và lý do Basel chuyển sang dùng nó. Với một mô hình vị trí-tỷ lệ (location-scale),

$\text{ES}_{\alpha}(t+1) = -\Bigl(\mu_{t+1} + \sigma_{t+1}\, \mathbb{E}\!\left[z \mid z \leq F_z^{-1}(1-\alpha)\right]\Bigr).$

Số hạng kỳ vọng đuôi có điều kiện $\mathbb{E}[z \mid z \le q]$ có dạng đóng cho các phân phối chuẩn. Với Gaussian, khi $q = \Phi^{-1}(1-\alpha)$,

$\mathbb{E}[z \mid z \le q] = -\frac{\phi(q)}{1-\alpha},$

trong đó $\phi$ là pdf chuẩn tắc. Với Student-$t$ chuẩn hóa có $\nu$ bậc tự do và $q = t_\nu^{-1}(1-\alpha)$ (trên thang chuẩn hóa), kỳ vọng đuôi là

$\mathbb{E}[z \mid z \le q] = -\frac{\nu + q^2}{\nu - 1}\cdot\frac{g_\nu(q)}{1-\alpha},$

trong đó $g_\nu$ là pdf của $t$ chuẩn hóa. Expected Shortfall của $t$ vượt của Gaussian nhiều hơn so với mức VaR vượt, bởi vì đuôi $t$ không chỉ nằm xa hơn — nó còn dày hơn, nên khoản lỗ trung bình vượt ngưỡng lớn một cách không tương xứng. Khoảng chênh lệch thêm đó là con số mà một mô hình Gaussian giấu bạn.

Tính VaR và ES từ một mô hình arch đã khớp

Các phân phối của arch cung cấp một phương thức ppf (phân vị), nên chúng ta có thể lấy phân vị chuẩn hóa trực tiếp và tránh phải suy dẫn lại bất cứ điều gì. Với ES, chúng ta tích phân bằng số, cách này vững chắc và hoạt động đồng nhất trên normal/t/skewt.

from scipy import integrate

def var_es_forecast(res, alpha=0.99):
    """
    One-step-ahead VaR and ES at level alpha, on the same (x100) scale
    as the returns fed to the model. Divide by 100 for fractional units.
    """
    fc = res.forecast(horizon=1, reindex=False)
    mu = fc.mean.iloc[-1, 0]
    sigma = np.sqrt(fc.variance.iloc[-1, 0])

    dist = res.model.distribution          # StudentsT / SkewStudent / Normal
    dp = [res.params[k] for k in res.params.index
          if k in ("nu", "eta", "lambda", "lam")]

    z_q = dist.ppf(1 - alpha, dp) if dp else dist.ppf(1 - alpha)
    z_q = float(np.atleast_1d(z_q)[0])

    def pdf(z):
        arr = np.atleast_1d(z).astype(float)
        lp = dist.loglikelihood(dp, arr, np.ones_like(arr), individual=True) \
             if dp else dist.loglikelihood([], arr, np.ones_like(arr),
                                           individual=True)
        return np.exp(lp)

    num, _ = integrate.quad(lambda z: z * pdf(z), -30, z_q, limit=200)
    es_z = num / (1 - alpha)               # E[z | z <= z_q]

    var = -(mu + sigma * z_q)
    es  = -(mu + sigma * es_z)
    return {"mu": mu, "sigma": sigma, "z_q": z_q,
            "VaR": var, "ES": es}

for name in ["GARCH-N", "GJR-t", "GJR-skewt"]:
    out = var_es_forecast(models[name], alpha=0.99)
    print(f"{name:11s}  sigma={out['sigma']:.2f}%  "
          f"z_q={out['z_q']:+.2f}  "
          f"VaR99={out['VaR']:.2f}%  ES99={out['ES']:.2f}%")

Cột z_q là toàn bộ câu chuyện gói trong một con số. Mô hình Gaussian dùng $z_q \approx -2.33$; $t$ với $\nu \approx 4.3$ dùng giá trị gần $-3.3$; skew-$t$ đẩy phân vị trái ra xa hơn nữa trong khi kéo phân vị phải vào trong. Cùng một $\sigma_{t+1}$, nhưng VaR lớn hơn đáng kể. Nếu bạn đang chạy VaR Gaussian trên crypto, đây chính là khoảng chênh lệch mà bạn đang âm thầm hứng chịu.

Một bước so với nhiều bước: Một lưu ý

Mọi thứ ở trên là một dự báo một ngày tới, và đó là nơi VaR GARCH sạch sẽ nhất. Hai điều làm phức tạp các đường chân trời dài hơn và bạn nên biết chúng trước khi ngoại suy.

Thứ nhất, các dự báo phương sai hồi quy về trung bình (mean-revert). Phương sai có điều kiện $h$ bước tới từ một GARCH dừng hội tụ về mức vô điều kiện $\bar\sigma^2$ khi $h$ tăng, và phương sai lũy kế $h$ ngày là tổng của các dự báo từng bước — nó không phải $h\times\sigma_{t+1}^2$ trừ khi biến động đang ở mức trung bình dài hạn. Phép tỷ lệ hóa "căn bậc hai của thời gian" ngây thơ $\text{VaR}(h) = \sqrt{h}\,\text{VaR}(1)$ bỏ qua sự hồi quy về trung bình này và sai chính xác ngay sau một cú sốc, khi bạn cần con số đó nhất. Hãy dùng đường phương sai nhiều bước của chính mô hình.

Thứ hai, phân phối của một lợi suất nhiều ngày không cùng dạng với cải tiến một ngày. Cộng nhiều cú sốc hằng ngày phân phối $t$ (thông qua phép truy hồi GARCH phi tuyến) không cho một phân phối $t$ tại đường chân trời $h$ ngày; không có dạng đóng sạch sẽ nào. Với VaR nhiều ngày, con đường trung thực là mô phỏng: lấy các đường cải tiến từ phân phối chuẩn hóa đã khớp, chạy chúng qua phép truy hồi GARCH để có các đường lợi suất mô phỏng, tổng hợp thành lợi suất $h$ ngày, và đọc phân vị thực nghiệm. Cách đó cũng xử lý một cách tự nhiên trường hợp skew-$t$, nơi hoàn toàn không tồn tại phân vị nhiều đường chân trời dạng giải tích. Các công thức giải tích một bước trong bài này là chính xác; hãy xem mọi lối tắt nhiều bước như một xấp xỉ cần kiểm chứng.

Backtest VaR: Kupiec và Christoffersen

Một dự báo VaR là một khẳng định xác suất: "khoản lỗ sẽ vượt ngưỡng này chỉ trong $(1-\alpha)$ số ngày." Bạn kiểm định nó bằng cách đếm số vi phạm (những ngày khoản lỗ thực hiện vượt VaR dự báo) qua một đánh giá walk-forward và kiểm tra hai điều. Thứ nhất, tỷ lệ vi phạm có đúng không? Thứ hai, các vi phạm có độc lập hay chúng phân cụm (nghĩa là mô hình thất bại đúng lúc quan trọng, trong các đợt biến động vọt lên)?

Gọi $I_t = \mathbf{1}\{\text{loss}_t > \text{VaR}_t\}$ là chuỗi vi phạm, $N = \sum I_t$ là số vi phạm trên $T$ ngày, và $\hat{\pi} = N/T$ là tỷ lệ quan sát được. Tỷ lệ mục tiêu $p = 1-\alpha$.

Kiểm định độ phủ vô điều kiện của Kupiec (1995) kiểm tra $\hat\pi \approx p$ thông qua một tỷ số hợp lý:

$LR_{uc} = -2\log\!\left[\frac{p^{N}(1-p)^{T-N}}{\hat\pi^{N}(1-\hat\pi)^{T-N}}\right] \;\sim\; \chi^2_1.$

Kiểm định độc lập của Christoffersen (1998) kiểm tra rằng một vi phạm hôm nay không được dự báo bởi một vi phạm hôm qua. Gọi $n_{ij}$ đếm các chuyển tiếp từ trạng thái $i$ sang trạng thái $j$ trong chuỗi vi phạm, $\pi_{01} = n_{01}/(n_{00}+n_{01})$, $\pi_{11} = n_{11}/(n_{10}+n_{11})$, và $\pi = (n_{01}+n_{11})/T$. Khi đó

$LR_{ind} = -2\log\!\left[\frac{(1-\pi)^{n_{00}+n_{10}}\,\pi^{\,n_{01}+n_{11}}}{(1-\pi_{01})^{n_{00}}\pi_{01}^{n_{01}}(1-\pi_{11})^{n_{10}}\pi_{11}^{n_{11}}}\right] \;\sim\; \chi^2_1.$

Hai kiểm định này kết hợp thành kiểm định độ phủ có điều kiện $LR_{cc} = LR_{uc} + LR_{ind} \sim \chi^2_2$, đồng thời kiểm tra tỷ lệ đúng và tính độc lập. Một mô hình có thể vượt qua Kupiec (số vi phạm đúng) nhưng lại trượt Christoffersen (chúng dồn hết vào một tuần sụp đổ) — đó là kiểu thất bại mà bạn muốn bắt được nhất, bởi vì các vi phạm phân cụm là những cái làm nổ tung một tài khoản.

from scipy.stats import chi2

def var_backtest(losses, var, p):
    """
    losses, var: 1D arrays, same length; loss>0 is a loss, var>0 threshold.
    p = 1 - alpha (target violation rate, e.g. 0.01 for 99% VaR).
    Returns Kupiec, Christoffersen-independence, and conditional-coverage LR.
    """
    losses = np.asarray(losses)
    var = np.asarray(var)
    I = (losses > var).astype(int)          # violation indicators
    T = len(I)
    N = int(I.sum())
    pi_hat = N / T

    eps = 1e-12
    ll_null = N * np.log(p + eps) + (T - N) * np.log(1 - p + eps)
    ll_alt  = N * np.log(pi_hat + eps) + (T - N) * np.log(1 - pi_hat + eps)
    LR_uc = -2 * (ll_null - ll_alt)
    p_uc = 1 - chi2.cdf(LR_uc, df=1)

    n00 = np.sum((I[:-1] == 0) & (I[1:] == 0))
    n01 = np.sum((I[:-1] == 0) & (I[1:] == 1))
    n10 = np.sum((I[:-1] == 1) & (I[1:] == 0))
    n11 = np.sum((I[:-1] == 1) & (I[1:] == 1))
    pi01 = n01 / (n00 + n01) if (n00 + n01) else 0.0
    pi11 = n11 / (n10 + n11) if (n10 + n11) else 0.0
    pi   = (n01 + n11) / (n00 + n01 + n10 + n11)

    def _ll(pi01, pi11, pi_pooled, use_pooled):
        if use_pooled:
            a = (n00 + n10) * np.log(1 - pi_pooled + eps)
            b = (n01 + n11) * np.log(pi_pooled + eps)
            return a + b
        return (n00 * np.log(1 - pi01 + eps) + n01 * np.log(pi01 + eps)
                + n10 * np.log(1 - pi11 + eps) + n11 * np.log(pi11 + eps))

    LR_ind = -2 * (_ll(pi01, pi11, pi, True) - _ll(pi01, pi11, pi, False))
    p_ind = 1 - chi2.cdf(LR_ind, df=1)

    LR_cc = LR_uc + LR_ind
    p_cc = 1 - chi2.cdf(LR_cc, df=2)

    return {
        "T": T, "violations": N,
        "obs_rate": pi_hat, "target_rate": p,
        "LR_uc": LR_uc, "p_uc": p_uc,
        "LR_ind": LR_ind, "p_ind": p_ind,
        "LR_cc": LR_cc, "p_cc": p_cc,
    }

Để tạo các đầu vào losses/var một cách trung thực, bạn khớp lại (hoặc ít nhất dự báo lại) trên một cửa sổ mở rộng hoặc lăn và ghi lại VaR một bước tới cho mỗi ngày ngoài mẫu, rồi so sánh nó với khoản lỗ thực hiện của ngày đó. Đừng bao giờ backtest VaR trong mẫu — một mô hình khớp trên chính đợt sụp đổ mà nó được yêu cầu dự báo sẽ trông tốt hơn thực tế rất nhiều. Đây là cùng kỷ luật với tính nhất quán giữa backtest và giao dịch thật: việc đánh giá chỉ được dùng thông tin có sẵn tại thời điểm ra quyết định.

def walk_forward_var(returns, alpha=0.99, start=750, refit_every=25,
                     vol="Garch", o=1, dist="t"):
    """
    Expanding-window one-step VaR. Refit every `refit_every` days
    (refitting daily is correct but slow; every ~25 days is a common
    compromise -- validate the shortcut on your data).
    """
    losses, vars_ = [], []
    res = None
    for t in range(start, len(returns)):
        if res is None or (t - start) % refit_every == 0:
            am = arch_model(returns.iloc[:t], mean="Constant",
                            vol=vol, p=1, o=o, q=1, dist=dist)
            res = am.fit(disp="off")
        fc = res.forecast(horizon=1, reindex=False)
        mu = fc.mean.iloc[-1, 0]
        sig = np.sqrt(fc.variance.iloc[-1, 0])
        dp = [res.params[k] for k in res.params.index
              if k in ("nu", "eta", "lambda", "lam")]
        z_q = float(np.atleast_1d(
            res.model.distribution.ppf(1 - alpha, dp) if dp
            else res.model.distribution.ppf(1 - alpha))[0])
        var_t = -(mu + sig * z_q)
        vars_.append(var_t)
        losses.append(-returns.iloc[t])     # realized loss for that day
    return np.array(losses), np.array(vars_)

losses, vars_ = walk_forward_var(r, alpha=0.99, dist="t")
print(var_backtest(losses, vars_, p=0.01))

Cách đọc: một VaR 99% được hiệu chỉnh tốt cho thấy tỷ lệ quan sát gần 1%, một Kupiec không có ý nghĩa (p_uc lớn), và một Christoffersen không có ý nghĩa (p_ind lớn) — không phân cụm. Trên thực tế, kết quả trung thực đối với crypto là GARCH-Normal trượt Kupiec (quá nhiều vi phạm, p_uc cực nhỏ) trong khi GJR-$t$ hoặc EGARCH-$t$ vượt qua hoặc suýt vượt. Sự tương phản đó là toàn bộ lập luận của bài viết này được diễn đạt dưới dạng một kiểm định giả thuyết. Nếu ngay cả mô hình $t$ cũng cho thấy các vi phạm phân cụm (p_ind nhỏ), thì động lực biến động của bạn vẫn bị đặc tả sai — thường là dấu hiệu bạn cần một bộ nhớ dài hơn (component/FIGARCH) hoặc một lớp chế độ (regime), điều này liên hệ tới phát hiện chế độ bằng HMM.

Xếp hạng mô hình theo tổn thất đuôi, không chỉ đạt/trượt

Kupiec và Christoffersen cho bạn một phán quyết nhị phân — mô hình có bị bác bỏ hay không. Điều đó cần thiết nhưng thô: hai mô hình đều có thể "đạt" trong khi một mô hình sắc nét hơn một cách đáng kể. Để xếp hạng các dự báo VaR cạnh tranh nhau, hãy chấm điểm chúng bằng một hàm tổn thất nhất quán chặt cho phân vị, tổn thất pinball (quantile loss):

$L_\tau(r, q) = \bigl(\tau - \mathbf{1}\{r < q\}\bigr)\,(r - q), \qquad \tau = 1-\alpha,$

trong đó $q$ là phân vị VaR (có dấu) và $r$ là lợi suất thực hiện. Lấy trung bình trên các ngày ngoài mẫu, một tổn thất pinball trung bình thấp hơn nghĩa là một phân vị được hiệu chỉnh tốt hơn và sắc nét hơn; bởi vì hàm tổn thất là nhất quán cho phân vị $\tau$, việc cực tiểu hóa nó không tưởng thưởng cho một mô hình chỉ vì nó rộng một cách lười biếng. Để so sánh hai mô hình một cách hình thức, hãy đưa các chênh lệch tổn thất theo từng ngày của chúng vào một kiểm định Diebold-Mariano.

def pinball_loss(returns, var, alpha=0.99):
    tau = 1 - alpha
    q = -np.asarray(var)               # VaR is a positive loss; quantile is negative
    r = np.asarray(returns)
    hit = (r < q).astype(float)
    return np.mean((tau - hit) * (r - q))

Riêng với Expected Shortfall, lưu ý rằng ES không thể tự nó khơi gợi được (not elicitable) (không tồn tại hàm tổn thất nào mà cực tiểu của nó chỉ là ES), đó là một nếp gấp lý thuyết thực sự: bạn đánh giá ES đồng thời với VaR bằng các quy tắc chấm điểm Fissler-Ziggel, hoặc bạn quay về thực hành đơn giản hơn là kiểm tra rằng độ lớn vi phạm trung bình khớp với ES mà mô hình dự báo. Một kiểm tra ES thô nhưng hữu ích: trong số các ngày vi phạm VaR, so sánh khoản lỗ thực hiện trung bình với ES dự báo trung bình vào những ngày đó — chúng nên gần nhau.

Khung quy định là cách tiếp cận đèn giao thông của Basel (Basel traffic-light): qua 250 ngày giao dịch, 0-4 vi phạm của một VaR 99% là "xanh" (chấp nhận được), 5-9 là "vàng" (cần soi xét), 10+ là "đỏ" (mô hình bị bác bỏ và các hệ số nhân vốn tăng lên). Đó là một người anh em họ thô hơn của Kupiec, nhưng nó là ngôn ngữ mà các ủy ban rủi ro thực sự nói, và đáng để báo cáo bên cạnh các thống kê LR.

Những cân nhắc thực tiễn

Khi các tham số bổ sung không đáng công

Mặc định trung thực là hoài nghi trước độ phức tạp. Mỗi tham số bạn thêm là một núm mà bộ tối ưu có thể overfit, và GARCH bất đối xứng đuôi dày có nhiều núm. Hướng dẫn cụ thể:

• Mẫu kém thanh khoản hoặc ngắn. Với vài trăm quan sát hằng ngày, sai số chuẩn trên $\gamma$ và $\lambda$ sẽ lớn, và bạn sẽ "phát hiện" các tính bất đối xứng vốn chỉ là nhiễu mẫu. Trên một altcoin mới hoặc mỏng, một GARCH-$t$ đối xứng thường là mô hình phức tạp nhất mà dữ liệu có thể chống đỡ. Khớp skew-$t$ EGARCH vào 200 ngày là tự lừa dối mình.
• Số hạng lệch thường không trang trải được chi phí của nó. Trên thực tế, chuyển Normal → $t$ là một cải thiện lớn, đáng tin cậy (đuôi dày là có thật và mạnh). Chuyển $t$ → skew-$t$ thường không đáng kể — một mức cải thiện BIC 1 hoặc 2, đôi khi âm. Chỉ thêm phần lệch khi dữ liệu rõ ràng đòi hỏi.
• EGARCH so với GJR thường ngang nhau trên dữ liệu hằng ngày. Chúng mã hóa cùng một câu chuyện định tính với các dạng hàm khác nhau. Hãy chọn theo backtest VaR ngoài mẫu, không theo cái nào có log-hợp lý đẹp hơn trong mẫu.
• Tần suất cao hơn làm thay đổi câu trả lời. Trên các thanh (bar) theo giờ hoặc theo phút, tính mùa vụ trong ngày và vi cấu trúc (microstructure) chiếm ưu thế, và một GARCH kiểu hằng ngày thuần túy bị đặc tả sai bất kể tính bất đối xứng. Bài toán khác, công cụ khác.

Đây là cùng bài học với đánh giá trung thực khi không có lợi thế vững chắc: một mô hình phức tạp hơn mà không sống sót qua kiểm định ngoài mẫu thì tệ hơn mô hình đơn giản mà nó thay thế, bởi vì nó mang theo ảo tưởng về độ chính xác. Hãy báo cáo kết quả âm — "phần lệch không giúp ích trên ETH" — như một phát hiện có thật, và dùng tối ưu hóa walk-forward làm trọng tài, chứ không phải AIC trong mẫu.

Đây là các phân phối biên mà mọi người khác xây dựng lên trên

Các mô hình ở đây không phải điểm cuối; chúng là khối xây dựng đơn biến (univariate) cho bộ máy rủi ro đồng thời. Bài viết các mô hình copula cho rủi ro crypto đồng thời dùng chính xác EGARCH/GJR-$t$ làm các phân phối biên GARCH-EVT trước khi khớp một vine copula — bạn khớp một GARCH bất đối xứng đuôi dày cho mỗi tài sản, trích xuất phần dư chuẩn hóa, và chỉ khi đó mới mô hình hóa sự phụ thuộc chéo giữa các tài sản. Nếu phân phối biên của bạn là một GARCH Gaussian đối xứng, thì copula thừa hưởng các sai số đuôi của nó dù mô hình phụ thuộc có tốt đến đâu. Phân phối biên rác, VaR đồng thời rác.

Với bài toán biến động đa biến — các tương quan thay đổi theo thời gian chứ không phải các phương sai theo từng tài sản — hãy xem Phần 3, DCC-GARCH, bài này xếp một mô hình tương quan động lên trên các khớp đơn biến này. Và để biến một dự báo biến động thành việc định cỡ vị thế và một backtest giao dịch, Phần 4 về nhắm mục tiêu biến động (volatility targeting) dùng các dự báo $\sigma_{t+1}$ từ chính các mô hình này để tỷ lệ hóa mức phơi nhiễm nghịch với rủi ro dự báo.

Một lựa chọn thay thế không phụ thuộc phân phối

Mọi thứ trong phần rủi ro đặt trên một giả định tham số: rằng phần dư chuẩn hóa tuân theo một $t$ hoặc skew-$t$. Giả định đó có thể kiểm định được và thường hợp lý, nhưng nó có thể thất bại. Nếu bạn muốn không cam kết vào một dạng đuôi nào cả, dự đoán conformal (conformal prediction) cho các khoảng dự đoán không phụ thuộc phân phối với các đảm bảo độ phủ mẫu hữu hạn — một triết lý thực sự khác biệt, không đưa ra khẳng định nào về phân phối cải tiến. Hai cách tiếp cận bổ trợ cho nhau: GARCH-$t$ tham số cho bạn một mật độ có điều kiện đầy đủ (và do đó ES, thứ mà các khoảng conformal không cung cấp trực tiếp), trong khi conformal cho bạn độ phủ vẫn đúng ngay cả khi mật độ của bạn sai. Trong môi trường thực địa (production), dùng cả hai làm kiểm tra chéo là một bảo hiểm rẻ tiền.

Vệ sinh số học và quy trình làm việc

• Tỷ lệ hóa lợi suất nhân 100. Các bộ tối ưu GARCH hội tụ đáng tin cậy hơn nhiều trên lợi suất phần trăm so với trên lợi suất phân số thô. Hãy nhớ bỏ tỷ lệ hóa VaR/ES nếu bạn báo cáo theo đơn vị phân số.
• Để mắt tới độ dai dẳng. Nếu $\alpha + \beta + \tfrac12\gamma$ ước lượng trên ~0.999, mô hình gần tích hợp (near-integrated, giống IGARCH); các dự báo hồi quy về trung bình cực kỳ chậm và các dự báo phương sai đường chân trời dài trở nên không đáng tin. Không nhất thiết là sai, nhưng hãy đánh dấu nó.
• Thất bại hội tụ trên các cửa sổ lăn. Dạng log của EGARCH tránh được các ràng buộc tính dương nhưng vẫn có thể không hội tụ trên một cửa sổ bệnh lý. Hãy bọc fit() trong một try/except và quay về các tham số của cửa sổ trước đó thay vì làm sập một backtest chạy trực tiếp.
• Mô hình trung bình. Chúng ta đã dùng một trung bình hằng số xuyên suốt. Với hầu hết crypto hằng ngày, trung bình có điều kiện gần bằng không và bị biến động lấn át; đừng tiêu độ phức tạp mô hình để cố dự báo nó trừ khi bạn có lý do thực sự.

Tóm tắt

• GARCH(1,1) thuần túy có hai khiếm khuyết cấu trúc: nó đối xứng (phản ứng với $+x\%$ và $-x\%$ giống hệt nhau vì các cú sốc đi vào dưới dạng $\varepsilon^2$) và nó giả định các cải tiến Gaussian (định giá thấp đuôi dày của crypto). Cả hai đều tốn tiền thật thông qua VaR lạc quan.
• GJR-GARCH thêm một số hạng ngưỡng $\gamma\,I(\varepsilon_{t-1}<0)\,\varepsilon_{t-1}^2$. Một $\gamma > 0$ có ý nghĩa chính là hiệu ứng đòn bẩy: tin xấu làm tăng biến động nhiều hơn. Tính dương cần $\alpha+\gamma\ge0$; độ dai dẳng là $\alpha+\beta+\tfrac12\gamma$.
• EGARCH mô hình hóa $\log\sigma_t^2$, nên không có ràng buộc tính dương và tính dừng chỉ là $|\beta|<1$. Tính bất đối xứng đi vào qua một số hạng có dấu $\gamma z_{t-1}$ (đòn bẩy là $\gamma<0$ trong quy ước này) tách biệt khỏi một số hạng độ lớn $|z_{t-1}|$.
• Đường cong tác động của tin tức — phương sai kỳ tiếp theo theo cú sốc gần nhất — làm cho tính bất đối xứng trở nên hữu hình và xác minh quy ước dấu của EGARCH chỉ trong một cái liếc.
• Các cải tiến Student-$t$ (dist='t') khắc phục đuôi thông qua một bậc tự do $\nu$ (thường là 3–6 cho crypto); skew-$t$ của Hansen (dist='skewt') thêm một độ lệch $\lambda$ cho một đuôi trái nặng hơn. Chuyển Normal → $t$ là một lợi ích lớn, đáng tin cậy; $t$ → skew-$t$ thường không đáng kể.
• VaR và ES suy ra từ phân phối có điều kiện đã khớp: $\text{VaR}_\alpha = -(\mu_{t+1} + \sigma_{t+1}F_z^{-1}(1-\alpha))$, với phân vị đuôi dày làm cho rủi ro lớn hơn một cách trung thực so với Gaussian. ES (nhất quán, $\approx$ CVaR) nắm bắt khoản lỗ trung bình vượt VaR.
• Backtest bằng Kupiec và Christoffersen. Kupiec kiểm tra tỷ lệ vi phạm; Christoffersen kiểm tra rằng các vi phạm không phân cụm. Một mô hình có thể vượt qua cái này và trượt cái kia — các vi phạm phân cụm là kiểu thất bại nguy hiểm. Hãy backtest nghiêm ngặt ngoài mẫu.
• Kỷ luật hơn độ phức tạp. Chỉ thêm tính bất đối xứng/độ lệch khi nó sống sót qua BIC và một backtest VaR ngoài mẫu. Trên các chuỗi ngắn hoặc kém thanh khoản, mô hình đơn giản hơn thường thắng.

References:

• Nelson, D. B. (1991). Conditional heteroskedasticity in asset returns: A new approach. Econometrica, 59(2), 347-370. DOI
• Glosten, L. R., Jagannathan, R., & Runkle, D. E. (1993). On the relation between the expected value and the volatility of the nominal excess return on stocks. Journal of Finance, 48(5), 1779-1801. DOI
• Bollerslev, T. (1987). A conditionally heteroskedastic time series model for speculative prices and rates of return. Review of Economics and Statistics, 69(3), 542-547. DOI
• Hansen, B. E. (1994). Autoregressive conditional density estimation. International Economic Review, 35(3), 705-730. DOI
• Engle, R. F., & Ng, V. K. (1993). Measuring and testing the impact of news on volatility. Journal of Finance, 48(5), 1749-1778. DOI
• Christoffersen, P. F. (1998). Evaluating interval forecasts. International Economic Review, 39(4), 841-862. DOI
• Kupiec, P. H. (1995). Techniques for verifying the accuracy of risk measurement models. Journal of Derivatives, 3(2), 73-84. DOI
• Black, F. (1976). Studies of stock price volatility changes. Proceedings of the American Statistical Association, Business and Economic Statistics Section, 177-181.
• McNeil, A. J., Frey, R., & Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques and Tools (2nd ed.). Princeton University Press.
• Sheppard, K. (2023). arch: Autoregressive conditional heteroskedasticity models in Python. GitHub.