Convergence issues

EDIT: the previous version of this post contained a small mistake (machine precision related I think). I rewrote it and updated the figures.

As I pointed out earlier, one need to check the code before any use in production. One practical aspect is deciding the number of steps to use. There is obviously a trade off between precision and calculation time.

One also needs to make sure that increasing the number of steps actually decreases the error. We have seen already that for the binomial tree for european option, it is not always the case. We have seen that the price oscillate.

One can see on the following image (taken from that older post):

Pricing without regularization using 15 steps is worse than using 10 steps. For american options, one cannot check against the Black and Scholes price but still convergence is important. In the case the dividend is late, we can see oscillations.

For the following graph, we have $S0=100$, $K=100$, $r=5%$, $\sigma=40%$ and $T=3$. We have a dividend at $T=2.9$ for 1 euro.

In the tree, the dividend is paid on a step but there will be a mismatch between the real dividend date and the step. This makes the price oscillate. As the step size goes toward 0, that mismatch will go toward 0 and the oscillations go away but for smaller number of steps this generates oscillations.

I try two methods to fix this. Once I know the mismatch between the step where the dividend will be paid and the real dividend date I can either:
1/increase the step size so that the step number multiplied by the step size falls right on the dividend date.
2/decrease the step size so that the dividend date falls right onto the next step.

In both cases if I keep the number of steps unchanged, the maturity of the option will be wrong. If I increase the step size I will need less steps for the same maturity, opposite if I decrease the size. And I will need a fractional step to reach exactly the maturity I want. That can pose problem for my recombining tree.

But remember that regularization is replacing the last step of the tree by a BS formula. I can use that and change the time used to have my "fraction" of step. One can see on the following graph the effect of this change. Note that I start plotting from 150 steps.

For a more global view, we see convergence of the 3 methods so, in theory, all methods works.

The adapted time step oscillate but much much less than the original.
Let's now see what happens for european option using such method. This graph shows a large point of view:

 

and this one is a close up

A few comments are in order:
It seems the adjustment of the time step creates a shift in the european call price. I don't see why at the moment. That difference is really small though. Increasing the time step size seem to introduce a downward bias so I would choose the technique decreasing the size of the time step.

I did a quick check and results are similar with ITM options.
Find enclosed the code for the decreased size version:

