Duan (1995) 使用 MATLAB 的 GARCH 期權定價模型
這是複制段在他的論文“GARCH 期權定價模型”中提出的期權定價模型的 MATLAB 程式碼。但是，文件中估計的參數與論文中提供的參數不匹配。我試圖修復它，但我仍然得到錯誤的參數值。
這是最大概似估計的 .m 文件：
function y = findGARCH_LLy(params,S,rf) % Finds log-likelihood for the GARCH option pricing model. alpha0 = params(1); alpha1 = params(2); beta1 = params(3); lambda = params(4); N = length(S); % Define the returns. r = log(S(2:N)./S(1:N-1)); r = [0; r]; % Initialize the GARCH log-likelihood at t=1. h(1) = var(r); e(1) = 0; LL(1) = -0.5*log(2*pi)-0.5*log(h(1))-0.5*e(1)^2/h(1); % Find the rest of the GARCH log-likelihood. for t=2:N h(t) = alpha0 + alpha1*(e(t-1) - lambda*sqrt(h(t-1)))^2 + beta1*h(t-1); e(t) = log(S(t)/S(t-1)) - rf + 0.5*h(t); LL(t) = 0.5*log(2*pi)-0.5*log(h(t))-0.5*e(t)^2/h(t); end % Return the negative log-likelihood. y = -sum(LL);這是主要的 .m 文件：
clc; clear; % Input the price levels and dates. [P, Dates] = xlsread('SP 100 Prices.xls','Sheet1'); Prices = P(:,2); % Risk free rate rf = 0; % Starting values for optimization. % Use estimates from Duan's paper as starting values. start = [0.000015 0.19 0.72 0.007]; A = [0 1 1 0]; b = 1; lb = [0 0 0 0]; ub = [+Inf +Inf +Inf +Inf]; params = fmincon(@(b) findGARCH_LLy(b,Prices,rf), start, A, b, [], [], lb, ub) alpha0 = params(1); alpha1 = params(2); beta1 = params(3); lambda = params(4); % Find the standard deviation (sigma) implied by the parameters. % Assume 365 days per year. variance = alpha0 / (1 - alpha1*(1+lambda^2) - beta1); sigma = sqrt(variance)*sqrt(365);為了資訊的完整性，論文中的估計參數值為 $ \alpha_0=1.524\times 10^{-5} $ , $ \alpha_1=0.1883 $ , $ \beta_1=0.7162 $ 和 $ \lambda=7.452\times 10^{-3} $ , 和標準差 $ \sigma=24.13% $ . 此外，該模型適用於 1986 年 1 月 2 日至 1989 年 12 月 15 日的標準普爾 100 每日指數。
有人可以幫我找出我的錯誤所在嗎？非常感謝您的幫助。
你的函式返回（減去）對數概似對我來說似乎很奇怪，我會去
function y = findGARCH_LLy(params,S,rf) % Finds log-likelihood for the GARCH option pricing model. alpha0 = params(1); alpha1 = params(2); beta1 = params(3); lambda = params(4); N = length(S); % Define the returns (pad first return with zero) r = [0, diff(log(S))]; % Infer the other conditional instantaneous variances [see EQ1] h(1) = var(r); % Initialize the conditional variance recursion for i=2:N h(i) = alpha0 + alpha1*(r(i-1) - rf + .5*h(i-1) - lambda*sqrt(h(i-1)))^2 + beta1*h(i-1)]; end % Return the negative log-likelihood (up to constant terms) [see EQ2] LL = log(h) + 1./h .* (r - rf + .5*h).^2; y = -sum(LL);話雖如此，如果一切正常但您的輸入價格錯誤，您將永遠無法獲得正確的參數，因此請從檢查此開始。
REM：我沒有嘗試執行它，所以我可能留下了一些錯別字，但你應該明白這一點。
編輯
我確實犯了一些錯別字（現在希望已更正），但我假設您正在使用物理測量下的以下動態 $$ \begin{align} \ln{\left(\frac{S_{t+1}}{S_t}\right)} = R_{t+1} &= r_F - \frac{1}{2}h_{t+1} + \sqrt{h_{t+1}}z_{t+1} \ h_{t+1} &= \alpha_0 + \alpha_1 (\sqrt{h_t} z_t - \sqrt{h_t} \lambda)^2 + \beta_1 h_t \end{align} $$
和 $ \mathbf{z}=(z_1,\dots,z_N) $ iid 高斯變數。
所以第一步是從數據序列中恢復隱藏的條件變異數。這是通過將上面的第一個等式插入到第二個等式中來完成的，從而產生以下遞歸： $$ h_{t+1} = \alpha_0 + \alpha_1 \left(R_t - r_f + \frac{1}{2}h_{t} - \sqrt{h_t} \lambda \right)^2 + \beta_1 h_t \tag{EQ1} $$ 從…開始 $ h_1 = \text{var}(\mathbf{r}) $ 和 $ \mathbf{r}=(r_1,\dots,r_N) $ 您的範例日誌返回。
從那裡開始，你就有了，有條件地 $ h_t $ $$ R_t \sim \mathcal{N}(r_F - \frac{1}{2}h_{t}, h_t) $$ 假設你有一個 iid 樣本 $ N $ 返回你然後形成對數概似為
$$ \begin{align} p_{R_1,R_2,…,R_N}(r_1,r_2,…,r_N; \Theta) &= \prod_{i=1}^N p_{R_i|\phi_{i-1}}(r_i; \Theta) \ &= \prod_{i=1}^N \frac{1}{\sqrt{2\pi h_i}} \exp \left( -\frac{1}{2h_i}(r_i - (r_F-\frac{1}{2}h_i))^2\right) \ &:= \mathcal{L}(\Theta) \end{align} $$
因此最後是負對數概似 $$ \begin{align} l(\Theta) &:= -\ln(\mathcal{L}(\Theta)) \nonumber\ &\propto \sum_{i=1}^N \ln(h_i) + \frac{1}{h_i} (r_i - r_F +\frac{1}{2}h_i))^2 \tag{EQ2} \end{align} $$
在哪裡， $ \mathbf{r}=(r_1,\dots,r_N) $ 是觀察到的對數回報和瞬時變異數 $ \mathbf{h}=(h_1,\dots,h_N) $ 如上所述獲得。
現在，請參見上述程式碼中的等式 (EQ1) 和 (EQ2)。