Option pricing using discrete fourier transform (python)

Question

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])

Answer 1

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.