Pricing Options with Fourier Series P3 — The Heston Model
Unlock the Power of Fourier Transform & The Heston Model to Master European Options Pricing — Dive into Part 3 of our Fourier Series journey. Missed the basics? Check out Part 1 & Part 2!
- 👨💻 Interactive Guide: Grab the Github!
- 📖 Based on cutting-edge “Stochastic Calculus” lectures by Prof. Michel Vellekoop at the University of Amsterdam. Also I included the articles which this project is based on.
Contents:
- Recap: Fourier Transformation, characteristic function
- Material used: The Heston Model, Characteristic Function of the Heston Model, Integration Limits (cumulants)
- Method of calculating IV, Risk Free Rate Model (Nelson Siegel Svensson)
- Collecting Option chain Data from yfinance
- 3D Option Chain Plots: (Price, IV) vs (Time) vs (Strike)
- Calibration of Heston Model Parameters. min IV(predicted) vs IV(observed)
- Error Analysis: Weighted RMSE
- Heston IV surface
- Discussion: Limitations of this method. Why Deep In/Out of the money options are hard to price correctly. (“The Fourier Trap”)
“As far as the laws of mathematics refer to reality, they are not certain, and as far as they are certain, they do not refer to reality.”
Albert Einstein
Fourier Transformation, characteristic function of Random Variable
Let’s quickly revisit the mathematical tool we’re using to price a European Put option. In Part 1, I delved into the Fourier Transformation, a technique that converts a function into a series of cosine waves. In Part 2, I discussed how to leverage the characteristic function of a random variable’s distribution to accurately approximate its PDF/CDF. This allows us to price the European Put option effectively.
What’s missing from our toolkit are the characteristic function of the distribution and the integration limits (cumulants), both of which were covered in Part 2. There, I demonstrated how to price a European Put option using Geometric Brownian Motion (GBM) and compared it to its analytical counterpart. The results showed that this method is highly accurate.
The Heston Model
Below are the equations that the Heston model uses to describe the rate of change of S(t) through a stochastic volatility model. I’m going to assume you’re already familiar with the Heston Model. Instead of rehashing its basics, let’s focus on some of its key properties that make Fourier Transform (FT) a valuable tool for this model.
If we employ Euler-Monte Carlo (MC) to approximate dS, we run into two issues:
- Time Efficiency: MC takes significantly longer than FT. When calibrating across an option chain containing approximately 6,000 European options, each set of parameters tested in MC took around 10 minutes. In contrast, FT took less than a second. This makes MC computationally infeasible when you have to test hundreds of parameter combinations.
- Negative Values in Vt: The term ζ×sqrt{dt×Vt}×Z can yield negative values for Vt. This is problematic because the square root of a negative number is undefined in the realm of real numbers. While you can modify the model to use |Vt| and a very small dt, this workaround makes it impossible to compare MC and FT results.
There are alternative approaches for MC, such as adding more Ito terms or explicitly calculating the Vt integral.
Source: Milan Mrázek and Jan Pospíšil [1]
Characteristic Function of The Heston Model
I highly recommend taking a look at “The Little Heston Trap”[2]. This source offers an in-depth explanation of the correct equation and dynamics for the characteristic function of the Heston model. Understanding this characteristic function is crucial for achieving stable approximations, regardless of the parameters you’re working with.
For this project, I relied on the characteristic function provided in my “Stochastic Calculus” course. The equation is as follows:
def heston_char(u):
t0 = 0.0 ; q = 0.0
m = log(S0) + (r - q)*(T-t0)
D = sqrt((rho*zeta*1j*u - kappa)**2 + zeta**2*(1j*u + u**2))
C = (kappa - rho*zeta*1j*u - D) / (kappa - rho*zeta*1j*u + D)
beta = ((kappa - rho*zeta*1j*u - D)*(1-exp(-D*(T-t0)))) / (zeta**2*(1-C*exp(-D*(T-t0))))
alpha = ((kappa*theta)/(zeta**2))*((kappa - rho*zeta*1j*u - D)*(T-t0) - 2*log((1-C*exp(-D*(T-t0))/(1-C))))
return exp(1j*u*m + alpha + beta*v0)🚨 Heads Up: Be cautious with the α and β values in the equation. They are not the same as the integration limits. (Don’t ask me how I know.)
Integration Limits, cumulants
Cumulants, as defined by Wikipedia, are essentially the exponentiated central moments of a distribution. The characteristic function and the cumulants are intrinsically linked, allowing us to use the second cumulant to calculate the variance. This, in turn, helps us define the integration limits [a,b].
In Part 2, I discussed the first and second cumulants for the GBM and how they relate to the integration limits. Similarly, the Heston Model has a cumulant-generating function from which c2 is derived.
Similarly, the Heston Model has a cumulant-generating function from which c2 is derived. Source: “A Generalized Fourier Transform Approach to Risk Measures” [3]
Defining these limits can be as challenging as the main model itself. The Heston PDF/CDF is notably skewed and exhibits higher kurtosis compared to GBM. While we previously used z=12, a higher zz value, such as z=24, yields better results here.
🔗 The Fourier Trap: The values of Z (integration limits) and N (number of terms) are interconnected. Increasing Z necessitates a higher N, but there’s a limit to this, on which I’ll elaborate later on as “The Fourier Trap.”
Python Code for Pricing a European Option
Below is the complete Python code that allows you to price a European Option using the Fourier Transformed Heston Model. A couple of key points to note:
- We’re pricing a put option, and the formula C=P+S0−K×exp(−r×T) is used to price a call option.
- This function was a labor of love, so if you find it useful, consider having a drink in my honor and maybe link back to this blog.
🍻 Cheers to Numba: Using a Just-In-Time (JIT) compiler like Numba can speed up the code execution by a whopping 164 times! We’re talking 26 ms (Numpy) vs. 159 µs (Numba). So, let’s raise a glass to Numba for making this process lightning-fast in Python.
import numpy as np
from numpy import sqrt, exp, pi, cos, sin, log, abs
from numba import njit
@njit
def Fourier_Heston_Put(S0, K, T, r,
# Heston Model Paramters
kappa, # Speed of the mean reversion
theta, # Long term mean
rho, # correlation between 2 random variables
zeta, # Volatility of volatility
v0, # Initial volatility
opt_type,
N = 1_012,
z = 24
):
def heston_char(u):
t0 = 0.0 ; q = 0.0
m = log(S0) + (r - q)*(T-t0)
D = sqrt((rho*zeta*1j*u - kappa)**2 + zeta**2*(1j*u + u**2))
C = (kappa - rho*zeta*1j*u - D) / (kappa - rho*zeta*1j*u + D)
beta = ((kappa - rho*zeta*1j*u - D)*(1-exp(-D*(T-t0)))) / (zeta**2*(1-C*exp(-D*(T-t0))))
alpha = ((kappa*theta)/(zeta**2))*((kappa - rho*zeta*1j*u - D)*(T-t0) - 2*log((1-C*exp(-D*(T-t0))/(1-C))))
return exp(1j*u*m + alpha + beta*v0)
# # Parameters for the Function to make sure the approximations are correct.
c1 = log(S0) + r*T - .5*theta*T
c2 = theta/(8*kappa**3)*(-zeta**2*exp(-2*kappa*T) + 4*zeta*exp(-kappa*T)*(zeta-2*kappa*rho)
+ 2*kappa*T*(4*kappa**2 + zeta**2 - 4*kappa*zeta*rho) + zeta*(8*kappa*rho - 3*zeta))
a = c1 - z*sqrt(abs(c2))
b = c1 + z*sqrt(abs(c2))
h = lambda n : (n*pi) / (b-a)
g_n = lambda n : (exp(a) - (K/h(n))*sin(h(n)*(a - log(K))) - K*cos(h(n)*(a - log(K)))) / (1 + h(n)**2)
g0 = K*(log(K) - a - 1) + exp(a)
F = g0
for n in range(1, N+1):
h_n = h(n)
F += 2*heston_char(h_n) * exp(-1j*a*h_n) * g_n(n)
F = exp(-r*T)/(b-a) * np.real(F)
F = F if opt_type == 'p' else F + S0 - K*exp(-r*T)
return F if F > 0 else 0
print('Speed Analysis of Fourier_Heston_Put: Per Option')
%timeit Fourier_Heston_Put(S0,K, T, r, kappa, theta, rho, zeta, v0, 'p', N = 1_012, z = 24)
""" Results (Macbook Air M1):
Speed Analysis of Fourier_Heston_Put: Per Option
159 µs ± 8.51 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
"""Implied Volatility (IV) and the Risk-Free Rate Model (Nelson Siegel Svensson Model)
Calculating Implied Volatility (IV)
When dealing with deep in-the-money or out-of-the-money options, calculating IV becomes a numerical game. This task can be as computationally intensive as the Fourier Transformation (FT). Typically, a Newton method-based approach is employed, which, while accurate, can be time-consuming due to its iterative nature.
🛠 Tool of Choice: For this reason, I opted for the py_vollib_vectorized.vectorized_implied_volatility library. It leverages Numba for faster calculations and is incredibly user-friendly.
import py_vollib_vectorized
price = 0.10 ; S = 95 ; K = 100 ; t = .2 ; r = .2 ; flag = 'c'
def implied_volatility(price, S, K, t, r, flag):
return py_vollib_vectorized.vectorized_implied_volatility(
price, S, K, t, r, flag, q=0.0, on_error='ignore', model='black_scholes_merton',return_as='numpy')
print('Speed Analysis of the implied volatility Function: Per Option')
%timeit implied_volatility(price, S, K, t, r, flag)
""" Result (Mackbook Air M1):
Speed Analysis of the implied volatility Function: Per Option
49 µs ± 2.79 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
"""Risk-Free Rate (RFR) Curve
To determine the RFR, I used U.S. Treasury bonds, which can be accessed through this link.
Nelson Siegel Svensson Method
To generate a yield curve, I employed the Nelson Siegel Svensson method. The equation for this method is as follows:
📈 Implementation: The process is straightforward. First, collect the yields from the U.S. Treasury. Then, fit these yields using the Nelson Siegel Svensson model with the Python library nelson_siegel_svensso(link).
from nelson_siegel_svensson.calibrate import calibrate_nss_ols
# Collected at : 06/01/2023
yield_maturities = np.array([1/12, 2/12, 3/12, 4/12, 6/12, 1, 2, 3, 5, 7, 10, 20, 30])
yields = np.array([5.30,5.39,5.50,5.50,5.44,5.11,4.33,3.98,3.70,3.66,3.61,3.98,3.84])
#NSS model calibrate
curve_fit, status = calibrate_nss_ols(yield_maturities,yields)
#Results: NelsonSiegelSvenssonCurve(beta0=1.2595121714700128, beta1=4.367186541206554, beta2=-1457.4165045532893, beta3=1458.9304400823748, tau1=6.83997027106895, tau2=6.877834685907713)
Collecting Option Chain Data from yfinance
Focus on ^SPX (S&P 500 Options)
We’ll be using the yfinance platform to access a complete option chain for ^SPX (S&P 500 Options). But before diving into the analysis, it’s crucial to clean and format the data.
- Retrieve the expiration dates and the last closing price of the S&P 500. This sets the stage for the rest of the analysis.
SPX = yf.Ticker("^SPX")
# Step 1: get last clossing price
S0 = SPX.history(perdiod='3mo')['Close'][-1] ; print(f'S0={S0}')
# Step 2: get Expiration Dates
expiration_dates = SPX.options2. Looping Through Expiration Dates: Iterate over all expiration dates to collect data on both Puts and Calls, giving us a comprehensive view of the option chain.
- Days Since Last Traded: Calculate the number of days since each option was last traded.
- Option Weight: Compute the weight of each option, which is essential for calibrating the Heston parameters.
- Option Price: Assume the price of each option to be the midpoint between the bid and ask prices. Note that this may not be accurate for illiquid options, but these will carry less weight in our analysis due to their spread.
- Maturity: Calculate the maturity date for each option based on its expiration date.
from tqdm import tqdm
from datetime import datetime
import yfinance as yf
from tqdm import tqdm
df = []
for expiration_date in tqdm(expiration_dates, desc='Collecting Option Data'):
# Collect all Calls
_df = SPX.option_chain(expiration_date).calls
_df['DaysSinceLastTraded'] = (pd.to_datetime(_df['lastTradeDate']).dt.tz_localize(None) - pd.to_datetime(datetime.now())).dt.days
_df = _df[_df['DaysSinceLastTraded'] > -14]
_df['Weight'] = 1 / (_df['bid'] - _df['ask'])**2
_df['price'] = (_df['bid'] + _df['ask'])/2
_df['maturity'] = (pd.to_datetime(expiration_date) - pd.to_datetime(datetime.now())).days / 365
_df['Type'] = 'c'
_df['S'] = S0
df.append(_df)
# Collect all Puts
_df = SPX.option_chain(expiration_date).puts
_df['DaysSinceLastTraded'] = (pd.to_datetime(_df['lastTradeDate']).dt.tz_localize(None) - pd.to_datetime(datetime.now())).dt.days
_df = _df[_df['DaysSinceLastTraded'] > -14]
_df['Weight'] = 1 / (_df['bid'] - _df['ask'])**2
_df['price'] = (_df['bid'] + _df['ask'])/2
_df['maturity'] = (pd.to_datetime(expiration_date) - pd.to_datetime(datetime.now())).days / 365
_df['Type'] = 'p'
_df['S'] = S0
df.append(_df)
df = pd.concat(df)
df = df[df['maturity'] > 0]3. Filter the options based on their last trade date to minimize fitting issues related to “irrelevant” options. Ideally, having precise details like the time at which an option was sold and the S&P 500’s position at that moment would aid in calculating the Implied Volatility (IV). However, this level of granularity is generally unavailable. Alternative methods for approximating option price, beyond simply taking the midpoint between bid and ask, could be explored.
4. Apply the Nelson-Siegel-Svensson model to estimate the risk-free rate corresponding to each option’s maturity. This model offers a robust way to account for the term structure of interest rates.
volSurfaceLong['rate'] = volSurfaceLong['maturity'].apply(curve_fit) / 100
volSurfaceLong = volSurfaceLong.reset_index()3D Option Chain Plots: Price, IV vs Time vs Strike
For European Options available on the S&P 500, I’ve employed Plotly to visualize the option chain, as it excels in rendering 3D plots.
We will utilize py_vollib_vectorized to compute the Implied Volatility (IV) for these options.
Key Observations:
The following points stand out from the plots:
- Smirk or Smile: A pronounced smirk or smile is evident, particularly in the IV plot.
- Symmetry in IV: For options with the same strike and maturity, the IVs of calls and puts are remarkably close.
- Low IV Near Expiry: Options that are both close to expiration and at-the-money exhibit the lowest IV.
- Skewed IV for Short Expirations: IV tends to be both higher and more skewed for options with shorter expirations.
Modeling Considerations:
These observations need to be encapsulated effectively by any model we decide to employ.
Note on IV Calculation: Given that our data collection and IV calculation methods aren’t flawless, there’s a potential for inaccuracies. For instance, if incorrect inputs are fed into py_vollib_vectorized, it defaults to zero. Given that IV is particularly sensitive for options near expiration, it's crucial to accurately input S0 and the transactional option price.
Calibration of Heston Model Parameters
Now that we have a tool for pricing European options based on specific parameters and inputs, the next step is to identify the optimal parameters of the Heston model that best align with market data. Let’s explore the challenges and nuances involved in this calibration:
- Optimization Metric: Option Price or Implied Volatility (% Annualized) — Critical. Implied Volatility (IV) treats all options equally, unlike option price, which gives more weight to deep in-the-money or out-of-the-money options.
- Optimization Method: Sequential Least Squares Programming optimizer (SLSQP), Simulated Annealing (SA), etc. — High. Your choice here impacts both the quality of the fit and the convergence time. I opted for SLSQP for its efficiency and reliability.
- Error Function: MAE, MSE, MAPE (Mean Absolute Percentage Error) — High. MSE might seem tempting due to its sensitivity to outliers and quicker convergence, but it’s not always the best choice.
- Weight: 1/(spread²),1/(|spread|), 1/(√spread) — Medium. The choice of weighting is discretionary, although, in my experience, options closer to expiration and at-the-money tend to have smaller spreads.
- Initial Heston Parameters and Search Space Limits: Refer to those used in the project — High. The initial parameters should be “reasonable” to enable effective convergence. That’s why others prefer methods like SA, but starting from a strong initial point makes SLSQP more advantageous.
- Tolerance (tol, ftol): I set this at 1×10−51×10−5, which is sufficient for our needs, given the % annualized IV like 30.05%. If the improvement is smaller, the algorithm exits — usually before hitting the max iterations.
- Handling Errors: When dealing with options close to expiration or with high strike prices, the Fourier Transform may struggle. You can either filter out these options or set their IV to zero. I chose the latter to maintain a robust set of parameters.
By leveraging both Fourier Transform (FT) and Numba, calibrating a full option chain of 4,200 calls and puts takes under 6 minutes on a single-core laptop. The results? A minimum weighted RMSE(IV) of 3.4009 and a MAE(IV) of 2.7232 (weighted) and 1.3806 (non-weighted). That’s an average annualized % volatility error of just 1%.
Example of Resulting Fit
Results:
v0=0.0228 ; kappa=0.0807 ; theta=0.0363 ; zeta=1.1760 ; rho=-0.3021
>>Zeros : 1
# -------------------------------- Error for Prices --------------------------------
MAE of Price : 13.6089
WMAE of Price : 1.2354
MAPE of Price : 54.966 %
WMAPE of Price : 77.121 %
Max APE of Price : 4961.760 %
RMSE of Price : 31.5974
WRMSE of Price : 3.9844
# -------------------------------- Error for IV --------------------------------
MAE of IV : 2.7232
WMAE of IV : 1.3806
MAPE of IV : 15.250 %
WMAPE of IV : 6.704 %
Max APE of IV : 230.815 %
RMSE of IV : 4.2373
-> WRMSE of IV: 3.4009
The Weighted Percent Error shows the weak spots in the algorithm’s fit, particularly for short-expiration and deep out/in-the-money options, which are inherently more challenging to price.
What We Expect from the Optimization Algorithm
- Converge to robust parameters.
- Execute efficiently, leveraging prior attempts.
- Avoid getting stuck by exploring various regions of the parameter space.
Code — Implementation
from scipy.optimize import minimize
# Define variables to be used in optimization
_K = volSurface['strike'].to_numpy()
_C = volSurface['price'].to_numpy()
_type = volSurface['Type'].to_numpy(dtype=str)
_T = volSurface['maturity'].to_numpy()
_r = volSurface['rate'].to_numpy()
_S = np.full_like(_K, S0)
_IV = volSurface['IV'].to_numpy()
_Weight = volSurface['Weight'].to_numpy()
# Parameters to search for in the simualation: --------------------------------
params = {
"v0": {"x0": 0.002874, "lbub": [1e-3, 1.2]},
"kappa": {"x0": 1.6891, "lbub": [1e-3, 10]},
"theta": {"x0": 0.0190, "lbub": [1e-3, 1.2]},
"zeta": {"x0": 3.7472, "lbub": [1e-2, 4]},
"rho": {"x0": -0.2731, "lbub": [-1, 1]},
}
x0 = [param["x0"] for key, param in params.items()]
bnds = [param["lbub"] for key, param in params.items()]
def SqErr(x):
v0, kappa, theta, zeta, rho = [param for param in x]
# Prices Using Heston Model
Price_Heston = get_resutls_array_Heston(volSurface, v0=v0, kappa=kappa, theta=theta, zeta=zeta, rho=rho, N=1_012, z=24)
# Inverse approximation of implied Volatilities
IV_Heston = implied_volatility(price=Price_Heston, S=_S, K=_K, t=_T, r=_r, flag=_type)
# Safty Feature: Exclude the IV caluclations which are not definable.
diff = IV_Heston - _IV
idx = np.isnan(diff) | np.isinf(diff)
diff[idx] = 0 - _IV[idx]
IV_Heston[idx] = 0
# RMSE
rmse = sqrt(np.mean((diff*100)**2 * _Weight))
zeros = int(np.where(IV_Heston==0, 1,0).sum())
wmae = np.mean(np.abs(diff*100) * _Weight)
print(f">>v0={v0:7.4}; kappa={kappa:5.4f}; theta={theta:5.4f}; zeta={zeta:5.4f}; rho={rho:6.4f} | WMAE(IV): {wmae:5e} | Nulls: {(idx).sum()}/{idx.shape[0]} | Zeros: {zeros}/{idx.shape[0]} | WRMSE(IV): {rmse:5e}")
return rmse
# Perform scipy minimize to find the best parameters in the Heston Model which match the observed values
result = minimize(SqErr, x0, tol = 1e-5, method='SLSQP', options={'maxiter': 80 , 'ftol': 1e-5, 'disp':True}, bounds=bnds, jac='3-point')
v0, kappa, theta, zeta, rho = [param for param in result.x]
""" Example output:
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606438e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789558e+00
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606397e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789579e+00
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606400e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789582e+00
>>v0=0.01459; kappa=1.0330; theta=0.0268; zeta=3.5362; rho=-0.2757 | WMAE(IV): 1.606403e+00 | Nulls: 0/4283 | Zeros: 6/4283 | WRMSE(IV): 3.789586e+00
...
Current function value: 3.4043308974510307
Iterations: 41
Function evaluations: 486
Gradient evaluations: 41
Time: 5:52minsVisualizing the Fit
Predicted prices diverge more for options with longer expirations, primarily because the algorithm focuses on fitting options closer to expiration due to their smaller spreads. The same logic applies when considering IV, making it a logical fit metric for the Heston Model.
Heston IV surface
The ultimate goal here is to construct a 3D surface plot for the Heston Implied Volatility (IV). By generating a 2D array containing data from 15,000 options, we can visualize the characteristics of this IV surface.
In a previous discussion about the Observed IV surface, the following properties were highlighted:
- Presence of a smirk or smile — This feature is apparent in our Heston model.
- Lowest IV for options that are both at-the-money and close to expiration — This holds true in the Heston model as well.
- Increased and skewed IV for short-term expirations — Again, the model corroborates this.
With parameters:
Heston: v0=0.001 ; kappa=1.6483 ; theta=0.0261 ; zeta=4.0000 ; rho=-0.0832
S(0)=4405.71 USD
Addressing Numerical Instabilities
One challenge that arises is the model’s failure to accurately calculate option values under certain parameters. This is particularly evident in the following plot, which focuses on options with maturities of less than one year.
“The Fourier Trap”: Understanding the Pitfalls of Option Pricing with FT
To shed light on these intricacies, I revisited the previous code to plot the Probability Density Function (PDF) of the Heston model. The parameters remain consistent with the implied volatility (IV) surface plot discussed earlier. The PDF corresponds to an option set to expire in one day, i.e., T=1/365.
It’s readily apparent that the plot with fewer terms (N=1,000) demonstrates more instability than its counterpart with more terms (N=10,000). So, what gives? Revisiting Part 1, we recall that as N approaches infinity, the Fourier series converges to the original function. The caveat, though, is that until that point, we’re working with approximations. Since our approximation leverages sinusoidal waves, an original function that’s close to zero becomes notoriously hard to approximate, especially for extreme z-scores.
def Fourier_cosine_char_function(char_function, x:np.array, a, b, N:int):
A = lambda n : 2/(b-a) * np.real(char_function((n*pi)/(b-a))*exp((-1j*n*pi*a)/(b-a)) )
fx = .5*A(0)
for n in np.arange(1, N+1):
fx += A(n) * cos(n*pi * (x-a)/(b-a))
return fxThe Fourier Trap
One might suspect that these inconsistencies have something to do with the Heston PDF’s peculiarities. However, that assumption is not corroborated.
It’s evident that regardless of the number of terms used (N), the approximation reaches a point where it ceases to converge to the original function. This pitfall is what I’ve coined as the “Fourier Trap”. My assumption is that a similar phenomenon occurs in the pricing of Put options.
Conclusion: The Fourier Trap
The plots above explains why the numerical instabilities occure in the left side of the calculation. The blue values are negative, which is impossible in a CDF. Similarly in the heston model the area which is missing is because there are negative values in the Heston-CDF which results to the price of the derivative to equal to 0.
Addressing the Problem
- Boosting the Term Count (N): Online sources generally skimp on terms while calculating options, focusing mainly on “convenient” maturities and near-the-money options. Given that adding terms in Fourier Series is computationally cheap, there’s no harm in using 1,024 terms. However past that point, benefits appear to plateau. For instance, upping the term count from 1,024 to 10,000 shows negligible improvements in reducing numerical instabilities but increases computation time by a factor of 10.
- Choosing Optimal Cumulants (c1, c2): While many articles overlook integral limits, Article [3] (2018) is an exception. Many methodologies, like in Article [4] (2008), seem to lack rigor. First off, the cumulants differ from those in [3]. Second, the proposed range of ±12*std falls short for short maturities or deep in/out-of-the-money options. Finally, they don’t thoroughly analyze the different Heston parameters. After extensive testing, using cumulants from Article [3] in conjunction with z=±24 yielded the lowest Mean Squared Error (MSE).
- Setting an Appropriate Range (z): Given the high kurtosis and skewness in the Heston Probability Density Function (PDF), especially when rho is significantly less than zero, it’s crucial to employ a broad enough range. A value of z=±24 seems optimal, based on the lowest MSE results
- Penalising Heston Parameters which result to numerical instabilities during fitting: During fitting we can artificially increase the MSE when the values cause numerical instabilies. That way we can ensure that we are going to converge to a set of values which are more stable. This has the drawback that the resulting fit might not be as good. But its hard to quantify “the penalty” of not being able to price some of the options in the chain. Without any penalty, less than 10/~6000 options are unable to be priced. (see point 5 and 6 bellow) The process of penalising the fit for numerical instabilities is another story its self, but in its simplest form you can test “hard to price” options and see how many numerical instabilities there are and penalize accordingly. Or a second option would be to increase MSE very much if there are any options which cannot be priced.
Let’s dive into some revealing plots:
Insights on the Fourier Trap in the Heston Model
- A larger z helps in pricing deep in/out-of-the-money options.
- The higherthe z value, the less stable the results.
- Increasing z demands a higher N for stable outcomes.
- A larger z enhances the value of the density at the tails, possibly indicating greater accuracy.
- The left tail tends to be more susceptible to breakdown, particularly when rho is negative (e.g., rho=-0.5).
- Numerical stability varies based on specific Heston parameters:
- Higher v0 relative to theta(θ) yields more stable results.
- A higher kappa(κ) improves stability.
- Elevated zeta(ζ) leads to instability.
- Rho(ρ) doesn’t noticeably influence result stability.
Further Work
- Objective 1: Develop a function to optimize the variables zz and NN based on the given Heston model parameters. The aim is to minimize numerical instabilities across a variety of option contracts.
- Objective 2: Dive deep into the underlying mathematical framework to identify potential modifications that could solve the Fourier Trap issue. An immediate area of investigation is the simultaneous use of sine and cosine functions, as their interaction may nullify each other, potentially yielding more stable, near-zero values.