European Spread Call Option with Stochastic Volatility
One of the main benefits of Monte Carlo simulations is to price options under multiple factors. By this I refer to multiple underlying asset prices or stochastic volatility or even changing interest rates.
In this tutorial we will explore the pricing of a European Spread Call Option on the difference between two stock indices \(S_1\) and \(S_2\) following a more general stochastic process. The SDE’s will be have stochastic volatility as described by the Heston Model (1993):
Underlying SDE’s under risk-neutral measure:
\(\large dS_{1,t} = (r-\delta_1) S_{1,t} dt + \sqrt{v_{1,t}} S_{1,t} dW^\mathbb{Q}_{1,t}\)
\(\large dS_{2,t} = (r-\delta_2) S_{2,t} dt + \sqrt{v_{2,t}} S_{2,t} dW^\mathbb{Q}_{2,t}\)
Variance Process:
\(\large dv_{1,t} = \kappa_1 (\theta_2 – v_{1,t})dt + \sigma_1 \sqrt{v_{1,t}} dW^\mathbb{Q}_{3,t}\)
\(\large dv_{2,t} = \kappa_2 (\theta_2 – v_{2,t})dt + \sigma_2 \sqrt{v_{2,t}} dW^\mathbb{Q}_{4,t}\)
The Monte Carlo procedure is exactly the same for a spread call option except the correlation matrix between Wiener processes is larger, as in we have four correlated normal variates to simulate the four processes.
Where:
\(\large \rho_{W^\mathbb{Q}} = \begin{vmatrix}1 & \rho_{12} & \rho_{13} & \rho_{14} \\ \rho_{12} & 1 & \rho_{23} & \rho_{24} \\ \rho_{13} & \rho_{23} & 1 & \rho_{34} \\ \rho_{14} & \rho_{24} & \rho_{34} & 1 \notag\end{vmatrix} \)
Notation:
– \(S_t\) Equity spot price, financial index
– \(v_t\) Variance.
– \(C\) European call option price.
– \(K\) Strike price.
– \(W_{1,2}\) Standard Brownian movements.
– \(r\) Interest rate.
– \(\delta\) discrete dividend payment
– \(\kappa\) Mean reversion rate.
– \(\theta\) Long run variance.
– \(v_0\) Initial variance.
– \(\sigma\) Volatility of variance.
– \(\rho\) Correlation parameter.
– \(t\) Current date.
– \(T\) Maturity date.
Example of Nasdaq vs S&P500 Index
A while back we discussed how to futures or options could be used to speculate on the divergence between the spread of the NASDAQ and S&P500 Indices. We follow on this example by incorperating stochastic volatility.
import numpy as np from time import time from scipy.linalg import cholesky # Parameters - as of 1-Mar-22 SP500 = 4373.94 NASD = 13751.40 div_SP500 = 0.0127 div_NASD = 0.0126 K = 9377 # current difference between index points T = 1 r = 0.01828 # 10yr US Treasury Bond Yield # Heston Model Parameters - made up for the example, # we have discussed how to complete this step in a previous video using market option prices theta1 = 0.02 theta2 = 0.03 kappa1 = 0.1 kappa2 = 0.12 sigma1 = 0.05 sigma2 = 0.06 # initial variances vt10 = 0.03 vt20 = 0.03 # Correlation matrix between wiener process under risk-neutral measure rho = np.array([[1,0.5,0.15,0.02], [0.5,1,0.01,0.25], [0.15,0.01,1,0.2], [0.02,0.25,0.2,1]])
Slow Implementation
Let’s first understand the implementation process
# Monte Carlo Specific Parameters N = 100 # discrete time steps M = 1000 # number of simulations # Start Timer start_time = time() # Precompute constants dt = T/N # log normal prices lnS1 = np.log(NASD) lnS2 = np.log(SP500) # Heston model adjustments for time steps kappa1dt = kappa1*dt kappa2dt = kappa2*dt sigma1sdt = sigma1*np.sqrt(dt) sigma2sdt = sigma2*np.sqrt(dt) vt1 = vt10 vt2 = vt20 # Perform (lower) cholesky decomposition lower_chol = cholesky(rho, lower=True) # Standard Error Placeholders sum_CT = 0 sum_CT2 = 0 # Monte Carlo Method for i in range(M): # for each simulation i in M lnSt1 = lnS1 lnSt2 = lnS2 for j in range(N): # for each time step j in N # Generate correlated Wiener variables Z = np.random.normal(0.0, 1.0, size=(4)) W = Z @ lower_chol # Simulate variance processes vt1 = vt1 + kappa1dt*(theta1 - vt1) + sigma1sdt*np.sqrt(vt1)*W[2] vt2 = vt2 + kappa2dt*(theta2 - vt2) + sigma2sdt*np.sqrt(vt2)*W[3] # Simulate log asset prices nu1dt = (r - div_NASD - 0.5*vt1)*dt nu2dt = (r - div_SP500 - 0.5*vt2)*dt lnSt1 = lnSt1 + nu1dt + np.sqrt(vt1*dt)*W[0] lnSt2 = lnSt2 + nu2dt + np.sqrt(vt2*dt)*W[1] St1 = np.exp(lnSt1) St2 = np.exp(lnSt2) CT = max(0, (St1 - St2) - K) sum_CT = sum_CT + CT sum_CT2 = sum_CT2 + CT*CT # Compute Expectation and SE C0 = np.exp(-r*T)*sum_CT/M sigma = np.sqrt( (sum_CT2 - sum_CT*sum_CT/M)*np.exp(-2*r*T) / (M-1) ) SE = sigma/np.sqrt(M) print("Call value is ${0} with SE +/- {1}".format(np.round(C0,2),np.round(SE,2))) print("Calculation time: {0} sec".format(round(time() - start_time,2)))
Fast Implementation
Vectorizing option pricing
# Monte Carlo Specific Parameters N = 100 # discrete time steps M = 1000 # number of simulations # Start Timer start_time = time() # Precompute constants dt = T/N # Heston model adjustments for time steps kappa1dt = kappa1*dt kappa2dt = kappa2*dt sigma1sdt = sigma1*np.sqrt(dt) sigma2sdt = sigma2*np.sqrt(dt) # Perform (lower) cholesky decomposition lower_chol = cholesky(rho, lower=True) # Generate correlated Wiener variables Z = np.random.normal(0.0, 1.0, size=(N,M,4)) W = Z @ lower_chol # arrays for storing prices and variances lnSt1 = np.full(shape=(N+1,M), fill_value=np.log(NASD)) lnSt2 = np.full(shape=(N+1,M), fill_value=np.log(SP500)) vt1 = np.full(shape=(N+1,M), fill_value=vt10) vt2 = np.full(shape=(N+1,M), fill_value=vt20) for j in range(1,N+1): # Simulate variance processes vt1[j] = vt1[j-1] + kappa1dt*(theta1 - vt1[j-1]) + sigma1sdt*np.sqrt(vt1[j-1])*W[j-1,:,2] vt2[j] = vt2[j-1] + kappa2dt*(theta2 - vt2[j-1]) + sigma2sdt*np.sqrt(vt2[j-1])*W[j-1,:,3] # Simulate log asset prices nu1dt = (r - div_NASD - 0.5*vt1[j])*dt nu2dt = (r - div_SP500 - 0.5*vt2[j])*dt lnSt1[j] = lnSt1[j-1] + nu1dt + np.sqrt(vt1[j]*dt)*W[j-1,:,0] lnSt2[j] = lnSt2[j-1] + nu2dt + np.sqrt(vt2[j]*dt)*W[j-1,:,1] St1 = np.exp(lnSt1) St2 = np.exp(lnSt2) # Compute Expectation and SE CT = np.maximum(0, (St1[-1] - St2[-1]) - K) C0 = np.exp(-r*T)*np.sum(CT)/M sigma = np.sqrt( np.sum( (np.exp(-r*T)*CT - C0)**2) / (M-1) ) SE= sigma/np.sqrt(M) print("Call value is ${0} with SE +/- {1}".format(np.round(C0,2),np.round(SE,2))) print("Calculation time: {0} sec".format(round(time() - start_time,2)))
Visualisation
import matplotlib.pyplot as plt t = np.linspace(0,1,len(St1)) fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(8,6)) ax[0,0].plot(t,St1[:,:4]) ax[0,0].set_title('$NASDAQ$') ax[0,1].plot(t,St2[:,:4]) ax[0,1].set_title('$S&P500$') ax[1,0].plot(t,vt1[:,:10]) ax[1,0].set_title('$v_1$') ax[1,1].plot(t,vt2[:,:10]) ax[1,1].set_title('$v_2$') fig.suptitle("Index Spread Option with Stochastic Volatility", fontsize=14) plt.show()