Cara menggunakan bivariate normal distribution python

I had to write an option model that was using a bivariate distribution in Python. However, I did not find a prebuilt function that was fast - some seem to be using the random scipy generator to emulate it with the multivariate function. BUT... if you really dig deep and see what the other financial packages are using, it's pretty much ALL code written by one guy, Alan Genz from the University of Washington. He pretty much writes everything in Fortan or MATLAB, and that's it. So if you look into packages that have the bivariate CDF you'll find his name and his code there (I found an old version in MATLAB actually). http://www.math.wsu.edu/faculty/genz/software/software.html

So why this isn't built into SciPy or NumPy already, I have no idea. But I rewrote it in several hours, using both MATLAB and Python to check the resulting code. He's solving with Gauss Legendre quadrature so this will always be much faster than a solution that uses a random number generator: https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature

#       Alan Genz - original MATLAB code converted by Matt Slezak to Python
#       http://www.math.wsu.edu/faculty/genz/software/matlab/bvnl.m
#       Department of Mathematics
#       Washington State University
#       This is the bivariate CDF distribution
#
#       dh 1st upper integration limit
#       dk 2nd upper integration limit
#       r   correlation coefficient

import numpy as np
from scipy.stats import norm

def bivariate_cdf(dh, dk, r):
    # dh and dk get signs flipped right away
    dh = float(-dh)
    dk = float(-dk)

    if (dh == np.inf) or (dk == np.inf): 
        return 0
    
    else:
        if dh == - np.inf:
            if dk == - np.inf:
                return 1
            else:
                return norm.cdf(- dk)
        if dk == - np.inf:
            return norm.cdf(- dh)
        else:
            if r == 0:
                return norm.cdf(- dh)*norm.cdf(- dk)
            else:
                    tp=2*np.pi
                    h=dh
                    k=dk
                    hk=h*k
                    bvn=0
                    if abs(r) < 0.3:
                        w=np.array([0.1713244923791705,0.3607615730481384,0.4679139345726904])
                        x=np.array([0.9324695142031522,0.6612093864662647,0.238619186083197])
                    else:
                        if abs(r) < 0.75:
                            w=np.array([0.04717533638651177,0.1069393259953183,0.1600783285433464,0.2031674267230659,0.2334925365383547,0.2491470458134029])
                            x=np.array([0.9815606342467191,0.904117256370475,0.769902674194305,0.5873179542866171,0.3678314989981802,0.1252334085114692])
                        else:
                            w=np.array([0.01761400713915212,0.04060142980038694,0.06267204833410905,0.08327674157670475,0.1019301198172404,0.1181945319615184,0.1316886384491766,0.1420961093183821,0.1491729864726037,0.1527533871307259])
                            x=np.array([0.9931285991850949, 0.9639719272779138, 0.9122344282513259,0.8391169718222188, 0.7463319064601508, 0.6360536807265150,0.5108670019508271,0.3737060887154196,0.2277858511416451,0.07652652113349733])
                    w = np.concatenate((w, w), axis=0)
                    x = np.concatenate((1 - x, 1 + x), axis=0)

                    if abs(r) < 0.925:
                        hs=(h*h+k*k) / 2
                        asr=np.arcsin(r) / 2
                        sn=np.sin(asr*x)
                        bvn=np.dot(np.exp((sn*hk-hs)/(1-sn**2)),w.T)
                        bvn=bvn*asr/tp + norm.cdf(-h)*norm.cdf(-k)

                    else:
                        if r < 0:
                            k=- k
                            hk=- hk
                            
                        if abs(r) < 1:
                            as1=1 - r ** 2
                            a=np.sqrt(as1)
                            bs=(h - k) ** 2
                            asr=- (bs / as1 + hk) / 2
                            c=(4 - hk) / 8
                            d=(12 - hk) / 80

                            if asr > - 100:
                                bvn= a*np.exp(asr)*(1-c*(bs-as1)*(1-d*bs)/3+c*d*as1**2)
                            if hk > - 100:
                                b=np.sqrt(bs)
                                sp=np.sqrt(tp)*norm.cdf(-b/a)
                                bvn=bvn - np.exp(-hk/2)*sp*b*( 1 - c*bs*(1-d*bs)/3 )
                            a=a / 2
                            xs=(a*x) ** 2
                            asr=- (bs / xs + hk) / 2
                            ix=asr > - 100
                            xs=xs[ix]
                            sp=( 1 + c*xs*(1+5*d*xs) )
                            rs=np.sqrt(1 - xs)
                            ep=np.exp(- (hk / 2)*xs / (1 + rs) ** 2) / rs
                            bvn=( a*( np.dot((np.exp(asr[ix])*(sp-ep)),w[ix].T) ) - bvn )/tp

                        if r > 0:
                            bvn=bvn + norm.cdf(- max(h,k))

                        else:
                            if h >= k:
                                bvn=- bvn

                            else:
                                if h < 0:
                                    L=norm.cdf(k) - norm.cdf(h)
                                else:
                                    L=norm.cdf(- h) - norm.cdf(- k)
                                bvn=L - bvn

                    return max(0,min(1,bvn))