Implementing a Fast Fourier Transform for Option P

Question

I'm in need of some tips regarding a small project I'm doing. My goal is an implementation of a Fast Fourier Transform algorithm (FFT) which can be applied to the pricing of options.

First concern: which FFT?

There are a lot of different FFT algorithms, the most famous one being Cooley-Tukey. My thoughts on this: I prefer the most simple one, since this is no thesis or big project, just a course on Algorithms. But it has to be compatible with option pricing (in contrast with the most - well in our general literature- referenced application of images/sound processing). So it depends on the form of input that is provided (on which I need some advice). I'm familiar with the several improvements, like a Fractional FFT, mixed-radix FFT etc. But these seem pretty complex and optimization/performance driven, which is not relevant for my project.

Second concern: which pricing model?

I Guess Black-Scholes (BS) is a bit too 'flat' and I am aware of the several models that emerged after BS. So, with the same objectives as stated above, I'd initially prefer the Heston model.

There are a lot of considerations, and the truth is that I just can't see the wood for the trees.

Some background info:

My background is a B.Sc in Mathematics (Theoretical), so I have some understanding of Fourier transforms.

The goal is a working FFT implementation for calculating option pricing. It does not have to be the fastest (no extreme optimization). The goals are trying to understand the chosen FFT and having a real-world working application.

So could you give some advice on the choices?

I've read a lot of papers on FFT + Option pricing, say all the decent hits on googles first few pages. But those studies were written with a much 'higher' cause.

Answer 1

I've been researching this topic - FFTs applied to options pricing - for a few weeks now. It turns out there's an extensive body of work on the subject, so it's hardly futile as Alexandre implies.

The most readable basic paper I've found is from Carr and Madan - www.math.nyu.edu/research/carrp/papers/pdf/jcfpub.pdf - but there are many others in varying levels of detail, which Google will find via "option pricing Fourier" search.

I may be coding this up in R in the near future; I'm trying to locate a decent source of options price data for testing.

Answer 2

An FFT is just an optimised implementation of the DFT. I suggest either using an existing FFT library, such as KissFFT, or if you really want to implement this from scratch then just implement the DFT rather than an FFT, since it is much simpler, and performance should not be an issue unless you have high data rates or large amounts of data.

Answer 3

If your goal is to make some use of the FFT, then your choices are poor: only affine models give you enough information to obtain the Fourier transform of the spot density. In practice, this means Black-Scholes, or Heston. Perhaps a few more, but none of the "useful" models.

Heston's model has peculiar features (pertaining to its implied vol dynamics) which makes it quite useless as a stochastic vol model. I suppose it is popular precisely because of the fact that you can price vanilla options in semi-closed form, through Fourier transforms. With modern technology, this is no longer a real asset.

If you are interested in option pricing, I'd therefore suggest you don't try too hard with FFT, and turn to PDE or Monte-Carlo methods: the range of models you can play with are much more interesting (and much more valuable on the job market, in case you're interested).

For the FFT part of your question, implementing Cooley-Tukey from scratch is not hard, and you can start there. Of course, in production code, you are better using a canned package (like FFTW).

Answer 4

I'm providing below the implementation of the radix-2 Decimation In Time Cooley-Tukey scheme in Matlab. The code is an iterative one and considers the scheme in the following figure:

A recursive approach is also possible.

As you will see, the implementation calculates also the number of performed multiplications and additions and compares it with the theoretical calculations reported in How many FLOPS for FFT?.

The code is obviously much slower than the highly optimized FFTW exploited by Matlab.

Note also that the twiddle factors omegaa^(interButterflyIndex * 2^(numStages - p)) can be calculated off-line and then restored from a lookup table, but this point is skipped in the code below.

% --- Radix-2 Decimation In Time - Iterative approach

clear all
close all
clc

N = 32;

x = randn(1, N);
xoriginal = x;
x = bitrevorder(x);
xhat = zeros(1, N);

numStages = log2(N);

omegaa = exp(-1i * 2 * pi / N);

mulCount = 0;
sumCount = 0;
tic
for p = 1 : numStages
    alpha = 2^(p - 1);
    butterflyStart = 1;
    while (butterflyStart <= (N - alpha))
        for interButterflyIndex = 0 : alpha - 1
            xhat(butterflyStart)          = x(butterflyStart) + x(butterflyStart + alpha) * omegaa^(interButterflyIndex * 2^(numStages - p)); 
            xhat(butterflyStart + alpha)  = x(butterflyStart) - x(butterflyStart + alpha) * omegaa^(interButterflyIndex * 2^(numStages - p));
            mulCount = mulCount + 4;
            sumCount = sumCount + 6;
            butterflyStart = butterflyStart + 1;
            if (interButterflyIndex == (alpha - 1))
                butterflyStart=butterflyStart + alpha;
            end;
        end;
    end;
    x = xhat;
end;
timeCooleyTukey = toc;

tic
xhatcheck = fft(xoriginal, N);
timeFFTW = toc;

rms = 100 * sqrt(sum(sum(abs(xhat - xhatcheck).^2)) / sum(sum(abs(xhat).^2)));

fprintf('Time Cooley-Tukey = %f; \t Time FFTW = %f\n\n', timeCooleyTukey, timeFFTW);
fprintf('Theoretical multiplications count \t = %i; \t Actual multiplications count \t = %i\n', ...
         2 * N * log2(N), mulCount);
fprintf('Theoretical additions count \t\t = %i; \t Actual additions count \t\t = %i\n\n', ...
         3 * N * log2(N), sumCount);
fprintf('Root mean square with FFTW implementation = %.10e\n', rms);