# Python Numpy FFT array size limit?

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]


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

• Good observation at the end, +1 – Bob Jansen May 18 '19 at 14:25