Convergence of MC

In the previous post we looked at the convergence of the binomial tree for american option. We noted that careful management of the times steps can really improve convergence.

Now we look at the MC pricer. Remember, it is based on the idea of approximating the distribution of possible stocks values at expiration by generating random variables according to the same law. As the number of random variable increases, it will converge toward the true distribution. We will then be able to price derivatives.

Let's take a very simple case: Strike and stock value of today set to 100, no dividends, rate 0.
What we should see is the same price for the call and the put, because of call the put parity. However our pricer computes one using paths below the strike and one using paths higher.

Since we are in the simple case of constant volatility, we know one step paths generated using the equation:
\[
S_t=S_0e^{(r-\sigma^2/2)t+\epsilon\sigma\sqrt{t}}
\]
are exact ($\epsilon$ is our gaussian random variable here), there is no approximation error. So we can just focus on the number of paths.
For a few paths with a pseudo random generator we get

The first graph is with nothing and the second using moment matching and antithetic paths. For reference in that case (time 3 years and vol 40%, a call is worth 27.1 euro so a 5 euro error as we can see on the first graph or even a 50 cents on the second graph is unacceptable.

If we use the Sobol sequence we get better results.

The first graph is plain Sobol and the second using moment matching (MM) and antithetic paths (AP). I tried using different technique to transform my Sobol numbers, uniformly distributed in [0,1]. It doesn't seem to make a difference in prices, I'll come back to that at the end of the post.

Sobol is doing better than pseudo random. I am surprised that MM and AP improve Sobol. Well to be honest, only MM improves, AP doesn't do much. But still one should consider it because it uses a Sobol sequence of half the size and get the other half by symmetry, it makes it faster than generating and turning into gaussian a full size Sobol sequence. We note a bias introduced by AP and MM but clearly a lower variance.

If we increase the number of paths we get these results

We note the same general trend, using AP and MM gives a bias but less error. I continue increasing the number of paths and I get:

Just to remind the reader the usefulness of Sobol, here is the same with pseudo random instead (AP and MM on):

So the point is: Sobol converges faster than random but even with Sobol moment matching and antithetic paths can help.
Antithetic paths can be a speed up as we use only half the length of the Sobol sequence, plus it makes the draw symmetric so all odd moments are 0 (mean, skewness,...). Using moment matching also ensure that the distribution has the proper variance. I don't know how to normalize a sequence to have a kurtosis of 3 (like the gaussian has) but I suspect this might help further.

While it did not change pricing results, I found out that using the following Sobol generator (that I found online):

1. class SobolSequence:  
2.     "From the Gnu Scientific Library"  
3.     def __init__(self,dim):  
4.         self.max_dim = 40  
5.         self.bit_count = 30  
6.         self.primitive_polynomials = [  
7.             1,     3,   7,  11,  13,  19,  25,  37,  59,  47,  
8.             61,   55,  41,  67,  97,  91, 109, 103, 115, 131,  
9.             193, 137, 145, 143, 241, 157, 185, 167, 229, 171,  
10.             213, 191, 253, 203, 211, 239, 247, 285, 369, 299  
11.             ]  
12.         self.degree_table = [  
13.             0, 1, 2, 3, 3, 4, 4, 5, 5, 5,  
14.             5, 5, 5, 6, 6, 6, 6, 6, 6, 7,  
15.             7, 7, 7, 7, 7, 7, 7, 7, 7, 7,   
16.             7, 7, 7, 7, 7, 7, 7, 8, 8, 8  
17.             ]  
18.         self.v_init = [  
19.             [0, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
20.              1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
21.              1, 1, 1, 1, 1, 1, 1, 1, 1, 1,  
22.              1, 1, 1, 1, 1, 1, 1, 1, 1, 1],  
23.             [0, 0, 1, 3, 1, 3, 1, 3, 3, 1,  
24.              3, 1, 3, 1, 3, 1, 1, 3, 1, 3,  
25.              1, 3, 1, 3, 3, 1, 3, 1, 3, 1,  
26.              3, 1, 1, 3, 1, 3, 1, 3, 1, 3],  
27.             [0, 0, 0, 7, 5, 1, 3, 3, 7, 5,  
28.              5, 7, 7, 1, 3, 3, 7, 5, 1, 1,  
29.              5, 3, 3, 1, 7, 5, 1, 3, 3, 7,  
30.              5, 1, 1, 5, 7, 7, 5, 1, 3, 3],  
31.             [0,  0,  0,  0,  0,  1,  7,  9, 13, 11,  
32.              1,  3,  7,  9,  5, 13, 13, 11,  3, 15,  
33.              5,  3, 15,  7,  9, 13,  9,  1, 11,  7,  
34.              5, 15,  1, 15, 11,  5,  3,  1,  7,  9],  
35.             [0,  0,  0,  0,  0,  0,  0,  9,  3, 27,  
36.              15, 29, 21, 23, 19, 11, 25,  7, 13, 17,  
37.              1, 25, 29,  3, 31, 11,  5, 23, 27, 19,  
38.              21,  5,  1, 17, 13,  7, 15,  9, 31,  9],  
39.             [0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  
40.              0,  0,  0, 37, 33,  7,  5, 11, 39, 63,  
41.              27, 17, 15, 23, 29,  3, 21, 13, 31, 25,  
42.              9, 49, 33, 19, 29, 11, 19, 27, 15, 25],  
43.             [0,   0,  0,  0,  0,  0,    0,  0,  0,   0,  
44.              0,   0,  0,  0,  0,  0,    0,  0,  0,  13,  
45.              33, 115, 41, 79, 17,  29, 119, 75, 73, 105,  
46.              7,  59, 65, 21,  3, 113,  61, 89, 45, 107],  
47.             [0, 0, 0, 0, 0, 0, 0, 0,  0,  0,  
48.              0, 0, 0, 0, 0, 0, 0, 0,  0,  0,  
49.              0, 0, 0, 0, 0, 0, 0, 0,  0,  0,  
50.              0, 0, 0, 0, 0, 0, 0, 7, 23, 39]  
51.             ] # end of self.v_init  
52.         self.dim = dim  
53.         #self.v_direction = [[0]*self.bit_count]*self.max_dim  
54.         #self.v_direction = [[0]*self.max_dim]*self.bit_count #??  
55.         self.v_direction = []  
56.         for i in range(self.bit_count):  
57.             array = []  
58.             for j in range(self.max_dim):  
59.                 array.append(0)  
60.             self.v_direction.append(array)  
61.           
62.   
63.         # Initialize direction table in dimension 0.  
64.         for k in range(self.bit_count): self.v_direction[k][0] = 1  
65.   
66.         # Initialize in remaining dimensions.  
67.         for i_dim in range(1,self.dim):  
68.             poly_index = i_dim  
69.             degree_i = self.degree_table[poly_index]  
70.             includ = [0]*8  
71.               
72.             # Expand the polynomial bit pattern to separate  
73.             # components of the logical array includ[].  
74.             p_i = self.primitive_polynomials[poly_index]  
75.             for k in range(degree_i-1,-1,-1):  
76.                 includ[k] = ((p_i % 2) == 1)  
77.                 p_i /= 2  
78.   
79.             # Leading elements for dimension i come from v_init[][].  
80.             for j in range(degree_i):  
81.                 self.v_direction[j][i_dim] = self.v_init[j][i_dim]  
82.   
83.             # Calculate remaining elements for this dimension,  
84.             # as explained in Bratley+Fox, section 2.  
85.             for j in range(degree_i, self.bit_count):  
86.                 newv = self.v_direction[j-degree_i][i_dim]  
87.                 ell = 1  
88.                 for k in range(degree_i):  
89.                     ell *= 2  
90.                     if includ[k]:  
91.                         newv = newv ^ (ell*self.v_direction[j-k-1][i_dim])  
92.                 self.v_direction[j][i_dim] = newv  
93.         # done with for i_dim in range(1,self.dim)  
94.   
95.         # Multiply columns of v by appropriate power of 2  
96.         ell = 1  
97.         for j in range(self.bit_count-2,-1,-1):  
98.             ell *= 2  
99.             for i_dim in range(self.dim):  
100.                 self.v_direction[j][i_dim] = self.v_direction[j][i_dim]*ell  
101.         # 1/(common denom of the elements in v_direction)  
102.   
103.         self.last_denominator_inv = 1.0 /(2.0 * ell)  
104.   
105.         # final setup  
106.         self.sequence_count = 0  
107.         self.last_numerator_vec = [0]*self.dim  
108.         return  
109.   
110.     def __call__(self):  
111.         ell = 0  
112.         c = self.sequence_count  
113.         while 1:  
114.             ell += 1  
115.             if c % 2 == 1:  
116.                 c/=2  
117.             else:  
118.                 break  
119.         if ell > self.bit_count:  
120.             raise "Sobol failed for %d" % self.sequence_count  
121.   
122.         v = [0]*self.dim  
123.         for i_dimension in range(self.dim):  
124.             direction_i = self.v_direction[ell-1][i_dimension]  
125.             old_numerator_i = self.last_numerator_vec[i_dimension]  
126.             new_numerator_i = old_numerator_i ^ direction_i  
127.             self.last_numerator_vec[i_dimension] = new_numerator_i  
128.             v[i_dimension] = new_numerator_i*self.last_denominator_inv  
129.         self.sequence_count += 1  
130.         return v  

Together with this python implementation of the algorithm proposed by Peter J. Acklam to transform my uniform draw in gaussian draws:

1. def ltqnorm( p ):  
2.     """ 
3.     Modified from the author's original perl code (original comments follow below) 
4.     by dfield@yahoo-inc.com.  May 3, 2004. 
5.  
6.     Lower tail quantile for standard normal distribution function. 
7.  
8.     This function returns an approximation of the inverse cumulative 
9.     standard normal distribution function.  I.e., given P, it returns 
10.     an approximation to the X satisfying P = Pr{Z <= X} where Z is a 
11.     random variable from the standard normal distribution. 
12.  
13.     The algorithm uses a minimax approximation by rational functions 
14.     and the result has a relative error whose absolute value is less 
15.     than 1.15e-9. 
16.  
17.     Author:      Peter J. Acklam 
18.     Time-stamp:  2000-07-19 18:26:14 
19.     E-mail:      pjacklam@online.no 
20.     WWW URL:     http://home.online.no/~pjacklam 
21.     """  
22.   
23.     if p <= 0 or p >= 1:  
24.         # The original perl code exits here, we'll throw an exception instead  
25.         raise ValueError( "Argument to ltqnorm %f must be in open interval (0,1)" % p )  
26.   
27.     # Coefficients in rational approximations.  
28.     a = (-3.969683028665376e+01,  2.209460984245205e+02, \  
29.          -2.759285104469687e+02,  1.383577518672690e+02, \  
30.          -3.066479806614716e+01,  2.506628277459239e+00)  
31.     b = (-5.447609879822406e+01,  1.615858368580409e+02, \  
32.          -1.556989798598866e+02,  6.680131188771972e+01, \  
33.          -1.328068155288572e+01 )  
34.     c = (-7.784894002430293e-03, -3.223964580411365e-01, \  
35.          -2.400758277161838e+00, -2.549732539343734e+00, \  
36.           4.374664141464968e+00,  2.938163982698783e+00)  
37.     d = ( 7.784695709041462e-03,  3.224671290700398e-01, \  
38.           2.445134137142996e+00,  3.754408661907416e+00)  
39.   
40.     # Define break-points.  
41.     plow  = 0.02425  
42.     phigh = 1 - plow  
43.     
44.     # Rational approximation for lower region:  
45.     if p < plow:  
46.           
47.         q  = np.sqrt(-2*np.log(p))  
48.         return (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \  
49.                ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1)  
50.   
51.     # Rational approximation for upper region:  
52.     if phigh < p:  
53.         q = np.sqrt(-2*np.log(1-p))  
54.         return -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / \  
55.                 ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1)  
56.   
57.     # Rational approximation for central region:  
58.     q = p - 0.5  
59.     r = q*q  
60.     return (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / \  
61.            (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1)  

gives much faster code. I then just modified my SobolRNG as follows:

1. def SobolRNG(I):  
2.         sseq=SobolSequence(1)       
3.         if AP == True :  
4.            ran = np.array([ltqnorm(sseq()[0]) for i in xrange(I/2)])  
5.            ran = np.concatenate (( ran ,- ran ))  
6.         else :  
7.            ran = np.array([ltqnorm(sseq()[0]) for i in xrange(I)])  
8.         if MM == True :  
9.            ran =ran - np.mean ( ran)  
10.            ran = ran / np.std ( ran )  
11.         return ran         

As you can see, you first initialize the Sobol object and then successive calls gives the sequence. So it is a scalar, that passes in a scalar function ltqnorm. I was expecting this to be slower than vector based construction but there might be some numpy optimization in the background because it is called for the creation of the array. The speed difference with the older code is significant.