American Fourier

The nice approach proposed by Surkov is built on the fact that instead of computing the convolution and get the expectation in one step, one can divide the time and perform a serie of convolutions. One obviously needs to adjust the characteristic function.

From a computational point of view this is slower but we can check early exercise at these steps. Globally speaking we go to the Fourier space, perform product, inverse the fourier transform. This gives us option prices for that time step, we calculate early exercise as for the pde or the tree.

Then we perform another step: fourier transform, convolution, inverse transform. We do so until we reach now (we go backward in time, one step at the time)

The code is straightforward:

1. def FourierST2(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=np.array([(dx*i+x_min) for i in range(N+1)])  
12.      
13.      
14.    #create vector in the fourier space  
15.    w_max=np.pi/dx;  
16.    dw=2*w_max/(N);  
17.    w1=np.array([(dw*i) for i in range(N/2+1)])  
18.    w2=np.array([(dw*i-w_max) for i in range(N/2)])  
19.    w=np.concatenate((w1,w2))  
20.      
21.   
22.    # Option payoff  
23.    s = S0*np.exp(x);  
24.    VC = np.maximum(s-K,0)  
25.    VP = np.maximum(K-s,0)  
26.      
27.    # FST method  
28.    char_exp_factor = np.exp((j*(r-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-r)*T/NTS)  
29.    for m in range(NTS):     
30.        VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
31.        VC=earlyCall(VC,s,K)  
32.          
33.        VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))  
34.        VP=earlyPut(VP,s,K)  
35.          
36.    #Interpolate option prices  
37.    tck=si.splrep(s,VC)  
38.    option.call.price=si.splev(S0,tck,der=0)     
39.    tck=si.splrep(s,VP)  
40.    option.put.price=si.splev(S0,tck,der=0)