Interface-aligned Chebyshev integral-tau methods for occupation-time and Parisian option pricing

A numerical analysis study has developed and tested a multidomain Chebyshev integral-tau method for pricing occupation-time and cumulative-Parisian options. The method aligns the monitoring level and strike with element interfaces, uses an exact treatment of occupation-time transport, and advances the remaining diffusion with a backward-Euler-damped, L-stable TR-BDF2 scheme. In experiments using a European call with identical market data and polynomial degrees, interface alignment reduced the price error from 1.12 for a single global element to 8.24 × 10^-8 at degree N = 48, a reported gain of 1.4 × 10^7.

The work concerns a single-asset mathematical-finance model with constant interest rate and volatility; it is a numerical-methods study, not an assessment of live trading performance or market outcomes. Occupation-time derivatives and Parisian options depend on how long an asset remains beyond a monitoring level, rather than only on whether the asset touches that level. In the cumulative-Parisian down-and-out call examined in the study, the claim pays (S_T − K)+ at maturity if the accumulated time spent at or below a barrier remains below a prescribed window; it becomes worthless when that occupation time reaches the window.

The valuation therefore requires a two-dimensional state consisting of the asset level and accumulated occupation time. The associated degenerate-parabolic equation diffuses in the log-price direction while transporting in the occupation-time direction, and the monitoring level switches the occupation clock on and off. Such contracts were introduced to remove the knock-out discontinuity of standard barrier options and make manipulation near the barrier harder, but their path dependence makes numerical pricing more demanding.

The main spatial difficulty is that the value function is non-smooth at two locations that are known before the calculation: the monitoring level, where the occupation source switches and the second derivative of the value can jump, and the strike, where the terminal payoff has a kink. Classical finite-difference schemes can handle such kinks but generally converge algebraically unless the mesh is refined locally. Global spectral methods can converge geometrically for analytic data, yet an interior kink or jump causes a Gibbs or corner defect and removes that advantage.

The proposed method instead divides the log-price domain so that both financial singularities lie on element faces. On each open element, the transformed data and the relevant right-hand side are assumed to extend analytically to a Bernstein ellipse, allowing Chebyshev coefficients to decay geometrically. The singularity is not removed; it is confined to interfaces and handled by a small number of algebraic tau conditions.

The method combines this interface alignment with an integration, or primitive-matrix, reformulation of spectral discretisation. Instead of building dense differentiation matrices for the second derivative, the formulation uses banded Chebyshev antiderivatives. An exact identity shows that the defect in the twice-integrated second-derivative operator acts only on the first two coefficient rows of each element.

Those rows are replaced by boundary and transmission conditions, while the remaining implicit spatial operator is block-pentadiagonal. Each element contributes two independent integration constants, and the resulting system is closed by exactly 2m tau rows, where m is the number of aligned elements. A direct argument shows that these rows are independent on the primitive nullspace: a piecewise affine function that satisfies zero external values and continuous value and flux at every internal face must be identically zero.

This independence establishes the well-posedness of the assembled tau system. Time stepping is designed to respect the special structure of the occupation clock. The transformed generator separates into diffusion in log price and transport in accumulated occupation time.

A Strang splitting alternates the two operations, and the transport step is exact because it reduces to an index shift through the occupation grid. Consequently, the clock transport introduces no separate time-discretisation error. The diffusion is advanced with TR-BDF2 after two backward-Euler startup steps that damp modes excited by the payoff kink.

The scheme is second order, A-acceptable and L-stable: as the argument moves deeply into the left half-plane, its amplification factor tends to zero, whereas the Crank–Nicolson factor approaches a magnitude of one. This difference matters for non-smooth data because L-stability removes stiff modes rather than leaving them weakly damped. In the reported time-convergence test, TR-BDF2 achieved an observed order of 2.00 at 320 time steps, while Crank–Nicolson remained first order.

The authors also analyse stability at continuous and discrete levels. An energy estimate shows that the coupled continuous operator is dissipative under the stated domain, boundary and transmission conditions, and therefore generates a contraction semigroup. At the discrete level, the theory gives unconditional power boundedness only under a stated hypothesis: the discrete generator must be diagonalisable, with its eigenvalues in the left half-plane, and its eigenvector conditioning must remain bounded for a uniform-in-degree conclusion.

