A B-Spline Galerkin Method for Option Pricing under the CGMY Tempered Stable Process

We develop and analyze a B-spline scaling-function Galerkin method for the tempered fractional partial differential equation governing European and American option prices under the CGMY process. The spatial generator couples exponentially tempered left- and right-sided Riemann–Liouville derivatives of order $Y\in (0,2)$, whose nonlocality and asymmetry present discretization challenges. We derive a corrected Fourier-symbol construction of the tempered connection-coefficient matrices whose interior blocks are Toeplitz, maintain exact fractional scaling, and enable $\mathcal{O}(N\log N)$ matrix–vector products via the fast Fourier transform. In the regime $Y\in (1,2)$, we prove full-line and finite-domain Gårding inequalities, coercive whenever the risk-free rate is positive, via zero-extension to the full-line symbol. On this foundation, we establish elliptic-projection, semidiscrete, and fully discrete error estimates with spatial rate $\mathcal{O}(h^{p-Y\mathrm{/}2})$, alongside an unconditional shifted-energy stability estimate for the forced Crank–Nicolson scheme with mesh-independent constants. For positive rates the admissible shift vanishes, governing the plain scheme exactly as implemented. A single backward solve prices an entire strike strip at fixed maturity with smooth, analytically differentiated Greeks. The same discretization prices American puts through a projected linear complementarity solve, with early-exercise premiums and exercise boundaries validated against a resolved deterministic reference and a Longstaff–Schwartz Monte Carlo benchmark. Calibration to S&P~500 option surfaces across low-volatility, crisis, and recovery regimes yields implied-volatility errors near two vol-points on well-identified surfaces, while an accuracy-matched comparison with a stabilized Grünwald–Letnikov baseline delineates precisely when the high-order variational structure pays off.