¿Límite de tamaño de la matriz FFT de Python Numpy?

Question

Estoy tratando de encontrar el precio de una Opción basado en la técnica fft dentro del modelo binomial y funciona bien hasta N>40000 donde empiezo a obtener valores negativos y convergentes extraños y no estoy seguro de si es un problema de codificación, una limitación de los arrays en numpy o un error computacional por ser N demasiado grande. Aquí está mi función para determinar el precio de la opción de compra en el momento 0.

def FFTBinCall(S0,R,sigma,K,T,N):
    pQ = 0.5
    qQ = 1- pQ
    Dt=T/N
    #Initialize Vectors of Proper Dimensions
    C_T=np.zeros(N+1)
    S_T=np.zeros(N+1)

    u=1+R*Dt+sigma*math.sqrt(Dt)
    d=1+R*Dt-sigma*math.sqrt(Dt)
    D=math.exp(-R*Dt)

    for i in (range(0,N+1)):
        S_T[i]=S0*u**(N-i)*d**(i)
        C_T[i]=max(S_T[i]-K,0)

    QDistr=np.concatenate([[pQ,qQ], np.zeros(N-1)])
    Discounted_QDistr=QDistr*D

    C_0=np.fft.fft(np.fft.ifft(C_T)*np.fft.fft(Discounted_QDistr)**N).real
    return C_0[0]

Answer 1

Voy a arriesgarme a decir que su problema es u**(N-i) . Los exponentes grandes son notoriamente de bajo rendimiento, yo buscaría primero reestructurar ese aspecto del código y luego aislar otras secciones de bajo rendimiento.

Por ejemplo, puede observar que:

S_T[i] = S0 * u**(N-i) * d**(i)

es equivalente a:

S_T[i] = S0 * u**N * (d/u)**i

entonces u**N puede ser extraído del bucle como una constante y te quedas con un iterador:

S[0] = 1
for i in range(1, N+1):
    S_T[i] = S_T[i-1] * (d/u)
S *= S0 * u**N

En general, la tolerancia de la máquina a los flotadores de 64 bits es de alrededor de 1e-15. Supongamos que (d/u)=0.999 entonces el número de multiplicaciones (su exponente) que se pueden realizar antes de perder la precisión en este caso es:

(d/u)**x = 1e-15
x = log(1e-15) / log(d/u)
x = 34521