MCPricer Class

Now it is the turn of the MC pricer to get a face lift. Note that I only have two versions of the computing code. One for solo, one for duo. The code is smart enough to directly jump to expiration if the option is european without dividend. If not the brownian bridge is used to generate the paths.

In case the option is american, LSM is used to check for early exercise. I included the class Sobol in the file but one still need the ltqnorm.py file that I provided earlier.
When one initialize the MC pricer, a value, depth, is used to generate the bridge so depth=10 means 2^10=1024 time steps.

1. import scipy.interpolate as si  
2.   
3.   
4. from OptionPricer import *  
5. from ltqnorm import *  
6.   
7. class SobolSequence:  
8.     "From the Gnu Scientific Library"  
9.     def __init__(self,dim):  
10.         self.max_dim = 40  
11.         self.bit_count = 30  
12.         self.primitive_polynomials = [  
13.             1,     3,   7,  11,  13,  19,  25,  37,  59,  47,  
14.             61,   55,  41,  67,  97,  91, 109, 103, 115, 131,  
15.             193, 137, 145, 143, 241, 157, 185, 167, 229, 171,  
16.             213, 191, 253, 203, 211, 239, 247, 285, 369, 299  
17.             ]  
18.         self.degree_table = [  
19.             0, 1, 2, 3, 3, 4, 4, 5, 5, 5,  
20.             5, 5, 5, 6, 6, 6, 6, 6, 6, 7,  
21.             7, 7, 7, 7, 7, 7, 7, 7, 7, 7,   
22.             7, 7, 7, 7, 7, 7, 7, 8, 8, 8  
23.             ]  
24.         self.v_init = [  
25.             [0, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
26.              1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
27.              1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
28.              1, 1, 1, 1, 1, 1, 1, 1, 1, 1],  
29.             [0, 0, 1, 3, 1, 3, 1, 3, 3, 1,  
30.              3, 1, 3, 1, 3, 1, 1, 3, 1, 3,  
31.              1, 3, 1, 3, 3, 1, 3, 1, 3, 1,  
32.              3, 1, 1, 3, 1, 3, 1, 3, 1, 3],  
33.             [0, 0, 0, 7, 5, 1, 3, 3, 7, 5,  
34.              5, 7, 7, 1, 3, 3, 7, 5, 1, 1,  
35.              5, 3, 3, 1, 7, 5, 1, 3, 3, 7,  
36.              5, 1, 1, 5, 7, 7, 5, 1, 3, 3],  
37.             [0,  0,  0,  0,  0,  1,  7,  9, 13, 11,  
38.              1,  3,  7,  9,  5, 13, 13, 11,  3, 15,  
39.              5,  3, 15,  7,  9, 13,  9,  1, 11,  7,  
40.              5, 15,  1, 15, 11,  5,  3,  1,  7,  9],  
41.             [0,  0,  0,  0,  0,  0,  0,  9,  3, 27,  
42.              15, 29, 21, 23, 19, 11, 25,  7, 13, 17,  
43.              1, 25, 29,  3, 31, 11,  5, 23, 27, 19,  
44.              21,  5,  1, 17, 13,  7, 15,  9, 31,  9],  
45.             [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  
46.              0,  0,  0, 37, 33,  7,  5, 11, 39, 63,  
47.              27, 17, 15, 23, 29,  3, 21, 13, 31, 25,  
48.              9, 49, 33, 19, 29, 11, 19, 27, 15, 25],  
49.             [0,   0,  0,  0,  0,  0,    0,  0,  0,   0,  
50.              0,   0,  0,  0,  0,  0,    0,  0,  0,  13,  
51.              33, 115, 41, 79, 17,  29, 119, 75, 73, 105,  
52.              7,  59, 65, 21,  3, 113,  61, 89, 45, 107],  
53.             [0, 0, 0, 0, 0, 0, 0, 0,  0,  0,  
54.              0, 0, 0, 0, 0, 0, 0, 0,  0,  0,  
55.              0, 0, 0, 0, 0, 0, 0, 0,  0,  0,  
56.              0, 0, 0, 0, 0, 0, 0, 7, 23, 39]  
57.             ] # end of self.v_init  
58.         self.dim = dim  
59.         #self.v_direction = [[0]*self.bit_count]*self.max_dim  
60.         #self.v_direction = [[0]*self.max_dim]*self.bit_count #??  
61.         self.v_direction = []  
62.         for i in range(self.bit_count):  
63.             array = []  
64.             for j in range(self.max_dim):  
65.                 array.append(0)  
66.             self.v_direction.append(array)  
67.           
68.   
69.         # Initialize direction table in dimension 0.  
70.         for k in range(self.bit_count): self.v_direction[k][0] = 1  
71.   
72.         # Initialize in remaining dimensions.  
73.         for i_dim in range(1,self.dim):  
74.             poly_index = i_dim  
75.             degree_i = self.degree_table[poly_index]  
76.             includ = [0]*8  
77.               
78.             # Expand the polynomial bit pattern to separate  
79.             # components of the logical array includ[].  
80.             p_i = self.primitive_polynomials[poly_index]  
81.             for k in range(degree_i-1,-1,-1):  
82.                 includ[k] = ((p_i % 2) == 1)  
83.                 p_i /= 2  
84.   
85.             # Leading elements for dimension i come from v_init[][].  
86.             for j in range(degree_i):  
87.                 self.v_direction[j][i_dim] = self.v_init[j][i_dim]  
88.   
89.             # Calculate remaining elements for this dimension,  
90.             # as explained in Bratley+Fox, section 2.  
91.             for j in range(degree_i, self.bit_count):  
92.                 newv = self.v_direction[j-degree_i][i_dim]  
93.                 ell = 1  
94.                 for k in range(degree_i):  
95.                     ell *= 2  
96.                     if includ[k]:  
97.                         newv = newv ^ (ell*self.v_direction[j-k-1][i_dim])  
98.                 self.v_direction[j][i_dim] = newv  
99.         # done with for i_dim in range(1,self.dim)  
100.   
101.         # Multiply columns of v by appropriate power of 2  
102.         ell = 1  
103.         for j in range(self.bit_count-2,-1,-1):  
104.             ell *= 2  
105.             for i_dim in range(self.dim):  
106.                 self.v_direction[j][i_dim] = self.v_direction[j][i_dim]*ell  
107.         # 1/(common denom of the elements in v_direction)  
108.   
109.         self.last_denominator_inv = 1.0 /(2.0 * ell)  
110.   
111.         # final setup  
112.         self.sequence_count = 0  
113.         self.last_numerator_vec = [0]*self.dim  
114.         return  
115.   
116.     def __call__(self):  
117.         ell = 0  
118.         c = self.sequence_count  
119.         while 1:  
120.             ell += 1  
121.             if c % 2 == 1:  
122.                 c/=2  
123.             else:  
124.                 break  
125.         if ell > self.bit_count:  
126.             raise "Sobol failed for %d" % self.sequence_count  
127.   
128.         v = [0]*self.dim  
129.         for i_dimension in range(self.dim):  
130.             direction_i = self.v_direction[ell-1][i_dimension]  
131.             old_numerator_i = self.last_numerator_vec[i_dimension]  
132.             new_numerator_i = old_numerator_i ^ direction_i  
133.             self.last_numerator_vec[i_dimension] = new_numerator_i  
134.             v[i_dimension] = new_numerator_i*self.last_denominator_inv  
135.         self.sequence_count += 1  
136.         return v  
137.          
138.   
139.   
140.   
141.   
142. class MCPricer(OptionPricer):  
143.     _NumberPath = 0  
144.     _StepAdapt = False  
145.     _Regul = False  
146.     _AP = False  
147.     _MM = True  
148.     _Depth = 1  
149.       
150.      
151.     @staticmethod  
152.     def RNG(I, ap, mm):  
153.             if ap:  
154.                 ran = np.random.standard_normal (I/2)  
155.                 ran = np.concatenate (( ran ,- ran ))  
156.             else :  
157.                 ran = np.random.standard_normal (I)  
158.             if mm:  
159.                 ran = ran - np.mean ( ran)  
160.                 ran = ran / np.std ( ran )  
161.             return ran  
162.   
163.     @staticmethod  
164.     def SobolRNG(I, ap, mm):  
165.              #we get a sobol sequence of proper size  
166.              #if we use it only for values at expiration, dimension1 is fine  
167.              ss = SobolSequence(1)       
168.              ss()  
169.              if ap:  
170.                 ran = np.array([ltqnorm(ss()[0]) for i in xrange(I/2)])  
171.                 ran = np.concatenate (( ran ,- ran ))  
172.              else :  
173.                 ran = np.array([ltqnorm(ss()[0]) for i in xrange(I)])  
174.              if mm:  
175.                 ran = ran - np.mean ( ran)  
176.                 ran = ran / np.std ( ran )  
177.              return ran  
178.                
179.                
180.     @staticmethod  
181.     def Bridge(x, depth, depthmax, I, sigma, ap, mm):  
182.     #we create path with the brownian bridge  
183.             if depth == depthmax:  
184.                 return x  
185.             xmid = np.average(x,0)+0.5*MCPricer.RNG(I, ap, mm)*np.sqrt(sigma/(2**(depth-1)))  
186.             y = MCPricer.Bridge([x[0][:], xmid[:]], depth+1, depthmax, I, sigma, ap, mm)  
187.             z = MCPricer.Bridge([xmid[:], x[-1][:]], depth+1, depthmax, I, sigma, ap, mm)  
188.             x = np.concatenate((y[:-1][:], z[:][:]))  
189.             return x  
190.          
191.     @staticmethod  
192.     def LSM(hC,VC,S,NP,df,t):   
193.         deg = 2  
194.         #get the ITM call  
195.         dummy = np.greater(hC[t,:], 0)  
196.         relevant = np.nonzero(dummy)  
197.         relevantS = np.compress(dummy, S[t,:])  
198.         p = len(relevantS)  
199.               
200.         #if any  
201.         if p==0:  
202.               contvalues = np.zeros((NP))  
203.         else:  
204.               relevantV = np.compress(dummy, VC[t+1,:]*df)  
205.               #rg is an approximation of the call price as a function of the stock                  
206.               #only in the ITM region  
207.               #to fit that function we discount call prices from thenext time step (we go backward in time)  
208.               rg = np.polyfit(relevantS[:], relevantV[:], deg)                
209.               contvalues = np.polyval(rg, relevantS[:])  
210.                       
211.         #for each call we check if exercise now is higher than approximate value of keeping it      
212.         erg = np.zeros((NP))  
213.         np.put(erg,relevant,contvalues)  
214.         return erg  
215.   
216.   
217.     @staticmethod  
218.     def MCDuo(pricingParameters, optionContract, depth, NP, ap, mm, stepAdapt,american):  
219.               
220.         #parameters  
221.         S0 = pricingParameters._S0  
222.         r2 = pricingParameters._r2  
223.         T = pricingParameters._T         
224.         sigma = optionContract._Sigma  
225.         divi = pricingParameters._div  
226.         NTS=2**depth  
227.         dt=float(T)/NTS  
228.         np.random.seed  
229.       
230.         dividends,dividendsStep=PricingParameters.preparedivs(divi,T,dt)  
231.         r,rArray,rateChangeStep,DFArray=PricingParameters.preparerates(r2,T,dt,NTS)  
232.       
233.         #we directly compute the final value of the brownian motion  
234.         ran=np.asarray(MCPricer.SobolRNG(NP, ap, mm))*sigma*np.sqrt(T)  
235.           
236.         #if we get the proper rate, we can directly get the european value  
237.         rate=pricingParameters._r  
238.           
239.         # generate stock price paths  
240.         SEuro= np.zeros ((NP),'d') # index level matrix  
241.         rr=r[1][1]  
242.         SEuro[:]=(S0)*np.exp((rate-0.5*(sigma**2))* T+ ran[:])  
243.           
244.         if (american or pricingParameters._pvdividends>0):  
245.             # generate stock price paths  
246.             S = np.zeros ((NTS,NP),'d') # index level matrix  
247.             S[0 ,:] = S0 # initial values  
248.       
249.             #we then proceed to generate intermediate values  
250.             ranBB = (np.concatenate(([ran-ran],[ran])))  
251.             ranBB = MCPricer.Bridge(ranBB,1,depth+1,ranBB.shape[1],sigma*sigma*T,ap,mm)  
252.             incBB = ranBB[1:][:]-ranBB[:-1][:]  
253.               
254.             for t in range (1,NTS): # loop over the days  
255.                 if (t in rateChangeStep[:]):  
256.                     rr = (r[1][np.nonzero(rateChangeStep==(t))[0]+1])  
257.       
258.                 S[t ,:] = S[t-1 ,:]* np.exp ((rr-0.5*sigma **2)* dt+incBB[t][:])  
259.                   
260.                 if (t in dividendsStep):  
261.                     div = dividends[1][np.nonzero(dividendsStep==(t))[0]]    
262.                     S[t,:] = np.maximum(S[t,:]-div,0)  
263.               
264.             if american:      
265.                 #apply LMS method  
266.                 VC = optionContract.earlyC(S,0.0)  
267.                 VP = optionContract.earlyP(S,0.0)  
268.                 hC = VC  
269.                 hP = VP  
270.                 for t in range(NTS-2,-1,-1):  
271.                     df = DFArray[t]  
272.                     erg = MCPricer.LSM(hC,VC,S,NP,df,t)  
273.                     VC[t,:] = np.where(hC[t,:]>erg,hC[t,:],VC[t+1,:]*df)  
274.                     erg = MCPricer.LSM(hP,VP,S,NP,df,t)  
275.                     VP[t,:] = np.where(hP[t,:]>erg,hP[t,:],VP[t+1,:]*df)  
276.                    
277.                 optionContract.call.price = sum(VC[0,:])/NP  
278.                 optionContract.put.price = sum(VP[0,:])/NP  
279.             else:  
280.                 df = np.exp(-rate*T)  
281.                 VC = optionContract.expiC(S[-1,:])  
282.                 optionContract.call.price = df*sum(VC)/NP  
283.                 VP = optionContract.expiP(S[-1,:])  
284.                 optionContract.put.price = df*sum(VP)/NP  
285.                       
286.         else:     
287.             #print 'euro pricing'  
288.             df = np.exp(-rate*T)  
289.             VC = optionContract.expiC(SEuro[:])  
290.             optionContract.call.price = df*sum(VC)/NP  
291.             VP = optionContract.expiP(SEuro[:])  
292.             optionContract.put.price = df*sum(VP)/NP  
293.           
294.   
295.     @staticmethod  
296.     def MCSolo(pricingParameters, optionContract, depth, NP, ap, mm, stepAdapt,american,option,payoff,early):  
297.               
298.         #parameters  
299.         S0 = pricingParameters._S0  
300.         r2 = pricingParameters._r2  
301.         T = pricingParameters._T         
302.         sigma = optionContract._Sigma  
303.         divi = pricingParameters._div  
304.         NTS=2**depth  
305.           
306.         dt=float(T)/NTS  
307.         np.random.seed  
308.       
309.         dividends,dividendsStep=PricingParameters.preparedivs(divi,T,dt)  
310.         r,rArray,rateChangeStep,DFArray=PricingParameters.preparerates(r2,T,dt,NTS)  
311.       
312.         #we directly compute the final value of the brownian motion  
313.         ran=np.asarray(MCPricer.SobolRNG(NP, ap, mm))*sigma*np.sqrt(T)  
314.           
315.         #if we get the proper rate, we can directly get the european value  
316.         rate=pricingParameters._r  
317.           
318.         # generate stock price paths  
319.         SEuro= np.zeros ((NP),'d') # index level matrix  
320.         rr=r[1][1]  
321.         SEuro[:]=(S0)*np.exp((rate-0.5*(sigma**2))* T+ ran[:])  
322.           
323.         if (american or pricingParameters._pvdividends>0):  
324.             # generate stock price paths  
325.             S = np.zeros ((NTS,NP),'d') # index level matrix  
326.             S[0 ,:] = S0 # initial values  
327.       
328.             #we then proceed to generate intermediate values  
329.             ranBB = (np.concatenate(([ran-ran],[ran])))  
330.             ranBB = MCPricer.Bridge(ranBB,1,depth+1,ranBB.shape[1],sigma*sigma*T,ap,mm)  
331.             incBB = ranBB[1:][:]-ranBB[:-1][:]  
332.               
333.             for t in range (1,NTS): # loop over the days  
334.                 if (t in rateChangeStep[:]):  
335.                     rr = (r[1][np.nonzero(rateChangeStep==(t))[0]+1])  
336.       
337.                 S[t ,:] = S[t-1 ,:]* np.exp ((rr-0.5*sigma **2)* dt+incBB[t][:])  
338.       
339.                 if (t in dividendsStep):  
340.                     div = dividends[1][np.nonzero(dividendsStep==(t))[0]]    
341.                     S[t,:] = np.maximum(S[t,:]-div,0)  
342.               
343.             if american:      
344.                 #apply LMS method  
345.                 VC = early(0.0)  
346.                 hC = VC  
347.                 for t in range(NTS-2,-1,-1):  
348.                     df = DFArray[t]  
349.                     erg = MCPricer.LSM(hC,VC,S,NP,df,t)  
350.                     VC[t,:] = np.where(hC[t,:]>erg,hC[t,:],VC[t+1,:]*df)  
351.                 option.price = sum(VC[0,:])/NP  
352.             else:  
353.                 df = np.exp(-rate*T)  
354.                 VC = payoff(S[-1,:])  
355.                 option.price = df*sum(VC)/NP  
356.         else:     
357.             #print 'euro pricing'  
358.             df = np.exp(-rate*T)  
359.             VC = payoff(SEuro[:])  
360.             option.price = df*sum(VC)/NP  
361.               
362.           
363.   
364.       
365.     def __init__(self):  
366.         self.SetPricerId("MC")    
367.         super(MCPricer, self).__init__()          
368.           
369.       
370.     def SetPricerParameters(self, NP, depth = 1, ap = False, mm = False, StepAdapt = False):  
371.         self._StepAdapt = StepAdapt  
372.         self._AP=ap  
373.         self._MM=mm  
374.         self._NumberPath=NP  
375.         self._Depth=depth  
376.         self._NumberSteps=2**depth  
377.           
378.   
379.     def PrintSelf(self):  
380.         print self._pricerId  
381.         print "Number of time steps: ", self._NumberSteps  
382.         print "Number of path: ", self._NumberPath  
383.         print "StepAdapt on: ", self._StepAdapt  
384.         print "MM on: ", self._MM  
385.         print "AP on: ", self._AP  
386.         self._PricingParameters.printSelf()  
387.           
388.         for optionContract in self._contractList:  
389.             optionContract.printSelf()  
390.   
391.     def CalculateContract(self, optionContract):  
392.         # Check input parameters  
393.          
394.         if isinstance(optionContract, OptionContractUS):  
395.             if isinstance(optionContract, OptionContractDuo):  
396.                 MCPricer.MCDuo(self._PricingParameters, optionContract, self._Depth, self._NumberPath, self._AP, self._MM, self._StepAdapt,True)  
397.             elif isinstance(optionContract, OptionContractCall):  
398.                 MCPricer.MCSolo(self._PricingParameters, optionContract, self._Depth, self._NumberPath, self._AP, self._MM, self._StepAdapt,True, optionContract.call, optionContract.expiC, optionContract.earlyC)  
399.             elif isinstance(optionContract, OptionContractPut):  
400.                 MCPricer.MCSolo(self._PricingParameters, optionContract, self._Depth, self._NumberPath, self._AP, self._MM, self._StepAdapt,True, optionContract.put, optionContract.expiP, optionContract.earlyP)  
401.         else:  
402.             if isinstance(optionContract, OptionContractDuo):  
403.                 MCPricer.MCDuo(self._PricingParameters, optionContract, self._Depth, self._NumberPath, self._AP, self._MM, self._StepAdapt,False)  
404.             elif isinstance(optionContract, OptionContractCall):  
405.                 MCPricer.MCSolo(self._PricingParameters, optionContract, self._Depth, self._NumberPath, self._AP, self._MM, self._StepAdapt,False, optionContract.call, optionContract.expiC,"")  
406.             elif isinstance(optionContract, OptionContractPut):  
407.                 MCPricer.MCSolo(self._PricingParameters, optionContract, self._Depth, self._NumberPath, self._AP, self._MM, self._StepAdapt,False, optionContract.put, optionContract.expiP,"")  
408.               
409.             else:  
410.                 raise PricerError("MCPricer", "Contract not yet supported by pricing algorithm")