1. def BinomialTree8c(option, params):    
2.     
3.     #we pass two objects: one for the option values and one for the tree parameters    
4.     N=params.N    
5.     [S0,K,sigma,T,r2,divi,american]=[option.S0,option.K,option.sigma,option.T,option.r,option.divi,option.american]    
6.         
7.     #calculate delta T        
8.     deltaT = float(T) / N    
9.     
10.      
11.         
12.     #we improve the previous tree by generalization of payoff    
13.     def expiCall(stock,strike):    
14.         return(np.maximum(stock[:]-strike[:][0],0))    
15.             
16.     def expiPut(stock,strike):    
17.         return(np.maximum(strike[:][0]-stock[:],0))    
18.         
19.     def earlyCall(Option,stock,strike):    
20.         return(np.maximum(stock[:]-strike[:][0],Option))    
21.         
22.     def earlyEuro(Option,stock,strike):    
23.         return(np.maximum(Option,0))    
24.         
25.     def earlyPut(Option,stock,strike):    
26.         return(np.maximum(strike[:][0]-stock[:],Option))        
27.             
28.          #regularisation will be usefull for delta calculations    
29.     def expiCallR(stock,strike,rr=0,vol=sigma,time=deltaT):    
30.         return(BlackScholes('C',stock,strike,rr,vol,time))    
31.         
32.     def expiPutR(stock,strike,rr=0,vol=sigma,time=deltaT):    
33.         return(BlackScholes('P',stock,strike,rr,vol,time))    
34.             
35.      
36.     if american=='false':    
37.         earlyCall=earlyEuro    
38.         earlyPut=earlyEuro    
39.     
40.        
41.     r=[[],[]]    
42.     dividends=[[],[]]            
43.     
44.     #2 columns vector, one for amount one for time.    
45.     #we make sure we don't take into account dividend happening after expiration    
46.     if (np.size(divi)>0 and divi[0][0]<T) :    
47.         lastdiv=np.nonzero(np.array(divi[0][:])<=T)[0][-1]    
48.         dividends[0]=divi[0][:lastdiv+1]            
49.         dividends[1]=divi[1][:lastdiv+1]     
50.       
51.   
52.     if (params.stepAdapt=='true'and np.size(divi)>0):      
53.         #in order to have the last dividend fall on a step we modify the length of the option  
54.         #we correct via regularisation with a smaller time step  
55.         #first we need the mismatch  
56.         #print deltaT,N,deltaT*N  
57.         TimeM=np.int(dividends[0][-1]/deltaT)  
58.         TimeMismatch=(dividends[0][-1]-TimeM*deltaT)-0.00001*deltaT  
59.         deltaT=deltaT-(deltaT-TimeMismatch)/TimeM  
60.         ExpMismatchInt=int((T-deltaT*N)/deltaT)   
61.         ExpMismatch=ExpMismatchInt-(T-deltaT*N)/deltaT  
62.         N=N+ExpMismatchInt  
63.         #print deltaT,N,deltaT*N,deltaT*N-ExpMismatch*deltaT  
64.           
65.     else:  
66.         ExpMismatch=0  
67.         ExpMismatchInt=0  
68.                   
69.     #BurnIn are 2 extra steps we add.    
70.     #that's used for calculation of delta,    
71.     if params.BurnIn=='true':    
72.         BI=2    
73.         N=N+BI    
74.     else:    
75.         BI=0        
76.     #similar structure for rates    
77.     # rates with date d1 is use until date d1 so we need the first date post maturity.       
78.     if np.size(r2)>0 :    
79.         lastrate=np.nonzero(np.array(r2[0][:])<=T)[0][-1]    
80.         r[0]=r2[0][:lastrate+2]            
81.         r[1]=r2[1][:lastrate+2]           
82.          
83.          
84.     #Transform the dividend date into a step    
85.     if np.size(dividends)>0:    
86.         dividendsStep=np.floor(np.multiply(dividends[0],1/deltaT))    
87.     else:    
88.         dividendsStep=[]        
89.         
90.     if np.size(r)>0:    
91.         rateChangeStep=np.floor(np.multiply(r[0],1/deltaT))    
92.     else:    
93.         rateChangeStep=[]          
94.         
95.     #we add the BurnIn to the dividend step for correction    
96.        
97.     dividendsStep=dividendsStep+np.asarray( [float(BI) for i in xrange(np.size(dividends[1][:]))])    
98.     rateChangeStep=rateChangeStep+np.asarray( [float(BI) for i in xrange(np.size(r[1][:]))])    
99.             
100.         
101.     #prepare the array for discount factors used for the discounting of dividends    
102.     DFArray=np.zeros((N+2,1),float)    
103.     rArray=np.zeros((N+2,1),float)       
104.         
105.     #get the steps where the rate changes    
106.     prevstep=0    
107.     prevrate=0    
108.     for step in range(len(rateChangeStep)-1):    
109.         DFArray[prevstep:rateChangeStep[step]]=np.exp(-prevrate*deltaT)    
110.         DFArray[rateChangeStep[step]:]=np.exp(-r[1][step]*deltaT)    
111.         rArray[prevstep:rateChangeStep[step],0]=prevrate    
112.         rArray[rateChangeStep[step]:,0]=r[1][step]    
113.         prevstep=rateChangeStep[step]    
114.         prevrate=r[1][step]    
115.             
116.             
117.     rArray[0]=rArray[1]       
118.     pvdividends=0    
119.         
120.     #get present value of the dividend        
121.     if np.size(dividends)>0:    
122.         pvdividends=np.sum(dividends[1])    
123.        #for index,dv in enumerate(list(dividends[1])):    
124.        #     pvdividends=pvdividends+dv*np.product(DFArray[0:np.int_(dividendsStep[index])])    
125.             
126.     else:    
127.         pvdividends=0    
128.         
129.     S0=S0-pvdividends    
130.     currentDividend=0    
131.         
132.     # up and down factor will be constant for the tree so we calculate outside the loop    
133.     u = np.exp(sigma * np.sqrt(deltaT))    
134.     d = 1.0 / u    
135.      
136.      
137.     #if we use regularisation, one need one less time step as we use BS for the last one.    
138.     if params.regul=='true':    
139.         N=N-1  
140.          
141.     #to work with vector we need to init the arrays using numpy    
142.     # we do call and put together     
143.     fsC =  np.asarray([0.0 for k in xrange(N + 1)])    
144.     fsP =  np.asarray([0.0 for k in xrange(N + 1)])    
145.             
146.     #we need the stock tree for calculations of expiration values    
147.     fs2 = np.asarray([(S0 * u**k * d**(N - k)) for k in xrange(N + 1)])    
148.         
149.         
150.     #we vectorize the strikes as well so the expiration check will be faster    
151.     fs3 =np.asarray( [float(K) for k in xrange(N + 1)])    
152.          
153.     # Compute the leaves, f_{N, j}    
154.     if params.regul=='true':    
155.         fsC = expiCallR(fs2,fs3,rArray[-1][0],sigma,(1.0-ExpMismatch)*deltaT)    
156.         fsP = expiPutR(fs2,fs3,rArray[-1][0],sigma,(1.0-ExpMismatch)*deltaT)    
157.     else:    
158.         fsC = expiCall(fs2,fs3)    
159.         fsP = expiPut(fs2,fs3)    
160.         
161.     #we get the rate    
162.     df=rArray[N-1]*deltaT    
163.     discountfactor=np.exp(-df)    
164.     p=(np.exp(df)-d)/(u-d)    
165.     divi=0    
166.        
167.     #calculate backward the option prices    
168.     #if we use BurnIn, we don't need to calculate the extra steps for the price     
169.     #but we need them for theta    
170.     for k in xrange(N-1, -1+BI, -1):    
171.            
172.           
173.        if (k in rateChangeStep[:]):    
174.               
175.           df= (r[1][np.nonzero(rateChangeStep==(k))[0]-1])*deltaT    
176.           discountfactor= np.exp(-df)          
177.           p=(np.exp(df)-d)/(u-d)    
178.              
179.        fsC[:-1]=discountfactor* (p * fsC[1:] + (1-p) * fsC[:-1])    
180.        fsP[:-1]=discountfactor* (p * fsP[1:] + (1-p) * fsP[:-1])    
181.        
182.        fs2[:]=fs2[:]*u    
183.           
184.          
185.        #we add back the dividends in current dividends as we move backward, useful    
186.        #only for american options    
187.        if (k in dividendsStep):    
188.                 
189.             divi=dividends[1][np.nonzero(dividendsStep==(k))[0]]      
190.             currentDividend=currentDividend+divi    
191.             #the VN method requires an interpolation               
192.             tckC=si.splrep(fs2[:],fsC[:],k=1)    
193.             tckP=si.splrep(fs2[:],fsP[:],k=1)    
194.             #stock tree is recalculated as if the dividends were in from the start    
195.             #carreful with indices we are in step i so we calculated N-i steps already                
196.             fs2[:]=np.asarray([((S0+currentDividend) * u**(j) * d**(- j+k)) for j in xrange(N + 1)])    
197.                 
198.             #option prices are interpolated on this tree and backward calculations continues    
199.             fsC[:]=si.splev(fs2[:]-divi,tckC,der=0)                
200.             fsP[:]=si.splev(fs2[:]-divi,tckP,der=0)    
201.                             
202.                     
203.        fsC[:]=earlyCall(fsC[:],fs2[:],fs3[:])    
204.        fsP[:]=earlyPut(fsP[:],fs2[:],fs3[:])    
205.    
206.        
207.     if params.BurnIn=='true':    
208.         option.call.price=fsC[BI/2]    
209.         option.put.price=fsP[BI/2]    
210.         option.call.delta=(fsC[BI/2+1]-fsC[BI/2-1])/(fs2[BI/2+1]-fs2[BI/2-1])    
211.         option.put.delta=(fsP[BI/2+1]-fsP[BI/2-1])/(fs2[BI/2+1]-fs2[BI/2-1])    
212.         s2a=fs2[BI/2+1]-fs2[BI/2]    
213.         s2b=fs2[BI/2]-fs2[BI/2-1]    
214.         option.call.gamma=((fsC[BI/2+1]-fsC[BI/2])/s2a +(-fsC[BI/2]+fsC[BI/2-1])/s2b)/((fs2[BI/2+1]-fs2[BI/2-1])/2)    
215.         option.put.gamma=((fsP[BI/2+1]-fsP[BI/2])/s2a +(-fsP[BI/2]+fsP[BI/2-1])/s2b)/((fs2[BI/2+1]-fs2[BI/2-1])/2)    
216.         
217.        
218.         #final steps are done outside the loop for theta calculation only,    
219.         #no american, rate or divs in there    
220.         for i in xrange(-1+BI, -1, -1):    
221.             
222.             fsC[:-1]=discountfactor* (p * fsC[1:] + (1-p) * fsC[:-1])    
223.             fsP[:-1]=discountfactor* (p * fsP[1:] + (1-p) * fsP[:-1])    
224.             
225.         option.call.theta=(option.call.price-fsC[0])/(BI*deltaT)    
226.         option.put.theta=(option.put.price-fsP[0])/(BI*deltaT)    
227.          
228.     else:                    
229.         option.call.price=fsC[0]    
230.         option.put.price=fsP[0]    
231.     return option   

In fact not much changes: the size time step is adapted in a specific if statement, the time used in the regularization as well. The rest stays all the same.

Note that if I have several late dividends, the method will adjust only to the latest one, one could still see oscillations. Here with dividends at 2.8 and 2.9