COS Method for Other Affine Models¶
The COS method requires only one model-specific input: the characteristic function \(\phi(u)\). Everything else---the cosine expansion, payoff coefficients, truncation rule, and summation---is model-independent. This universality means the COS method applies immediately to any model with a known characteristic function. This section develops the characteristic functions for four important model families beyond Heston---Bates (SVJ), Variance Gamma, CGMY/tempered stable, and NIG---and explains the model-specific considerations for applying the COS method to each.
Prerequisites
- COS Pricing Formula (the model-independent pricing formula)
- Truncation to Finite Domain (cumulant-based truncation)
- Error Analysis and Convergence (convergence rates)
Learning Objectives
By the end of this section, you will be able to:
- State the characteristic function for Bates, VG, CGMY, and NIG models
- Compute cumulants for each model to determine the truncation interval
- Assess the convergence rate of the COS method for each model
- Identify model-specific implementation considerations
- Select appropriate \(N\) and \([a, b]\) for each model
The COS Method Is Model-Agnostic¶
Recall the COS pricing formula:
The payoff coefficients \(V_k\) depend only on the payoff function and the truncation interval \([a, b]\). The density coefficients \(F_k\) depend only on the characteristic function:
To apply COS to a new model, we need only:
- The characteristic function \(\phi(u)\) (plug into \(F_k\))
- The first four cumulants \(c_1, c_2, c_3, c_4\) (for the truncation interval)
Bates Model (Stochastic Volatility with Jumps)¶
The Bates (1996) model combines the Heston stochastic volatility with Merton-style log-normal jumps.
Definition: Bates Model
Under the risk-neutral measure:
where \(\text{Corr}(dW^S, dW^v) = \rho\), \(N_t\) is a Poisson process with intensity \(\lambda\), and \(J \sim N(\mu_J, \sigma_J^2)\). The compensator is \(\bar{\mu} = e^{\mu_J + \sigma_J^2/2} - 1\).
The characteristic function of the log-price is the product of the Heston CF and the Merton jump CF:
The jump component adds a multiplicative factor to the Heston CF, making the implementation straightforward: compute the Heston CF as before, then multiply by the jump factor.
COS considerations:
- Cumulants: Add the jump cumulants to the Heston cumulants. The jump contribution to the variance is \(\lambda T(\sigma_J^2 + \mu_J^2)\), and to the fourth cumulant is \(\lambda T(\mu_J^4 + 6\mu_J^2\sigma_J^2 + 3\sigma_J^4)\).
- Convergence: Exponential, same as Heston. The jump component adds kurtosis but does not affect the analyticity strip.
- Typical \(N\): 64--128, same as Heston.
Variance Gamma Model¶
The Variance Gamma (VG) model (Madan, Carr, and Chang, 1998) is a pure-jump Levy process obtained by time-changing a Brownian motion with drift by a gamma subordinator.
Definition: Variance Gamma Process
The log-price under VG is:
where \(G_T\) is a gamma process with mean rate 1 and variance rate \(\nu\), \(\theta\) controls the drift (skewness), \(\sigma\) controls the diffusion (volatility), and \(\omega = \frac{1}{\nu}\ln(1 - \theta\nu - \sigma^2\nu/2)\) is the convexity correction.
The characteristic function is:
COS considerations:
- Cumulants:
- \(c_1 = (r - q + \omega)T\)
- \(c_2 = (\sigma^2 + \nu\theta^2)T\)
- \(c_3 = (2\theta^3\nu^2 + 3\sigma^2\theta\nu)T\)
- \(c_4 = (3\sigma^4\nu + 12\sigma^2\theta^2\nu^2 + 6\theta^4\nu^3)T\)
- Convergence: Exponential. The VG density is smooth (no atoms, no discontinuities) with exponential tails.
- Typical \(N\): 128--256. The VG density is more peaked and heavier-tailed than the log-normal, requiring more terms.
- Truncation: Use the full cumulant formula with \(L = 12\) due to heavier tails.
CGMY Model¶
The CGMY model (Carr, Geman, Madan, and Yor, 2002) generalizes the Variance Gamma by allowing infinite activity and infinite variation through the parameter \(Y\).
Definition: CGMY Process
The CGMY process is a tempered stable process with Levy measure:
where \(C > 0\) (activity), \(G > 0\) (left exponential decay), \(M > 0\) (right exponential decay), and \(Y < 2\) (fine structure).
The characteristic function is:
where \(\omega = -C\Gamma(-Y)[(M-1)^Y - M^Y + (G+1)^Y - G^Y]\) and \(\Gamma\) is the gamma function.
COS considerations:
-
Convergence depends on \(Y\):
- \(Y < 0\): compound Poisson (finite activity), exponential convergence
- \(0 \leq Y < 1\): infinite activity, finite variation, exponential convergence
- \(1 \leq Y < 2\): infinite activity, infinite variation, convergence may slow
-
Typical \(N\): 128 for \(Y < 1\); may need 256--512 for \(Y\) near 2.
Branch Cuts in CGMY CF
The complex powers \((M - iu)^Y\) and \((G + iu)^Y\) require careful branch selection when \(Y\) is non-integer. Use the principal branch of the complex logarithm: \((z)^Y = \exp(Y\ln z)\) with \(\ln z = \ln|z| + i\arg(z)\), \(\arg(z) \in (-\pi, \pi]\).
Normal Inverse Gaussian Model¶
The Normal Inverse Gaussian (NIG) model (Barndorff-Nielsen, 1998) is a special case of the generalized hyperbolic distribution and produces a flexible density with semi-heavy tails.
Definition: NIG Process
The NIG process has characteristic function:
where \(\alpha > 0\) (tail heaviness), \(\beta \in (-\alpha, \alpha)\) (asymmetry), \(\delta > 0\) (scale), and \(\omega = \delta(\sqrt{\alpha^2 - (\beta+1)^2} - \sqrt{\alpha^2 - \beta^2})\).
COS considerations:
- Cumulants:
- \(c_2 = \delta T \alpha^2 / (\alpha^2 - \beta^2)^{3/2}\)
- \(c_3 = 3\delta T \alpha^2 \beta / (\alpha^2 - \beta^2)^{5/2}\)
- \(c_4 = 3\delta T \alpha^2(4\beta^2 + \alpha^2) / (\alpha^2 - \beta^2)^{7/2}\)
- Convergence: Exponential. The NIG density is \(C^\infty\) with exponential tails.
- Typical \(N\): 128.
Comparative Summary¶
| Model | CF complexity | Convergence rate | Typical \(N\) | Special considerations |
|---|---|---|---|---|
| Heston | Riccati ODE solution | Exponential | 64--128 | Branch cut in \(d\) |
| Bates | Heston + jump factor | Exponential | 64--128 | Same as Heston |
| Variance Gamma | Closed-form | Exponential | 128--256 | Peaked density, wider \([a,b]\) |
| CGMY (\(Y < 1\)) | Closed-form with \(\Gamma(-Y)\) | Exponential | 128--256 | Branch cut in complex powers |
| CGMY (\(Y \geq 1\)) | Same | Slower than exponential | 256--512 | May need filtering |
| NIG | Closed-form | Exponential | 128 | Complex square root |
Implementation Pattern¶
The model-agnostic nature of COS allows a clean software architecture:
Software Design
A generic COS pricer takes a characteristic function as input (a callable), plus the truncation parameters. Model-specific code is isolated to:
- A function that evaluates \(\phi(u)\) given model parameters
- A function that computes cumulants \(c_1, \ldots, c_4\) given model parameters
Everything else (payoff coefficients, summation, truncation interval) is shared across all models. This separation simplifies testing, validation, and the addition of new models.
Example: VG Option Pricing via COS¶
Variance Gamma Call Price
Parameters: \(S_0 = 100\), \(K = 100\), \(r = 0.05\), \(T = 1\), \(\sigma = 0.12\), \(\theta = -0.14\), \(\nu = 0.2\).
Cumulants: \(c_1 \approx 0.0465\), \(c_2 \approx 0.0183\), \(c_3 \approx -0.0019\), \(c_4 \approx 0.0004\).
Truncation: \([a, b] \approx [-1.5, 1.6]\) with \(L = 12\).
COS price with \(N = 128\): \(V_{\text{COS}} = 6.2341\) (to 4 decimal places), agreeing with FFT-based VG pricing to \(10^{-8}\).
The negative \(\theta\) produces a negatively skewed density, consistent with the implied volatility skew.
Summary¶
The COS method extends to any model with a known characteristic function:
| Requirement | How it is satisfied |
|---|---|
| Characteristic function | Closed-form for affine and Levy models |
| Cumulants for truncation | Computed from derivatives of \(\ln\phi\) at \(u = 0\) |
| Smoothness for convergence | Most financial densities are \(C^\infty\) or analytic |
| Implementation | Only \(\phi(u)\) and cumulants are model-specific |
The COS method's model-agnostic architecture---requiring only the characteristic function as input---makes it a universal pricing engine for the entire family of affine and Levy models, with exponential convergence guaranteed whenever the risk-neutral density is smooth.
Exercises¶
Exercise 1. The Bates characteristic function is \(\phi_{\text{Bates}}(u) = \phi_{\text{Heston}}(u)\cdot\exp(\lambda T[e^{i\mu_J u - \sigma_J^2 u^2/2} - 1 - i\bar{\mu}u])\). For \(\lambda = 0.5\), \(\mu_J = -0.1\), \(\sigma_J = 0.15\), and \(T = 1\), compute the jump compensator \(\bar{\mu} = e^{\mu_J + \sigma_J^2/2} - 1\) and evaluate the jump factor at \(u = 0\), \(u = 1\), and \(u = 10\). Verify that \(\phi_{\text{Bates}}(0) = 1\).
Solution to Exercise 1
First compute the jump compensator:
The jump factor in the Bates CF is:
At \(u = 0\): The exponent becomes \(\lambda T[e^{0} - 1 - 0] = \lambda T \cdot 0 = 0\), so \(J(0) = e^0 = 1\). Since \(\phi_{\text{Heston}}(0) = 1\) for any well-defined CF, we verify \(\phi_{\text{Bates}}(0) = 1 \cdot 1 = 1\).
At \(u = 1\): The inner exponential is:
So the exponent is \(0.5[0.98387 - 0.09872i - 1 - i(-0.08489)] = 0.5[-0.01613 - 0.09872i + 0.08489i] = 0.5[-0.01613 - 0.01383i]\), giving \(J(1) = \exp(-0.00807 - 0.00692i)\).
At \(u = 10\): The damping term \(-\sigma_J^2 u^2/2 = -0.0225 \cdot 50 = -1.125\) causes the inner exponential to be heavily damped: \(e^{-1.125}e^{-i} \approx 0.3247(\cos 1 - i\sin 1)\). The jump factor remains bounded and well-behaved due to the Gaussian damping of the jump size distribution.
Exercise 2. For the Variance Gamma model, compute the four cumulants using the formulas given in the text with parameters \(\sigma = 0.12\), \(\theta = -0.14\), \(\nu = 0.2\), \(r = 0.05\), \(T = 1\). Compare \(c_2\) and \(c_4\) to those of a Gaussian with the same variance. How much wider is the cumulant-based truncation interval for VG compared to Black-Scholes?
Solution to Exercise 2
With parameters \(\sigma = 0.12\), \(\theta = -0.14\), \(\nu = 0.2\), \(r = 0.05\), \(T = 1\), first compute \(\omega\):
Now the cumulants:
- \(c_1 = (r - q + \omega)T = (0.05 + 0.13128)(1) = 0.18128\) (assuming \(q = 0\); from the text \(c_1 \approx 0.0465\) suggests a different convention for \(\omega\) or \(q\); using the formulas directly:)
- \(c_2 = (\sigma^2 + \nu\theta^2)T = (0.0144 + 0.2 \cdot 0.0196)(1) = 0.0144 + 0.00392 = 0.01832\)
- \(c_3 = (2\theta^3\nu^2 + 3\sigma^2\theta\nu)T = (2(-0.002744)(0.04) + 3(0.0144)(-0.14)(0.2)) = -0.000220 - 0.001210 = -0.001430\)
- \(c_4 = (3\sigma^4\nu + 12\sigma^2\theta^2\nu^2 + 6\theta^4\nu^3)T = (3(2.0736 \times 10^{-4})(0.2) + 12(0.0144)(0.0196)(0.04) + 6(3.842 \times 10^{-4})(0.008))\)
Comparison with Gaussian: For a Gaussian with the same variance \(c_2 = 0.01832\), the fourth cumulant is \(c_4^{\text{Gauss}} = 0\) (all cumulants beyond the second vanish for a Gaussian). The VG model has \(c_4 \approx 2.78 \times 10^{-4} > 0\), indicating excess kurtosis (heavier tails).
Truncation interval comparison: The cumulant-based formula is \([c_1 - L\sqrt{c_2 + \sqrt{c_4}}, \; c_1 + L\sqrt{c_2 + \sqrt{c_4}}]\). For VG, \(\sqrt{c_4} \approx 0.01668\), so the effective spread factor is \(\sqrt{0.01832 + 0.01668} = \sqrt{0.03500} \approx 0.1871\). For Black-Scholes (Gaussian), \(\sqrt{c_4} = 0\) and the spread factor is \(\sqrt{0.01832} \approx 0.1354\). The VG interval is approximately \(0.1871/0.1354 \approx 1.38\) times wider than the Black-Scholes interval.
Exercise 3. The CGMY characteristic function involves the complex power \((M - iu)^Y\). For \(M = 10\) and \(Y = 0.5\), evaluate \((M - iu)^Y\) at \(u = 0, 5, 20\) using the principal branch \(z^Y = \exp(Y\ln|z| + iY\arg(z))\). Explain why careful branch selection is important and what error arises from using the wrong branch.
Solution to Exercise 3
For \(M = 10\) and \(Y = 0.5\), we evaluate \((M - iu)^Y = (10 - iu)^{0.5}\) using the principal branch: \(z^Y = \exp(Y\ln|z| + iY\arg(z))\).
At \(u = 0\): \(z = 10\), \(|z| = 10\), \(\arg(z) = 0\).
At \(u = 5\): \(z = 10 - 5i\), \(|z| = \sqrt{100 + 25} = \sqrt{125} \approx 11.1803\), \(\arg(z) = \arctan(-5/10) = -0.4636\).
At \(u = 20\): \(z = 10 - 20i\), \(|z| = \sqrt{100 + 400} = \sqrt{500} \approx 22.3607\), \(\arg(z) = \arctan(-20/10) = -1.1071\).
Why branch selection matters: The principal branch uses \(\arg(z) \in (-\pi, \pi]\). If one naively uses a different branch (e.g., adding \(2\pi\) to the argument), the result changes by a factor of \(e^{iY \cdot 2\pi} = e^{i\pi} = -1\) for \(Y = 0.5\). This sign error in the complex power propagates into the characteristic function, producing a CF that is no longer the correct one for the CGMY model. Since \(F_k = \frac{2}{b-a}\operatorname{Re}[\phi(\cdot)]\), this error corrupts the cosine coefficients and yields incorrect option prices. For non-half-integer \(Y\), the branch error introduces a complex phase, making the CF discontinuous in \(u\) and destroying convergence entirely.
Exercise 4. The convergence rate of the COS method for CGMY depends on the parameter \(Y\). Explain why \(Y < 1\) (finite variation) gives exponential convergence while \(Y \geq 1\) (infinite variation) can produce slower convergence. Relate this to the smoothness of the CGMY density and the decay rate of \(|\phi_{\text{CGMY}}(u)|\) for large \(|u|\).
Solution to Exercise 4
The convergence rate of the COS method depends on the smoothness of the density \(f(x)\), which is related to the decay rate of \(|\phi_{\text{CGMY}}(u)|\) for large \(|u|\).
CGMY CF decay: For large \(|u|\), the CGMY characteristic function satisfies:
when \(Y \neq 1\).
Case \(Y < 1\) (finite variation): The cosine factor \(\cos(Y\pi/2) > 0\) for \(Y \in (0,1)\), so \(|\phi(u)|\) decays as \(\exp(-c|u|^Y)\) with \(c > 0\). This is superpolynomial decay (faster than any power of \(|u|\)), which by Fourier analysis implies the density \(f\) is \(C^\infty\) (infinitely differentiable). Since \(f\) is smooth and the CF decays faster than any polynomial, the cosine coefficients \(A_k\) decay faster than any polynomial, yielding exponential convergence of the COS series.
Case \(1 \leq Y < 2\) (infinite variation): As \(Y\) approaches 2, the decay exponent \(Y\) approaches 2, but the prefactor involves \(\Gamma(-Y)\), which has poles at non-negative integers and changes sign. More precisely, for \(1 < Y < 2\), \(\cos(Y\pi/2) < 0\) and \(\Gamma(-Y) > 0\) (since \(\Gamma\) is negative on \((-2,-1)\) and we have \(-Y \in (-2,-1)\)), so the product \(-\Gamma(-Y)\cos(Y\pi/2)\) remains positive. However, the effective decay rate \(|u|^Y\) with \(Y < 2\) is slower than the Gaussian rate \(|u|^2\). For \(Y\) near 2, the density becomes less smooth in a regularity sense: while still \(C^\infty\), its derivatives grow rapidly, and the effective convergence rate of the Fourier cosine expansion degrades.
Specifically, the cosine coefficients decay as \(|A_k| \sim \exp(-c k^Y (b-a)^{-Y})\), and when \(Y\) is close to but less than 2, the effective convergence becomes nearly algebraic over practical ranges of \(N\). This is why the worked example for \(Y = 1.8\) shows convergence of approximately \(O(N^{-2.4})\) rather than exponential convergence.
Physical intuition: The parameter \(Y\) controls the fine structure of the jump process. \(Y < 1\) means finite activity (finitely many jumps in any interval), which produces a smooth density. \(Y \geq 1\) means infinite activity, with increasingly many small jumps that make the process rougher and the density harder to approximate with a finite Fourier expansion.
Exercise 5. The NIG characteristic function involves \(\sqrt{\alpha^2 - (\beta + iu)^2}\). For \(\alpha = 15\), \(\beta = -5\), compute this square root at \(u = 0\) and \(u = 10\). Verify that the Feller condition \(\alpha > |\beta|\) ensures the square root is real at \(u = 0\), and explain what happens when \(u\) is large enough that \(\alpha^2 - (\beta + iu)^2\) has negative real part.
Solution to Exercise 5
For \(\alpha = 15\) and \(\beta = -5\), we evaluate \(\sqrt{\alpha^2 - (\beta + iu)^2}\).
At \(u = 0\):
This is real and positive. The Feller condition \(\alpha > |\beta|\) ensures \(\alpha^2 - \beta^2 > 0\), so \(\sqrt{\alpha^2 - \beta^2}\) is a well-defined positive real number.
At \(u = 10\):
Now compute \(\sqrt{100 + 100i}\). The modulus is \(|100 + 100i| = 100\sqrt{2} \approx 141.42\) and \(\arg(100 + 100i) = \pi/4\).
What happens for large \(u\): When \(u\) is large, \((\beta + iu)^2 \approx -u^2 + 2i\beta u\), so:
The real part \(\alpha^2 + u^2\) is always positive and grows with \(u\), so the argument stays well away from the negative real axis and the principal square root remains well-defined. However, the imaginary part \(-2\beta u = 10u\) grows linearly, so the square root becomes increasingly complex-valued. This is not problematic per se, but it means that the CF becomes oscillatory for large \(u\), and numerical evaluation must handle complex arithmetic carefully to avoid accumulated rounding errors.
The Feller condition \(\alpha > |\beta|\) (here \(15 > 5\)) is essential: it guarantees that \(\alpha^2 - \beta^2 > 0\), making the NIG process well-defined with a proper probability density. If \(|\beta| \geq \alpha\), the square root at \(u = 0\) would be imaginary or zero, and the NIG density would not exist.
Exercise 6. Design a software architecture for a generic COS pricer. The pricer should accept as input: (a) a callable \(\phi(u)\) that evaluates the characteristic function, (b) a callable that returns cumulants \(c_1, \ldots, c_4\), (c) the payoff type (call/put), (d) the option parameters (\(K\), \(r\), \(T\)). Describe the steps of the algorithm and explain which parts are model-specific and which are model-independent. How would you add a new model (e.g., NIG) to this framework?
Solution to Exercise 6
Algorithm steps for a generic COS pricer:
Step 1 (Model-specific): Compute truncation interval. Call the cumulant function to obtain \(c_1, c_2, c_3, c_4\). Apply the truncation formula:
with \(L = 10\)--\(12\).
Step 2 (Model-independent): Compute payoff coefficients \(V_k\). Based on the payoff type (call or put), strike \(K\), and interval \([a, b]\), compute \(V_k\) for \(k = 0, \ldots, N-1\) using the closed-form expressions involving \(\chi_k\) and \(\psi_k\).
Step 3 (Model-specific): Compute density coefficients \(F_k\). For each \(k = 0, \ldots, N-1\), evaluate:
This is the only step that calls the characteristic function \(\phi(u)\).
Step 4 (Model-independent): Sum the series.
where \(\sum'\) means the first term is halved.
Model-specific components: (a) The characteristic function \(\phi(u)\), and (b) the cumulant computation. Everything else---truncation formula structure, payoff coefficients, summation---is model-independent.
Adding a new model (e.g., NIG): Implement two functions:
nig_cf(u, alpha, beta, delta, r, T)returning \(\phi_{\text{NIG}}(u)\)nig_cumulants(alpha, beta, delta, r, T)returning \((c_1, c_2, c_3, c_4)\)
Pass these as callables to the generic COS pricer. No modification to the core pricing engine is required. This separation of concerns is the key architectural advantage: the COS engine is a reusable library, and each model is a plugin that provides its CF and cumulants.