Fractional-Parabolic Deformations with Sinh-Acceleration

Abstract

We develop a general simple methodology for very fast and accurate evaluation of integrals of functions that admit analytic continuation into wide regions of the complex plane or into appropriate Riemann surfaces. We use a family of fractional-parabolic deformations of the contours of integration to appropriate Riemann surfaces, make the corresponding conformal changes of variables and use the simplified trapezoid rule. An additional change variables leads to numerical schemes of complexity of the order of $A\pi^{-2}E^2$, where $E=\ln(1/\eps)$, $\eps$ is the error tolerance, and $A$ is, typically, in the range $[0.5,3]$, or even $A\pi^{-2}E \ln E$ (the Laplace and Fourier inversion, the Laplace transform and other integrals with exponentially decreasing or oscillating factors). Simple efficient error bounds make adaptive procedures unnecessary. As applications, we consider special functions of several classes, and outline applications to the Wiener-Hopf factorization, calculation of various probability distributions and pricing vanillas and exotics in the L'evy models, Heston model, and more general stochastic volatility models and affine models. The method can be applied to various boundary problems for partial differential and pseudo-differential equations, and in essentially all situations when the saddle point method or reduction to a cut are applied. Simple approximate recommendations for the choice of parameters of the scheme lead to programs that are accurate and competitive with Mathematica and Python, even without parallelization. In a number of situations, for the error tolerance $10^{-15}$, the resulting programs are dozens, hundreds and even thousands times faster than functions implemented in MATLAB, 8-15 times faster than SymPy. For $ _1F_1$ and $ _2F_1$, in some regions in the parameter space, the new programs are faster than functions implemented in Mathematica. Pricing of vanillas in the Heston model requires the summation of 6-15 terms for options of maturities T ≥ 0.004, in a wide region of strikes; the number of terms is a decreasing function of $T$.