Pricing Options with Fourier Series P2

The how to go from knowing what a Fourier Series is to being able to price Options under Black-.

• Please check out the Part 1 if you want to see more about the fundamentals of the Fourier Series.
• Accompanied Notebook: Github
• This article is based on notes from my course in Stochastic Calculus (2023) at the University of Amsterdam under Professor Michel Vellekoop.
• Check out part 3: Link

Contents:

• Fourier cosine expansion
• The density function of S(T)
• Euler's formula for cosine
• The characteristic function — The characteristic function of the Geometric Brownian motion.
• Option Pricing formulas
• Convergence over N terms and the importance of picking good limits
• Error analysis over the analytical solution

“There is nothing impossible to him who will try”
Alexander the Great

Fourier cosine expansion

Similar to the Sine-Cosine Expansion, we can do the same thing with only the cos part of the equation:

Going from the Cosine-Sine transformation to the Cosine Transformation

You can see from the following function which we also did in Part 1, that the cosine is able to approximate very well the same functions

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

mean = 100
std  = .1 *sqrt(5)*100 
f    = lambda x : (1/(std*sqrt(2*pi))) * exp(-(x-mean)**2/(2*std**2)) #*12/0.0175

a = mean - 12 * std
b = mean + 12 * std

Square Function — Source:Notebook

Line Function — Source: Notebook

Normal Distribution — Source: Notebook

Avg Residuals
N   Square f   Line f      Normal Distribution
8   1.311711   0.447389    1.593045e-04
16  0.784683   0.264635    1.214211e-07
32  0.440387   0.153540    7.100142e-18
64  0.268449   0.088745    3.628299e-14
128 0.154604   0.052147    2.992953e-12

Remarks

• The conclusions are the same as Part 1, that the normal distribution fits very well and we need to take an adequate range in the values [a, b]
• Since only the cosine part is used. It simplifies a lot of the next steps we have to do and speeds up the process as we have to calculate fewer things.

Density Function of S(T)

To start gaining intuition on the problem we have to solve it’s good to start by making a function and its Fourier transform, which given some initial values, S0, T, r, sigma it calculates the Probability of S_T = X.

Method to calculate the density function of S(T)

It's important to note that an interval of [-12 sigma, 12 sigma] covers 99. plus another 33 9s! Someone could say it's extreme but it's important to use a high sigma so the density function covers deep out-of-the-money options. (I will come back later to that)

Density function of S(T) with S0=100, T=5Y, r=5%, sigma=10% — Source: Notebook

N    Avg Residual
8    1.897794e-02
16   7.159815e-03
32   5.616416e-04
64   2.326956e-06
128  1.384069e-10

Pricing Options with Fourier Series

Step 1: Euler’s Formula to transform cos(x)

We gonna start our journey of deriving a function to price an option with Fourier Series by using Euler’s Formula to transform the cos(x). This will help with the next steps we have to do.

Step 2: Characteristic Function of a random variable

In probability theory, a characteristic function is a complex-valued function defined for a random variable. It captures the moments and distribution properties of the random variable, allowing for the analysis of various statistical properties using techniques such as Fourier transforms. For example for a random variable following a normal distribution or a Poisson distribution:

Characteristic Function of a Random Variable

Step 3: Analytical solution to the integral using the characteristic function

In our use case, we can transform the density function of i.e. log(S_T). Continuing the transformation of f(x) from Step 1 of An in combination with Step 2 and knowing that f(x) is the density of a random variable with a known characteristic function:

Knowing that we can transform the integral to the characteristic function allows us to explicitly calculate the integral which I did numerically with the quad from scipy before. The expectation is assumed for a large enough interval of [a, b] since f(z) is the density function of a random variable.

Step 4: Tying it all together to calculate explicitly the An coefficients

Where En is part of the integral outside of [a,b] which is very small if for probability mass f(x) we pick a,b where f(x) is close to zero (±12sigma, for the normal). A similar concept can be applied if f(x) is the payoff of an option and the payoff is very small for a,b

def Fourier_cosine_char_function(char_function, x:np.array, a, b, N:int):
    A  = lambda n : 2/(b-a) * np.real( char_function((n*pi)/(b-a))*exp((-1j*n*pi*a)/(b-a))  )
    fx = .5*A(0)
    for n in range(1, N+1):
        fx += A(n) * cos(n*pi * (x-a)/(b-a))
    
    return fx

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

# For normal distribution 
mu      = 100
sigma   = 25
char_function = lambda u: exp(1j*u*mu -.5*sigma**2*u**2)
f             = lambda x: norm.pdf(x, loc=mu, scale=sigma)

a   = 100 - 12 * sigma
b   = 100 + 12 * sigma

The density of a normal distribution using an explicit solution with the characteristic function - Source: Notebook

N   Avg. Residual    Execution Time
8   1.182498e-03     0.000059
16  1.588727e-04     0.000092
32  1.160981e-07     0.000180
64  2.983220e-18     0.000358
128 2.986019e-18     0.000767

Remarks

• As you can see we are doing very well. Execution time is low, and we get a very quick convergence with only 32 terms
• Using the characteristic function increases the accuracy of the Fourier Series approximation and makes the calculations simpler.
• Using the imaginary numbers transformation for cos allows us to transform the integral in the coefficients to the characteristic function of a density function.

Option pricing using the Fourier series

Step 1: Recap

For a function f where it represents the density function of a random variable z, we have the Fourier cosine expansion:

Density Function of random variable z with known characteristic function

Step 2: Characteristic function of the density of future asset prices

