Carr-Madan FFT Option Pricing: Formula and Python

In 1999 Peter Carr and Dilip Madan published Option Valuation Using the Fast Fourier Transform. Their input was a characteristic function. Their output was a grid of European call prices obtained in $O(N\log N)$ operations.

The FFT is not the central mathematical idea. Damping is. A call price is not integrable over all log-strikes, so its ordinary Fourier transform does not exist. Carr and Madan multiplied it by an exponential, transformed the result and only then used the FFT as a fast quadrature engine.

This method is useful for the Heston model, Merton jump diffusion and other models with tractable characteristic functions.

Transform convention

Let $X_T=\log S_T$ and use

$$\phi_T(u)=\mathbb E^{\mathbb Q}[e^{iuX_T}].$$

Write $k=\log K$ and let $C_T(k)$ be the discounted call price,

$$C_T(k)=e^{-rT}\mathbb E[(S_T-e^k)^+].$$

As $k\to-\infty$, this tends to the discounted forward value rather than zero. Define the damped price

$$c_T(k)=e^{\alpha k}C_T(k),\qquad \alpha>0.$$

The Fourier transform $\psi_T(v)=\int_{-\infty}^{\infty}e^{ivk}c_T(k)\,dk$ is

$$\boxed{ \psi_T(v)= \frac{e^{-rT}\phi_T(v-i(\alpha+1))} {\alpha^2+\alpha-v^2+i(2\alpha+1)v} }.$$

With this sign convention, inversion gives

$$C_T(k)=\frac{e^{-\alpha k}}{\pi} \int_0^\infty \operatorname{Re}\!\left[e^{-ivk}\psi_T(v)\right]dv.$$

Mixing a return characteristic function with a log-price characteristic function shifts every strike. Mixing discounted and undiscounted calls adds or removes $e^{-rT}$. Both bugs produce plausible curves, so unit tests matter.

The denominator comes from integrating the payoff over log-strike before taking the expectation. For a fixed terminal log-price $x$, the damped payoff contributes only for $k<x$:

$$\int_{-\infty}^{x} e^{(\alpha+iv)k}(e^x-e^k)\,dk =\frac{e^{(\alpha+1+iv)x}} {(\alpha+iv)(\alpha+1+iv)}.$$

Expanding the two factors in the denominator gives

$$(\alpha+iv)(\alpha+1+iv) =\alpha^2+\alpha-v^2+i(2\alpha+1)v.$$

Taking the discounted expectation of the numerator evaluates the characteristic function at the complex argument $v-i(\alpha+1)$. This derivation makes the convention checks mechanical. If the transform uses $e^{-ivk}$ instead, both the inversion phase and complex shift change. Altering only one sign reflects or shifts the price grid.

The damping condition

$\alpha>0$ controls the low-strike tail, but positivity alone is not sufficient. The complex characteristic function must exist at $v-i(\alpha+1)$, which requires the moment

$$\mathbb E[S_T^{\alpha+1}]<\infty.$$

Black-Scholes admits every positive moment. Lévy and stochastic-volatility models usually have a finite moment strip. Choosing $\alpha$ outside it makes the transform undefined even if floating-point code returns a number.

There is no universal best $\alpha$. Values around $1.5$ are common, not guaranteed. Price at several admissible values and require agreement. If results move materially, investigate the moment strip, integration range and characteristic-function branch before tuning $\alpha$ further.

FFT grid

Set

$$v_j=j\eta,\qquad k_n=-b+n\lambda, \qquad j,n=0,\ldots,N-1,$$

with

$$\lambda\eta=\frac{2\pi}{N}, \qquad b=\frac{N\lambda}{2}.$$

The grid trade-off is coupled. Increasing $\eta$ extends the Fourier cut-off and makes the log-strike spacing $\lambda$ finer, but coarsens the integration step. Increasing $N$ is the cleanest way to enlarge both resolutions, at extra cost.

For the code below, $N=4096$ and $\eta=0.25$ imply $\lambda=0.00613592$ and $b=12.56637$. Because $k$ is absolute log-strike, the resulting strike interval is approximately $[e^{-b},e^b]=[3.49\times10^{-6},286751]$. That range is far wider than a typical desk needs. It is a reminder that a default grid can spend most output points on irrelevant strikes while the integration cut-off $v_{\max}\approx1023.75$ may be unnecessarily large for a smooth characteristic function. Grid design should start from the required strike interval and accuracy, not from a familiar tuple.

