Fourier

To sum up what we did so far for pricing option (with constant volatility)

-We can solve the expectation of the option value at expiration. It can be done explicitly under the risk neutral measure, that's the Black-Scholes formula. It works for european options.

-We can simulate replicating the option via a dynamic stock trading strategy. It has been shown that the costs of this replication converges toward the BS price for european options under some assumptions. This is binomial pricing, we can use it for european and american options.

-We can discretize the differential equation for the stock and the option and solve it step by step via PDE approach. This can be done for european and american options.

-We can approximate the expectation using simulation (MC) methods. We know that if we use enough random (quasi or pseudo) draws the measured average will converge toward the expectation. This can also be used for american options.

-We will turn to another technique: the expectation of the option is given by the convolution between the log normal distribution and the payoff of the option.

Basically, what we want to know is how the option will be worth at expiration, given that we know that the distribution of the stock at expiration. More formally (taking $r=0$):
\[
E_t^Q \left[\Phi(S(T))\right]=\int^{\infty}_{-\infty}\Phi(S(t)e^y)f_x(y)dy=\Phi(S(t)e)*f_x
\]
where $\Phi(S)$ is the payoff of the function, $f_x$ is the stock price density function and $E_t^Q$ denotes the expectation under the risk neutral measure (also called Q measure) and $*$ denotes the convolution. Calculating this can be done using Fourier transform. If I note $F(a)$ the fourier transform of $a$ and $a.b$ the usual product of $a$ and $b$.
\[
F(a*b)=F(a).F(b)
\]
So if somehow we can get
\[
F^{-1}(  F(\Phi(S(t)e).F(f_x))
\]
in an fast way, that would be a valid pricing method.
The characteristic function is actually the Fourier transform of the density function. So we just need to compute the Fourier transform of the payoff multiplied by exp, mutiply this by the characteristic function, get an inverse Fourier transform of the result et voila.
The characteristic function for lognormal density is given by:
\[
e^{i(r-\sigma^2/2)\omega -\frac{\sigma^2\omega^2}{2}}
\]

This is excatly what this program does. It is a rewriting in python of some matlab code of Dr Surkov (find his PhD thesis there).

1. def FourierST(option,params):  
2.    N=params.N  
3.    [S0,K,sigma,T,r,divi,american]=[option.S0,option.K,option.sigma,option.T,option.r,option.divi,option.american]  
4.   
5.    j=complex(0,1)  
6.    #create vector in the real space  
7.    x_min = -7.5  
8.    x_max = 7.5  
9.    dx=(x_max-x_min)/(N-1)  
10.    x=np.array([(dx*i+x_min) for i in range(N+1)])  
11.      
12.      
13.    #create vector in the fourier space  
14.    w_max=np.pi/dx;  
15.    dw=2*w_max/(N);  
16.    w1=np.array([(dw*i) for i in range(N/2+1)])  
17.    w2=np.array([(dw*i-w_max) for i in range(N/2)])  
18.    w=np.concatenate((w1,w2))  
19.      
20.   
21.    # Option payoff  
22.    s = S0*np.exp(x);  
23.    v_call = np.maximum(s-K,0)  
24.    v_put = np.maximum(K-s,0)  
25.      
26.    # FST method  
27.    char_exp_factor = np.exp((j*(r-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-r)*T)  
28.    VC = np.real(np.fft.ifft(np.fft.fft(v_call)*char_exp_factor))  
29.    VP = np.real(np.fft.ifft(np.fft.fft(v_put)*char_exp_factor))  
30.   
31.    #Interpolate option prices  
32.    tck=si.splrep(s,VC)  
33.    option.call.price=si.splev(S0,tck,der=0)     
34.    tck=si.splrep(s,VP)  
35.    option.put.price=si.splev(S0,tck,der=0)   

To get the complex number, I imported scypi.
To use that code, one could use the following piece:

1. #test  
2. S0 = 100.0  
3. K = 100.0  
4. r=0.0  
5. sigma = .3  
6. T = 3  
7. Otype='C'  
8. div=[[0.9,1.9,2.9],[0,0,0]]  
9.   
10. myOption=OptionClass(S0,K,r,sigma,T,'false')  
11. myOption.divi=div  
12.           
13. myStrike=strikeClass(S0,K,r,sigma,T,'false')  
14. div=[[0.090,0.095,0.099],[0,0,0]]  
15. myStrike.divi=[[],[]]  
16.   
17.       
18. myParams=ParamsTree(2**15)  
19.   
20. print "S0\tstock price at time 0:", S0  
21. print "K\tstrike price:", K  
22. print "r\tcontinuously compounded risk-free rate:", r  
23. print "sigma\tvolatility of the stock price per year:", sigma  
24. print "T\ttime to maturity in trading years:", T  
25.   
26.   
27. t=time.time()  
28. c_BS = BlackScholes(Otype,S0, K, r, sigma, T)  
29. elapsed=time.time()-t  
30. print "c_BS\tBlack-Scholes price:", c_BS, elapsed  
31.   
32. t=time.time()  
33. S3 = FourierST(myStrike,myParams)  
34. elapsed=time.time()-t  
35. print "c_BT\tFourier:",myStrike.call.price,elapsed  

Note that I use a ParamsTree structure as I just need a number of steps (used to discretize the space of of possible stock values and via nyquist sampling also used for the discretization of the fourier space $\omega$) but in a next version I'll make a proper structure for that pricer.
The number of step must be dividable by 2.