Pricing Options with Fast Fourier Transform (FFT) — Heston, Variance-Gamma and Black-Scholes (with code)

In this article, I’ll show how you can value options with Fast Fourier Transform (FFT). The concepts and part of the code, that I’ll present here is based on the course “Financial Engineering and Risk Management” by Columbia University. It is offered through Coursera and I highly recommend it.

The FFT is a mathematical technique that efficiently computes the discrete Fourier transform and has been adapted for option pricing due to its computational advantages. By employing FFT in option valuation, you can quickly and accurately determine the prices of various financial derivatives.

First, I’ll show the general concept and how to perform the computations. Then I’ll show 3 pricing models and how to compute option prices based on their characteristic functions. We'll work with 3 models: Black-Scholes, Heston, and Variance Gamma.

FFT technique can be applied also to caps, floors, and swaptions, which can be expressed as simple put and call options, thus enabling fast valuation of these instruments.

Foundations

Even though one could price the option by integration, this is not an efficient method. The better method is FFT pricing.

Of course, we could write our explicit FFT functions to price the options, however, the NumPy module np.fft.fft is much faster, so we’ll use it instead.

There is a set of rules we have to follow to make our calculations correct.
1) log Strike — quite straightforward, you just have to write e.g. Strike = 100 as log(100)
2) For FFT to be efficient you should always use 2^n (vectors, e.g. 2³ = 8)
3) Logarithmic of a stock price process
4) Defined parameters: eta, alpha, n and beta

S0 = 100
K = 80
k = np.log(K)
r = 0.055
q = 0.03

These are example parameters

# Parameters for FFT 
n = 12
N = 2**n
#step-size
eta = 0.25
# damping factor
alpha = 1

lambda_ = (2*np.pi/N)/eta;
beta = np.log(K)

These are the basics. Now, let’s review the models will be working with and implement the FFT pricing.

Black-Scholes

The Black-Scholes model is a fundamental mathematical model used in finance to calculate the theoretical price of European-style options.

• Efficient Markets: The model assumes that financial markets are efficient and there are no arbitrage opportunities. It assumes continuous trading, no transaction costs, and that all information is reflected in the current market price.
• Log-Normal Distribution: The model assumes that the underlying asset’s price follows a log-normal distribution. This assumption, combined with continuous trading, leads to a geometric Brownian motion for the stock price.
• Constant Volatility: It assumes that the volatility of the underlying asset is constant over the option’s life. This assumption is one of the primary limitations of the Black-Scholes model, as it contradicts empirical observations where volatility tends to change over time.
• Risk-Free Interest Rate: The risk-free interest rate is assumed to be known and constant throughout the option’s life.

The Black-Scholes model is a foundational tool in options pricing, providing a theoretical framework to estimate the fair price of European options under certain assumptions. Despite its simplicity and various assumptions, it laid the groundwork for further advancements in financial derivatives and option pricing theory. However, its limitations in capturing changing volatility and other market dynamics have led to the development of more sophisticated models like the Heston model and the Variance Gamma model mentioned later.

The following equations represent the BS process and its characteristic function necessary for FFT calculations:

Heston

The Heston model is a widely used stochastic volatility model in finance that addresses some limitations of the Black-Scholes model by introducing stochastic volatility. The model aims to capture the volatility smile observed in options markets and the tendency for volatility to fluctuate over time.

• Stochastic Volatility: Unlike the Black-Scholes model, where volatility is assumed to be constant, the Heston model incorporates stochasticity in volatility. It assumes that the variance itself follows a stochastic process, typically modeled as a mean-reverting square-root diffusion process (CIR process).
• Mean-Reversion: Heston’s stochastic volatility model includes mean-reverting behavior, where volatility tends to move towards a long-term average level. This mean-reverting property helps to capture the observed behavior of volatility in financial markets.
• Volatility Smile: One of the primary motivations for the Heston model is its ability to produce the volatility smile, an empirical phenomenon where implied volatility varies with strike price and maturity in options markets.
• Correlation between Asset Price and Volatility: The Heston model allows for a correlation parameter (typically denoted by ρ) between the asset price process and the volatility process. This correlation captures the relationship between asset prices and their corresponding volatilities. If the correlation is zero, it implies independence between the two processes.

