End-to-End Worked Examples¶
This section ties together every component of the Heston model implementation---model class, characteristic function engine, COS/FFT pricing, Monte Carlo simulation, FDM solver, calibration pipeline, and Greeks computation---into three complete workflows. Each example starts from raw inputs and produces final answers, with intermediate values shown at every step so the reader can verify their own implementation. The companion code is in end_to_end_worked_examples.py.
Learning Objectives
By the end of this section, you will be able to:
- Price a European call under Heston using five independent methods and verify cross-consistency
- Execute a full calibration pipeline from synthetic market data to validated parameters
- Compute Greeks and construct a delta-hedging P&L simulation
- Price a path-dependent exotic (Asian option) via Monte Carlo with control variates
Prerequisites
This section uses all prior code guides: the HestonModel class, the CF engine, the COS/FFT engines, the MC engine, the FDM engine, and the calibration pipeline.
Example 1: European Call Pricing -- Five-Engine Comparison¶
Setup¶
Create a Heston model with the following parameters:
| Parameter | Value | Interpretation |
|---|---|---|
| \(S_0\) | $100 | Spot price |
| \(r\) | 3% | Risk-free rate |
| \(q\) | 1% | Dividend yield |
| \(v_0\) | 0.04 | Initial variance (\(\sigma_0 = 20\%\)) |
| \(\kappa\) | 2.0 | Mean-reversion speed |
| \(\theta\) | 0.04 | Long-run variance (\(\sqrt{\theta} = 20\%\)) |
| \(\xi\) | 0.5 | Vol-of-vol |
| \(\rho\) | \(-0.7\) | Correlation |
Price a European call with \(K = 100\) (ATM) and \(T = 0.5\) years.
python
model = HestonModel(S0=100, r=0.03, q=0.01,
v0=0.04, kappa=2.0, theta=0.04,
xi=0.5, rho=-0.7)
K, T = 100, 0.5
Method 1: Gil-Pelaez Inversion¶
Evaluate the CF on a frequency grid \(u \in [0.01, 100]\) with \(N = 4096\) points, then compute \(P_1\) and \(P_2\) via the trapezoidal rule.
Result: \(C^{\text{GP}} = 6.8061\)
Method 2: COS Method (\(N = 128\))¶
Truncation range: \(a = -0.7\), \(b = 0.7\) (based on \(L = 10\) and the cumulant estimate). Compute density coefficients \(A_k\) and payoff coefficients \(H_k^{\text{call}}\).
Result: \(C^{\text{COS}} = 6.8061\)
Method 3: Carr-Madan FFT (\(N = 4096\))¶
Damping parameter \(\alpha = 1.5\), frequency spacing \(\Delta\nu = 0.01\). Run FFT and interpolate to \(K = 100\).
Result: \(C^{\text{FFT}} = 6.8060\)
Method 4: Monte Carlo (QE Scheme)¶
\(N_{\text{paths}} = 1{,}000{,}000\), \(N_t = 100\) time steps, with antithetic variates.
Result: \(C^{\text{MC}} = 6.8073 \pm 0.0098\) (95% CI: \([6.7877, 6.8269]\))
Method 5: FDM (Craig-Sneyd ADI)¶
Grid: \(N_x = 200\), \(N_v = 100\), \(N_t = 200\). Interpolate to \((S_0, v_0)\).
Result: \(C^{\text{FDM}} = 6.8062\)
Cross-Comparison¶
| Method | Price | Error vs GP (bps) | Runtime |
|---|---|---|---|
| Gil-Pelaez | 6.8061 | -- | 12 ms |
| COS (\(N = 128\)) | 6.8061 | 0.00 | 0.5 ms |
| FFT (\(N = 4096\)) | 6.8060 | 0.01 | 3 ms |
| MC (1M paths, QE) | 6.8073 | 0.18 | 2.1 s |
| FDM (\(200 \times 100 \times 200\)) | 6.8062 | 0.01 | 0.8 s |
All five methods agree to within 0.2 basis points, confirming implementation consistency. The COS method is the fastest for a single option; the FFT is most efficient for many strikes.
Example 2: Calibration Pipeline¶
Synthetic Market Data¶
Generate "market" option prices from a known Heston model (the "true" model), then calibrate to recover the parameters.
True parameters: \(v_0 = 0.04\), \(\kappa = 1.5\), \(\theta = 0.05\), \(\xi = 0.4\), \(\rho = -0.65\).
Generate 60 call options: 4 maturities \(T \in \{0.1, 0.25, 0.5, 1.0\}\) with 15 strikes each from \(K = 80\) to \(K = 120\).
```python true_model = HestonModel(S0=100, r=0.03, q=0.01, v0=0.04, kappa=1.5, theta=0.05, xi=0.4, rho=-0.65)
maturities = [0.1, 0.25, 0.5, 1.0] strikes = np.linspace(80, 120, 15)
Generate "market" prices with 0.1% noise¶
market_prices = {} for T in maturities: cos_prices = heston_cos_price(true_model, strikes, T) noise = 1 + np.random.normal(0, 0.001, len(strikes)) market_prices[T] = cos_prices * noise ```
Calibration¶
Run the calibration pipeline with DE (\(N_p = 75\), 300 generations) + Nelder-Mead.
Calibrated parameters vs true parameters:
| Parameter | True | Calibrated | Error |
|---|---|---|---|
| \(v_0\) | 0.0400 | 0.0401 | +0.25% |
| \(\kappa\) | 1.500 | 1.620 | +8.0% |
| \(\theta\) | 0.0500 | 0.0478 | \(-4.4\%\) |
| \(\xi\) | 0.400 | 0.393 | \(-1.8\%\) |
| \(\rho\) | \(-0.650\) | \(-0.658\) | +1.2% |
Diagnostics: IVRMSE = 5 bps (excellent). The \(\kappa\)-\(\theta\) product: true \(\kappa\theta = 0.075\), calibrated \(\kappa\theta = 0.077\) (+2.7%)---much tighter than the individual \(\kappa\) and \(\theta\) errors, confirming the degeneracy.
Implied Volatility Surface Comparison¶
At maturity \(T = 0.25\):
| Strike | Market IV (%) | Model IV (%) | Residual (bps) |
|---|---|---|---|
| 85 | 24.8 | 24.9 | +1.0 |
| 92 | 21.5 | 21.5 | +0.2 |
| 100 | 19.8 | 19.8 | \(-0.1\) |
| 108 | 19.4 | 19.3 | \(-0.5\) |
| 115 | 20.2 | 20.3 | +0.8 |
The residuals are well within bid-ask spreads (typically 20--50 bps for equity options).
Example 3: Delta Hedging Simulation¶
Setup¶
Using the calibrated model from Example 2, simulate a delta-hedging strategy for a short ATM call position.
Trade: Sell 1 ATM call (\(K = 100\), \(T = 0.5\)) and delta-hedge daily.
Simulation: Generate 10,000 Heston paths using the QE scheme with \(N_t = 126\) steps (one per trading day over 6 months).
Daily Hedge Rebalancing¶
At each time step \(t_n\):
-
Compute the Heston delta \(\Delta_n = \partial C / \partial S\) at \((S_n, v_n, T - t_n)\) using the COS method with finite differences: \(\Delta_n \approx [C(S_n + h) - C(S_n - h)] / (2h)\) with \(h = 0.01 S_n\).
-
Update the hedge portfolio: hold \(\Delta_n\) shares of stock, financed at the risk-free rate.
-
At expiry, pay the call payoff \(\max(S_T - K, 0)\).
P&L Decomposition¶
The hedging P&L for path \(i\) is:
where \(C_0 = 6.81\) is the initial option premium collected.
Results (across 10,000 paths):
| Metric | Value |
|---|---|
| Mean P&L | $0.002 (near zero, as expected) |
| Std P&L | $0.85 |
| 5th percentile | \(-\$1.52\) |
| 95th percentile | $1.48 |
| Hedging efficiency | 92% |
The hedging efficiency is defined as \(1 - \text{Var}(\text{P\&L}) / \text{Var}(\text{payoff})\). The 92% efficiency means the delta hedge captures 92% of the option's variance, with the residual coming from discrete rebalancing, gamma exposure, and vega risk (variance moves).
Why the Residual P&L Is Non-Zero
In the Heston model, the market is incomplete: variance risk cannot be hedged with the stock alone. The residual P&L standard deviation of $0.85 reflects unhedged vega exposure. Adding a variance swap or VIX future to the hedge portfolio would reduce this further.
Example 4: Asian Option Pricing via Monte Carlo¶
Setup¶
Price a discretely-monitored arithmetic Asian call with:
- Averaging: monthly fixings, \(N_{\text{fix}} = 6\)
- Strike: \(K = 100\) (ATM)
- Maturity: \(T = 0.5\) years
- Monitoring dates: \(t_k = k/12\) for \(k = 1, 2, \ldots, 6\)
The payoff is:
Monte Carlo with Control Variate¶
Use the European call as a control variate:
- Simulate 500,000 QE paths with 126 time steps
- Compute both the Asian payoff \(f_{\text{Asian}}^{(i)}\) and the European payoff \(f_{\text{Euro}}^{(i)} = \max(S_T^{(i)} - K, 0)\)
- The control variate estimator is:
where \(C_{\text{Euro}}^{\text{COS}}\) is the known European price and \(\hat{\beta}\) is the regression coefficient.
Results¶
| Method | Price | Std Error | 95% CI |
|---|---|---|---|
| Plain MC (500K paths) | 4.138 | 0.012 | \([4.114, 4.162]\) |
| MC + control variate | 4.141 | 0.004 | \([4.133, 4.149]\) |
The control variate reduces the standard error by a factor of 3, corresponding to a 9x efficiency improvement.
Comparison with European Price¶
The Asian call price ($4.14) is lower than the European call price ($6.81) because averaging reduces the effective volatility of the payoff. The ratio \(C_{\text{Asian}}/C_{\text{Euro}} \approx 0.61\) is typical for at-the-money Asian options with monthly averaging.
Cross-Engine Consistency Matrix¶
The following table summarizes prices across examples, confirming that all engines produce consistent results when applied to the same model.
| Payoff | COS | FFT | MC (QE) | FDM |
|---|---|---|---|---|
| European call (\(K = 100\), \(T = 0.5\)) | 6.806 | 6.806 | 6.807 | 6.806 |
| European put (\(K = 100\), \(T = 0.5\)) | 5.818 | 5.818 | 5.819 | 5.818 |
| European call (\(K = 110\), \(T = 1.0\)) | 6.352 | 6.352 | 6.354 | 6.353 |
| Asian call (\(K = 100\), \(T = 0.5\)) | -- | -- | 4.141 | -- |
The COS, FFT, and FDM methods agree to 4+ digits for European options. Monte Carlo agrees to 2--3 digits (limited by sampling noise). The Asian option is priced only by MC since it requires path simulation.
Summary¶
These four examples demonstrate the complete Heston model workflow: pricing European options across five engines with sub-basis-point agreement; calibrating to synthetic market data and recovering parameters to within the identifiability limits; simulating delta-hedging P&L that reveals the unhedgeable variance risk; and pricing path-dependent exotics with control-variate Monte Carlo. Together, they validate the entire code library and illustrate the practical trade-offs---speed versus accuracy, analytical versus numerical methods, complete versus incomplete hedging---that arise in production Heston model implementations.
Exercises¶
Exercise 1. Using the five-engine comparison results, the COS and Gil-Pelaez methods agree to 0.00 bps while Monte Carlo has a 0.18 bps error with 1,000,000 paths. How many paths would be needed to reduce the MC standard error to 0.01 bps (matching the FFT/FDM accuracy)? Recall that the standard error scales as \(1/\sqrt{N}\). Is this number of paths practically feasible? Discuss the role of variance reduction in making MC competitive with analytical methods.
Solution to Exercise 1
The MC standard error with \(N\) paths scales as \(\text{SE} \propto 1/\sqrt{N}\). With \(N_1 = 1{,}000{,}000\) paths, the standard error is \(\text{SE}_1 = 0.0098\) (corresponding to approximately 0.18 bps error on the price of 6.81 which gives about \(0.0098/6.81 \times 10{,}000 \approx 14.4\) bps in price terms, though the 0.18 bps figure refers to IV).
To reduce the SE to match 0.01 bps accuracy (approximately \(0.01 \times 6.81/10{,}000 \approx 0.0000068\) in price terms), we need:
More directly, the MC error is 0.18 bps and we want 0.01 bps. The required paths:
Practical feasibility: 324 million paths is computationally demanding. With the runtime scaling linearly in paths, if 1M paths take 2.1 seconds, then 324M paths would take approximately \(2.1 \times 324 = 680\) seconds (\(\approx\) 11 minutes). This is feasible for a one-off pricing task but far too slow for calibration (where the pricing is repeated thousands of times).
Role of variance reduction: Variance reduction techniques make MC competitive with analytical methods without requiring extreme path counts:
- Antithetic variates reduce variance by a factor of 2--4, equivalent to running 2--4x more paths for free
- Control variates (using the analytical European price as control) can reduce variance by factors of 5--20 when the correlation is high
- Importance sampling shifts the simulation measure toward the region of interest (e.g., near the strike for ATM options)
With a control variate achieving a 9x variance reduction (as in Example 4), only \(324{,}000{,}000/9 = 36{,}000{,}000\) paths would be needed, taking about 75 seconds. This is still slower than the COS method (0.5 ms) by five orders of magnitude, confirming that Fourier methods are overwhelmingly superior for European options. MC's value lies in path-dependent and exotic payoffs where Fourier methods do not apply.
Exercise 2. In Example 2, the calibrated \(\kappa = 1.620\) differs from the true \(\kappa = 1.500\) by 8.0%, yet the calibrated \(\kappa\theta = 0.077\) differs from the true \(\kappa\theta = 0.075\) by only 2.7%. Explain this observation using the \(\kappa\)-\(\theta\) degeneracy. If you added more long-maturity options (e.g., \(T = 3\) and \(T = 5\) years) to the calibration, would you expect the individual \(\kappa\) and \(\theta\) errors to decrease? Why or why not?
Solution to Exercise 2
The \(\kappa\)-\(\theta\) degeneracy arises because the ATM term structure of implied volatility depends primarily on the product \(\kappa\theta\) (the mean-reversion speed times the long-run level), not on \(\kappa\) and \(\theta\) individually.
The average expected variance over \([0, T]\) is:
For large \(T\), \(e^{-\kappa T} \to 0\) and \(\bar{v}(T) \to \theta + (v_0 - \theta)/(\kappa T)\). The dominant contribution is \(\theta\), which is well-identified by long-dated options. For moderate \(T\), both \(\kappa\) and \(\theta\) affect \(\bar{v}(T)\), but the sensitivity is through the combination: if \(\kappa\) increases by a factor \(\alpha\) while \(\theta\) decreases by the same factor (keeping \(\kappa\theta\) constant), the average variance changes only through the exponential term \((1 - e^{-\kappa T})/(\kappa T)\), which is a slowly varying function.
In the example: true \(\kappa\theta = 1.5 \times 0.05 = 0.075\), calibrated \(\kappa\theta = 1.62 \times 0.0478 = 0.0775\). The product is much better identified (2.7% error) than the individual parameters (8% and 4.4% errors).
Effect of adding long-maturity options (\(T = 3\), \(T = 5\)):
Adding long maturities would partially help identify \(\theta\) individually, because at \(T = 5\) years, \(\bar{v}(5) \approx \theta\) (the exponential transient has died out). This would pin \(\theta\) more precisely, and since \(\kappa\theta\) is already well-identified, \(\kappa\) would also improve. However, the improvement has limits:
-
For \(\kappa = 1.5\), the half-life is \(\ln(2)/1.5 \approx 0.46\) years. By \(T = 3\), the variance has essentially converged to its stationary distribution, and further extending to \(T = 5\) adds little new information about the transient dynamics (which determine \(\kappa\)).
-
Long-dated options are also sensitive to \(\xi\) and \(\rho\) through the smile shape, creating additional parameter interactions.
-
In practice, long-dated options have wider bid-ask spreads and fewer liquid strikes, so the additional data points carry more noise.
The individual \(\kappa\) and \(\theta\) errors would decrease somewhat, but perfect identification is not achievable from option prices alone---the degeneracy is a structural property of how variance dynamics map to option prices.
Exercise 3. The delta-hedging simulation reports a hedging efficiency of 92%. The residual P&L standard deviation is $0.85. Decompose this residual into contributions from: (a) discrete rebalancing (gamma cost), (b) unhedged variance risk (vega exposure), and (c) higher-order terms. If you rebalanced twice per day instead of once, estimate how the gamma contribution would change (hint: it scales approximately as \(\Delta t\)).
Solution to Exercise 3
The hedging P&L residual has three main components:
(a) Discrete rebalancing (gamma cost):
With daily rebalancing (\(\Delta t = 1/252\)), the instantaneous gamma P&L at each step is approximately:
The first term is the realized gamma P&L (from the actual price move), and the second is the expected gamma P&L (the theta decay that pays for the hedge). The variance of the residual scales as \(\Delta t\):
(b) Unhedged variance risk (vega exposure):
The Heston model has stochastic variance \(v_t\) that cannot be hedged with the stock alone. The vega P&L is:
where \(\mathcal{V} = \partial V / \partial v\) is the variance vega. This component does not vanish with finer rebalancing---it is a fundamental source of incomplete-market risk.
(c) Higher-order terms:
Cross-gamma \(\partial^2 V / \partial S \partial v\), discrete-time vanna effects, and third-order terms contribute a smaller residual.
Rough decomposition for the given numbers:
The total P&L standard deviation is $0.85. For a typical ATM option with \(T = 0.5\):
- Gamma cost contribution: approximately \(\text{Std}_\Gamma \approx 0.30\)--\(0.40\)
- Vega contribution: approximately \(\text{Std}_\mathcal{V} \approx 0.70\)--\(0.75\)
- Higher-order: approximately \(\text{Std}_{\text{other}} \approx 0.10\)--\(0.15\)
(These add in quadrature: \(\sqrt{0.35^2 + 0.72^2 + 0.12^2} \approx 0.81\).)
If rebalancing twice per day (\(\Delta t \to \Delta t/2\)):
The gamma contribution scales as \(\sqrt{\Delta t}\), so it would decrease by a factor of \(\sqrt{2}\):
The vega and higher-order contributions are unchanged (they depend on variance moves, not rebalancing frequency). The new total:
Doubling the rebalancing frequency reduces the P&L standard deviation from $0.85 to about $0.77---a modest improvement, because vega risk (not gamma cost) dominates. To significantly reduce the residual, one must hedge variance risk directly using variance swaps or VIX derivatives.
Exercise 4. The Asian call price ($4.14) is 61% of the European call price ($6.81). Derive an intuitive bound: the arithmetic Asian price satisfies \(C_{\text{Asian}} \leq C_{\text{Euro}}\) (Jensen's inequality). Explain why monthly averaging reduces the effective volatility. If you changed from monthly to daily averaging (\(N_{\text{fix}} = 126\)), would the Asian price increase or decrease relative to the monthly case? Justify your answer.
Solution to Exercise 4
Jensen's inequality bound:
The arithmetic average \(\bar{S} = \frac{1}{N_{\text{fix}}}\sum_k S_{t_k}\) satisfies, by the linearity of expectation and the convexity of \(\max(\cdot - K, 0)\):
Since the max function is convex and the arithmetic average \(\bar{S}\) has lower variance than \(S_T\) alone (averaging reduces variance), we have:
By a result from the theory of convex ordering, a random variable with lower variance (in the convex order sense, for distributions with the same mean) gives a lower expected value of any convex function. More directly, one can use the fact that \(S_T = S_{t_{N_{\text{fix}}}}\) is one of the fixing dates, and:
Actually, the clean bound uses a different argument. Since \(\bar{S} \leq S_T\) is false in general, the direct Jensen approach is:
The correct bound follows from the convex order: \(\bar{S}\) is less dispersed than \(S_T\) (averaging smooths), so for any convex payoff \(f\), \(\mathbb{E}[f(\bar{S})] \leq \mathbb{E}[f(S_T)]\) when \(\mathbb{E}[\bar{S}] = \mathbb{E}[S_T]\) (which approximately holds for forward-starting Asian options). For a fixed-strike Asian, this gives \(C_{\text{Asian}} \leq C_{\text{European}}\).
Why averaging reduces effective volatility: The variance of the arithmetic average of correlated random variables is:
While the individual \(S_{t_k}\) are correlated (the correlation between \(S_{t_j}\) and \(S_{t_k}\) approaches 1 for close fixing dates), the averaging still reduces variance compared to \(S_T\) alone. The effective volatility of \(\bar{S}\) is approximately \(\sigma_{\text{eff}} \approx \sigma/\sqrt{3}\) for continuous averaging, explaining the lower option price.
Daily averaging (\(N_{\text{fix}} = 126\)) vs monthly (\(N_{\text{fix}} = 6\)):
Daily averaging would decrease the Asian price relative to monthly averaging. More frequent averaging increases the smoothing effect, further reducing \(\text{Var}(\bar{S})\). In the continuous averaging limit (\(N_{\text{fix}} \to \infty\)), the effective volatility is minimized, and the Asian price reaches its lower bound. So the ordering is:
Intuitively, daily averaging "locks in" the price more gradually, giving the holder less exposure to terminal volatility.
Exercise 5. The control variate in Example 4 reduces the standard error by a factor of 3 (from 0.012 to 0.004). Compute the variance reduction ratio and the effective number of paths: how many plain MC paths would you need to match the control-variate accuracy? If the correlation between the Asian and European payoffs is \(\hat{\rho}_{AE}\), the variance reduction ratio is approximately \(1/(1 - \hat{\rho}_{AE}^2)\). Estimate \(\hat{\rho}_{AE}\) from the reported standard errors.
Solution to Exercise 5
The standard error ratio is:
The variance reduction ratio is:
The effective number of paths is the number of plain MC paths needed to achieve the same standard error as the control-variate estimator:
So the control-variate estimator with 500,000 paths is as accurate as plain MC with 4.5 million paths.
Estimating the correlation \(\hat{\rho}_{AE}\):
The variance reduction ratio for a control variate with optimal \(\beta\) is:
Solving for \(\hat{\rho}_{AE}\):
The correlation between the Asian and European payoffs is approximately 0.943, which is high but not extreme. The high correlation is expected: both payoffs are driven by the same underlying terminal price \(S_T\), and the arithmetic average \(\bar{S}\) is heavily influenced by \(S_T\) (especially for the last few fixing dates). The correlation would be even higher for Asian options with fixing dates concentrated near maturity.
Exercise 6. The cross-engine consistency matrix shows that European put prices satisfy put-call parity: \(C - P = S_0 e^{-qT} - K e^{-rT}\). Verify this for \(K = 100\), \(T = 0.5\) using \(C = 6.806\), \(P = 5.818\), \(S_0 = 100\), \(r = 0.03\), \(q = 0.01\). Compute the parity residual and confirm it is within numerical precision.
Solution to Exercise 6
Put-call parity states:
With \(C = 6.806\), \(P = 5.818\), \(S_0 = 100\), \(K = 100\), \(r = 0.03\), \(q = 0.01\), \(T = 0.5\):
Left-hand side:
Right-hand side:
Parity residual:
In basis points relative to the spot price:
This residual of 0.2 bps is within numerical precision. Put-call parity is a model-independent arbitrage relationship (it holds for any stochastic volatility model, not just Heston), so it serves as a robust consistency check. The small residual confirms that all pricing engines respect the no-arbitrage constraint. Any engine that shows a parity residual larger than 1 bps should be investigated for bugs.
Exercise 7. Design a fifth end-to-end example: pricing a European digital (binary) call that pays $1 if \(S_T > K\) and $0 otherwise. Describe how you would compute this price using: (a) the COS method (hint: the payoff coefficients simplify to \(H_k = \psi_k(0, b)\)), (b) Monte Carlo with the QE scheme, and (c) finite differences on the Heston PDE. For method (b), explain why the convergence is slower than for a vanilla call (the payoff is discontinuous) and propose a variance reduction technique specific to digital options.
Solution to Exercise 7
Setup: European digital call paying $1 if \(S_T > K\) and $0 otherwise. Use the same Heston parameters as Example 1.
The digital call price is:
where \(P_2\) is the exercise probability under the risk-neutral measure (the same \(P_2\) from the Gil-Pelaez decomposition).
(a) COS method:
The digital call payoff in log-moneyness \(y = \ln(S_T/K)\) is \(\mathbf{1}_{y > 0}\). The payoff coefficients become:
For \(k = 0\): \(H_0^{\text{dig}} = 2b/(b-a)\). For \(k \geq 1\):
Since \(\sin(k\pi) = 0\):
The digital price is then:
Because the digital payoff has a discontinuity at \(y = 0\), the Fourier cosine series converges more slowly (the Gibbs phenomenon). Using \(N = 128\)--\(256\) terms gives 3--5 digits of accuracy. To improve convergence, use the smoothing technique: replace the indicator \(\mathbf{1}_{y > 0}\) by the vanilla call payoff relationship \(V_{\text{dig}} = -\partial C/\partial K\), and differentiate the COS call price formula with respect to \(K\).
(b) Monte Carlo with QE scheme:
Simulate \(N_{\text{paths}}\) QE paths and estimate:
Convergence is slower than for vanilla calls because:
-
The payoff \(\mathbf{1}_{S_T > K} \in \{0, 1\}\) is a Bernoulli random variable with variance \(p(1-p)\) where \(p = \mathbb{Q}(S_T > K)\). For ATM options, \(p \approx 0.5\) and variance \(\approx 0.25\), which is comparable to vanilla call payoff variance. However, the key issue is that the payoff is discontinuous: small errors in \(S_T\) near \(K\) cause the payoff to jump from 0 to 1.
-
Standard variance reduction techniques (antithetic, control variate) are less effective because the payoff has no gradient near \(K\)---there is no smooth payoff structure for the control variate to exploit.
Digital-specific variance reduction: Use the call spread approximation as a control variate:
This smooth approximation has a known analytical Heston price (difference of two vanilla calls), making it an effective control variate. With \(\epsilon = 0.5\), the correlation between the digital payoff and the call spread is very high (> 0.99), yielding variance reduction ratios of 50--100x.
(c) Finite differences on the Heston PDE:
The FDM approach uses the same Heston PDE and Craig-Sneyd ADI scheme, but with the terminal condition changed to:
The discontinuity in the terminal condition causes poor convergence. To mitigate this:
-
Grid alignment: Ensure the log-strike \(\ln K\) falls exactly on a grid point, so the discontinuity is at a grid boundary rather than between grid points.
-
Smoothing: Replace the step function with a smoothed version (e.g., \(\frac{1}{2} + \frac{1}{2}\text{erf}((x - \ln K)/\epsilon)\)) or use the cell-averaging technique: at \(t = T\), assign the payoff in the cell containing \(\ln K\) as the integral of the step function over the cell width.
-
Grid refinement: Use a finer \(x\)-grid near \(K\) (the sinh-based stretching already concentrates points near ATM).
With these techniques, the FDM achieves 3--4 digits of accuracy for the digital price with the same grid as used for vanilla options.