Convergence FT

We now look at the convergence of the Fourier based pricer for american and european options. As I expected, one can see on the following graph the sawtooth pattern that existed with the binomial tree. It is present for both american and european options. First the american option:

and the the european one:

As you can guess, I fixed that with an adaptive step (curves with AS employ this technique). The very small oscillations we see for the american options with Fourier AS can be removed if we test early exercise only after the dividend is paid instead of at each step. This would be a slight speed up as well. However if we introduce stochastic volatility or jumps, it might not be the case that a call should only be exercised early before a dividend is paid so I prefer to remain general.

1. def FourierST6(option,params):  
2.    N=params.NAS  
3.    NTS=params.NTS  
4.    [S0,K,sigma,T,r2,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.    dt=float(T/NTS)  
23.   
24.    dividends,dividendsStep=preparedivs(divi,float(T),dt)  
25.    r,rArray,rateChangeStep,DFArray=preparerates(r2,T,dt,NTS)  
26.    
27.    rr=rArray[-1][0]    
28.   
29.    # FST method  
30.    char_exp_factor = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dt)  
31.    for k in range(NTS,0,-1):   
32.          
33.        if (k in rateChangeStep[:]):  
34.           rr= (r[1][np.nonzero(rateChangeStep==(k))[0]])  
35.           char_exp_factor = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dt)  
36.         
37.        #if there is a dividend, we interpolate the grid  
38.        if (k in dividendsStep):  
39.             if params.stepAdapt==True:  
40.                 #first get the fraction of dt  
41.                 fraction=dividends[0][np.nonzero(dividendsStep==(k))[0]]-k*dt  
42.                 dtf=dt-fraction  
43.                 #perform fraction of step  
44.                 char_exp_factorf = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dtf)  
45.                 VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factorf))  
46.                 VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factorf))  
47.                   
48.                 #handle dividends  
49.                 div=dividends[1][np.nonzero(dividendsStep==(k))[0]]    
50.                   
51.                 #a call is an always increasing function of the strike  
52.                 mid=int(N/2)  
53.                 delta=VC[1:mid]-VC[:mid-1]  
54.                 BoS=np.where(delta<0)[0][-1]  
55.                 VC[:BoS]=0  
56.                   
57.                 #The put always decreasing  
58.                 delta=VP[1:mid]-VP[:mid-1]  
59.                 BoS=np.where(delta<0)[0][0]  
60.                 VP[:BoS]=VP[BoS]  
61.         
62.                 tckC=si.splrep(s,VC)  
63.                 tckP=si.splrep(s,VP)  
64.                 #option prices are interpolated on the grid  
65.                 VC=si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der=0)              
66.                 VP=si.splev(np.maximum(s-div,np.exp(x_min)) ,tckP,der=0)   
67.                   
68.                 if american==True:       
69.                       VC=earlyCall(VC,s,K)  
70.                       VP=earlyPut(VP,s,K)  
71.                   
72.                 #finish the step  
73.                 dtf=fraction   
74.                 char_exp_factorf = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dtf)  
75.                 VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factorf))  
76.                 VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factorf))  
77.                  
78.             else:    
79.                 VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
80.                 VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))  
81.                   
82.                 div=dividends[1][np.nonzero(dividendsStep==(k))[0]]    
83.                   
84.                 #a call is an always increasing function of the strike  
85.                 mid=int(N/2)  
86.                 delta=VC[1:mid]-VC[:mid-1]  
87.                 BoS=np.where(delta<0)[0][-1]  
88.                 VC[:BoS]=0  
89.                   
90.                 #The put always decreasing  
91.                 delta=VP[1:mid]-VP[:mid-1]  
92.                 BoS=np.where(delta<0)[0][0]  
93.                 VP[:BoS]=VP[BoS]  
94.                  
95.                 tckC=si.splrep(s,VC)  
96.                 tckP=si.splrep(s,VP)  
97.                 #option prices are interpolated on the grid  
98.                 VC=si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der=0)              
99.                 VP=si.splev(np.maximum(s-div,np.exp(x_min)) ,tckP,der=0)   
100.                 if american==True:       
101.                       VC=earlyCall(VC,s,K)  
102.                       VP=earlyPut(VP,s,K)  
103.        else:  
104.            #no dividend at this step  
105.            VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
106.            VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))      
107.                   
108.            if american==True:       
109.                VP=earlyPut(VP,s,K)  
110.                VC=earlyCall(VC,s,K)  
111.                  
112.          
113.    #Interpolate option prices  
114.    tck=si.splrep(s,VC)  
115.    option.call.price=si.splev(S0,tck,der=0)     
116.    tck=si.splrep(s,VP)  
117.    option.put.price=si.splev(S0,tck,der=0)   

Compared to the adaptive stepping for the tree where we try to make the step match a dividend, here we cut the step where the dividend falls into two sub steps of appropriate size to make sure that the dividend payment falls right where it should. The good thing is that we can do that for an arbitrary number of dividends.

Assume now we have two late (2.8 and 2.9 years for a 3 year option) and large (5€ each) dividends. Compare the AS Fourier and AS tree on the following picture.

We still have the little waves but there is no comparison with the tree. I don't reproduce figures but similar sawtooth pattern can be seen with the puts (much smaller though) and AS improves the situation here as well.

One should note that I also included the time structure of rates as well. We could do some adaptive stepping for the rates as well but I suppose the error coming from here is very small in our case.