The Heston model’s introduction of stochastic volatility addresses some of the shortcomings of the Black-Scholes model by providing a more dynamic representation of market behavior. Its ability to capture the volatility smile and account for the tendency of volatility to fluctuate over time makes it a popular choice in options pricing and risk management applications within financial markets.

Variance-Gamma

The Variance Gamma (VG) model is an extension of the Black-Scholes model that incorporates both jumps in asset prices and stochastic volatility

• Incorporation of Jumps: The VG model introduces jump components into the price dynamics, allowing for more realistic behavior in financial markets. These jumps help in capturing sudden and large movements that occur occasionally in asset prices.
• Flexibility in Skewness and Kurtosis: VG allows for greater flexibility in modeling skewness (asymmetry) and kurtosis (fat-tailedness) compared to simpler models like Black-Scholes.
• Stochastic Volatility: Similar to the Heston model, VG incorporates stochastic volatility, allowing volatility itself to vary stochastically over time. This feature provides a better representation of the empirical behavior of volatility in financial markets.

When both skewness and kurtosis parameters are set to zero, the VG model converges to the Black-Scholes model. Additionally, the VG model reduces to the Merton jump-diffusion model when the stochastic volatility term is set to zero.

Its ability to encompass both jumps and stochastic volatility allows for a more realistic representation of market dynamics and helps in addressing some of the limitations of simpler models.

VG characteristic function of the Logarithm of Stock Price:

Code implementation

After reviewing the models, we’ll be working with let’s see the code implementation. The first function “compute_characteristic_function” is just the direct implementation of the formulas above.

The function takes input parameters specific to each model and generates the corresponding characteristic function. The phi variable holds the resulting characteristic function, which is then returned.

def compute_characteristic_function(u, params, S0, r, q, T, model):
    if model == 'BS':
        volatility = params[0]
        drift = np.log(S0) + (r - q - volatility ** 2 / 2) * T
        diffusion = volatility * np.sqrt(T)
        phi = np.exp(1j * drift * u - (diffusion * u) ** 2 / 2)
    elif model == 'Heston':
        kappa = params[0]
        theta = params[1]
        sigma = params[2]
        rho = params[3]
        v0 = params[4]
        
        temp1 = (kappa - 1j * rho * sigma * u)
        g = np.sqrt((sigma ** 2) * (u ** 2 + 1j * u) + temp1 ** 2)
        power_term = 2 * kappa * theta / (sigma ** 2)
        numerator = (kappa * theta * T * temp1) / (sigma ** 2) + 1j * u * T * r + 1j * u * np.log(S0)
        log_denominator = power_term * np.log(np.cosh(g * T / 2) + (temp1 / g) * np.sinh(g * T / 2))
        temp2 = ((u * u + 1j * u) * v0) / (g / np.tanh(g * T / 2) + temp1)
        log_phi = numerator - log_denominator - temp2
        phi = np.exp(log_phi)
    elif model == 'VG':
        sigma = params[0]
        nu = params[1]
        theta = params[2]
        
        if nu == 0:
            mu = np.log(S0) + (r - q - theta - 0.5 * sigma ** 2) * T
            phi = np.exp(1j * u * mu) * np.exp((1j * theta * u - 0.5 * sigma ** 2 * u ** 2) * T)
        else:
            mu = np.log(S0) + (r - q + np.log(1 - theta * nu - 0.5 * sigma ** 2 * nu) / nu) * T
            phi = np.exp(1j * u * mu) * ((1 - 1j * nu * theta * u + 0.5 * nu * sigma ** 2 * u ** 2) ** (-T / nu))
    return phi

The function “genericFFT” essentially uses FFT to efficiently compute option prices by transforming characteristic function values and strike prices into a suitable space for pricing calculation, enabling faster and more accurate option valuation.

