PDEPricer class

The PDE pricer is also now available as a separate class. As for the tree, the class is large because I have 4 variations of the same code, for handling duo and being a bit faster with european options by not checking exercise.

If the goal is readability, one pricing code with a few if/then would work as well.

1. import scipy.interpolate as si  
2.   
3.   
4. from OptionPricer import *  
5.   
6. class PDEPricer(OptionPricer):  
7.     _NumberSteps = 0  
8.     _NumberAssetSteps = 0  
9.       
10.     _StepAdapt = False  
11.     _Regul = False  
12.   
13.     @staticmethod  
14.     def preparedivs(divi, T, dt):  
15.         dividends = [[], []]   
16.         #2 columns vector, one for amount one for time.  
17.         #we make sure we don't take into account dividend happening after expiration  
18.      
19.         if (np.size(divi)>0 and divi[0][0]<T) :  
20.             lastdiv = np.nonzero(np.array(divi[0][:])<= T)[0][-1]  
21.             dividends[0] = divi[0][:lastdiv+1]          
22.             dividends[1] = divi[1][:lastdiv+1]   
23.       
24.         #Transform the dividend date into a step  
25.         if np.size(dividends)>0:  
26.             dividendsStep = np.floor(np.multiply(dividends[0], 1/dt))  
27.         else:  
28.             dividendsStep = []   
29.         return dividends, dividendsStep  
30.   
31.     @staticmethod  
32.     def preparerates(r2, T, dt, NTS):  
33.         r = [[], []]  
34.           
35.         #similar structure for rates  
36.         # rates with date d1 is use until date d1 so we need the first date post maturity.     
37.         if np.size(r2)>0 :  
38.             lastrate = np.nonzero(np.array(r2[0][:])>= T)[0]+1  
39.             #i take the one just after expiration  
40.             r[0] = r2[0][:lastrate]          
41.             r[1] = r2[1][:lastrate]         
42.         if np.size(r)>0:  
43.             rateChangeStep = np.floor(np.multiply(r[0], 1/dt))  
44.         else:  
45.             rateChangeStep = []        
46.           
47.         #prepare the array for rates used for the PDE  
48.         rArray = np.zeros((NTS+2, 1), float)     
49.         DFArray = np.zeros((NTS+2, 1), float)  
50.       
51.         #get the steps where the rate changes  
52.         prevstep = 0  
53.           
54.         for step in range(len(rateChangeStep)):  
55.               
56.             DFArray[prevstep:min(rateChangeStep[step], NTS+2)] = np.exp(-r[1][step]*dt)  
57.             rArray[prevstep:min(rateChangeStep[step], NTS+2)] = r[1][step]  
58.             prevstep = rateChangeStep[step]  
59.   
60.         return r, rArray, rateChangeStep, DFArray  
61.   
62.     @staticmethod  
63.     def PDEEUDuo(pricingParameters, optionContract, NTS, NAS, regul, stepAdapt):  
64.               
65.         #parameters  
66.         S0 = pricingParameters._S0  
67.         r2 = pricingParameters._r2  
68.         T = pricingParameters._T         
69.         sigma = optionContract._Sigma  
70.         divi = pricingParameters._div  
71.         K = optionContract._K  
72.           
73.         #calculate delta T      
74.         deltaT = float(T) / NTS  
75.    
76.         Zmin = -7.5  
77.         Zmax = 7.5  
78.         
79.         dZ = float((Zmax-Zmin)/NAS)     
80.           
81.         #time step is constraint by stability to be a function  of dZ  
82.         #but if the user wants a smaller time step, that is fine too  
83.         dt = min(dZ**2/(3*sigma*sigma), float(T)/NTS)  
84.           
85.         dividends, dividendsStep = PDEPricer.preparedivs(divi, T, deltaT)  
86.         r, rArray, rateChangeStep, DFArray = PDEPricer.preparerates(r2, T, dt, NTS)  
87.        
88.               
89.         #we get the rate  
90.         rr = rArray[-1][0]  
91.         NTS = int(T/dt)+1  
92.           
93.         a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
94.         b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
95.         c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
96.              
97.          
98.         VC = np.asarray( [0 for i in xrange(int(NAS)+1)])  
99.         VP = np.asarray( [0 for i in xrange(int(NAS)+1)])  
100.         Z = np.array([(Zmin+i*dZ) for i in xrange(int(NAS)+1)])  
101.           
102.           
103.         if regul:  
104.             VC = optionContract.expiCR(np.exp(Z), rr, sigma, dt)  
105.             VP = optionContract.expiPR(np.exp(Z), rr, sigma, dt)  
106.             NTS = NTS-1  
107.         else:  
108.             VC = optionContract.expiC(np.exp(Z))  
109.             VP = optionContract.expiP(np.exp(Z))  
110.       
111.         a1 = (Z[NAS]-Z[NAS-1])  
112.         b1 = (Z[NAS-1]-Z[NAS-2])  
113.               
114.       
115.         for k in range(NTS, -1, -1):  
116.               
117.             #we need to change the a, b, c only when rate changes  
118.             if (k in rateChangeStep[:]):  
119.                 rr = (r[1][np.nonzero(rateChangeStep == (k))[0]])  
120.                 a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
121.                 b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
122.                 c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
123.               
124.             if (k in dividendsStep):  
125.                 if stepAdapt:  
126.                       
127.                     #we perform a fraction of step for the pde, place the dividend and finish the step  
128.                     #so no matter the step size, the div is at the right place.  
129.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
130.                     fraction = -(k*dt-dividends[0][np.nonzero(dividendsStep == (k))[0]])  
131.                     dtf = dt-fraction  
132.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
133.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
134.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
135.               
136.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
137.                     VP[1:-1] = af*VP[:-2]+bf*VP[1:-1]+cf*VP[2:]  
138.                  
139.                     VC[0] = 0  
140.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
141.                     VP[0] = K  
142.                     VP[NAS] = 0  
143.                       
144.                     tckC = si.splrep(Z, VC, k = 1)  
145.                     tckP = si.splrep(Z, VP, k = 1)  
146.                     #option prices are interpolated on the grid  
147.                     VC = si.splev(np.log(np.maximum(np.exp(Z)-div, np.exp(Zmin))), tckC, der = 0)              
148.                     VP = si.splev(np.log(np.maximum(np.exp(Z)-div, np.exp(Zmin))), tckP, der = 0)             
149.                       
150.                                    
151.                     VC[0] = 0  
152.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
153.                     VP[0] = K  
154.                     VP[NAS] = 0  
155.                       
156.                     #we finish the step  
157.                     dtf = fraction  
158.                      
159.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
160.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
161.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
162.               
163.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
164.                     VP[1:-1] = af*VP[:-2]+bf*VP[1:-1]+cf*VP[2:]  
165.                  
166.                 else:  
167.                     #otherwise perform the step and then place the dividend  
168.                      #call  
169.                     VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
170.                     VC[0] = VC[0]*(1-rr*dt)  
171.                       
172.                     #put  
173.                     VP[1:-1] = a*VP[:-2]+b*VP[1:-1]+c*VP[2:]  
174.                       
175.                     VC[0] = 0  
176.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
177.                     VP[0] = K  
178.                     VP[NAS] = 0  
179.                       
180.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
181.                       
182.                     tckC = si.splrep(np.exp(Z), VC, k = 1)  
183.                     tckP = si.splrep(np.exp(Z), VP, k = 1)  
184.                     #option prices are interpolated on the grid  
185.                     VC = si.splev(np.maximum(np.exp(Z)-div, np.exp(Zmin)), tckC, der = 0)              
186.                     VP = si.splev(np.maximum(np.exp(Z)-div, np.exp(Zmin)), tckP, der = 0)    
187.                      
188.                       
189.             else:  
190.                 #there is no dividends  
191.                 #call  
192.                 VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
193.                #put  
194.                 VP[1:-1] = a*VP[:-2]+b*VP[1:-1]+c*VP[2:]  
195.                    
196.             VC[0] = 0  
197.             VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
198.             VP[0] = K  
199.             VP[NAS] = 0  
200.             #check early exercise  
201.             #VC = earlyCall(VC, np.exp(Z), Kk)  
202.             #VP = earlyPut(VP, np.exp(Z), Kk)  
203.               
204.         #get the prices and the grid from the vector      
205.         tck = si.splrep(Z, VC)  
206.         optionContract.call.price = si.splev(np.log(S0), tck, der = 0)   
207.         p1 = si.splev(np.log(S0+1), tck, der = 0)  
208.         p2 = si.splev(np.log(S0-1), tck, der = 0)  
209.         optionContract.call.delta = (p1-p2)/2  
210.         optionContract.call.gamma = p1+p2-2*optionContract.call.price  
211.         VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
212.         tck = si.splrep(Z, VC)  
213.         optionContract.call.theta = (optionContract.call.price-si.splev(np.log(S0), tck, der = 0) )/dt  
214.           
215.           
216.         tck = si.splrep(Z, VP)  
217.         optionContract.put.price = si.splev(np.log(S0), tck, der = 0)   
218.         p1 = si.splev(np.log(S0+1), tck, der = 0)  
219.         p2 = si.splev(np.log(S0-1), tck, der = 0)  
220.         optionContract.put.delta = (p1-p2)/2  
221.         optionContract.put.gamma = p1+p2-2*optionContract.put.price  
222.         VP[1:-1] = a*VP[:-2]+b*VP[1:-1]+c*VP[2:]  
223.         tck = si.splrep(Z, VP)  
224.         optionContract.put.theta = (optionContract.put.price-si.splev(np.log(S0), tck, der = 0) )/dt  
225.           
226.           
227.         return VC        
228.           
229.     @staticmethod  
230.     def PDEEUSolo(pricingParameters, optionContract, NTS, NAS, regul, stepAdapt, option, payoff):  
231.               
232.         #parameters  
233.         S0 = pricingParameters._S0  
234.         r2 = pricingParameters._r2  
235.         T = pricingParameters._T         
236.         sigma = optionContract._Sigma  
237.         divi = pricingParameters._div  
238.           
239.         #calculate delta T      
240.         deltaT = float(T) / NTS  
241.    
242.         Zmin = -7.5  
243.         Zmax = 7.5  
244.         
245.         dZ = float((Zmax-Zmin)/NAS)     
246.           
247.         #time step is constraint by stability to be a function  of dZ  
248.         #but if the user wants a smaller time step, that is fine too  
249.         dt = min(dZ**2/(3*sigma*sigma), float(T)/NTS)  
250.           
251.         dividends, dividendsStep = PDEPricer.preparedivs(divi, T, deltaT)  
252.         r, rArray, rateChangeStep, DFArray = PDEPricer.preparerates(r2, T, dt, NTS)  
253.        
254.               
255.         #we get the rate  
256.         rr = rArray[-1][0]  
257.         NTS = int(T/dt)+1  
258.           
259.         a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
260.         b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
261.         c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
262.              
263.          
264.         VC = np.asarray( [0 for i in xrange(int(NAS)+1)])  
265.         Z = np.array([(Zmin+i*dZ) for i in xrange(int(NAS)+1)])  
266.           
267.           
268.         if regul:  
269.             VC = payoff(np.exp(Z), rr, sigma, dt)  
270.             NTS = NTS-1  
271.         else:  
272.             VC = payoff(np.exp(Z))  
273.               
274.         a1 = (Z[NAS]-Z[NAS-1])  
275.         b1 = (Z[NAS-1]-Z[NAS-2])  
276.               
277.       
278.         for k in range(NTS, -1, -1):  
279.               
280.             #we need to change the a, b, c only when rate changes  
281.             if (k in rateChangeStep[:]):  
282.                 rr = (r[1][np.nonzero(rateChangeStep == (k))[0]])  
283.                 a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
284.                 b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
285.                 c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
286.               
287.             if (k in dividendsStep):  
288.                 if stepAdapt:  
289.                       
290.                     #we perform a fraction of step for the pde, place the dividend and finish the step  
291.                     #so no matter the step size, the div is at the right place.  
292.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
293.                     fraction = -(k*dt-dividends[0][np.nonzero(dividendsStep == (k))[0]])  
294.                     dtf = dt-fraction  
295.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
296.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
297.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
298.               
299.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
300.                       
301.                     VC[0] = VC[0] = VC[0]*(1-rr*dt)  
302.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
303.                       
304.                     tckC = si.splrep(Z, VC, k = 1)  
305.                     #option prices are interpolated on the grid  
306.                     VC = si.splev(np.log(np.maximum(np.exp(Z)-div, np.exp(Zmin))), tckC, der = 0)              
307.                       
308.                                    
309.                     VC[0] = VC[0]*(1-rr*dt)  
310.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
311.                       
312.                     #check for exercise  
313.                     #VC = earlyCall(VC, np.exp(Z), Kk)  
314.                       
315.                     #we finish the step  
316.                     dtf = fraction  
317.                      
318.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
319.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
320.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
321.               
322.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
323.                       
324.                  
325.                 else:  
326.                     #otherwise perform the step and then place the dividend  
327.                      #call  
328.                     VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
329.                     VC[0] = VC[0]*(1-rr*dt)  
330.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
331.                       
332.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
333.                       
334.                     tckC = si.splrep(np.exp(Z), VC, k = 1)  
335.                     #option prices are interpolated on the grid  
336.                     VC = si.splev(np.maximum(np.exp(Z)-div, np.exp(Zmin)), tckC, der = 0)              
337.                       
338.             else:  
339.                 #there is no dividends  
340.                 #call  
341.                 VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
342.                  
343.                    
344.             VC[0] = VC[0]*(1-rr*dt)  
345.             VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
346.               
347.             #check early exercise  
348.             #VC = earlyCall(VC, np.exp(Z), Kk)  
349.               
350.               
351.         #get the prices and the grid from the vector      
352.         tck = si.splrep(Z, VC)  
353.         option.price = si.splev(np.log(S0), tck, der = 0)   
354.         p1 = si.splev(np.log(S0+1), tck, der = 0)  
355.         p2 = si.splev(np.log(S0-1), tck, der = 0)  
356.         option.delta = (p1-p2)/2  
357.         option.gamma = p1+p2-2*option.price  
358.         VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
359.         tck = si.splrep(Z, VC)  
360.         option.theta = (option.price-si.splev(np.log(S0), tck, der = 0) )/dt  
361.                  
362.           
363.         return VC        
364.     @staticmethod  
365.     def PDEUSDuo(pricingParameters, optionContract, NTS, NAS, regul, stepAdapt):  
366.               
367.         #parameters  
368.         S0 = pricingParameters._S0  
369.         r2 = pricingParameters._r2  
370.         T = pricingParameters._T         
371.         sigma = optionContract._Sigma  
372.         K = optionContract._K  
373.         divi = pricingParameters._div  
374.           
375.         #calculate delta T      
376.         deltaT = float(T) / NTS  
377.    
378.         Zmin = -7.5  
379.         Zmax = 7.5  
380.         
381.         dZ = float((Zmax-Zmin)/NAS)     
382.           
383.         #time step is constraint by stability to be a function  of dZ  
384.         #but if the user wants a smaller time step, that is fine too  
385.         dt = min(dZ**2/(3*sigma*sigma), float(T)/NTS)  
386.           
387.         dividends, dividendsStep = PDEPricer.preparedivs(divi, T, deltaT)  
388.         r, rArray, rateChangeStep, DFArray = PDEPricer.preparerates(r2, T, dt, NTS)  
389.        
390.               
391.         #we get the rate  
392.         rr = rArray[-1][0]  
393.         NTS = int(T/dt)+1  
394.           
395.         a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
396.         b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
397.         c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
398.              
399.          
400.         VC = np.asarray( [0 for i in xrange(int(NAS)+1)])  
401.         VP = np.asarray( [0 for i in xrange(int(NAS)+1)])  
402.         Z = np.array([(Zmin+i*dZ) for i in xrange(int(NAS)+1)])  
403.           
404.           
405.         if regul:  
406.             VC = optionContract.expiCR(np.exp(Z), rr, sigma, dt)  
407.             VP = optionContract.expiPR(np.exp(Z), rr, sigma, dt)  
408.             NTS = NTS-1  
409.         else:  
410.             VC = optionContract.expiC(np.exp(Z))  
411.             VP = optionContract.expiP(np.exp(Z))  
412.       
413.         a1 = (Z[NAS]-Z[NAS-1])  
414.         b1 = (Z[NAS-1]-Z[NAS-2])  
415.               
416.       
417.         for k in range(NTS, -1, -1):  
418.               
419.             #we need to change the a, b, c only when rate changes  
420.             if (k in rateChangeStep[:]):  
421.                 rr = (r[1][np.nonzero(rateChangeStep == (k))[0]])  
422.                 a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
423.                 b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
424.                 c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
425.               
426.             if (k in dividendsStep):  
427.                 if stepAdapt:  
428.                       
429.                     #we perform a fraction of step for the pde, place the dividend and finish the step  
430.                     #so no matter the step size, the div is at the right place.  
431.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
432.                     fraction = -(k*dt-dividends[0][np.nonzero(dividendsStep == (k))[0]])  
433.                     dtf = dt-fraction  
434.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
435.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
436.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
437.               
438.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
439.                     VP[1:-1] = af*VP[:-2]+bf*VP[1:-1]+cf*VP[2:]  
440.                  
441.                     VC[0] = 0  
442.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
443.                     VP[0] = K  
444.                     VP[NAS] = 0  
445.                       
446.                     tckC = si.splrep(Z, VC, k = 1)  
447.                     tckP = si.splrep(Z, VP, k = 1)  
448.                     #option prices are interpolated on the grid  
449.                     VC = si.splev(np.log(np.maximum(np.exp(Z)-div, np.exp(Zmin))), tckC, der = 0)              
450.                     VP = si.splev(np.log(np.maximum(np.exp(Z)-div, np.exp(Zmin))), tckP, der = 0)             
451.                       
452.                                    
453.                     VC[0] = 0  
454.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
455.                     VP[0] = K  
456.                     VP[NAS] = 0  
457.                       
458.                     #check for exercise  
459.                     VC = optionContract.earlyC(VC, np.exp(Z))  
460.                     VP = optionContract.earlyP(VP, np.exp(Z))  
461.                       
462.                     #we finish the step  
463.                     dtf = fraction  
464.                      
465.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
466.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
467.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
468.               
469.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
470.                     VP[1:-1] = af*VP[:-2]+bf*VP[1:-1]+cf*VP[2:]  
471.                  
472.                 else:  
473.                     #otherwise perform the step and then place the dividend  
474.                      #call  
475.                     VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
476.                     VC[0] = VC[0]*(1-rr*dt)  
477.                       
478.                     #put  
479.                     VP[1:-1] = a*VP[:-2]+b*VP[1:-1]+c*VP[2:]  
480.                       
481.                     VC[0] = 0  
482.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
483.                     VP[0] = K  
484.                     VP[NAS] = 0  
485.                       
486.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
487.                       
488.                     tckC = si.splrep(np.exp(Z), VC, k = 1)  
489.                     tckP = si.splrep(np.exp(Z), VP, k = 1)  
490.                     #option prices are interpolated on the grid  
491.                     VC = si.splev(np.maximum(np.exp(Z)-div, np.exp(Zmin)), tckC, der = 0)              
492.                     VP = si.splev(np.maximum(np.exp(Z)-div, np.exp(Zmin)), tckP, der = 0)    
493.                      
494.                       
495.             else:  
496.                 #there is no dividends  
497.                 #call  
498.                 VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
499.                #put  
500.                 VP[1:-1] = a*VP[:-2]+b*VP[1:-1]+c*VP[2:]  
501.                    
502.             VC[0] = 0  
503.             VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
504.             VP[0] = K  
505.             VP[NAS] = 0  
506.             #check early exercise  
507.             VC = optionContract.earlyC(VC, np.exp(Z))  
508.             VP = optionContract.earlyP(VP, np.exp(Z))  
509.               
510.         #get the prices and the grid from the vector      
511.         tck = si.splrep(Z, VC)  
512.         optionContract.call.price = si.splev(np.log(S0), tck, der = 0)   
513.         p1 = si.splev(np.log(S0+1), tck, der = 0)  
514.         p2 = si.splev(np.log(S0-1), tck, der = 0)  
515.         optionContract.call.delta = (p1-p2)/2  
516.         optionContract.call.gamma = p1+p2-2*optionContract.call.price  
517.         VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
518.         tck = si.splrep(Z, VC)  
519.         optionContract.call.theta = (optionContract.call.price-si.splev(np.log(S0), tck, der = 0) )/dt  
520.           
521.           
522.         tck = si.splrep(Z, VP)  
523.         optionContract.put.price = si.splev(np.log(S0), tck, der = 0)   
524.         p1 = si.splev(np.log(S0+1), tck, der = 0)  
525.         p2 = si.splev(np.log(S0-1), tck, der = 0)  
526.         optionContract.put.delta = (p1-p2)/2  
527.         optionContract.put.gamma = p1+p2-2*optionContract.put.price  
528.         VP[1:-1] = a*VP[:-2]+b*VP[1:-1]+c*VP[2:]  
529.         tck = si.splrep(Z, VP)  
530.         optionContract.put.theta = (optionContract.put.price-si.splev(np.log(S0), tck, der = 0) )/dt  
531.           
532.           
533.         return VC                
534.     @staticmethod  
535.     def PDEUSSolo(pricingParameters, optionContract, NTS, NAS, regul, stepAdapt, option, payoff, early):  
536.               
537.         #parameters  
538.         S0 = pricingParameters._S0  
539.         r2 = pricingParameters._r2  
540.         T = pricingParameters._T         
541.         sigma = optionContract._Sigma  
542.         divi = pricingParameters._div  
543.           
544.         #calculate delta T      
545.         deltaT = float(T) / NTS  
546.    
547.         Zmin = -7.5  
548.         Zmax = 7.5  
549.         
550.         dZ = float((Zmax-Zmin)/NAS)     
551.           
552.         #time step is constraint by stability to be a function  of dZ  
553.         #but if the user wants a smaller time step, that is fine too  
554.         dt = min(dZ**2/(3*sigma*sigma), float(T)/NTS)  
555.           
556.         dividends, dividendsStep = PDEPricer.preparedivs(divi, T, deltaT)  
557.         r, rArray, rateChangeStep, DFArray = PDEPricer.preparerates(r2, T, dt, NTS)  
558.        
559.               
560.         #we get the rate  
561.         rr = rArray[-1][0]  
562.         NTS = int(T/dt)+1  
563.           
564.         a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
565.         b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
566.         c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
567.              
568.          
569.         VC = np.asarray( [0 for i in xrange(int(NAS)+1)])  
570.         Z = np.array([(Zmin+i*dZ) for i in xrange(int(NAS)+1)])  
571.           
572.           
573.         if regul:  
574.             VC = payoff(np.exp(Z), rr, sigma, dt)  
575.             NTS = NTS-1  
576.         else:  
577.             VC = payoff(np.exp(Z))  
578.               
579.         a1 = (Z[NAS]-Z[NAS-1])  
580.         b1 = (Z[NAS-1]-Z[NAS-2])  
581.               
582.       
583.         for k in range(NTS, -1, -1):  
584.               
585.             #we need to change the a, b, c only when rate changes  
586.             if (k in rateChangeStep[:]):  
587.                 rr = (r[1][np.nonzero(rateChangeStep == (k))[0]])  
588.                 a = np.asarray( [(1/(1+rr*dt))*(-(rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
589.                 b = np.asarray( [(1/(1+rr*dt))*(1-(sigma**2)*dt/(dZ*dZ)) for i in xrange(int(NAS-1))])  
590.                 c = np.asarray( [(1/(1+rr*dt))*((rr-sigma**2/2)*dt/(2*dZ)+(sigma**2)*dt/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
591.               
592.             if (k in dividendsStep):  
593.                 if stepAdapt:  
594.                       
595.                     #we perform a fraction of step for the pde, place the dividend and finish the step  
596.                     #so no matter the step size, the div is at the right place.  
597.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
598.                     fraction = -(k*dt-dividends[0][np.nonzero(dividendsStep == (k))[0]])  
599.                     dtf = dt-fraction  
600.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
601.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
602.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
603.               
604.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
605.                       
606.                     VC[0] = VC[0] = VC[0]*(1-rr*dt)  
607.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
608.                       
609.                     tckC = si.splrep(Z, VC, k = 1)  
610.                     #option prices are interpolated on the grid  
611.                     VC = si.splev(np.log(np.maximum(np.exp(Z)-div, np.exp(Zmin))), tckC, der = 0)              
612.                       
613.                                    
614.                     VC[0] = VC[0]*(1-rr*dt)  
615.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
616.                       
617.                     #check for exercise  
618.                     VC = early(VC, np.exp(Z))  
619.                       
620.                     #we finish the step  
621.                     dtf = fraction  
622.                      
623.                     af = np.asarray( [(1/(1+rr*dtf))*(-(rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
624.                     bf = np.asarray( [(1/(1+rr*dtf))*(1-(sigma**2)*dtf/(dZ*dZ)) for i in xrange(int(NAS-1))])  
625.                     cf = np.asarray( [(1/(1+rr*dtf))*((rr-sigma**2/2)*dtf/(2*dZ)+(sigma**2)*dtf/(2*dZ*dZ)) for i in xrange(int(NAS-1))])  
626.               
627.                     VC[1:-1] = af*VC[:-2]+bf*VC[1:-1]+cf*VC[2:]  
628.                       
629.                  
630.                 else:  
631.                     #otherwise perform the step and then place the dividend  
632.                      #call  
633.                     VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
634.                     VC[0] = VC[0]*(1-rr*dt)  
635.                     VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
636.                       
637.                     div = dividends[1][np.nonzero(dividendsStep == (k))[0]]   
638.                       
639.                     tckC = si.splrep(np.exp(Z), VC, k = 1)  
640.                     #option prices are interpolated on the grid  
641.                     VC = si.splev(np.maximum(np.exp(Z)-div, np.exp(Zmin)), tckC, der = 0)              
642.                      
643.                       
644.             else:  
645.                 #there is no dividends  
646.                 #call  
647.                 VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
648.                  
649.                    
650.             VC[0] = VC[0]*(1-rr*dt)  
651.             VC[NAS] = (VC[NAS-1]/a1+VC[NAS-1]/b1 -VC[NAS-2]/b1)*a1  
652.               
653.             #check early exercise  
654.             VC = early(VC, np.exp(Z))  
655.               
656.               
657.         #get the prices and the grid from the vector      
658.         tck = si.splrep(Z, VC)  
659.         option.price = si.splev(np.log(S0), tck, der = 0)   
660.         p1 = si.splev(np.log(S0+1), tck, der = 0)  
661.         p2 = si.splev(np.log(S0-1), tck, der = 0)  
662.         option.delta = (p1-p2)/2  
663.         option.gamma = p1+p2-2*option.price  
664.         VC[1:-1] = a*VC[:-2]+b*VC[1:-1]+c*VC[2:]  
665.         tck = si.splrep(Z, VC)  
666.         option.theta = (option.price-si.splev(np.log(S0), tck, der = 0) )/dt  
667.                      
668.           
669.         return VC   
670.     def __init__(self):  
671.         self.SetPricerId("PDE")    
672.         super(PDEPricer, self).__init__()          
673.           
674.       
675.     def SetPricerParameters(self, N, NAS = 1, Regul = False, StepAdapt = False):  
676.         self._NumberAssetSteps = NAS  
677.         self._NumberSteps = N  
678.         self._Regul = Regul  
679.         self._StepAdapt = StepAdapt  
680.   
681.     def PrintSelf(self):  
682.         print self._pricerId  
683.         print "Number of time steps: ", self._NumberSteps  
684.         print "Number of asset steps: ", self._NumberAssetSteps  
685.         print "StepAdapt on: ", self._StepAdapt  
686.         print "Regul on: ", self._Regul  
687.         self._PricingParameters.printSelf()  
688.           
689.         for optionContract in self._contractList:  
690.             optionContract.printSelf()  
691.   
692.     def CalculateContract(self, optionContract):  
693.         # Check input parameters  
694.          
695.         if isinstance(optionContract, OptionContractUS):  
696.               
697.             if isinstance(optionContract, OptionContractDuo):  
698.                 PDEPricer.PDEUSDuo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, self._StepAdapt)  
699.             elif isinstance(optionContract, OptionContractCall):  
700.                 if self._Regul:  
701.                     PDEPricer.PDEUSSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, self._StepAdapt, optionContract.call, optionContract.expiCR, optionContract.earlyC)  
702.                 else:  
703.                     PDEPricer.PDEUSSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, False, optionContract.call, optionContract.expiC, optionContract.earlyC)  
704.             elif isinstance(optionContract, OptionContractPut):  
705.                 if self._Regul:  
706.                     PDEPricer.PDEUSSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, self._StepAdapt, optionContract.put, optionContract.expiPR, optionContract.earlyP)  
707.                 else:  
708.                     PDEPricer.PDEUSSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, False, optionContract.put, optionContract.expiP, optionContract.earlyP)  
709.         else:  
710.               
711.             if isinstance(optionContract, OptionContractDuo):  
712.                 PDEPricer.PDEEUDuo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, self._StepAdapt)  
713.             elif isinstance(optionContract, OptionContractCall):  
714.                 if self._Regul:  
715.                     PDEPricer.PDEEUSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, self._StepAdapt, optionContract.call, optionContract.expiCR)  
716.                 else:  
717.                     PDEPricer.PDEEUSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, False, optionContract.call, optionContract.expiC)  
718.             elif isinstance(optionContract, OptionContractPut):  
719.                 if self._Regul:  
720.                     PDEPricer.PDEEUSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, self._StepAdapt, optionContract.put, optionContract.expiPR)  
721.                 else:  
722.                     PDEPricer.PDEEUSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._Regul, False, optionContract.put, optionContract.expiP)  
723.             else:  
724.                 raise PricerError("PDEPricer", "Contract not yet supported by pricing algorithm")  

Note there there is not specific handling of the cash or nothing options. They use proper payoff function and it works. This is because I have pretty generic boundary conditions. It will not be the case for all possible options.