The numerical tests found that the TR-BDF2 propagator eigenvalues lay in the unit disk and that its spectral radius remained 0.999996 across the tested polynomial degrees, from N = 8 to N = 128. Those results certify power boundedness of the propagator in the reported tests, not by themselves the left-half-plane location of the generator spectrum. A direct eigenvalue computation at N = 32 found every listed generator eigenvalue to have a real part no greater than −3 × 10^-3, while the paper also presents a discrete energy identity for the conforming space.

The computational efficiency comes from a Sherman–Morrison–Woodbury correction applied to the bordered-banded system. For a fixed number of interfaces, the banded core has fixed bandwidth and the finite-rank correction has size independent of the polynomial degree, giving a theoretical operation count linear in the number of coefficients. Timing experiments showed that the advantage increased as the system grew: at 2050 degrees of freedom, the bordered-banded solve took 2.97 × 10^-4 seconds, compared with 2.73 × 10^-2 seconds for dense LU factorization, a reported speed-up of 92.

The trade-off is conditioning. Nevertheless, the tested pricing degrees retained six to eight significant digits because, for analytic data, the most sensitive high modes carried exponentially small weight. Checks against independent calculations supported the method’s accuracy.

For the occupation-time down-and-out call with S_0 = 100, K = 100, B = 90, r = 0.05, σ = 0.2, maturity T = 1 and occupation window D = 0.1, the aligned tau method produced a price of 9.93994 at N = 48 and M = 800. A second-order finite-difference reference gave 9.94901, while a Monte Carlo estimate based on two million paths gave 9.93347 ± 0.02047. The aligned result, the finite-difference result and the Monte Carlo estimate were consistent within the reported Monte Carlo uncertainty.

Changing the computational truncation width from 7 to 12 altered the price by less than 10^-3, supporting the domain-truncation estimate. A Richardson limit of 9.94081 was also reported. The sensitivity results also improved: at a coarse time step, maximum Gamma error was 2.22 × 10^-5 for TR-BDF2 compared with 2.62 × 10^-3 for Crank–Nicolson, and the aligned Black–Scholes Gamma at S = 110 fell from 4.7 × 10^-4 at degree 12 to 3.7 × 10^-10 at degree 48.

The paper identifies limitations in the scope of the analysis. The pricing model uses one asset, constant coefficients and a single occupation grid tied to the time step through exact-shift splitting. The authors discuss local or stochastic volatility, adaptive tracking of moving interfaces and implicit occupation transport as possible extensions, but these are proposals rather than results established by the study.

The authors also caution that an equilibrated or hybrid preconditioner may be needed beyond the degrees used for the reported pricing calculations. Within the tested contract and model scope, the findings show that known financial singularities can be handled as interface conditions, that exact clock transport can separate path dependence from time-stepping error, and that a banded integral formulation can preserve practical accuracy while reducing the per-step cost relative to dense factorization. The study’s central evidence is numerical and theoretical within that defined model, rather than a claim about the value of any particular derivative contract in financial markets.

Subject of Research: Numerical methods for pricing single-asset occupation-time and cumulative-Parisian options with known payoff and monitoring singularities.

Article Title: Interface-aligned Chebyshev integral-tau methods for occupation-time and Parisian option pricing

Article References: Pindza, E., & Bambe Moutsinga, C. R. (2026). Interface-aligned Chebyshev integral-tau methods for occupation-time and Parisian option pricing. Frontiers in Applied Mathematics and Statistics, 12, 1944707. https://doi.org/10.3389/fams.2026.1944707

Image Credits: AI Generated

DOI: 10.3389/fams.2026.1944707

Keywords: mathematical finance, occupation-time options, Parisian options, Chebyshev spectral methods, integral-tau discretization, TR-BDF2 stability, path-dependent PDEs, Sherman–Morrison–Woodbury solver, option-pricing accuracy, Interface-aligned, Chebyshev, integral-tau