如何在施瓦茨均值回歸模型下執行蒙特卡羅模擬來為遠期合約定價？
目標：（1）實施歐拉顯式方法求解施瓦茨均值回歸模型下期權價格的偏微分方程。(2) 與蒙特卡羅模擬比較。
我堅持第 1 點（請參閱此問題），因此同時我正在研究第 2 點。基本上，我有 200 個商品（石油）的月現貨價格，遵循 Schwartz 均值回歸 SDE $ dS = \alpha(\mu-\log S)Sdt + \sigma S dW $ . 該商品是遠期合約的標的資產，交割時價格 $ F(S,T) $ ， 在哪裡 $ S $ 是時間 0 的現貨價格和 $ T $ 是到期時間。這樣的契約相當於有一個選擇權 $ \text{payoff}(S_T)=S_T-K $ ， 在哪裡 $ K=F(S,T) $ 是執行價格。
讓 $ \tau=T-t $ 為到期時間，遠期合約價格的公式為
$$ \tag{1}F(S,\tau)=\exp\Big(e^{-\alpha\tau}\log S +(\mu-\frac{\sigma^2}{2\alpha}-\frac{\mu-r}{\alpha})(1-e^{-\alpha\tau})+\frac{\sigma^2}{4\alpha}(1-e^{-2\alpha\tau})\Big) $$
等價期權的價值是 $ V(S,t)=e^{-r(T-t)}F(S,t) $ . ( $ \alpha, \sigma, \mu $ 根據數據估算）
知道這些資訊後，如何執行（我使用的是 matlab，但其他語言也可以）蒙特卡羅模擬來為選項定價？
我讀了這些是步驟
- 為標的資產生成大量隨機價格路徑
- 為每條路徑計算期權的收益
- 計算所有收益的平均值
- 折扣到今天的平均值
但由於我使用的是石油現貨價格的真實數據，我有點困惑：我們如何通過生成隨機現貨價格路徑來獲得具有已知現貨價格的標的商品的期權價值？
此外，如何使用上述公式執行後續步驟？
編輯：第一次嘗試
我編寫的程式碼計算我認為（你能確認嗎？）是期權定價的精確解（
V_exact），然後通過蒙地卡羅模擬（V_Monte_Carlo）計算近似解。MC 模擬使用第一個現貨價格（真實數據）和估計參數“隨機”計算下一個現貨價格。計算現貨價格，而不是使用等式 $ dS $ 我在上面寫的它更容易使用從 $ dS $ 通過讓 $ X=\log S $ .$$ \tag{2}X_t = X_0 e^{-\alpha t} + \mu^*(1-e^{-\alpha t}) + \underbrace{\sigma e^{-\alpha t} \int_0^t e^{\alpha s} dW_s}_Z $$
在哪裡 $ \mu^* = \mu-\dfrac{\sigma^2}{2\alpha} $ 和 $ Z \sim \sigma\sqrt{\dfrac{1-e^{-2\alpha t}}{2\alpha}} \cdot N(0,1) $
這些是蒙特卡羅模擬的步驟：
- 遠期合約的第一個價格
F- 以及期權的第一個價格V- 設置為 0，因為簽訂遠期合約沒有成本。計算第一個值 $ X $ 取第一個現貨價格的對數- 計算下一個值 $ X $ 使用公式（2）
- 通過取指數計算現貨價格的下一個值 $ X $
F使用公式（1）計算遠期合約的下一個價格V使用計算期權的下一個價格 $ V(S,t)=e^{-r(T-t)}F(S,t) $V_Monte_Carlo計算期權價格的平均值- 重複步驟 2-6，直到計算出所有值
但是，如下圖所示，通過 MC 模擬獲得的期權價格曲線是平滑的，而我所期望的是一條類似於精確解的曲線。所以我在某個地方犯了一些錯誤
% S = oil spot prices S = [ 22.93 15.45 12.61 12.84 15.38 13.43 11.58 15.10 14.87 14.90 15.22 16.11 18.65 17.75 18.30 18.68 19.44 20.07 21.34 20.31 19.53 19.86 18.85 17.27 17.13 16.80 16.20 17.86 17.42 16.53 15.50 15.52 14.54 13.77 14.14 16.38 18.02 17.94 19.48 21.07 20.12 20.05 19.78 18.58 19.59 20.10 19.86 21.10 22.86 22.11 20.39 18.43 18.20 16.70 18.45 27.31 33.51 36.04 32.33 27.28 25.23 20.48 19.90 20.83 21.23 20.19 21.40 21.69 21.89 23.23 22.46 19.50 18.79 19.01 18.92 20.23 20.98 22.38 21.78 21.34 21.88 21.69 20.34 19.41 19.03 20.09 20.32 20.25 19.95 19.09 17.89 18.01 17.50 18.15 16.61 14.51 15.03 14.78 14.68 16.42 17.89 19.06 19.65 18.38 17.45 17.72 18.07 17.16 18.04 18.57 18.54 19.90 19.74 18.45 17.33 18.02 18.23 17.43 17.99 19.03 18.85 19.09 21.33 23.50 21.17 20.42 21.30 21.90 23.97 24.88 23.71 25.23 25.13 22.18 20.97 19.70 20.82 19.26 19.66 19.95 19.80 21.33 20.19 18.33 16.72 16.06 15.12 15.35 14.91 13.72 14.17 13.47 15.03 14.46 13.00 11.35 12.51 12.01 14.68 17.31 17.72 17.92 20.10 21.28 23.80 22.69 25.00 26.10 27.26 29.37 29.84 25.72 28.79 31.82 29.70 31.26 33.88 33.11 34.42 28.44 29.59 29.61 27.24 27.49 28.63 27.60 26.42 27.37 26.20 22.17 19.64 19.39 19.71 20.72 24.53 26.18 27.04 25.52 26.97 28.39 ]; r = .1; % yearly instantaneous interest rate T = 1/2; % expiry time alpha = 0.0692; sigma = 0.0876; % values estimated from data mu = 3.0582; t = linspace(0,T,numel(S)); tau = T-t; lambdahat = (mu-r)/alpha; muhat = mu-sigma^2/2/alpha-lambdahat; %% Is this V the exact solution? % this F is the formula (1) on stackexchange F = exp( exp(-alpha*tau).*log(S) + muhat*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) ); % K = F(end); % shouldn't we use also the strike price? V_exact = exp( -r*tau ) .* F; % V contains the prices of the option plot(t,S) hold on plot(t,V_exact) %% Monte Carlo simulation n = numel(S); N = 10000; % number of samples mustar = mu-sigma^2/2/alpha; X = zeros(N,n); X(:,1) = log(S(1)); S = zeros(N,n); S(:,1) = S(1); F = zeros(N,n); V = zeros(N,n); V_Monte_Carlo = zeros(1,n); dt = 1; % which values for dt has to be pick? for k = 2:n tau = T-t(k); % this X is the formula (2) on stackexchange X(:,k) = X(k-1)*exp(-alpha*dt) + mustar*(1-exp(-alpha*dt)) + sigma*sqrt(1-exp(-2*alpha*dt))/sqrt(2*alpha)*randn(N,1); S(:,k) = exp(X(:,k)); F(:,k) = exp( exp(-alpha*tau).*log(S(:,k)) + muhat*(1-exp(-alpha*tau)) + sigma^2/4/alpha*(1-exp(-2*alpha*tau)) ); V(:,k) = exp(-r*tau)*F(:,k); V_Monte_Carlo(k) = mean(V(:,k)); end plot(t,V_Monte_Carlo,'g') legend('Spot prices (S)','Option prices (V) from exact solution','Option prices from Monte-Carlo simulation')
讓我們來看看。
對於該度量下的股票價格，您有以下 SDE $ P $ :
$$ dS(t) = \alpha \cdot (\mu - \log S(t)) \cdot S(t) \cdot dt + \sigma \cdot S \cdot dW^P(t), $$
有初始條件 $ S(0) = S_0 $ . 此外，定義 $ X(t) = \log S(t) $ , 假設風險的市場價格不變 $ \lambda = \mu - r $ 並執行更改措施，您將獲得流程 $ X(t) $ 在風險中性措施下 $ Q $ 作為：
$$ dX(t) = \alpha \cdot \left( \mu - \frac{\sigma^2}{2\alpha} - \frac{(\mu - r) \cdot \sigma}{\alpha} - X(t) \right) \cdot dt + \sigma\cdot dW^Q(t). $$
讓我們定義漂移和擴散函式：
$$ \begin{aligned} f(X_t, t) &= \alpha \cdot \left( \mu - \frac{\sigma^2}{2\alpha} - \frac{(\mu - r) \cdot \sigma}{\alpha} - X_t \right),\ g(X_t, t) &= \sigma. \end{aligned} $$
現在，只需應用 EulerMaruyama 方案：
$$ X(t_i) = X(t_{i-1}) + f(X(t_{i-1}), t_{i-1}) \cdot \Delta t_i + g(X(t_{i-1}), t_{i-1}) \cdot \left[ W(t_i) - W(t_{i-1}) \right], $$
在哪裡 $ W(t_i) - W(t_{i-1}) = \mathcal{N}(0, \Delta t_i) $ .
我在 Python 中也沒有在 Matlab 中編寫程式碼，但你可以在維基百科中找到這些實現，這裡是。
我建議使用 Julia 語言程式來模擬隨機微分方程。
您想計算以下衍生品價格：
$$ \begin{aligned} V(0) &= E_t^Q \left[D(0, T) \cdot \left( S(T) - K \right) \right]\ V(0) &= E_t^Q \left[D(0, T) \cdot \left( \exp \left( X(T) \right) - K \right) \right]\ \end{aligned} $$
由於折扣因子是確定性的，並且 $ K $ 是恆定的：
$$ \begin{aligned} V(0) &= D(0, T) \cdot \left( E_t^Q \left[ \exp \left( X(T) \right) \right] - K \right) \end{aligned} $$
現在，通過蒙地卡羅方法來近似期望。我們必須模擬這個過程 $ X(t) $ 為了 $ N $ 路徑，及時評估 $ T $ 並應用指數函式。
$$ E_t^Q \left[ \exp \left( X(T) \right) \right] = \frac{1}{N} \cdot \sum_{n=1}^N \exp \left( X_n(T) \right) $$
衍生品價格為
$$ V(0) = D(0, T) \cdot \left[ \frac{1}{N} \cdot \sum_{n=1}^N \exp \left( X_n(T) \right) - K \right] $$
對於任何 $ K $ . 特別是，如果您使用 $ K = F(S_0, T) $ 你會得到衍生價格 $ V(0) = 0 $ ，這就是現實中發生的事情，其中 $ K $ 假定取該值。請注意您的功能 $ F(S_0, T) $ 可能有錯誤，有一個 $ \sigma $ 缺少乘法 $ \mu - r $ ，請參閱下面的函式
F定義中的程式碼。無論如何，刪除它 $ \sigma $ 在兩個 SDE 中 $ X(t) $ 和 $ F(S_0, T) $ 屈服於正確的價格。但我想強調我在本文第 2.1 節看到的這種差異。就是這樣了。這裡有 Julia 中的程式碼：
using StochasticDiffEq using Parameters using Statistics using Printf using Plots function drift(u, p, t) @unpack α, μ, σ, r = p return α * ((μ - σ^2 / (2α) - (μ - r) / α) - u) end function diffusion(u, p, t) @unpack σ = p return σ end function F(S0, τ, p) @unpack α, μ, σ, r = p return exp( exp(-α * τ) * log(S0) + (μ - σ^2 / (2α) - (μ - r) / α) * (1 - exp(-α * τ)) + σ^2 / (4 * α) * (1 - exp(-2 * α * τ)) ) end let params = ( α = 0.0692, σ = 0.0876, μ = 0.25, r = 0.10 ) S0 = 22.93 X0 = log(S0) T = 1/2 trajectories = 100_000 prob = SDEProblem(drift, diffusion, X0, (0., T), params) sde = EnsembleProblem(prob) sol = solve(sde, SRIW1(), trajectories = trajectories) evs = zeros(trajectories) for n in 1:trajectories evs[n] = exp(sol[n](T)) end # fair value and standard deviation r = params.r K = sort(push!(collect(15:0.5:25), F(S0, T, params))) x = [] y = [] for k in K μ = mean(exp(-r * T) .* (evs .- k)) σ = stdm(evs, μ; corrected = true) / sqrt(trajectories) mcapprox = μ analytic = exp(-r * T) * (F(S0, T, params) - k) push!(x, mcapprox) push!(y, analytic) # @printf("%e\t%e\n", mcapprox, analytic) end plot(K, x) plot!(K, y) end我得到一個等於零的衍生品價格 $ K = F(S_0, T) $ . 此外，我使用蒙地卡羅方法得到的價格與使用解析解時的價格相同 $ K $ 是別的什麼！
在這裡，您可以查看帶有結果的圖，其中近似值與分析結果相匹配，並且 $ V(0) $ 等於 $ 0 $ 預期時：
當然價格可以是正數或負數，這取決於罷工的價值 $ K $ . 如果 $ K < F(S_0, T) $ , 期權價格是正的，因為為了公平起見，持有人（將收到 $ S(T) $ ）， 必須付款 $ V(0) $ 提前。如果發生相反的情況 $ K > F(S_0, T) $ .
