Matlab <> python

When I coded the Fourier pricing technique for this blog, it was the first time I used such technique.
Speed wasn't that impressive but without any reference, it's difficult to judge.

Recalling from my digital signal processing course, we know that FFT stands for fast Fourier transform and it relies on symmetries to achieve its full speed for signals of length $2^N$.

I looked further and two things happened:
1/I mis translated matlab code
"x=x_min:dx:x_max"
into
"x=np.array([(dx*i+x_min) for i in range(N+1)]) "
and made a similar translation for the vector of frequency w. Price was converging it seemed but the length of the produced vector is odd.
2/numpy FFT behaved too nicely, instead of complaining about the vector size not being in the form $2^N$ it ran other optimized Fourier transforms. Such transforms however aren't nearly as fast as FFT.

To give an idea, with only 1024 steps, the previous version was 4 cents off, and PC parity off by less than 2 cents. When we keep in mind the slow convergence of the MC, we could accept such values. But the correct version is less than 0.05 cent from the theoretical price. Call put parity is perfect and my time measure says 0. Using 2^15 steps only takes 0.05 second and give an accuracy of 0.00005 cents. Find below the corrected code for the european option pricer and the one for the american pricer.

1. def FourierST3(option,params):  
2.    N=params.NAS  
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=linspace(x_min,x_max,N)  
11.      
12.    #create vector in the fourier space  
13.    w_max=np.pi/dx;  
14.    dw=2.0*w_max/(N);  
15.    w=np.concatenate((linspace(0,w_max,N/2+1),linspace(-w_max+dw,-dw,N/2-1)))  
16.   
17.    # Option payoff  
18.    s = S0*np.exp(x);  
19.    v_call = np.maximum(s-K,0)  
20.    v_put = np.maximum(K-s,0)  
21.      
22.    # FST method  
23.    char_exp_factor = np.exp((j*(r-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-r)*T)  
24.    VC = np.real(np.fft.ifft(np.fft.fft(v_call)*char_exp_factor))  
25.    VP = np.real(np.fft.ifft(np.fft.fft(v_put)*char_exp_factor))  
26.   
27.    #Interpolate option prices  
28.    tck=si.splrep(s,VC)  
29.    option.call.price=si.splev(S0,tck,der=0)     
30.    tck=si.splrep(s,VP)  
31.    option.put.price=si.splev(S0,tck,der=0)    
32.     

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

I chose to keep the older code for reference but clearly the pricers FourierST3 and FourierST4 should ALWAYS be preferred to FourierST and FourierST2.