Fourier with cash dividends

This one took me a bit of time to figure out. I wanted to apply the same technique we used with trees and PDE's to price options with cash dividends to Fourier pricing. It kinda of work but not always.

That is it doesn't give plain wrong values except for very large dividends. The problem is linked to the Fourier transform itself. See, the Fourier theory usually expects signal going to $-\infty$ to $\infty$. In our case, our log stock cannot cover this range and we have to cut it. However cutting introduce some border effects. Such effect are amplified as we repeatedly use Ft's and inverse Ft's  for each time step.

In practice that means that my call value, as a function of the stock, is monotonically increasing except for very very large stock values and very very small ones. For European options without dividends, that doesn't matter much as it is extreme values.
For cash dividends though we interpolate the option value to ensure the following condition:
\[
V(S,t_d^-)=V(S-D,t_d^+)
\]
with $V(S,t_d^-)$ the value of the option right before the dividend is paid and  $V(S-D,t_d^+)$ right after. $D$ is the dividend. Indeed it must be equal: the value of the derivative is not going to change when the dividend gets paid (unless we have early exercise possible of course).

Here comes the trouble: the very very low stock values for which we have a border effect have a large effect on this interpolation. The high ones do not matter much because $D$ is always positive.

For the call, the border effect makes the call price increase, then decrease. Then it increases, as it should, with the stock. Before doing the interpolation, I find the part where it stops decreasing and replace everything up to that point with 0. For the upside part, after a point, it starts decreasing. I could replace it by the intrinsic value but it doesn't seem to make a difference.
If we price an american option it will be replaced anyway.

For the put, the value of the derivative starts by increasing then decreases, as it should. I replace this increasing part by a flat part. If we would price american options this part would be exercised early and the problem would go away.
I just do this replacement right before the interpolation for the dividend. As the main problem is the effect of these borders on the interpolation.

Find the code below:

1. def FourierST5(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.    dt=float(T/NTS)  
23.    dividends,dividendsStep=preparedivs(divi,float(T),dt)  
24.        
25.   
26.    # FST method  
27.    char_exp_factor = np.exp((j*(r-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-r)*dt)  
28.    for k in range(NTS,0,-1):     
29.        VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
30.        VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))  
31.          
32.         
33.        #if there is a dividend, we interpolate the grid  
34.        if (k in dividendsStep):  
35.             div=dividends[1][np.nonzero(dividendsStep==(k))[0]]    
36.               
37.             #a call is an always increasing function of the strike  
38.             mid=int(N/2)  
39.             delta=VC[1:mid]-VC[:mid-1]  
40.             BoS=np.where(delta<0)[0][-1]  
41.             VC[:BoS]=0  
42.               
43.             #The put always decreasing  
44.             delta=VP[1:mid]-VP[:mid-1]  
45.             BoS=np.where(delta<0)[0][0]  
46.             VP[:BoS]=VP[BoS]  
47.               
48.             tckC=si.splrep(s,VC)  
49.             tckP=si.splrep(s,VP)  
50.             #option prices are interpolated on the grid  
51.             VC=si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der=0)              
52.             VP=si.splev(np.maximum(s-div,np.exp(x_min)) ,tckP,der=0)   
53.                           
54.        if american==True:       
55.            VC=earlyCall(VC,s,K)  
56.            VP=earlyPut(VP,s,K)  
57.          
58.    #Interpolate option prices  
59.    tck=si.splrep(s,VC)  
60.    option.call.price=si.splev(S0,tck,der=0)     
61.    tck=si.splrep(s,VP)  
62.    option.put.price=si.splev(S0,tck,der=0)   

A couple of comments:
I changed the loop to show the backward time marching. It is consistent with the other models.
I still need to add a time structure of rates, that shouldn't be a problem.
For the Fourier pricer and the PDE, we don't need the step to fall right on the dividend. If, say the real dividend payment is half between two steps, we can do a convolution with $dt/2$, manage the dividend and then do another $dt/2$ convolution. For the PDE it is the same idea. It should reduce the oscillations we saw for the binomial tree. I did not check if that's needed yet though.

In terms of performance, we look at a 3 year option with $\sigma$=40%, $r$=5%, $K$=100 and $S=100$.
If I price a european options with a 10 euros dividend to be paid in 2.5 years, I find $c=28.657$ and $p=23.553$ for the binomial tree, $c=28.655$ and $p=23.545$ for the PDE and $c=28.649$ and $p=23.544$ for Fourier.
If I place the dividend in 0.1 instead, I find $c=25.995$ and $p=22.017$ for the binomial tree, $c=25.995$ and $p=22.011$ for the PDE and $c=25.988$ and $p=22.009$ for Fourier.
The values are not the same, but quite close, that's a good sign for such a large dividend.
(Fourier use 2^10 for NAS and 10k for NTS, PDE NAS=1000 and the tree 1k steps).