Pricing Options with Fourier Series P1

6 min readJun 19, 2023

Theory and intuition on how to use Fourier Series to price derivatives.
Part 1 out of 3.

Notebook: Github
Please check out the Part 2 and Part 3 if you want to see more.

This is the first article where I am going to help you gain an intuition behind how to use this new powerful tool to solve semi-analytical problems in finance and replace your Monte-Carlo Methods.

Introduction

We all know and love the monte carlo methods to numerically integrate but what if I told you that you can use imaginary numbers and Fourier Series to replace Monte Carlo?
The main benefit is speed which in option pricing is very important. This is very important as the Heston model for pricing stock options needs to numerically integrate which with Monte Carlo would take around 100ms and with the Fourier Series a dozen of μs as you can see in Part 3 of my article series.

Section 1: But what is a Fourier Series?

For any function f and for an interval a,b we can approximate f(x) as the infinite sum of cosines and sines, with L = b-a.

Section 2: Apply the math formulas to Python

def get_fourier_approx(f, x:np.array, a:float, b:float, N:int):
    fa = lambda x, n : f(x) * cos((2*pi*n*x)/(b - a))
    fb = lambda x, n : f(x) * sin((2*pi*n*x)/(b - a))

    A0 = 1/(b - a) * quad(f, a, b, limit=200)[0]

    Cosine_Sine_Sum = np.zeros_like(x)
    for n in range(1, N+1):
        A = 2/(b - a) * quad(fa, a, b, args=(n), limit=200)[0]
        B = 2/(b - a) * quad(fb, a, b, args=(n), limit=200)[0]
        Cosine_Sine_Sum += A*cos((2*pi*n*x)/(b - a)) + B*sin((2*pi*n*x)/(b - a))

    fx = A0 + Cosine_Sine_Sum
    return fx

a = -6
b = 6
x = np.linspace(a, b, 1_000)
y = f(x)

fig, (ax1, ax2) = plt.subplots(2, figsize=(20,12))
blue_shades = ['#0000FF', '#3399FF', '#66B2FF', '#99CCFF', '#CCE5FF']

avg_residuals = []
Ns = [8, 16, 32, 64, 128]
for i, N in enumerate(Ns):
    fx = get_fourier_approx(f=f, x=x, a=a, b=b, N=N)
    ax1.plot(x,fx, blue_shades[i], label=f'N = {N}')
    ax2.plot(x,y-fx, blue_shades[i], label=f'N = {N}')
    avg_residuals.append(np.abs(y-fx).mean())

ax1.set_title('Fourier Transform of f(x)')
ax1.plot(x,y,'tab:red', linestyle='--')
ax2.set_title('Residuals')
plt.tight_layout() ; ax1.legend();ax2.legend() ; plt.show()

pd.Series(avg_residuals, index=Ns, name='Avg Residual')

Square Function:

N      Avg. Residual
--------------------
8      1.311711
16     0.784683
32     0.440387
64     0.268449
128    0.154604

Line Function:

N      Avg. Residual
--------------------
8      0.447389
16     0.264635
32     0.153540
64     0.088745
128    0.052147

Normal Distribution

Scaled so y in [0, 12] with:
- mean = 100
- std = 0.1 *sqrt(5)*100
- a = 100 -12 * std
- b = 100 +12 * std

N      Avg. Residual
--------------------
8      1.092374e-01
16     8.326020e-05
32     6.878539e-14
64     5.721031e-14
128    5.170898e-14

Remarks

All the distributions are scaled such that y ranges from [0,12], so we can compare the magnitude of the residuals.
You can see from the plots and the residuals, that the more curvy a function is, the faster the Fourier series converges to the correct value. We leverage this property as the normal and the log-normal doesn't need a lot of terms to calculate with sufficient accuracy in our approximation.
There is a significantly higher error at the beginning and the end of the data. Therefore it's a good practice to include higher limits than what is intended to be used. For example, calculate ±4std when you need ±3. This makes deep out-of-the-money options trickier to calculate.

Section 3: Log-normal Distribution of S(T)

S_T follows a simple GBM under Q, where we can derive the probability density of S_T using the following equations:

Now we can define f(S_T) in Python with the following functions and define the lower limit as ( 0, S_0*exp(r*T) + 12 * sigma*sqrt(T)*S_0 )

S0      = 100
r       = 0.05
sigma   = 0.1
T       = 5.0

Z = lambda St : np.where(St > 0, ((log(St/S0) - (r - .5*sigma)*T)/(sqrt(T)*sigma)), -np.inf)
f = lambda x : norm.pdf(Z(x))

a   = S0*exp(r*T) - 12 * sigma*sqrt(T)*S0
b   = S0*exp(r*T) + 12 * sigma*sqrt(T)*S0

N       Avg. (scaled) Residual      Avg. Residual       Execution Time (sec)
----------------------------------------------------------------------------
8       0.176429                    5.880975e-03        0.112720
16      0.004235                    1.411566e-04        0.246473
32      0.000030                    9.855127e-07        0.624209
64      0.000027                    8.918504e-07        1.936948
128     0.000026                    8.530034e-07        6.741019

Remarks:

I have included both scaled and non-scaled residuals. Scaled Residuals correspond to the (incorrectly) scaled probability such that max{y}=12, wherein the (correct) non-scaled, max{y}=0.4. This was done to compare the fit of the log-normal distribution to the rest of the functions plotted above.
We can deduce that the log-normal distribution is harder to fit that the normal distribution due to the nonsymmetrical shape.
We can see that with 64 terms in the non-scaled version, the expected error in calculating P(S_T=x) is very small, less than 0.0001%.
It's very important to center the distribution around the expected value at T, with 12stds symmetrically left and right. I made a version where a, and b are not symmetrical, and the residuals are not evenly spread out. Intuitively, you would take a=0+ but it will not generate desirable outcomes. *The residuals of S_T are calculated for values of S_T>0 since it's not possible and we don't care about values less than 0.

The density of S_T with a=1e-8, showing undesirable properties at the residuals. Source: Notebook

Limitations — Drawbacks — Improvements:

There is no benefit in approximating any of the previously mentioned functions as a Fourier series and using the numerical integration as a means to calculate An and Bn.
It is very important to calculate the An and Bn coefficients analytically so the only numerical part is calculating the series.
The benefit lies elsewhere. It is when f(x) does not have an explicit form and numerical integration is needed, and we can solve the problem semi-analytically with the characteristic function and Fourier series.
If we were to price an Option with standard BS in Python with scipy.norm it takes around 0.06 milliseconds. So it doesn't make sense to use Fourier for the standard BS.
However, if we analytically solve the integrals A0, An, Bn, and use the version with complex numbers we get around 0.6miliseconds, which is comparable. We will use this to our benefit in the Heston model which is the industry standard for pricing such options, in part 3.

Next Article: Part 2

Medium: Link

In the next article, we see how to use this knowledge in combination with the characteristic function, (the use of imaginary numbers) to calculate the value of a European Call option using the standard Black-Scholes model.