FourierPricer class

And the last class (to date) is the Fourier pricer. Like the MC class, I only have to computing methods, one for solo and one for duo. For european options without dividends, the pricer only compute in one step.
If we price european options with cash dividends, we need to provide time steps.
The value provided for the number of asset steps is a power of 2 so 10 gives 2^10=1024 asset steps.

1. import scipy.interpolate as si  
2.   
3.   
4. from OptionPricer import *  
5. from scipy import *  
6.   
7. class FourierPricer(OptionPricer):  
8.     _NumberAssetStep = 0  
9.     _StepAdapt = False  
10.     _Depth =1  
11.       
12.        
13.     @staticmethod  
14.     def FourierDuo(pricingParameters, optionContract, nts, N, stepAdapt,american):  
15.         #parameters  
16.         S0 = pricingParameters._S0  
17.         r2 = pricingParameters._r2  
18.         T = pricingParameters._T         
19.         sigma = optionContract._Sigma  
20.         divi = pricingParameters._div  
21.         NTS = nts  
22.         dt = float(T)/NTS  
23.            
24.         dividends,dividendsStep = PricingParameters.preparedivs(divi,T,dt)  
25.         r,rArray,rateChangeStep,DFArray = PricingParameters.preparerates(r2,T,dt,NTS)  
26.               
27.         j = complex(0,1)  
28.         #create vector in the real space  
29.         x_min = -7.5  
30.         x_max = 7.5  
31.         dx = (x_max-x_min)/(N-1)  
32.         x = linspace(x_min,x_max,N)  
33.          
34.         #create vector in the fourier space  
35.         w_max = np.pi/dx;  
36.         dw = 2*w_max/(N);  
37.         w = np.concatenate((linspace(0,w_max,N/2+1),linspace(-w_max+dw,-dw,N/2-1)))  
38.       
39.         # Option payoff  
40.         s = S0*np.exp(x);  
41.         VC = optionContract.expiC(s)  
42.         VP = optionContract.expiP(s)  
43.          
44.         # FST method  
45.         if (american or pricingParameters._pvdividends>0):  
46.           
47.             rr = rArray[-1][0]   
48.             char_exp_factor = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dt)  
49.               
50.             for k in range(NTS,0,-1):   
51.              
52.                if (k in rateChangeStep[:]):  
53.                    rr =  (r[1][np.nonzero(rateChangeStep == (k))[0]])  
54.                    char_exp_factor = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dt)  
55.             
56.                #if there is a dividend, we adjust the path  
57.                if (k in dividendsStep):  
58.                      
59.                    if stepAdapt:  
60.                         #first get the fraction of dt  
61.                         fraction = dividends[0][np.nonzero(dividendsStep == (k))[0]]-k*dt  
62.                         dtf = dt-fraction  
63.                         #perform fraction of step  
64.                         char_exp_factorf = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dtf)  
65.                         VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factorf))  
66.                         VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factorf))  
67.                           
68.                         #handle dividends  
69.                         div = dividends[1][np.nonzero(dividendsStep == (k))[0]]    
70.                           
71.                         #a call is an always increasing function of the strike  
72.                         mid = int(N/2)  
73.                         delta = VC[1:mid]-VC[:mid-1]  
74.                         BoS = np.where(delta<0)[0][-1]  
75.                         VC[:BoS] = 0  
76.                           
77.                         #The put always decreasing  
78.                         delta = VP[1:mid]-VP[:mid-1]  
79.                         BoS = np.where(delta<0)[0][0]  
80.                         VP[:BoS] = VP[BoS]  
81.                 
82.                         tckC = si.splrep(s,VC)  
83.                         tckP = si.splrep(s,VP)  
84.                         #option prices are interpolated on the grid  
85.                         VC = si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der = 0)              
86.                         VP = si.splev(np.maximum(s-div,np.exp(x_min)) ,tckP,der = 0)   
87.                           
88.                         if american:       
89.                               VC = optionContract.earlyC(VC,s)  
90.                               VP = optionContract.earlyP(VP,s)  
91.                           
92.                         #finish the step  
93.                         dtf = fraction   
94.                         char_exp_factorf = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dtf)  
95.                         VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factorf))  
96.                         VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factorf))  
97.                    else:    
98.                         VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
99.                         VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))  
100.                       
101.                         div = dividends[1][np.nonzero(dividendsStep==(k))[0]]    
102.                           
103.                         #a call is an always increasing function of the strike  
104.                         mid = int(N/2)  
105.                         delta = VC[1:mid]-VC[:mid-1]  
106.                         BoS = np.where(delta<0)[0][-1]  
107.                         VC[:BoS] = 0  
108.                           
109.                         #The put always decreasing  
110.                         delta = VP[1:mid]-VP[:mid-1]  
111.                         BoS = np.where(delta<0)[0][0]  
112.                         VP[:BoS] = VP[BoS]  
113.                          
114.                         tckC = si.splrep(s,VC)  
115.                         tckP = si.splrep(s,VP)  
116.                         #option prices are interpolated on the grid  
117.                         VC = si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der = 0)              
118.                         VP = si.splev(np.maximum(s-div,np.exp(x_min)) ,tckP,der = 0)   
119.                         if american:       
120.                               VC = optionContract.earlyC(VC,s)  
121.                               VP = optionContract.earlyP(VP,s)  
122.                else:  
123.                    #no dividend at this step  
124.                    VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
125.                    VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))      
126.                           
127.                    if american:       
128.                         VC = optionContract.earlyC(VC,s)  
129.                         VP = optionContract.earlyP(VP,s)  
130.         else:  
131.           char_exp_factor = np.exp((j*(pricingParameters._r-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-pricingParameters._r)*T)   
132.           VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
133.           VP = np.real(np.fft.ifft(np.fft.fft(VP)*char_exp_factor))      
134.                  
135.          
136.         #Interpolate option prices  
137.         tck = si.splrep(s,VC)  
138.         optionContract.call.price = si.splev(S0,tck,der = 0)     
139.         tck = si.splrep(s,VP)  
140.         optionContract.put.price = si.splev(S0,tck,der = 0)       
141.   
142.   
143.     @staticmethod  
144.     def FourierSolo(pricingParameters, optionContract, nts, N, stepAdapt,american,option,payoff,early):  
145.         #parameters  
146.         S0 = pricingParameters._S0  
147.         r2 = pricingParameters._r2  
148.         T = pricingParameters._T         
149.         sigma = optionContract._Sigma  
150.         divi = pricingParameters._div  
151.         NTS = nts  
152.         dt = float(T)/NTS  
153.   
154.         dividends,dividendsStep = PricingParameters.preparedivs(divi,T,dt)  
155.         r,rArray,rateChangeStep,DFArray = PricingParameters.preparerates(r2,T,dt,NTS)  
156.               
157.         j = complex(0,1)  
158.         #create vector in the real space  
159.         x_min = -7.5  
160.         x_max = 7.5  
161.         dx = (x_max-x_min)/(N-1)  
162.         x = linspace(x_min,x_max,N)  
163.          
164.         #create vector in the fourier space  
165.         w_max = np.pi/dx;  
166.         dw = 2*w_max/(N);  
167.         w = np.concatenate((linspace(0,w_max,N/2+1),linspace(-w_max+dw,-dw,N/2-1)))  
168.       
169.         # Option payoff  
170.         s = S0*np.exp(x);  
171.         VC = payoff(s)  
172.           
173.         # FST method  
174.         if (american or pricingParameters._pvdividends>0):  
175.             rr = rArray[-1][0]   
176.             char_exp_factor = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dt)  
177.             for k in range(NTS,0,-1):   
178.              
179.                if (k in rateChangeStep[:]):  
180.                    rr =  (r[1][np.nonzero(rateChangeStep == (k))[0]])  
181.                    char_exp_factor = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dt)  
182.             
183.                #if there is a dividend, we adjust the path  
184.                if (k in dividendsStep):  
185.                    if stepAdapt:  
186.                         #first get the fraction of dt  
187.                         fraction = dividends[0][np.nonzero(dividendsStep == (k))[0]]-k*dt  
188.                         dtf = dt-fraction  
189.                         #perform fraction of step  
190.                         char_exp_factorf = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dtf)  
191.                         VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factorf))  
192.                           
193.                         #handle dividends  
194.                         div = dividends[1][np.nonzero(dividendsStep == (k))[0]]    
195.                           
196.                         if isinstance(optionContract,OptionContractCall):  
197.                             #a call is an always increasing function of the strike  
198.                             mid = int(N/2)  
199.                             delta = VC[1:mid]-VC[:mid-1]  
200.                             BoS = np.where(delta<0)[0][-1]  
201.                             VC[:BoS] = 0  
202.                           
203.                         if isinstance(optionContract,OptionContractPut):  
204.                             #The put always decreasing  
205.                             mid = int(N/2)  
206.                             delta = VC[1:mid]-VC[:mid-1]  
207.                             BoS = np.where(delta<0)[0][0]  
208.                             VC[:BoS] = VC[BoS]  
209.                 
210.                         tckC = si.splrep(s,VC)  
211.                         #option prices are interpolated on the grid  
212.                         VC = si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der = 0)              
213.                           
214.                         if american:       
215.                               VC = earlyC(VC,s)  
216.                                 
217.                         #finish the step  
218.                         dtf = fraction   
219.                         char_exp_factorf = np.exp((j*(rr-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-rr)*dtf)  
220.                         VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factorf))  
221.                    else:    
222.                         VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
223.                           
224.                         div = dividends[1][np.nonzero(dividendsStep == (k))[0]]    
225.                           
226.                         if isinstance(optionContract,OptionContractCall):  
227.                             #a call is an always increasing function of the strike  
228.                             mid = int(N/2)  
229.                             delta = VC[1:mid]-VC[:mid-1]  
230.                             BoS = np.where(delta<0)[0][-1]  
231.                             VC[:BoS] = 0  
232.                           
233.                         if isinstance(optionContract,OptionContractPut):  
234.                             #The put always decreasing  
235.                             mid = int(N/2)  
236.                             delta = VC[1:mid]-VC[:mid-1]  
237.                             BoS = np.where(delta<0)[0][0]  
238.                             VC[:BoS] = VC[BoS]  
239.                          
240.                         tckC = si.splrep(s,VC)  
241.                         #option prices are interpolated on the grid  
242.                         VC = si.splev(np.maximum(s-div,np.exp(x_min)) ,tckC,der = 0)              
243.                           
244.                         if american:       
245.                               VC = early(VC,s)  
246.                                 
247.                else:  
248.                    #no dividend at this step  
249.                    VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
250.                           
251.                    if american:       
252.                         VC = early(VC,s)  
253.                       
254.         else:  
255.           char_exp_factor = np.exp((j*(pricingParameters._r-0.5*sigma**2)*w - 0.5*sigma**2*(w**2)-pricingParameters._r)*T)   
256.           VC = np.real(np.fft.ifft(np.fft.fft(VC)*char_exp_factor))  
257.          
258.         #Interpolate option prices  
259.         tck = si.splrep(s,VC)  
260.         option.price = si.splev(S0,tck,der = 0)     
261.       
262.     def __init__(self):  
263.         self.SetPricerId("FourierPricer")    
264.         super(FourierPricer, self).__init__()          
265.       
266.     def SetPricerParameters(self, NTS, depth = 1, StepAdapt = False):  
267.         self._StepAdapt = StepAdapt  
268.         self._Depth = depth  
269.         self._NumberAssetSteps = 2**depth  
270.         self._NumberSteps=NTS  
271.     
272.     def PrintSelf(self):  
273.         print self._pricerId  
274.         print "Number of time steps: ", self._NumberSteps  
275.         print "Number of asset steps: ", self._NumberAssetSteps  
276.         print "StepAdapt on: ", self._StepAdapt  
277.         self._PricingParameters.printSelf()  
278.           
279.         for optionContract in self._contractList:  
280.             optionContract.printSelf()  
281.   
282.     def CalculateContract(self, optionContract):  
283.         # Check input parameters  
284.         if isinstance(optionContract, OptionContractUS):  
285.             if isinstance(optionContract, OptionContractDuo):  
286.                 FourierPricer.FourierDuo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._StepAdapt, True)  
287.             elif isinstance(optionContract, OptionContractCall):  
288.                 FourierPricer.FourierSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._StepAdapt, True, optionContract.call, optionContract.expiC, optionContract.earlyC)  
289.             elif isinstance(optionContract, OptionContractPut):  
290.                 FourierPricer.FourierSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._StepAdapt, True, optionContract.put, optionContract.expiP, optionContract.earlyP)  
291.         else:  
292.             if isinstance(optionContract, OptionContractDuo):  
293.                 FourierPricer.FourierDuo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._StepAdapt, False)  
294.             elif isinstance(optionContract, OptionContractCall):  
295.                 FourierPricer.FourierSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._StepAdapt, False, optionContract.call, optionContract.expiC,"")  
296.             elif isinstance(optionContract, OptionContractPut):  
297.                 FourierPricer.FourierSolo(self._PricingParameters, optionContract, self._NumberSteps, self._NumberAssetSteps, self._StepAdapt, False, optionContract.put, optionContract.expiP,"")  
298.               
299.             else:  
300.                 raise PricerError("FourierPricer", "Contract not yet supported by pricing algorithm")