Search code examples
pythonmatlabcythonnumerical-integrationquad

Python vs. MATLAB computing an integral to infinity with different results, alternative (i.e. expand Gauss-Legendre quadrature to -x-> Infinity)?

I am getting inconsistent results between MATLAB's quadgk and Python's quad routine for an integral from (-x or 0) -> infinity. I believe the MATLAB version is correct (based on a sense check of switching the flag parameter from 1 to -1) whereas the Python version gives erroneous results, in this case 0. MATLAB produces 0.1022. The integrands are identical and I have tied out each step of the way, even inserting the x values generated by MATLAB's quadgk into Python (which results in the Python version generating the same value as MATLAB, just passing them to the integrand function). At this point I'm looking to use another routine rather than SciPy such as Gauss-Legendre quadrature here https://sourceforge.net/projects/fastgausslegendrequadrature/ but I am not sure how to extend it from its a/b range into -a->infinity (I've seen these methods that only go to a finite number:

https://upload.wikimedia.org/math/a/7/f/ Different intervals for Gauss-Legendre quadrature in numpy whereas b=np.Inf results in NaN. Also not sure how to setup the integration from the returned nodes and weights, although I've been reading the transformation but only for a and b finite ranges: https://pomax.github.io/bezierinfo/legendre-gauss.html Either that or if someone knows a Python library that can handle this - I don't really like the fact quad isn't vectorized though and will likely write the code in Cython, as I have to integrate 600,000 functions quickly (i.e. link to the C++ library link above). What's really odd here is that I've managed to get identical results by moving up the vol input anywhere >= 0.39, below that Python's results collapse at 0. Very confusing. Any help is appreciated, it has been years since Calculus... here is the Python code:

from scipy.stats import norm, lognorm
from scipy.integrate import quad
import numpy as np

def integrand(x, flag, F, K, vol, T2, T1):
    d1 = (np.log(x / (x+K)) + 0.5 * (vol**2) * (T2-T1)) / (vol * np.sqrt(T2 - T1))
    d2 = d1 - vol*np.sqrt(T2 - T1)
    mu = np.log(F) - 0.5 *vol **2 * T1
    sigma = vol * np.sqrt(T1)
    value = lognorm.pdf(x, scale=np.exp(mu), s=sigma) * (flag * x*norm.cdf(flag * d1) - flag * (x+K)*norm.cdf(flag * d2))
    return value

if __name__ == '__main__':
        flag = 1
        F = 54.31
        K = 1.1967
        vol = 0.1328
        T2 = 0.0411
        T1 = 0.0137
        quad(integrand, 0, np.Inf, args=(flag, F, K, vol, T2, T1), epsabs=1e-12)[0]

And here is the MATLAB code (have to save the integrand as a .M then can enter the script in a command window):

function value = integrand(x, flag, F,K,vol,T2,T1)
d1 = (log(x ./ (x+K)) + 0.5 .* (vol.^2) .* (T2-T1)) ./ (vol .* sqrt(T2 - T1));
d2 = d1 - vol.*sqrt(T2 - T1);
mu = log(F) - 0.5 .*vol .^2 .* T1;
sigma = vol .* sqrt(T1);
value = lognpdf(x, mu, sigma) .* (flag .* x.*normcdf(flag .* d1) - flag .* (x+K).*normcdf(flag .* d2));
end

% script part

flag = 1
F = 54.31
K = 1.1967
vol = 0.1328
T2 = 0.0411
T1 = 0.0137
quadgk(@(x) integrand(x,flag, F, K, vol, T2, T1), 0, Inf, 'AbsTol',1e-12)

I should note that MATLAB and Python generate the same results using quad when these inputs are passed instead (transposed above variables):

 current_opt = [  -1.0000    1.2075    0.1251    0.4300    0.0685    0.0411     
                1.0000    1.2075    0.0512    0.5600    0.0685    0.0411]
Solution

  • Okay this is interesting. The integration falls apart unless the epsabs variable is set ridiculously high. I have managed to replicate results between MATLAB and Python using epsabs=-1e1000. Slow as can be but at least it works.