Carr and Madan use modified Simpson coefficients. In zero-based Python indexing they are

$$w_j=\frac{\eta}{3} \left(3-(-1)^j-\mathbf 1_{\{j=0\}}\right).$$

This starts $\eta/3,4\eta/3,2\eta/3,\ldots$. Writing $3+(-1)^j$ is an indexing error: it starts with weight $\eta$ and materially overprices calls.

The discrete inversion is

$$C_T(k_n)\approx \frac{e^{-\alpha k_n}}{\pi} \operatorname{Re}\!\left[ \sum_{j=0}^{N-1} e^{-2\pi ijn/N} e^{ibv_j}\psi_T(v_j)w_j \right].$$

Simpson's rule improves quadrature order for a smooth integrand. It does not remove Fourier truncation, periodic aliasing, interpolation error or cancellation in deep wings.

Working Python

This implementation accepts any log-price characteristic function and includes a Black-Scholes test.

import numpy as np
from scipy.stats import norm

def carr_madan_fft(cf, T, r, alpha=1.5, N=4096, eta=0.25):
    if N & (N - 1):
        raise ValueError("N must be a power of two")
    dk = 2*np.pi / (N*eta)
    b = 0.5*N*dk
    j, v = np.arange(N), eta*np.arange(N)
    k = -b + dk*np.arange(N)
    psi = (np.exp(-r*T) * cf(v - 1j*(alpha + 1))
           / (alpha**2 + alpha - v**2 + 1j*(2*alpha + 1)*v))
    simpson = (3.0 - (-1.0)**j - (j == 0)) / 3.0
    fft_input = np.exp(1j*b*v) * psi * eta * simpson
    calls = np.exp(-alpha*k) * np.fft.fft(fft_input).real / np.pi
    return k, calls

def bs_cf(u, S0, T, r, q, vol):
    mean = np.log(S0) + (r - q - 0.5*vol**2)*T
    return np.exp(1j*u*mean - 0.5*vol**2*u**2*T)

S0, T, r, q, vol = 100, 0.5, 0.03, 0.0, 0.20
cf = lambda u: bs_cf(u, S0, T, r, q, vol)
log_strikes, calls = carr_madan_fft(cf, T, r)
call_at_100 = np.interp(np.log(100), log_strikes, calls)

The output grid is in log-strike, not strike. Interpolate locally and include interpolation in convergence tests. A shape-preserving method is safer than an unconstrained cubic spline if the grid is not already fine.

For the parameters above, the FFT gives $6.37209$ at $K=100$ against the Black-Scholes value $6.37103$. At $K=90$ it gives $12.79979$ against $12.79930$. Changing $N$, $\eta$ and $\alpha$ should drive those gaps down. The Black-Scholes closed form is the first regression test every characteristic-function engine should pass.

Error budget and grid diagnostics

An FFT price contains several distinct numerical errors. Quadrature error comes from replacing the integral with spaced samples. Truncation error comes from stopping at $v_{\max}=(N-1)\eta$. Aliasing comes from the periodic extension implicit in the discrete Fourier transform. Strike interpolation adds another error when the requested $K$ is not on the log grid. Floating-point cancellation matters for very small option values after multiplication by $e^{-\alpha k}$.

These errors react differently to parameters. Reducing $\eta$ improves low-frequency quadrature but, with fixed $N$, reduces $v_{\max}$ and increases log-strike spacing. Increasing $N$ while holding $\eta$ fixed extends the integration range and refines strike spacing together. Changing $\alpha$ alters integrand decay and numerical scale but does not fix a characteristic-function branch error. A convergence study should vary one quantity at a time, then vary coupled choices that preserve either $v_{\max}$ or the desired strike interval.

A useful production harness prices a small panel by both FFT and adaptive one-dimensional quadrature using the same $\psi_T$. Include near-ATM, moderately out-of-the-money and deep-wing strikes at short and long maturities. Agreement at the centre alone is weak evidence because interpolation and cancellation are easiest there. For each panel point, record absolute price error, implied-volatility error and error relative to bid-ask spread.

