FFT and Carr-Madan Pricing

Setup

Motivation

For many models of practical importance — Heston, Variance Gamma, CGMY, NIG — the dynamics of $\ln S_T$ under the risk-neutral measure are not log-normal, and closed-form option prices do not exist. However, the characteristic function of the log-price is often available in closed form.

The Carr-Madan (1999) method prices European options by Fourier-inverting the characteristic function. Its efficiency comes from evaluating option prices at all strikes simultaneously using a single Fast Fourier Transform (FFT) pass — pricing at $N$ strikes for the cost of a single $O(N \log N)$ FFT, compared to $O(N)$ separate numerical integrations.

Characteristic Function

Let $x_T = \ln(S_T / S_0)$ be the log-return. Under the risk-neutral measure $\mathbb{Q}$:

$\varphi_T(u) = \mathbb{E}^{\mathbb{Q}}\!\left[e^{iu x_T}\right], \qquad u \in \mathbb{R}.$

For GBM: $\varphi_T(u) = \exp\!\left(iu(r - \frac{1}{2}\sigma^2)T - \frac{1}{2}u^2\sigma^2 T\right)$ (a known Gaussian characteristic function). For Heston: see the preceding module. For Lévy models (VG, NIG, CGMY): the characteristic function is an exponential of the Lévy exponent, available analytically.

Key property. $\varphi_T(u)$ is a bounded, uniformly continuous function of $u$. It decays as $|u| \to \infty$ whenever $S_T$ has a sufficiently smooth density.

Notation

• $k = \ln K$: log-strike.
• $F = S_0 e^{rT}$: forward price.
• $C_T(k) = e^{-rT}\mathbb{E}^{\mathbb{Q}}[(S_T - e^k)^+]$: call price as a function of log-strike.
• $\varphi_T(u)$: characteristic function of $x_T = \ln(S_T/S_0)$, so $\mathbb{E}^{\mathbb{Q}}[e^{iux_T}] = \varphi_T(u)$.

The Carr-Madan Formula

Problem with the Naive Approach

The undiscounted call payoff $(S_T - K)^+ = (e^{x_T + \ln S_0} - e^k)^+$ is not square-integrable in $k$: as $k \to -\infty$, $C_T(k) \to S_0 - e^{k-rT} \to S_0$ (a constant), so $C_T \notin L^2(\mathbb{R})$. The Fourier transform does not exist in the classical sense.

Dampening

Multiply $C_T(k)$ by an exponential damping factor $e^{\alpha k}$ for some $\alpha > 0$:

$z_T(k) = e^{\alpha k} C_T(k).$

For $\alpha > 0$ and sufficiently small (specifically $\alpha < \bar{\alpha}$ where $\mathbb{E}[S_T^{1+\alpha}] < \infty$): as $k \to -\infty$, $z_T(k) = e^{\alpha k} \cdot (S_0 - e^{k-rT} + O(e^{2k})) \to 0$. As $k \to +\infty$, $z_T(k) = e^{\alpha k} \cdot O(e^{-k}) = O(e^{(\alpha-1)k}) \to 0$ for $\alpha < 1$. So $z_T \in L^2(\mathbb{R})$ and its Fourier transform exists.

Fourier Transform of the Modified Price

Define $\Psi_T(v) = \int_{-\infty}^{\infty} e^{ivk} z_T(k)\, dk$. Substituting $z_T(k) = e^{\alpha k + ivk} e^{-rT}\mathbb{E}^{\mathbb{Q}}[(e^{x_T + \ln S_0} - e^k)^+]$ and exchanging expectation and integration (Fubini, justified under the moment condition):

$\Psi_T(v) = e^{-rT} \int_{-\infty}^{\infty} e^{(\alpha + iv)k} \mathbb{E}^{\mathbb{Q}}\!\left[(S_0 e^{x_T} - e^k)^+\right] dk.$

Switching the order of integration (integrating over $k$ from $-\infty$ to $x_T + \ln S_0$ where the payoff is positive):

$\Psi_T(v) = \frac{e^{-rT} \varphi_T(v - (\alpha + 1)i)}{\alpha^2 + \alpha - v^2 + i(2\alpha + 1)v}.$

Derivation. Inside the expectation, $\int_{-\infty}^{x_T + \ln S_0} e^{(\alpha+iv)k}(S_0 e^{x_T} - e^k)\,dk$ evaluates to a ratio involving $e^{(\alpha+1+iv)(x_T + \ln S_0)}$. Taking the expectation introduces $\mathbb{E}[e^{(\alpha+1+iv)x_T}] = \varphi_T(-i(\alpha+1) + v)$. Collecting terms yields the formula above.

Inversion

The call price is recovered by Fourier inversion:

$C_T(k) = e^{-\alpha k} z_T(k) = \frac{e^{-\alpha k}}{\pi} \int_0^\infty \mathrm{Re}\!\left[e^{-ivk} \Psi_T(v)\right] dv,$

using the symmetry of the real part of the integrand (since the imaginary part integrates to zero for real $C_T$):

$\boxed{C_T(k) = \frac{e^{-\alpha k}}{\pi} \int_0^{\infty} e^{-ivk} \Psi_T(v)\, dv \quad \text{(real part implied)}.}$