def calculate_fft_values(params, S0, K, r, q, T, alpha, eta, n, model):
    """
    S0: Initial asset price.
    K: Strike price of the option.
    r: Risk-free interest rate.
    q: Continuous dividend yield.
    T: Time to maturity.
    alpha: Constant used in the Carr-Madan method for damping.
    eta: Step size for numerical integration in FFT.
    n: Exponential factor to determine the number of steps in the FFT.
    """
    N = 2 ** n
    delta = (2 * np.pi / N) / eta
    beta = np.log(K)
    km_values = np.zeros(N)
    x_values = np.zeros(N)
    # discount factor
    discount_factor = np.exp(-r * T)

    nuJ = np.arange(N) * eta
    psi_nuJ = generic_characteristic_function(nuJ - (alpha + 1) * 1j, params, S0, r, q, T, model) / (
                (alpha + 1j * nuJ) * (alpha + 1 + 1j * nuJ))

    km_values = beta + delta * np.arange(N)
    w_values = eta * np.ones(N)
    w_values[0] = eta / 2

    x_values = np.exp(-1j * beta * nuJ) * discount_factor * psi_nuJ * w_values
    y_values = np.fft.fft(x_values)

    cT_km_values = np.zeros(N)
    multiplier = np.exp(-alpha * km_values) / np.pi

    cT_km_values = multiplier * np.real(y_values)

    return km_values, cT_km_values

And now code for pricing puts. It initializes a matrix put_matrix to store results. Each row represents a unique combination of parameters, and columns represent eta, n, alpha, and put price.

def calculate_put_prices(params, S0, K, r, q, T, model, alpha_values, eta_values, n_values):
    num_prices = len(eta_values) * len(n_values) * len(alpha_values)
    # Columns correspond to eta, n, alpha, and put price
    put_matrix = np.zeros([num_prices, 4])
    i = 0
    for eta in eta_values:
        for n in n_values:
            for alpha in alpha_values:
                # Pricing via FFT
                put = 0 
                k_values, option_values = genericFFT(params, S0, K, r, q, T, alpha, eta, n, model)
                put = option_values[0]  # Considering only one strike
                put_matrix[i] = np.array([eta, n, alpha, put])
                i += 1
    return put_matrix

Now, let’s use the code and price some options! We’ll start with Black Scholes

mod = 'BS'
sig = 0.3
params = [sig]

def print_put_prices(params, S0, K, r, q, T, model, alpha_values, eta_values, n_values):
    put_matrix = calculate_put_prices(params, S0, K, r, q, T, model, alpha_values, eta_values, n_values)
    num_prices = put_matrix.shape[0]
    
    print('Model = ' + model)
    print('eta\tN\talpha\tput')
    
    for i in range(num_prices):
        print('%.2f\t2^%i\t%.2f\t%.4f' % (put_matrix[i, 0], put_matrix[i, 1], put_matrix[i, 2], put_matrix[i, 3]))

The output:
— — — — — — — — — — —
Model = BS
eta N alpha put
0.10 2⁶ -1.01 89.7431
0.10 2⁶ -1.25 2.7396
0.10 2⁶ -1.50 2.7569
0.10 2⁶ -1.75 2.7701
0.10 2⁶ -2.00 2.7789
0.10 2⁶ -5.00 2.6727
0.10 2¹⁰ -1.01 89.7316
0.10 2¹⁰ -1.25 2.7080
0.10 2¹⁰ -1.50 2.7080
0.10 2¹⁰ -1.75 2.7080
0.10 2¹⁰ -2.00 2.7080
0.10 2¹⁰ -5.00 2.7080
0.25 2⁶ -1.01 269.0367
0.25 2⁶ -1.25 2.8504
0.25 2⁶ -1.50 2.7083
0.25 2⁶ -1.75 2.7080
0.25 2⁶ -2.00 2.7080
0.25 2⁶ -5.00 2.7080
0.25 2¹⁰ -1.01 269.0367
0.25 2¹⁰ -1.25 2.8504
0.25 2¹⁰ -1.50 2.7083
0.25 2¹⁰ -1.75 2.7080
0.25 2¹⁰ -2.00 2.7080
0.25 2¹⁰ -5.00 2.7080

Heston:

mod = 'Heston'
kappa = 2.
theta = 0.05
lda = 0.3
rho = -0.7

v0 = 0.04
params = [kappa, theta, lda, rho, v0]
heston_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)

print('Model = ' + mod)
print('eta\tN\talpha\tput')
for i in range(num_prices):
    print('%.2f\t2^%i\t%.2f\t%.4f' % (heston_matrix[i,0], heston_matrix[i,1], heston_matrix[i,2], heston_matrix[i,3]))

