程式
Heston Model python MC仿真
我有這個練習。 $ \\ $ 尋找參數的實際值併計算到期的歐式看漲期權的價格 $ T = 0.5 $ 和 $ S_0 = 1 $ 罷工價值 $ K = 0.5,0.6, ……, 1.5 $ . 顯示波動微笑。
我在 python 中發布了我的程式碼，但我不明白為什麼不顯示波動性笑臉。
import numpy as np import seaborn as sns import matplotlib.pyplot as plt from py_vollib_vectorized import vectorized_implied_volatility as implied_vol # Parameters # simulation dependent S0 = 1 # asset price T = 0.5 # time in years r = 0.00003 # risk-free rate N = 125 # number of time steps in simulation M = 10000 # number of simulations # Heston dependent parameters kappa = 1.1 # rate of mean reversion of variance under risk-neutral dynamics theta = 0.02 # long-term mean of variance under risk-neutral dynamics v0 = 0.025 # initial variance under risk-neutral dynamics rho = -0.7 # correlation between returns and variances under risk-neutral dynamics sigma = 0.1 # volatility of volatility def heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M): """ Inputs: - S0, v0: initial parameters for asset and variance - rho : correlation between asset returns and variance - kappa : rate of mean reversion in variance process - theta : long-term mean of variance process - sigma : vol of vol / volatility of variance process - T : time of simulation - N : number of time steps - M : number of scenarios / simulations Outputs: - asset prices over time (numpy array) - variance over time (numpy array) """ # initialise other parameters dt = T/N mu = np.array([0,0]) cov = np.array([[1,rho], [rho,1]]) # arrays for storing prices and variances S = np.full(shape=(N+1,M), fill_value=S0) v = np.full(shape=(N+1,M), fill_value=v0) # sampling correlated brownian motions under risk-neutral measure Z = np.random.multivariate_normal(mu, cov, (N,M)) for i in range(1,N+1): S[i] = S[i-1] * np.exp( (r - 0.5*v[i-1])*dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0] ) v[i] = np.maximum(v[i-1] + kappa*(theta-v[i-1])*dt + sigma*np.sqrt(v[i-1]*dt)*Z[i-1,:,1],0) return S, v rho_p = 0.98 rho_n = -0.98 S_p,v_p = heston_model_sim(S0, v0, rho_p, kappa, theta, sigma,T, N, M) S_n,v_n = heston_model_sim(S0, v0, rho_n, kappa, theta, sigma,T, N, M) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5)) time = np.linspace(0,T,N+1) ax1.plot(time,S_p) ax1.set_title('Heston Model Asset Prices') ax1.set_xlabel('Time') ax1.set_ylabel('Asset Prices') ax2.plot(time,v_p) ax2.set_title('Volatility Monte Carlo Simulation') ax2.set_xlabel('Time') ax2.set_ylabel('Variance') plt.show() # simulate gbm process at time T gbm = S0*np.exp( (r - theta**2/2)*T + np.sqrt(theta)*np.sqrt(T)*np.random.normal(0,1,M) ) fig, ax = plt.subplots() ax = sns.kdeplot(S_p[-1], label=r"$\rho= 0.98$", ax=ax) ax = sns.kdeplot(S_n[-1], label=r"$\rho= -0.98$", ax=ax) ax = sns.kdeplot(gbm, label="GBM", ax=ax) plt.title(r'Asset Price Density under Heston Model') plt.xlim([0.5, 1.5]) plt.xlabel('$S_T$') plt.ylabel('Density') plt.legend() plt.show() rho = -0.7 S,v = heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M) # Set strikes and complete MC option price for different strikes K = np.arange(0.5,1.5) calls = np.array([np.exp(-r*T)*np.mean(np.maximum(S-k,0)) for k in K]) call_ivs = implied_vol(calls, S0, K, T, r, flag='c', q=0, return_as='numpy', on_error='ignore') plt.plot(K, call_ivs, label=r'Call') plt.ylabel('Implied Volatility') plt.xlabel('Strike') plt.title('Implied Volatility Smile from Heston Model') plt.legend() plt.show()感謝您的支持。
發現問題：
- 主要問題：您定義int
S0 = 1類型Python為int。您的股票價格軌跡是建立在基礎之上的S0，因此也是IntS類型（更具體地說，是 Int32）。這意味著您估計中的任何小數都將被截斷。要解決這個問題，只需定義，它現在是float類型。S0 = 1.0- **次要問題：**在定義您的一組罷工時 $ K $ 您需要在 numpy arange 函式中定義步長（自動設置為 1）
K = np.arange(0.1,1.6,0.1)，. 這將為您提供指定的一組罷工。- 在測試程式碼時，我發現該
implied_vol函式可能會遇到一些優化問題，這會產生奇怪的微笑。經過幾次重新執行該功能後，我得到了一個體面的微笑（如下所示）。可能還有其他一些小問題，但解決上述問題給了我一個可行的解決方案。
確認：
我冒昧地修復了上述問題並執行了程式碼。但是，我使用了在Youtube 教程中找到的原始 Heston 參數，您清楚地遵循了這些參數來生成模擬、密度和微笑。這導致以下圖表：
以及下面的微笑：
我將把圖表的解釋留給你。我希望這有幫助。
附錄：Python 程式碼
# Parameters # simulation dependent S0 = 1.0 # asset price T = 0.5 # time in years r = 0.00003 # risk-free rate #r=0.02 N = 252 # number of time steps in simulation M = 10000 # number of simulations # Heston dependent parameters #kappa = 1.1 # rate of mean reversion of variance under risk-neutral dynamics #theta = 0.02 # long-term mean of variance under risk-neutral dynamics #v0 = 0.025 # initial variance under risk-neutral dynamics #rho = -0.7 # correlation between returns and variances under risk-neutral dynamics #sigma = 0.1 # volatility of volatility kappa = 3 theta = 0.20**2 v0=0.25**2 #rho = 0.7 sigma = 0.6 def heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M): """ Inputs: - S0, v0: initial parameters for asset and variance - rho : correlation between asset returns and variance - kappa : rate of mean reversion in variance process - theta : long-term mean of variance process - sigma : vol of vol / volatility of variance process - T : time of simulation - N : number of time steps - M : number of scenarios / simulations Outputs: - asset prices over time (numpy array) - variance over time (numpy array) """ # initialise other parameters dt = T/N mu = np.array([0,0]) cov = np.array([[1,rho], [rho,1]]) # arrays for storing prices and variances S = np.full(shape=(N+1,M), fill_value=S0) v = np.full(shape=(N+1,M), fill_value=v0) # sampling correlated brownian motions under risk-neutral measure Z = np.random.multivariate_normal(mu, cov, (N,M)) for i in range(1,N+1): S[i] = S[i-1] * np.exp((r - 0.5*v[i-1])*dt + np.sqrt(v[i-1] * dt) * Z[i-1,:,0]) v[i] = np.maximum(v[i-1] + kappa*(theta-v[i-1])*dt + sigma*np.sqrt(v[i-1]*dt)*Z[i-1,:,1],0) return S, v rho_p = 0.98 rho_n = -0.98 S_p,v_p = heston_model_sim(S0, v0, rho_p, kappa, theta, sigma,T, N, M) S_n,v_n = heston_model_sim(S0, v0, rho_n, kappa, theta, sigma,T, N, M) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5)) time = np.linspace(0,T,N+1) ax1.plot(time,S_p) ax1.set_title('Heston Model Asset Prices') ax1.set_xlabel('Time') ax1.set_ylabel('Asset Prices') ax2.plot(time,v_p) ax2.set_title('Volatility Monte Carlo Simulation') ax2.set_xlabel('Time') ax2.set_ylabel('Variance') plt.show() # simulate gbm process at time T gbm = S0*np.exp( (r - theta**2/2)*T + np.sqrt(theta)*np.sqrt(T)*np.random.normal(0,1,M) ) fig, ax = plt.subplots() ax = sns.kdeplot(S_p[-1], label=r"$\rho= 0.98$", ax=ax) ax = sns.kdeplot(S_n[-1], label=r"$\rho= -0.98$", ax=ax) ax = sns.kdeplot(gbm, label="GBM", ax=ax) plt.title(r'Asset Price Density under Heston Model') plt.xlim([0.5, 1.5]) plt.xlabel('$S_T$') plt.ylabel('Density') plt.legend() plt.show() rho = -0.7 S,v = heston_model_sim(S0, v0, rho, kappa, theta, sigma,T, N, M) # Set strikes and complete MC option price for different strikes K = np.arange(0.1,1.6,0.1) calls = np.array([np.exp(-r*T)*np.mean(np.maximum(S-k,0)) for k in K]) call_ivs = implied_vol(calls, S0, K, T, r, flag='c', q=0, return_as='numpy', on_error='ignore') plt.plot(K, call_ivs, label=r'Call') plt.ylabel('Implied Volatility') plt.xlabel('Strike') plt.title('Implied Volatility Smile from Heston Model') plt.legend() plt.show()