This is exact for any model with a known characteristic function $\varphi_T$, provided the dampening condition $\mathbb{E}^{\mathbb{Q}}[S_T^{1+\alpha}] < \infty$ holds.

FFT Implementation

Discretisation of the Integral

Approximate the integral by a truncated sum over frequencies $v_j = \eta j$ for $j = 0, 1, \ldots, N-1$, with spacing $\eta > 0$ and upper limit $v_{\max} = \eta(N-1)$:

$C_T(k) \approx \frac{e^{-\alpha k}}{\pi} \sum_{j=0}^{N-1} e^{-iv_j k} \Psi_T(v_j)\, \eta \cdot w_j,$

where $w_j$ are quadrature weights (e.g., trapezoidal: $w_0 = 1/2$, $w_j = 1$ for $j > 0$; or Simpson's rule alternating $1/3, 4/3, 2/3, 4/3, \ldots$).

Evaluate at a grid of log-strikes $k_u = -b + \lambda u$ for $u = 0, 1, \ldots, N-1$, where $b = N\lambda/2$ centres the grid near the money:

$C_T(k_u) \approx \frac{e^{-\alpha k_u}}{\pi} \sum_{j=0}^{N-1} e^{-iv_j k_u} \Psi_T(v_j)\, \eta\, w_j.$

Substituting $v_j = \eta j$ and $k_u = -b + \lambda u$:

$e^{-iv_j k_u} = e^{-i\eta j(-b + \lambda u)} = e^{i\eta j b} \cdot e^{-i\eta\lambda j u}.$

The FFT Condition

For the sum $\sum_{j=0}^{N-1} e^{-i\eta\lambda ju} f_j$ to be computable by FFT, we need:

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

With $N$ a power of 2, the Cooley-Tukey FFT computes all $N$ values simultaneously in $O(N \log N)$ operations. This is the Nyquist-Shannon link between the frequency spacing $\eta$ and the log-strike spacing $\lambda$.

The full algorithm:

1. Choose $N = 2^n$, upper truncation $v_{\max} = \eta N$, frequency spacing $\eta$.
2. Set log-strike spacing $\lambda = 2\pi / (\eta N)$ and grid centre $b = N\lambda/2$.
3. Compute: $x_j = e^{i v_j b} \Psi_T(v_j) \eta w_j$ for $j = 0, \ldots, N-1$.
4. Apply FFT: $y_u = \sum_{j=0}^{N-1} e^{-i(2\pi/N)ju} x_j$.
5. Recover prices: $C_T(k_u) = e^{-\alpha k_u} \cdot \mathrm{Re}(y_u) / \pi$.

Total cost: one characteristic function evaluation at $N$ frequencies + one FFT = $O(N \log N)$. Compare to $N$ separate numerical integrations at cost $O(N \cdot M_{\mathrm{quad}})$ where $M_{\mathrm{quad}}$ is the number of quadrature points per integration.

Grid Design and Parameter Choice

Frequency Grid ($\eta$)

$\eta$ controls the upper truncation $v_{\max} = \eta N$ and the log-strike spacing. Increasing $\eta$ extends the integration to higher frequencies (important for models with heavy tails or discontinuities) but coarsens the log-strike grid. A typical value: $\eta = 0.25$ with $N = 1024$ gives $\lambda = 2\pi/(0.25 \times 1024) \approx 0.024$ (log-strike spacing of 2.4%, i.e., roughly every 2.4% in relative strike).

Dampening Parameter $\alpha$

The dampening condition requires $\mathbb{E}^{\mathbb{Q}}[S_T^{1+\alpha}] < \infty$. For GBM this is always satisfied. For Heston, the moment condition can fail for large $\alpha$ when the Feller condition is violated (large $\xi$). The standard choice $\alpha = 1.5$ works well for most equity models. The integrand $\Psi_T(v)$ must decay as $|v| \to \infty$ — if $\alpha$ is too large, the denominator factor $\alpha^2 + \alpha - v^2$ changes sign at $v = \sqrt{\alpha^2 + \alpha}$, creating a pole in $\Psi_T$ that corrupts the integration. Verify that $\Psi_T(v)$ is smooth on $[0, v_{\max}]$.

Aliasing

FFT implicitly assumes the integrand is periodic with period $N\lambda = 2\pi/\eta$. Aliasing occurs when the characteristic function has significant energy at frequencies above $v_{\max} = \eta N$. Signs of aliasing: prices wrap around from deep OTM options at high strikes to low strikes; option prices become non-monotone or negative. Mitigation: increase $N$ or $\eta$ (larger $v_{\max}$), or use the fractional FFT (FRFT) which breaks the link $\eta\lambda = 2\pi/N$, allowing independent control of $\eta$ and $\lambda$.

Limitations

Accuracy for deep OTM options. The Fourier inversion is most accurate near the money. For deep OTM options, $C_T(k)$ is very small (machine-epsilon territory), and the max(Re(y_u), 0) applied to remove negative prices (from numerical errors) can introduce relative errors of 100%. For deep OTM options, compute the out-of-the-money put and use put-call parity.

Single maturity. The FFT as described prices all strikes at a single maturity $T$ in one pass. For multiple maturities (calibration to a surface), one FFT pass per maturity is needed. This is still faster than option-by-option integration.

Dampening sensitivity. The choice of $\alpha$ materially affects accuracy for some models. The optimal $\alpha$ depends on model parameters and maturity. In practice, test multiple values of $\alpha$ and cross-validate against analytic prices (when available, e.g., Black-Scholes).

Model restriction. Carr-Madan applies to models where the characteristic function is available in closed form. For local-stochastic-vol models (e.g., SABR-Heston), the characteristic function is not available analytically; Monte Carlo or PDE methods are required.

Fractional FFT (FRFT). The standard FFT imposes $\eta\lambda = 2\pi/N$. The FRFT (Bailey-Swarztrauber) removes this constraint at the cost of doubling the FFT length, enabling fine log-strike grids without large $\eta$. This is the method of choice when fine strike resolution is needed (e.g., calibration to illiquid OTM strikes).

Interview Angle

L1. What is the characteristic function of $\ln S_T$ under Black-Scholes? Why is it useful for option pricing even though Black-Scholes has a closed-form formula?

GBM characteristic function. Under $\mathbb{Q}$: $\ln S_T \sim \mathcal{N}(\ln S_0 + (r-\sigma^2/2)T, \sigma^2 T)$, so $\varphi_T(u) = \exp(iu(\ln S_0 + (r-\sigma^2/2)T) - u^2\sigma^2 T/2)$. It is useful as a template: Heston, VG, and NIG all modify this baseline by replacing the diffusion exponent with a more complex function. The Fourier method works for all such models uniformly — the only change is the characteristic function plugged into the same numerical recipe.

L2. Derive the Carr-Madan formula for $\Psi_T(v)$. State the dampening condition and explain why $\alpha > 0$ is needed. Show how the FFT condition $\eta\lambda = 2\pi/N$ arises.

Derivation outline. Start from $z_T(k) = e^{\alpha k}C_T(k)$, take the Fourier transform, swap expectation and integral via Fubini (valid under the moment condition), evaluate the inner integral over $k$:

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

Combining: $\Psi_T(v) = e^{-rT}\varphi_T(v-(\alpha+1)i)/(\alpha^2+\alpha-v^2+i(2\alpha+1)v)$.

FFT condition: We want to compute $\sum_{j=0}^{N-1} e^{-i\eta\lambda ju} f_j$ using FFT. The FFT computes sums of the form $\sum_j e^{-i(2\pi/N)ju} f_j$. Matching: $\eta\lambda = 2\pi/N$. This forces a trade-off: finer frequency grid ($\eta$ small) → coarser strike grid ($\lambda$ large), and vice versa.

L3. Explain the aliasing phenomenon in Carr-Madan. Derive the FRFT extension that decouples $\eta$ and $\lambda$. For Heston with large $\xi$, when does the moment condition $\mathbb{E}[S_T^{1+\alpha}] < \infty$ fail, and what is the practical implication for $\Psi_T(v)$?

Aliasing. The discrete Fourier sum implicitly wraps the integrand with period $N\lambda$ (the total width of the log-strike grid). If the call price function $z_T(k)$ has significant support outside $[-N\lambda/2, N\lambda/2]$ (i.e., at very high or very low log-strikes), the periodic wrap causes the contribution from outside the grid to be added to points inside — aliasing error. Increase $N$ or $\eta$ (larger $v_{\max}$) to reduce the tails of $\Psi_T(v)$ and the aliasing.

FRFT. Define a fractional DFT: $y_u = \sum_{j=0}^{N-1} e^{-i\delta ju}x_j$ for arbitrary $\delta$ (not necessarily $2\pi/N$). This is computable as a product of three standard FFTs of length $2N$ (using the chirp-z transform). Setting $\delta = \eta\lambda$ (arbitrary, no longer constrained to $2\pi/N$), FRFT allows $\eta$ and $\lambda$ to be chosen independently.

Heston moment failure. The moment generating function $\mathbb{E}[S_T^{1+\alpha}] = \varphi_T(-i(1+\alpha))$ requires evaluating the characteristic function at an imaginary argument. For Heston, $\varphi_T(u)$ is defined via the Riccati ODE for $B(u,T)$; for $u = -i(1+\alpha)$, the discriminant $d(u)$ can become real and large. If $2\kappa\theta < \xi^2(1+\alpha)\alpha$ (a condition on the Feller parameter relative to $\alpha$), the CIR process variance can dominate and $\varphi_T(-i(1+\alpha)) = +\infty$. Practically: the denominator $\alpha^2+\alpha-v^2+i(2\alpha+1)v$ in $\Psi_T(v)$ has a real zero at $v_0 = \sqrt{\alpha^2+\alpha}$. If $\varphi_T$ also has a singularity near this $v_0$, the integrand blows up and numerical integration fails. Solution: reduce $\alpha$ until the moment condition is satisfied, or use the Lewis formula (contour at $\mathrm{Im}(u) = -1/2$) which avoids the pole entirely.