If oscillations appear across strikes, inspect the integrand before increasing $N$. Plot or sample the real and imaginary parts of $\psi_T(v)$ and their magnitudes. Failure to decay indicates an invalid damping moment or a transform formula error. Smooth decay followed by a hard cut-off suggests truncation. Alternating noise near the cut-off indicates inadequate resolution. A discontinuity in complex phase often points to a logarithm branch, especially in Heston implementations.

FFT output should not be repaired silently. Clipping negative calls to zero can conceal oscillation, and enforcing monotonicity after pricing can distort calibration gradients. Reject a grid that violates material no-arbitrage bounds in the strike region used for fitting. Tiny violations far outside the requested range may be excluded, but the exclusion and tolerance should be explicit.

Practitioners also need to choose the right transform method for the workload. During calibration, one FFT can provide many strikes for one maturity, but quoted strikes are irregular and the optimiser may need only a few dozen values. Direct quadrature can then be competitive because it avoids generating and interpolating thousands of unused prices. Fractional FFT variants relax the reciprocal grid constraint. The COS method can converge rapidly when reliable cumulant-based density bounds are available. Benchmark total calibration time and accuracy rather than assuming the asymptotic $O(N\log N)$ label settles the decision.

Production checks

Characteristic-function convention. Verify $\phi_T(0)=1$, $\phi_T(-i)=\mathbb E[S_T]=S_0e^{(r-q)T}$ and conjugate symmetry on the real axis. The $-i$ check also catches a missing martingale correction.

Static bounds. Calls should satisfy

$$\max(S_0e^{-qT}-Ke^{-rT},0) \leq C(K,T)\leq S_0e^{-qT},$$

decrease with strike and remain convex. Small negative wing prices are numerical error, not a model feature.

Heston branches. A naive complex logarithm can jump branches across the integration grid. Albrecher, Mayer, Schoutens and Tistaert's “little Heston trap” formulation reduces this problem; continuity across maturities still deserves a test.

Greeks. Bumping and rerunning an FFT is simple, but interpolation noise contaminates small bumps. Analytic differentiation of the characteristic function, complex-step methods where valid or automatic differentiation can be cleaner.

For one or two strikes, direct adaptive quadrature can be faster and easier to control. FFT wins when a calibration needs many strikes on a regular grid. The COS method of Fang and Oosterlee is another strong option, particularly when density truncation bounds are well controlled.

Jump models can also be solved by partial integro-differential equations. Saying that PDE methods are inapplicable to Lévy processes is too strong; the non-local integral term simply changes the numerical problem.

Limits of the method

The transform prices European payoffs at one maturity. It does not perform early-exercise backwards induction. Bermudan contracts need a recursive Fourier method or simulation such as Longstaff-Schwartz.

Very short maturities and far out-of-the-money prices are demanding because the density is concentrated and the option value is obtained through cancellation. No fixed tuple $(N,\eta,\alpha)$ is reliable for every model, maturity and strike range.

Frequently Asked Questions

What does Carr-Madan FFT price?

It prices a grid of European calls from the risk-neutral characteristic function of terminal log-price, provided the required complex moment exists.

Why is exponential damping needed?

The call price does not vanish as log-strike tends to minus infinity. Multiplication by $e^{\alpha k}$ makes it integrable when $\alpha>0$ and the $(\alpha+1)$ moment is finite.

What are the correct Simpson weights in Python?

With $j=0,\ldots,N-1$, use $\eta[3-(-1)^j-\mathbf1_{\{j=0\}}]/3$. The opposite sign on $(-1)^j$ belongs to a different one-based indexing convention.

How should alpha be selected?

Keep $\alpha+1$ inside the model's finite moment strip, then check price stability over several admissible values. A conventional value such as $1.5$ is only a starting point.

Can Carr-Madan price American options?

Not through this one-period inversion. American and Bermudan exercise requires backwards recursion, a tree or PDE/PIDE, or an optimal-stopping simulation method.

What is the original source?

Use Carr and Madan's 1999 paper, Journal of Computational Finance 2(4), pages 61 to 73. The FFT is fast only after the transform conventions are right.