When we price assets we have a density function of the log of the value of an asset which accepts parameters: t, T, S(t), r, sigma, and outputs the probability of S(T). It is straightforward in stochastic calculus to calculate the characteristic function of a random variable, based on Step 2 of the previous section. We are going to use this new tool to price a derivative F with a defined payoff F(T, S_T):

Formula to calculate the price of a derivative with a known payoff at T

The steps that I did above are:

• Discount the value of the expected payoff of the derivative at T. (line 1)
• Substitute to the density function the Fourier Series Transformation we derived above. (line 2)
• Separate the terms which contain z and package them to the coefficient gn

Step 3: Solving g(n) analytically

Since the integral in gn is for values at T we can analytically calculate the integral if the payoff function is explicit. In this case, we will try to price a European Put Option with payoff max{K-S(T), 0}. Therefore:

Analytical solution for Gn coefficients

To obtain the solutions shown here you can use an online definite integral calculator, and for g0 you are going to use the limit of gn as n approaches 0 to get a solution as substituting n=0 means g0 cannot be defined.

Step 4: GBM Characteristic function

Now that we have our “weapon” we only miss the bullets. In this case, the random variable z is assumed to follow a geometric Brownian motion under Q. With a characteristic function the “bullets”:

The characteristic function of Geometric Brownian motion

def Fourier_BS_Put(S0, K, T, sigma, r, N):
    _sigma, _T = max(.1, sigma), max(1, T)
    a = log(S0) + r * _T - 12 * _sigma * np.sqrt(_T) # Upper limit -12 std from the expected value
    b = log(S0) + r * _T + 12 * _sigma * np.sqrt(_T) # Lower limit -12 std from the expected value

    h       = lambda n : (n*pi) / (b-a) 
    g       = lambda n : (exp(a) - (K/h(n))*sin(h(n)*(a - log(K))) - K*cos(h(n)*(a - log(K)))) / (1 + h(n)**2)
    g0      = K*(log(K) - a - 1) + exp(a)
    
    m       = log(S0) + (r - .5*sigma**2)*(T)
    gbm_char= lambda u : exp(1j * u * m -0.5*(u**2 * sigma**2)*T)

    F = g0 
    for n in range(1, N+1):
        h_n = h(n)
        F += 2*bs_char(h_n) * exp(-1j*a*h_n) * g(n)

    F = exp((-r*T))/(b-a) * np.real(F)
    return F

Price of an option over different strike prices. S0=100, T=5Y, r=5%, sigma=10% — Source: Notebook

N  avg_residuals Ex. Time
8   5.075343e-01 0.000185
16  2.556765e-02 0.000179
32  5.810656e-06 0.000265
64  2.319065e-14 0.000461
128 2.319066e-14 0.000787

Remarks

• This method is proving to be a very fast and efficient way to calculate the value of a European Put Option
• The only things that we control to influence the error in the approximation are N, the number of terms, and the interval [a, b] of the values of S(T) which we look at.

Picking the correct limits [a, b]

The graph below illustrates 2 very important things someone needs to understand when using the Fourier series.
1. That the limits [a=S_T(-z_score), b=S_T(z_score)] need to cover all possible strike prices
2. As sigma increases so does the error in the approximation

Illustration of selecting correctly the limits to the distribution and how the error in the approximation increases with larger sigmas — Source: Notebook

In the first subplot, you can see that if we had a strike price at S(T) = 600 we would not be able to price it correctly since we simulated for values only S(-12) to S(12). Therefore you can see in the second plot that the value of such a strike price would be very wrong!
In the third plot, you can see the sigmas covering a wide range of strike prices. However, due to the shape of the payoff it's harder to approximate the function accurately with the same number of terms. Therefore as sigma increases you may need to increase the number of terms N.

The average error of the first subplot is: 7.676e-12
The average error of the first subplot is: 2.830e-02
The average error of the first subplot is: 1.405e-03
*for N = 128

Error in the approximation over a different number of terms N

Error in approximating the Analytical solution of a European Put option — Source: Notebook

I conducted a simple experiment to see how the error is reduced over a different number of terms. You can see that we have an exponential convergence while the computational cost is linear. The dream of every computational scientist. This means that when we replace the characteristic function with the Heston model’s we don't have to pay the penalty of the Monte Carlo numerical integration, in which the error decreases with the square root of N.

Error analysis over all possible input values for a European Put

For a Stock of value 100 EUR. I simulated the conditions of pricing a Put where:
K = [0, 500]
T = [0, 30Y]
r = [-3%, 20%]
sigma = [0.5%, 100%]

Visualization of average error over different values of the parameters — Source: Notebook

Remarks

• It appears that for options at expiration T=0, the approximation has increased error but this is attributed to the known payoff which is a line function that I showed at the beginning of the article, which is more challenging to approximate than a curvy function. To combat this, a higher number of terms can be used.
• Sigma has a linear relationship with the error in the approximation as it changes the shape of the function to something that is harder to approximate. (As I showed earlier) Again, to combat this, a higher number of terms can be used.

The Total Average Error: 0.00289 EUR over all the given ranges of values K, T, r, and sigma I defined earlier.

Conclusion

Step by step I explained how to go from the general Fourier transformation to custom-made functions for calculating the payoff of a European Put option.
However, if you want to price derivatives where the payoff does not depend on fixed intervals, you have to use more complex numerical techniques. Like Finite Difference, or Monte-Carlo.

Despite this being a ‘useless’ tool so far we laid the groundwork for pricing options using the Heston model which we can't do analytically. Hopeful you can use the same principles for pricing options with

I hope you found this article useful I hope to see you in Part 3/3 which we are going to see how to use the Heston model with the Fourier Series and how to calibrate it over on the Option chain.

Link to part 3: Link