# Option pricing using discrete fourier transform (python)

I am trying to implement the pricing formula for a European (call) option given in Ales Cerny's paper "Introduction to Fast Fourier Transform in Finance" (paper can be found here), as follows:

My python code below does not return the correct answer, and in particular if I significantly increase the number of steps then I get a much larger answer. Where have I gone wrong?

import numpy as np
from numpy.fft import fft, ifft

def price_vanilla_option(s: float,
k: float,
r: float,
ro: float,
t: float) -> float:
"""
price vanilla option using Fast Fourier Transform
"""

steps = 1023  # 2^n - 1 for efficient fft
d_t = t / steps
discount = 1/(1 + r * d_t)

# use CRR probabilities
u = np.exp(ro * np.sqrt(d_t))
d = np.exp(-ro * np.sqrt(d_t))
p = (np.exp(r * d_t) - d)/(u - d)

# set up terminal vector and prob vector
c_n = np.zeros(steps + 1)
c_n[0] = s * (d ** steps)
for i in range(1, steps + 1):
c_n[i] = c_n[i - 1] * u / d
c_n = np.maximum(c_n - k, 0)
p_vec = np.pad([p, 1 - p], (0, steps - 1))

# fast fourier transform
c_0 = fft(ifft(c_n) * np.power(fft(p_vec) * discount, steps))
return np.real(c_0[0])

• Hi and welcome! I don't think that it is driving the problem here, but shouldn't you multiply by discount at the second-to-last line of your code? Feb 2 at 7:27
• ... and another question: are you sure that your functions in the second to last line are able to handle complex numbers? Feb 2 at 7:43
• Thanks v much @Kermittfrog - you're completely correct re the discount, I have updated accordingly. On your second point, I think they can do - I have tried both fft and ifft for simple examples with complex numbers and they appear to come to the expected result. Feb 2 at 18:52
• Ok, even the numpy.power function, yes? Feb 2 at 19:14
• Good point @Kermittfrog - I had not in fact checked that, but I have done so now and np.power() also returns the expected result for arrays of complex numbers Feb 2 at 20:08

## 1 Answer

In case this is useful for anyone else who comes across this, the issue was that I had set up my vector of terminal values the wrong way round (ie from smallest to largest rather than largest to smallest). The code to set up the terminal vector should read as follows:

# set up terminal vector and prob vector
c_n = np.zeros(steps + 1)
c_n[0] = s * (u ** steps)
for i in range(1, steps + 1):
c_n[i] = c_n[i - 1] * d / u
c_n = np.maximum(c_n - k, 0)



With this change, the code then produces the expected results.