Semi-Closed-Form via Fourier Inversion¶
The Heston model's characteristic function enables a "semi-closed-form" pricing formula for European options: the price is expressed as a deterministic integral involving the characteristic function, reducing option pricing to numerical quadrature. This formula decomposes the call price into two risk-adjusted probabilities \(P_1\) and \(P_2\), each computed by a one-dimensional Fourier inversion integral. The result is exact up to the numerical integration accuracy, providing a benchmark against which all other Heston pricing methods (COS, FFT, Monte Carlo) are validated.
Learning Objectives
By the end of this section, you will be able to:
- Derive the semi-closed-form European call price from the Heston characteristic function
- Understand the \(P_1\)/\(P_2\) probability decomposition and its measure-theoretic meaning
- Evaluate the Fourier inversion integrals using standard quadrature methods
- Assess the accuracy and computational cost of semi-closed-form pricing
Intuition¶
In the Black-Scholes model, the call price \(C = S_0 N(d_1) - K e^{-rT} N(d_2)\) involves two normal CDFs. Replacing the lognormal distribution with the Heston distribution means replacing \(N(d_j)\) with integrals that invert the characteristic function to recover the corresponding exercise probabilities. The characteristic function is known in closed form, so the only numerical step is one-dimensional integration --- making this approach "semi-closed-form." Each integral converges rapidly because the CF decays fast for large arguments, and standard quadrature (Gauss-Laguerre, Simpson, adaptive) achieves 10--15 digits of accuracy with modest computational effort.
European Call Price Decomposition¶
The risk-neutral pricing formula for a European call with strike \(K\) and maturity \(T\) is
where \(\tau = T - t\). Conditioning on the event \(\{S_T > K\}\):
The first term can be rewritten using the stock-price numeraire:
Theorem (Semi-Closed-Form Call Price)
The European call price under the Heston model is
where
are the exercise probabilities under the stock-price numeraire and money-market numeraire respectively.
The European put price follows from put-call parity:
Fourier Inversion of Exercise Probabilities¶
Each probability \(P_j\) is recovered from the corresponding characteristic function by Fourier inversion.
Theorem (Gil-Pelaez Formula for P1 and P2)
The exercise probabilities are
where \(\varphi_1(u)\) is the CF of \(\log S_T\) under \(\mathbb{Q}^S\) and \(\varphi_2(u)\) is the CF under \(\mathbb{Q}\).
Characteristic functions. The \(\mathbb{Q}\)-CF is the standard Heston CF:
with \(C, D\) solving the Heston Riccati system. The \(\mathbb{Q}^S\)-CF is obtained by the shift formula:
Properties of the Integrands¶
The integrands \(f_j(u) = \text{Re}[e^{-iu\log K}\varphi_j(u)/(iu)]\) have important properties that guide numerical integration.
Behavior at \(u = 0\). Both integrands have a removable singularity at \(u = 0\). Using L'Hopital's rule:
In practice, the limit is finite and can be computed from the CF derivatives.
Decay for large \(u\). For the Heston model, \(|\varphi_j(u)| \leq C e^{-\alpha |u|}\) for some \(\alpha > 0\). This exponential decay ensures:
- The integral converges absolutely
- Truncation at a finite upper limit introduces exponentially small error
- Standard quadrature rules achieve rapid convergence
Oscillatory behavior. The factor \(e^{-iu\log K}\) introduces oscillation at frequency proportional to \(|\log K|\). For deep in-the-money or out-of-the-money options (\(K\) far from \(S_0\)), the oscillation is rapid, requiring finer quadrature grids.
Deep OTM Options
For very deep out-of-the-money options (\(K \gg S_0\) or \(K \ll S_0\)), the integrands oscillate rapidly while the option value is small. This causes cancellation errors in the Fourier integral. For such cases, using the in-the-money counterpart via put-call parity and then subtracting is more numerically stable.
Numerical Integration Methods¶
Several quadrature methods are suitable for evaluating the Fourier integrals.
Method 1: Gauss-Laguerre Quadrature¶
The semi-infinite integral \(\int_0^\infty f(u) \, du\) maps naturally to Gauss-Laguerre quadrature after the substitution \(u = t\):
where \((u_n, w_n)\) are the Gauss-Laguerre nodes and weights. With \(N = 32\) to 64 nodes, this typically achieves 10--12 digits of accuracy.
Method 2: Truncated Simpson's Rule¶
Truncate the integral at \(u_{\max}\) and apply composite Simpson's rule on \([0, u_{\max}]\) with \(M\) subintervals:
where \(h = u_{\max}/(2M)\) and \(c_n = (1, 4, 2, 4, 2, \ldots, 4, 1)\) are Simpson weights.
Choosing \(u_{\max} = 50\) and \(M = 500\) (so \(h = 0.05\)) gives approximately 8 digits of accuracy.
Method 3: Adaptive Quadrature¶
Use a library routine (e.g., MATLAB's integral, Python's scipy.integrate.quad) that automatically adjusts the step size based on local error estimates. This is the most robust approach and is recommended for benchmarking.
| Method | CF evaluations | Typical accuracy | Speed |
|---|---|---|---|
| Gauss-Laguerre 64 | 64 | \(10^{-12}\) | 0.1 ms |
| Simpson \(M=500\) | 1001 | \(10^{-8}\) | 1 ms |
| Adaptive quadrature | 200--500 | \(10^{-10}\) | 0.5 ms |
Putting It Together: Complete Algorithm¶
Input: Heston parameters \((\kappa, \theta, \xi, \rho, v_0)\), market data \((S_0, K, r, q, \tau)\).
Step 1. Choose the Albrecher stable formulation for the Heston CF.
Step 2. For \(j = 2\) (\(\mathbb{Q}\) measure): evaluate the integrand
at the quadrature nodes and compute \(P_2 = \frac{1}{2} + \frac{1}{\pi}\int_0^\infty f_2(u) \, du\).
Step 3. For \(j = 1\) (\(\mathbb{Q}^S\) measure): evaluate
using \(\varphi_1(u) = \varphi_2(u-i)/(S_0 e^{(r-q)\tau})\), and compute \(P_1\).
Step 4. Compute the call price:
Numerical Example¶
Parameters: \(S_0 = 100\), \(K = 100\), \(r = 0.05\), \(q = 0\), \(v_0 = 0.04\), \(\kappa = 1.5\), \(\theta = 0.04\), \(\xi = 0.3\), \(\rho = -0.7\), \(\tau = 1\).
Step 1: Compute \(P_2\).
Step 2: Compute \(P_1\).
Step 3: Call price.
Implied volatility: The Black-Scholes implied volatility corresponding to \(C = 8.628\) is approximately \(\sigma_{\text{imp}} = 20.1\%\), close to \(\sqrt{v_0} = 20\%\) since the option is near-ATM and the Heston parameters are close to the lognormal case.
Sensitivity to Quadrature Method
Computing the same price with different methods:
| Method | Call Price | Abs Error vs GL-128 |
|---|---|---|
| Gauss-Laguerre 32 | 8.62812345 | \(2 \times 10^{-8}\) |
| Gauss-Laguerre 64 | 8.62812344 | \(< 10^{-11}\) |
| Gauss-Laguerre 128 | 8.62812344 | reference |
| Simpson \(M=200\) | 8.62812 | \(3 \times 10^{-6}\) |
| Simpson \(M=1000\) | 8.6281234 | \(4 \times 10^{-8}\) |
| Adaptive (tol=\(10^{-10}\)) | 8.62812344 | \(< 10^{-10}\) |
Gauss-Laguerre with 64 nodes achieves near-machine precision, making it the most efficient choice for single-strike pricing.
Advantages and Limitations¶
Advantages:
- Accuracy: 10--15 digits achievable with appropriate quadrature
- Simplicity: requires only a one-dimensional integral
- Benchmark quality: serves as the reference for validating all other methods
- Direct probability extraction: \(P_1\) and \(P_2\) are useful for delta computation and digital pricing
Limitations:
- Single-strike: each strike requires a separate integration (unlike FFT which prices all strikes simultaneously)
- Speed: slower than COS for individual strikes (64 vs 128 CF evaluations, but COS per-evaluation cost is lower)
- Oscillatory integrand: deep OTM options require careful treatment
Role as Benchmark
The semi-closed-form approach is the gold standard for Heston option pricing accuracy. Every implementation of COS, FFT, or Monte Carlo should be validated against the semi-closed-form result with high-order Gauss-Laguerre quadrature. Agreement to 8+ significant digits provides confidence that the implementation is correct.
Summary¶
The semi-closed-form Heston pricing formula expresses the European call as \(C = S_0 e^{-q\tau} P_1 - K e^{-r\tau} P_2\), where \(P_1\) and \(P_2\) are exercise probabilities computed by Gil-Pelaez Fourier inversion of the Heston characteristic function under the stock-price numeraire and risk-neutral measures respectively. The integrals converge rapidly due to the exponential decay of the CF and can be evaluated to 10--15 digit accuracy using Gauss-Laguerre quadrature with 32--64 nodes. While not the fastest method for production pricing (COS and FFT are faster for multiple strikes), the semi-closed-form approach provides the definitive benchmark for validating all other Heston pricing implementations.
Exercises¶
Exercise 1. Starting from the call price \(C = S_0 e^{-q\tau} P_1 - K e^{-r\tau} P_2\), derive the put price using put-call parity. Express the result in terms of \(P_1\) and \(P_2\). Verify that for the worked example (\(S_0 = 100\), \(K = 100\), \(r = 5\%\), \(q = 0\), \(\tau = 1\), \(C = 8.628\)), the put price is \(P = 8.628 - 100 + 100e^{-0.05} = 3.755\).
Solution to Exercise 1
Put-call parity states:
Solving for the put:
Substituting \(C = S_0 e^{-q\tau} P_1 - K e^{-r\tau} P_2\):
This is the Heston European put price, expressed in terms of the probabilities of the option expiring out of the money under each measure: \((1 - P_2) = \mathbb{Q}(S_T \leq K)\) and \((1 - P_1) = \mathbb{Q}^S(S_T \leq K)\).
Verification with the worked example: \(S_0 = 100\), \(K = 100\), \(r = 0.05\), \(q = 0\), \(\tau = 1\), \(C = 8.628\).
Alternatively, using the \(P_1\), \(P_2\) formula directly:
The slight discrepancy (\(3.751\) vs. \(3.758\)) is due to rounding in the stated values of \(P_1\), \(P_2\), and \(C\). To higher precision, both computations yield \(P \approx 3.755\).
Exercise 2. The exercise probability \(P_2 = \mathbb{Q}(S_T > K)\) under the risk-neutral measure, and \(P_1 = \mathbb{Q}^S(S_T > K)\) under the stock-price numeraire. Explain the economic interpretation of each probability. Why is \(P_1 > P_2\) in general for positive-drift assets? Show that in the Black-Scholes model, \(P_1 = N(d_1)\) and \(P_2 = N(d_2)\) with \(d_1 - d_2 = \sigma\sqrt{T}\).
Solution to Exercise 2
Economic interpretation of \(P_2 = \mathbb{Q}(S_T > K)\): This is the risk-neutral probability that the option finishes in the money. It represents the probability of exercise from the perspective of the money-market numeraire. Under \(\mathbb{Q}\), the discounted stock price \(e^{-rt}S_t\) is a martingale, so \(P_2\) measures the likelihood of exercise after adjusting for the risk-free rate of return.
Economic interpretation of \(P_1 = \mathbb{Q}^S(S_T > K)\): This is the exercise probability under the stock-price numeraire (the "share measure" \(\mathbb{Q}^S\)). Under \(\mathbb{Q}^S\), the process \(K e^{-r(T-t)} / S_t\) is a martingale, meaning probabilities are weighted by the stock price. States where \(S_T\) is large receive more weight, making \(P_1 > P_2\) in general.
Why \(P_1 > P_2\): Under \(\mathbb{Q}^S\), the Radon-Nikodym derivative is \(d\mathbb{Q}^S / d\mathbb{Q} = S_T / (S_0 e^{(r-q)\tau})\). This tilts the probability distribution toward high-\(S_T\) states. Since \(\{S_T > K\}\) is precisely the event where \(S_T\) is large, this event receives greater weight under \(\mathbb{Q}^S\) than under \(\mathbb{Q}\), so \(P_1 \geq P_2\) with equality only when \(S_T\) is deterministic.
Black-Scholes verification: In the Black-Scholes model, \(\log S_T \sim N(\log S_0 + (r - q - \sigma^2/2)\tau, \, \sigma^2 \tau)\) under \(\mathbb{Q}\). Then:
where
Under \(\mathbb{Q}^S\), the drift of \(\log S_T\) shifts from \((r - q - \sigma^2/2)\) to \((r - q + \sigma^2/2)\) (by Girsanov's theorem, the Brownian motion gains a drift of \(\sigma\)). Therefore:
The difference is:
Since \(\sigma > 0\) and \(\tau > 0\), we have \(d_1 > d_2\) and thus \(P_1 = N(d_1) > N(d_2) = P_2\).
Exercise 3. The semi-closed-form approach computes \(P_j\) via:
For \(P_2\), the CF is \(\varphi_2(u) = \varphi(u)\) (the standard Heston CF). For \(P_1\), the CF is \(\varphi_1(u) = \varphi(u - i)/\varphi(-i)\). Verify that \(\varphi(-i) = \mathbb{E}[S_T/S_0] \cdot e^{iu\ln S_0}|_{u=-i}\) equals the forward price ratio \(e^{(r-q)\tau}\) (up to the log-spot factor).
Solution to Exercise 3
The \(\mathbb{Q}^S\)-characteristic function is defined as:
where \(\varphi_2(u) = \mathbb{E}^{\mathbb{Q}}[e^{iu \log S_T}]\) is the standard Heston CF under \(\mathbb{Q}\). We need to verify that \(\varphi_2(-i) = S_0 e^{(r-q)\tau}\).
Evaluate \(\varphi_2(-i)\):
Under the risk-neutral measure \(\mathbb{Q}\), the discounted stock price \(e^{-(r-q)t}S_t\) is a martingale (accounting for continuous dividend yield \(q\)). Therefore:
This confirms \(\varphi_2(-i) = S_0 e^{(r-q)\tau}\).
Now verify the measure change. Under \(\mathbb{Q}^S\), the CF of \(\log S_T\) is:
Since \(S_T = e^{\log S_T}\):
This completes the verification: the ratio \(\varphi_2(u-i)/\varphi_2(-i)\) correctly implements the change of numeraire from the money-market account to the stock price, with \(\varphi_2(-i)\) equal to the forward price ratio \(S_0 e^{(r-q)\tau}\).
Exercise 4. Using Gauss-Laguerre quadrature with \(N = 32\) nodes, estimate the total number of CF evaluations needed to price a single European call (accounting for both \(P_1\) and \(P_2\) integrals). If each CF evaluation costs \(T_{\text{cf}} = 0.5\) microseconds, estimate the pricing time. Compare with the COS method using \(N_{\text{cos}} = 128\) terms.
Solution to Exercise 4
Each exercise probability \(P_j\) requires one Fourier inversion integral. Using Gauss-Laguerre with \(N = 32\) nodes, each integral requires 32 CF evaluations. Pricing a single European call requires both \(P_1\) and \(P_2\), so:
Note: in practice, \(\varphi_1(u) = \varphi_2(u-i)/\varphi_2(-i)\) can reuse parts of the \(\varphi_2\) computation, but each call to \(\varphi_1\) at a different node requires a separate CF evaluation at a shifted argument, so we count 64 total.
Pricing time estimate: At \(T_{\text{cf}} = 0.5\) microseconds per CF evaluation:
Adding overhead for quadrature weights, summation, and the final price computation (negligible compared to CF evaluations), the total is approximately \(0.03\)--\(0.05\) ms.
Comparison with COS method: The COS method with \(N_{\text{cos}} = 128\) terms requires 128 CF evaluations (one per cosine expansion coefficient). The pricing time is:
However, the COS method has two advantages: (1) each CF evaluation is at an integer multiple of a base frequency, allowing potential vectorization, and (2) after computing the cosine coefficients once, additional strikes can be priced with only \(O(N_{\text{cos}})\) multiply-add operations (no additional CF evaluations). For a single strike, the semi-closed-form approach with Gauss-Laguerre 32 is faster (64 vs. 128 CF evaluations). For multiple strikes with the same maturity, COS becomes more efficient after the initial 128 CF evaluations are amortized.
Exercise 5. The sensitivity table shows that Simpson's rule with \(M = 200\) points achieves an error of \(3 \times 10^{-6}\), while Gauss-Laguerre with only 32 nodes achieves \(2 \times 10^{-8}\). Explain why Gauss-Laguerre is so much more efficient. Hint: the Heston CF decays as \(e^{-cu^2}\) for large \(u\), and Gauss-Laguerre quadrature is exact for integrands of the form \(p(u)e^{-u}\) where \(p\) is a polynomial of degree \(\leq 2N - 1\).
Solution to Exercise 5
Gauss-Laguerre quadrature is dramatically more efficient than Simpson's rule for the Gil-Pelaez integral because of a fundamental match between the quadrature's design and the integrand's structure.
Simpson's rule is a polynomial-based method that approximates the integrand as piecewise cubic on each subinterval. Its error is \(\mathcal{O}(h^4 f^{(4)})\), giving algebraic (polynomial) convergence. With \(M\) subintervals on \([0, u_{\max}]\), \(h = u_{\max}/M\), and the error decreases as \(M^{-4}\). Achieving \(10^{-8}\) accuracy requires approximately 1000 nodes.
Gauss-Laguerre quadrature with \(N\) nodes is exact for integrands of the form \(p(u)e^{-u}\) where \(p\) is a polynomial of degree \(\leq 2N - 1\). For smooth integrands that can be well-approximated by such products, the convergence is exponential (spectral): the error decreases as \(O(C^{-N})\) for some \(C > 1\).
The Heston CF decays as \(|\varphi(u)| \sim e^{-cu^2}\) for large \(u\), which decays even faster than \(e^{-u}\). After the substitution \(u = \alpha t\) and multiplication by \(e^t\) (to extract the Laguerre weight), the remaining integrand \(g(\alpha t) e^t\) is an entire function of \(t\) that grows at most polynomially times \(e^{t - c\alpha^2 t^2}\), which decreases super-exponentially. This means \(g(\alpha t) e^t\) is extremely well-approximated by polynomials, and Gauss-Laguerre quadrature converges exponentially fast.
In contrast, Simpson's rule does not exploit the smoothness of the integrand globally --- it only uses local cubic approximations. Even though the integrand is entire and rapidly decaying, Simpson's rule converges only algebraically.
Quantitatively: Gauss-Laguerre with 32 nodes achieves \(2 \times 10^{-8}\) accuracy (as shown in the table), while Simpson's rule with 200 subintervals (401 nodes) achieves only \(3 \times 10^{-6}\). Gauss-Laguerre uses 12.5 times fewer function evaluations yet achieves 150 times better accuracy. This efficiency gap widens further at higher accuracy targets.
Exercise 6. A deep OTM put with \(K = 70\), \(S_0 = 100\), \(\tau = 0.25\) has a very small price. Compute the probabilities \(P_1\) and \(P_2\) by noting that \(\ln(K/S_0) = \ln(0.7) = -0.357\). The factor \(e^{-iu\ln K}\) causes the integrand to oscillate with period \(2\pi / |\ln K| \approx 1.77\). Estimate how many Gauss-Laguerre nodes are needed to resolve these oscillations accurately. Would adaptive quadrature be preferable in this case?
Solution to Exercise 6
For the deep OTM put with \(K = 70\), \(S_0 = 100\), \(\tau = 0.25\):
The log-strike is \(\ln K = \ln 70 = 4.248\). The oscillation in the integrand comes from the factor \(e^{-iu \ln K}\), which oscillates with period:
However, the relevant oscillation frequency for distinguishing the option from ATM behavior is driven by \(\ln(K/S_0) = -0.357\), giving period \(2\pi / 0.357 \approx 17.6\). The absolute \(\ln K = 4.248\) determines the raw oscillation, but the \(\ln S_0\) component cancels against the CF's \(e^{iu \ln S_0}\) factor.
Using \(|\ln K| = 4.248\), the oscillation period is approximately 1.48. The Heston CF decays with effective support up to \(u_{\text{eff}} \approx \sqrt{40 / (v_0 \tau)} = \sqrt{40 / (0.04 \times 0.25)} = \sqrt{4000} \approx 63\). The number of oscillation cycles in \([0, u_{\text{eff}}]\) is:
To accurately resolve 43 oscillation cycles with Gauss-Laguerre, we need at least 6--8 nodes per cycle (for reliable integration of oscillatory functions), giving:
Standard Gauss-Laguerre with 64 nodes provides only \(64 / 43 \approx 1.5\) nodes per cycle, which is grossly insufficient for accurate quadrature of this oscillatory integrand.
Would adaptive quadrature be preferable? Yes, adaptive quadrature is strongly preferred in this case. An adaptive scheme (e.g., Gauss-Kronrod) automatically places more nodes in the region \([0, \sim 20]\) where the integrand oscillates significantly and fewer nodes in the tail where it has decayed. Typically 300--500 function evaluations suffice for 10-digit accuracy, compared to the 300+ Gauss-Laguerre nodes that would be needed.
An alternative is to use the put-call parity approach: price the deep ITM call (\(K = 70\) call is 30% ITM) where the log-moneyness \(\ln(S_0/K) = 0.357\) is moderate, then convert to the put price. The ITM call integrand oscillates more slowly (same frequency but the large deterministic component \(S_0 e^{-q\tau} - Ke^{-r\tau}\) does not require the Fourier integral), making standard quadrature more reliable.
Exercise 7. The semi-closed-form formula directly provides \(P_1\), which is the Heston delta: \(\Delta = e^{-q\tau} P_1\). Compute the delta for the worked example (\(P_1 = 0.57481\), \(q = 0\), \(\tau = 1\)) and compare with the Black-Scholes delta \(N(d_1)\) at the implied volatility \(\sigma_{\text{imp}} = 20.1\%\). Explain why the Heston delta differs from the Black-Scholes delta even when the implied volatility is similar, and discuss the implications for hedging.
Solution to Exercise 7
Heston delta computation: The Heston delta for a European call is:
With \(P_1 = 0.57481\), \(q = 0\), \(\tau = 1\):
Black-Scholes delta: At implied volatility \(\sigma_{\text{imp}} = 0.201\), with \(S_0 = K = 100\), \(r = 0.05\), \(q = 0\), \(\tau = 1\):
Therefore \(\Delta_{\text{Heston}} = 0.575\) versus \(\Delta_{\text{BS}} = 0.636\). The Heston delta is lower by about \(0.061\), a substantial difference of roughly 6 percentage points of notional.
Why the deltas differ despite similar implied volatilities:
-
Skew effect on delta: The Heston model with \(\rho = -0.7\) generates a pronounced negative volatility skew (implied volatility decreases with strike). When \(S_0\) increases, the option moves further ITM, and the relevant implied volatility shifts to a lower-strike point on the skew, where implied volatility is higher. The Black-Scholes delta ignores this strike-dependence of implied volatility, while the Heston delta implicitly accounts for it.
-
Stochastic volatility adjustment: In the Heston model, the delta includes the effect of the correlation between stock returns and variance changes. With \(\rho < 0\), a stock price increase is associated with a variance decrease, which partially offsets the increase in option value. This makes the Heston delta smaller than the BS delta.
-
The "sticky strike" vs. "sticky delta" distinction: The BS delta \(N(d_1)\) assumes the implied volatility surface does not move when \(S_0\) changes ("sticky strike"). The Heston delta incorporates the dynamic response of the entire volatility surface to stock price changes, producing a different hedge ratio.
Hedging implications: Using the BS delta (\(0.636\)) instead of the Heston delta (\(0.575\)) would lead to over-hedging: the portfolio would hold too many shares, resulting in a hedge P&L that is systematically exposed to volatility risk. For a negative-correlation stochastic volatility model, the correct delta is always below the BS delta for ATM calls. The difference (\(\sim 6\%\) of notional in this case) is economically significant and grows larger for longer maturities and more negative correlation.