風險中性密度的 Breeden-Litzenberger 公式
基於此主題： 如何從 BS 波動率推導出隱含機率分佈？
我正在嘗試實施 Breeden-Litzenberger 公式來計算標準普爾 500 指數在某些報價日期的市場隱含風險中性密度。我採取的步驟如下：
第 1 步：提取給定期限 T 的 call_strikes
c_strikes和相應的市場價格css。第 2 步：一旦我有了行使價和市場價格，我就通過函式計算隱含波動率，
ImplieVolatilities.m我 100% 確定該函式有效。*第 3 步：*接下來，我使用 Matlab 命令對隱含波動率曲線進行插值
interp1。現在，向量xq是罷工網格，向量vq是相應的插值隱含波動率。第 4 步：由於我需要波動率曲線，我可以計算市場隱含密度：
[Math Processing Error] 在哪裡是罷工網格大小。我只展示了 Matlab 程序的核心部分：
c_strikes = Call_r_strikes(find(imp_vols > 0)); css = Call_r_prices(find(imp_vols > 0)); impvolss = ImpliedVolatilities(S,c_strikes,r,q,Time,0,css); %implied volatilities xq = (min(c_strikes):0.2:max(c_strikes)); %grid of strikes vq = interp1(c_strikes,impvolss,xq); %interpolatd values %interpolated implied volatilities f = zeros(1,length(xq)); %risk-neutral density function [f] = secondDerivativeNonUniformMesh(x, y) dx = diff(x); %grid size (uniform) dxp = dx(2:end); d2k = dxp.^2; %squared grid size f = (1./d2k).*(y(1:end-2)-2*y(2:end-1)+y(3:end)); end figure(2) plot(xq,f);注意：
blsprice是一個實習生 Matlab 函式，因此絕對可以工作。這是內部點的風險中性 pdf：
數據下表顯示期權數據（左欄：行使價，中：市場價格，右欄：隱含波動率）：
650 387.5 0.337024 700 346.45 0.325662 750 306.8 0.313846 800 268.95 0.302428 850 232.85 0.290759 900 199.15 0.279979 950 168.1 0.270041 975 153.7 0.265506 1000 139.9 0.260885 1025 126.15 0.255063 1050 112.9 0.248928 1075 100.7 0.243514 1100 89.45 0.238585 1125 79.25 0.234326 1150 69.7 0.229944 1175 60.8 0.22543 1200 52.8 0.221322 1225 45.8 0.21791 1250 39.6 0.21486 1275 33.9 0.21154 1300 28.85 0.208372 1325 24.4 0.205328 1350 20.6 0.202698 1375 17.3 0.200198 1400 13.35 0.193611 1450 9.25 0.190194 1500 6.55 0.188604 1550 4.2 0.184105 1600 2.65 0.180255 1650 1.675 0.177356 1700 1.125 0.176491 1800 0.575 0.177828是成熟的時間，無風險利率和股息收益率。現貨價格.
通過看跌期權平價計算無風險利率和股息收益率 [rMath Processing Error]是到期日對應的無風險利率[TMath Processing Error]和[qMath Processing Error]是股息收益率。在數值實現中，我們推導出[rMath Processing Error]和[qMath Processing Error]再次利用看跌期權平價。為此，我們假設以下線性關係：
在哪裡和， 和. 常數[αMath Processing Error]和[βMath Processing Error]然後通過執行線性回歸來計算。因此，無風險利率[rMath Processing Error]和股息收益率[qMath Processing Error]然後由 [數學處理錯誤]
我假設為了近似看漲價格的二階導數在罷工域的邊界（請參閱程式碼的最後一個 for 循環的前 2 個“if”案例），您嘗試設置邊界條件。
在右邊界，您的近似值可能有意義或者至少足夠大。
然而，在左邊界，不能近似為零，因為看漲價格幾乎只是低執行價的內在價值。您可以使用遠期價格減去行使價的貼現差值而不是零，， 如果或者至少足夠小。
這解釋了為什麼您在 LHS 上觀察到奇怪的行為（帶有負 pdf），而我猜如果您放大或刪除該異常值，您應該沒問題。
作為第一步，您可以簡單地計算內部點 xq(2:end-1) 的 pdf（即 for 循環的最後一種情況）並檢查它給出的結果。
仔細觀察，您的輸入價格似乎允許靜態套利機會，更具體地說是蝴蝶套利。因此，您不能期望獲得合理的 pdf。
雖然它遠非完美，但也許下面的腳本會有所幫助。在其中，我嘗試尋找導致套利的罷工，並丟棄相應的輸入以獲得所謂的“乾淨”pdf。這是一種非參數方法。正如評論中所討論的，使用參數化方法（例如，使用無套利表示擬合 vol 微笑）將產生更平滑的 pdf。
function [ ] = someBreedenLitzenbergerCode( ) % Get data and interpolate [K, C, ~] = getData(); dk = 0.5; k = min(K):dk:max(K); c = interp1(K, C, k, 'spline'); % Notifying vertical arbitrage opportunities dC = diff(C); if any(dC > 0) warning('Input call prices allow for arbtirage (vert. spread)'); end % Notifying and dealing with butterfly arbitrage opportunities f = secondDerivativeNonUniformMesh(k, c); idxkArb = find(f<0)+1; if ~isempty(idxkArb) warning('Input call prices allow for arbitrage (butterfly)'); end % Finding indices of input strikes/prices leading to arb idxKArb = []; for i = 1:length(idxkArb) idxHi = find(K>=k(idxkArb(i)),1,'first'); idxLo = find(K<=k(idxkArb(i)),1,'last'); idxKArb = [idxKArb, idxHi, idxLo]; end idxKArb = unique(idxKArb); disp('Strikes leading to butterfly arbitrages'); disp(K(idxKArb)); % Cleaning inputs and calculating pdf K_clean = K; K_clean(idxKArb) = []; C_clean = C; C_clean(idxKArb) = []; k_clean = min(K_clean):dk:max(K_clean); c_clean = interp1(K_clean, C_clean, k_clean, 'spline'); f_clean = secondDerivativeNonUniformMesh(k_clean, c_clean); % Plotting results if ishandle(1), clf(1); end; figure(1); subplot(221); plot(K, C, 'x'); hold on; plot(k, c, 'r-'); plot(k(idxkArb), c(idxkArb), 'rd'); grid on; title('Original call prices, C(K)'); subplot(222); plot(k(2:end-1), f); hold on; plot(k(idxkArb),f(idxkArb-1), 'rd'); grid on; title('Pdf, \phi(K)'); subplot(223); plot(k_clean, c_clean); hold on; lh.Box = 'off'; grid on; title('Clean call prices, C(K)'); subplot(224); plot(k_clean(2:end-1), f_clean); hold on; a = ksr(k_clean(2:end-1), f_clean, 60); % Available at Matlab Central plot(a.x, a.f, 'k-'); lh = legend('FD','FD + Smoothing'); lh.Box = 'Off'; grid on; title('Pdf, \phi(K)'); end function d2y_dx2 = secondDerivativeNonUniformMesh(x, y) dx = diff(x); dxp = dx(2:end); dxm = dx(1:end-1); d2k = dxp.*dxm; d2kp = dxp.*(dxm+dxp); d2km = dxm.*(dxm+dxp); d2y_dx2 = 2./d2km .* y(1:end-2) - 2./d2k .*y(2:end-1) + 2./d2kp .*y(3:end); end function [strikes, prices, volatilities] = getData() mat = [650 387.5 0.33702 700 346.45 0.32566 750 306.8 0.31385 800 268.95 0.30243 850 232.85 0.29076 900 199.15 0.27998 950 168.1 0.27004 975 153.7 0.26551 1000 139.9 0.26088 1025 126.15 0.25506 1050 112.9 0.24893 1075 100.7 0.24351 1100 89.45 0.23858 1125 79.25 0.23433 1150 69.7 0.22994 1175 60.8 0.22543 1200 52.8 0.22132 1225 45.8 0.21791 1250 39.6 0.21486 1275 33.9 0.21154 1300 28.85 0.20837 1325 24.4 0.20533 1350 20.6 0.2027 1375 17.3 0.2002 1400 13.35 0.19361 1450 9.25 0.19019 1500 6.55 0.1886 1550 4.2 0.1841 1600 2.65 0.18025 1650 1.675 0.17736 1700 1.125 0.17649 1800 0.575 0.17783]; strikes = mat(:,1); prices = mat(:,2); volatilities = mat(:,3); end
