Python Numpy FFT array size limit?

Question

I am trying to find the price of an Option based on the fft technique within the binomial model and it works fine until N>40000 where I start getting negative values and weird convergene and I am not sure whether it's a coding problem, a limitation of the arrays in numpy or computational error from N being too large. Here's my function to determine the price of the Call Option at time 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

Im going to hazard a guess that your problem is u**(N-i). Large exponents are notoriously poor performers, I would first look to restructure that aspect of the code and then isolate other poorly performant sections afterwards.

For example you might observe that:

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

is equivalent to:

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

then u**N can be extracted out of the loop as a constant and you are left with an iterator:

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

Broadly the machine tolerance of 64bit floats is around 1e-15. Suppose that (d/u)=0.999 then the number of multiplications (your exponent) that can be performed before precision is lost in this case is:

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