Reducing Numerical Error in FFT Binomial Option Pricing
Summary
A question reports unstable or negative option values when the step count in an FFT-based binomial calculation becomes large. The accepted response suspects numerical precision loss in the terminal stock-price calculation, which raises up and down factors to large powers. It proposes rewriting the stock price as a constant factor times successive powers of the ratio of the down and up factors. That ratio can be applied iteratively, while the shared large-power factor is computed outside the loop.
The answer explains that repeated floating-point operations can lose precision and illustrates how repeated multiplication of a value close to one can eventually reach the scale of machine precision. This is a plausible diagnostic and a way to restructure the computation, not a demonstrated diagnosis of the specific code. The thread does not compare alternative pricing algorithms or establish that this change alone resolves the observed convergence behavior; other numerical or modeling issues may also matter.
Key ideas
- Large powers in binomial terminal-price calculations can create numerical precision problems.
- Factor the common up-factor power outside the loop and iterate using the down-to-up ratio.
- Repeated floating-point multiplication can accumulate error when the ratio is near one.
- The proposed cause is a hypothesis and is not validated with a corrected implementation or benchmark.
Tags
Full text
# Python Numpy FFT array size limit?
# 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]
```
## Answer by Attack68 (score 3, accepted)
https://quant.stackexchange.com/a/45672
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
```Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.