The output:
— — — — — — — — — — -
Model = Heston
eta N alpha put
0.10 2⁶ -1.01 88.0744
0.10 2⁶ -1.25 1.1102
0.10 2⁶ -1.50 1.1666
0.10 2⁶ -1.75 1.2167
0.10 2⁶ -2.00 1.2605
0.10 2⁶ -5.00 1.3979
0.10 2¹⁰ -1.01 88.3648
0.10 2¹⁰ -1.25 1.3412
0.10 2¹⁰ -1.50 1.3412
0.10 2¹⁰ -1.75 1.3412
0.10 2¹⁰ -2.00 1.3412
0.10 2¹⁰ -5.00 1.3412
0.25 2⁶ -1.01 267.6738
0.25 2⁶ -1.25 1.4873
0.25 2⁶ -1.50 1.3450
0.25 2⁶ -1.75 1.3445
0.25 2⁶ -2.00 1.3442
0.25 2⁶ -5.00 1.3415
0.25 2¹⁰ -1.01 267.6698
0.25 2¹⁰ -1.25 1.4835
0.25 2¹⁰ -1.50 1.3414
0.25 2¹⁰ -1.75 1.3412
0.25 2¹⁰ -2.00 1.3412
0.25 2¹⁰ -5.00 1.3412

And Variance-Gamma:

# model-specific parameters
mod = 'VG'
sigma = 0.3
nu = 0.5
theta = -0.4
params = [sigma, nu, theta]

vg_matrix = price_all_puts(params, S0, K, r, q, T, mod, alpha_vec, eta_vec, n_vec)
print('Model = ' + mod)
print('eta\tN\talpha\tput')
for i in range(num_prices):
    print('%.2f\t2^%i\t%.2f\t%.4f' % (vg_matrix[i,0], vg_matrix[i,1], vg_matrix[i,2], vg_matrix[i,3]))

The output:
— — — — — — — — — — -
Model = VG
eta N alpha put
0.10 2⁶ -1.01 92.2103
0.10 2⁶ -1.25 5.2047
0.10 2⁶ -1.50 5.2230
0.10 2⁶ -1.75 5.2401
0.10 2⁶ -2.00 5.2556
0.10 2⁶ -5.00 0.0030
0.10 2¹⁰ -1.01 92.3373
0.10 2¹⁰ -1.25 5.3137
0.10 2¹⁰ -1.50 5.3137
0.10 2¹⁰ -1.75 5.3137
0.10 2¹⁰ -2.00 5.3137
0.10 2¹⁰ -5.00 0.0000
0.25 2⁶ -1.01 271.6399
0.25 2⁶ -1.25 5.4539
0.25 2⁶ -1.50 5.3121
0.25 2⁶ -1.75 5.3122
0.25 2⁶ -2.00 5.3124
0.25 2⁶ -5.00 0.0019
0.25 2¹⁰ -1.01 271.6423
0.25 2¹⁰ -1.25 5.4560
0.25 2¹⁰ -1.50 5.3139
0.25 2¹⁰ -1.75 5.3137
0.25 2¹⁰ -2.00 5.3137
0.25 2¹⁰ -5.00 0.0020

Interpretation of the results

• eta: Represents the step size for numerical integration in the FFT method. Smaller values typically provide more precise results but require more computational resources.
• n: Indicates the exponential factor used in determining the number of steps in the FFT. A larger n implies more precise calculations but increased computational complexity.
• alpha: This is a constant used in the Carr-Madan method for damping. It affects the smoothness of the FFT-generated option price curve.
• put: Denotes the computed put option price.

By observing the rows for different combinations of eta, n, and alpha, you can see how changing these parameters influences the resulting put option prices.

The results should not be interpreted without context! The reason for comparing the values is to compare them with the analytical solutions and determine good parameter choice for further pricing.

Conclusion

In summary, the presented analysis offers a comprehensive perspective on the influence of diverse parameters on put option prices within the VG model. This scrutiny underscores the model’s sensitivity to varying parameter inputs, crucial for ensuring precision and stability in pricing.

The exploration into FFT-based option pricing provides a foundational comprehension, highlighting the pivotal role of parameter sensitivity. This knowledge stands to aid quantitative finance practitioners in efficiently valuing derivatives. However, it also underscores the necessity for exacting parameter selections to ensure pricing accuracy in financial valuation practices.

References:

• An Introduction to the Mathematics of Financial Derivatives — Ali Hirsa, Salih A. Neftci
• Financial Engineering and Risk Management Specialization; Computational Methods in Pricing and Model Calibration