Monte Carlo simulations and option pricing by Bingqian Lu Undergraduate Mathematics Department Pennsylvania State University University Park, PA 16802 Project Supervisor: Professor Anna Mazzucato July, 2011

Abstract Monte Carlo simulation is a legitimate and widely used technique for dealing with uncertainty in many aspects of business operations. The purpose of this report is to explore the application of this technique to the stock volality and to test its accuracy by comparing the result computed by Monte Carlo Estimate with the result of Black-Schole model and the Variance Reduction by Antitheric Variattes. The mathematical computer softwear application that we use to compute and test the relationship between the sample size and the accuracy of Monte Carlo Simulation is itshapeMathematica. It also provides numerical and geometrical evidence for our conclusion.

0.1

Introduction to Monte Carlo Simulaion

Monte Carlo Option Price is a method often used in Mathematical finance to calculate the value of an option with multiple sources of uncertainties and random features, such as changing interest rates, stock prices or exchange rates, etc.. This method is called Monte Carlo simulation, naming after the city of Monte Carlo, which is noted for its casinos. In my project, I use Mathematica, a mathematics computer software, we can easily create a sequence of random number indicating the uncertainties that we might have for the stock prices for example.

0.2

Pricing Financial Options by Flipping a Coin

A distcrete model for change in price of a stock over a time interval [0,T] is √ Sn+1 = Sn + µSn ∆t + σSn εn+1 ∆t, S0 = s (1) where Sn = Stn is the stock price at time tn = n∆t, n = 0, 1, ..., N − 1, ∆t = T /N , µ is the annual growth rate of the stock, and σ is a measure of the stocks annual price volatility or tendency to fluctuate. Highly volatile stocks have large values of σ. Each term in sequence ε1 , ε2 ... takes on the value of 1 or -1 depending on the outcoming value of a coin tossing experiment, heads or tails respectively. In other words, for each n=1,2,... εn = 1 with probability = 1/2 −1 with probability = 1/2 (2)

By using Mathematica, it is very easy to create a sequence of random number. With this sequence, the equation (1) can then be used to simulate a sample path or trajectory of stock prices, {s, S1 , S2 , ..., SN }. For our purpose here, it has been shown as a relatively accurate method of pricing options and very useful for options that depend on paths.

0.3
0.3.1

Proof of highly volatile stocks have large values for σ σ = 0.01

Let us simulate several sample trajectories of (1) for the following parameter values and plot the trajectoris: µ = 0.12, σ = 0.01, T = 1, s = $40, N = 254. The following figures are the graphs that we got for σ = 0.01 Since we have ε to be the value generated by flipping a coin, it gives us arbitrary values and thus, we have different graphs for parameter σ = 0.01, µ = 0.12, T = 1, s = $40, N = 254 This is another possibile graph. 1

Figure 1: Let M be the value of S254 of different trajectories, k is the number of trajectories

0.3.2

σ = 0.7

Then we repeated the experiment using the value of σ = 0.7 for the volality and other parameters remain the same. Similarly, it should have a number of different graphs due to the arbitrary value of ε we generated by Mathematica. The following figures are the graphs that we got for σ = 0.7 Conclusion: From the two experiments above with the large different σ and constant other parameters, we can tell that the larger the σ, the greater degree of variability in their behavior forthe ε ’s it is permissible to use random number generator that creates normally distributed random numbers with mean zero and variance one. Recall that the standard normal distribution has the bell-shape with a standard deviation of 1.0 and standard normal random variable has a mean of zero.

0.4
0.4.1

Monte Carlo Method vs. Black-Scholes Model
Monte Carlo Method and its computing

Monte Carlo Method √ In the formular (1), the random terms Sn εn+1 ∆t on the right-hand side can be consider as shocks or distrubances that model functuations in the stock price. After repeatedly simulating stock price trajectories, as we did in the previous chapter, and computing appropriate averages, it is possible to obtain estimates of the price of a European call option, a type of

2

45

44

43

42

41

0.2

0.4

0.6

0.8

1.0

Figure 2: S254 = 45.415, σ = 0.01 financial derivative. A statistical simulation algorithm of this type is what we known as ”Monte Carlo method” A European call option is a contract between two parties, a holder and a writer, whereby, for a premium paid to the writer, the holder can purchase the stock at a future date T (the expiration date) at a price K (the strike price) agreed upon in the contract. If the buyer elect to exercise the option on the expiration date, the writer is obligated to sell the inderlying stock to the buyer at the price K, the strike price. Thus, the option has a payoff function f (S) = max(S − K, 0) (3) where S = S(T ) is the price of the underlying stock at the time T when the option expires. This equation (3) produces one possible option value at expiration and after computing this thousands of times in order to obtain a feel for the possible error in estimating the price. Equation(3) is also known as the value of the option at time T since if S(T ) > K, the holder can purchase, at price K, stock with market value S(T) and thereby make a profit equal toS(T ) − K not counting the option premium. However, on the other hand, if S(T ) < K, the holder will simply let the option expire since there would be no reason to purchase stock at a price that exceeds the market value. In other words, the option valuation problem is determine the correct and fair price of the option at the time that the holder and writer enter into the contract. In order to estimate the price call of a call option using a Monte Carlo method, an ensemble SN = S (k) (T ), k = 1, ...M 3
(k)

(4)

45

44

43

42

41

0.2

0.4

0.6

0.8

1.0

Figure 3: S254 =45.1316, σ = 0.01 of M stock orices at expiration is generated using the difference equation √ (k) (k) (k) (k) (k) (k) S0 = s (5) Sn+1 = Sn + rSn ∆t + σSn εn+1 ∆t, Equation (5) is identical to equation (1) for eachk = 1, ..., M , except the growth rate µ is replaces by the annual interest r that it costs the writer to borrow money. Option pricing theory requires that the average value of (k) the payoffs f (SN 0, k = 1, ..., M be equal to the compounded total return ˆ obtained by investing the option premium, C(s), at rate r over the life of option, M 1 (k) ˆ f (sN ) = (1 + r∆t)N C(s). (6) M k=1 ˆ Solving(6) for C(s) yields the Monte Carlo estimate ˆ C(s) = (1 + r∆t)−N 1 M
M

f (sN ) k=1 (k)

(7)

ˆ for the option price. So, the Monte Carlo estimateC(s) is the present value of the average of the payoffs computed using rules of compound interest.

0.4.2

Computing Monte Carlo Estimate

We use equation (7) to compute a Monte Carlo estimate of the value of a five 5 month call option, in other word T = 12 years, for the following parameter values: r = 0.06, σ = 0.2, N = 254, andK = $50. N is the number of times of steps for each trajectories.

4

45

44

43

42

41

0.2

0.4

0.6

0.8

1.0

Figure 4: S254 =45.1881, σ = 0.01

0.5

Comparing to the Exact Black-Scholes Formular

Monte Carlo has been used to price standard European options, but as we known that Black-Scholes model is the correct method of pricing these options, so it is not necessary to use Monte Carlo simulation. Here is the formular for exact Black-Scholes model: √ s K −d1 C(s) = erf c( √ ) − e−rT erf c(f rac−d2 2) (8) 2 2 2 where √ s σ2 1 d1 = √ [ln( ) + (r + )T ], d2 = d1 − σ T k 2 σ T and erfc(x) is the complementary error function, 2 erf c(x) = √ π
∞ x

(9)

e−t dt

2

(10)

Now we insert all data we have to the Black-Schole formula to check the accuracy of our results by comparing the Monte Carlo approximation with the value computed from exact Black-Schole formula. We generated BlackScholes Model with parameterr = 0.06, σ = 0.2, K = $50, k = 1, ..., M (whereT = N ∆t), N = 200. And we got: C(40) = 1.01189, C(45) = 2.71716 and C(50) = 5.49477. The error of the Monte Carlo Estimate seems to be very large. Thus we repeat the previous procedure and increased our sample size (M) to 50,000 and 100,000. Then, I made a chart to check if the accuracy of Monte Carlo Simulation increases by the increasing of the sample size.

5

45

44

43

42

41

0.2

0.4

0.6

0.8

1.0

Figure 5: S254 =45.5859, σ = 0.01 Comparison of the accuracy of the Monte Carlo Estimate to the BlackSchole Model

MCE(1,000) ˆ C(40) = 1.67263 ˆ C(45) = 4.6343 ˆ C(50) = 8.1409

MCE (10,000) ˆ C(40) = 1.61767 ˆ C(45) = 4.24609 ˆ C(50) = 7.92441

MCE (50,000) ˆ C(40) = 1.76289 ˆ C(45) = 4.3014 ˆ C(50) = 7.7889

MCE (100,000) ˆ C(40) = 1.7496 ˆ C(45) = 4.2097 ˆ C(50) = 7.6864

Black-Schole Model BS(40) = 1.71179 BS(45) = 4.11716 BS(50) = 7.49477

By the comparing Monte Carlo Estimate and the Black-Schole modle, we can tell quite easily that the error gets smaller as the sample size increases. Here we calculated the relative error by the equation M onteCarloEstimate − Black − Scholemodle Black − ScholeM odle The results are error(1,000) E(40) = 0.60593 E(45) = 0.32330 E(50) = 0.24272 error (10,000) E(40) = 0.30232 E(45) = 0.19466 E(50) = 0.07819 error(50,000) E(40) = 0.14037 E(45) = 0.13862 E(50) = 0.05254 error (100,000) E(40) = 0.14037 E(45) = 0.08833 E(50) = 0.02157 Black-Schole Model BS(40) = 1.71179 BS(45) = 4.11716 BS(50) = 7.49477 (11)

Conclusion: Monte Carlo Simulation gives the option price is a sample average, thus according to the most elementary principle of statistics, its standard deviation is the standard deviation of the sample divided by the square root of the sample size. So, the error reduces at the rate of 1 over the square root of the sample size. To sum up, the accuracy of Monte Carlo Simulation is increasing by increasing the size of the sample. 6

Figure 6: replaced previous value of sigma with 0.7

0.6

Comparing to the Variance Reduction by Antitheric Variattes

Variance Reduction by Antithetic Variates is a simple and more widely used way to increase the accuracy of the Monte Carlo Simulation. It is the technique used in some certain situations with an additional increase in computational complexity is the method of antithetic variates. In order to achieve greater accuracy, one method of doing so is simple and automatically doubles the sample size with only a minimum increase in computational time. This is called the antithetic variate method. Because we are generating obervations of a standard normal random variable which is distributed with a mean of zero, a variance of 1.0 and symmetric, there is an equally likely chance of having drawn the observed value times−1. Thus, for each arbitrary ε we draw, there should be an artificially observed companion observation of −ε that can be legitimately created by us. This is the antithetic variate. For each k=1,...,M use the sequence ε1 , ..., εN −1
(k) (k)

(12)
(k) (k)

k+1 in equation (5) to simulate a payoff f (SN )and also use the sequence −ε1 , ..., −εN −1 k− in equation (5) to simulate an associated payoff f (SN ). Now the payoffs k+ k− are simulated inpairs f (SN ), f (SN ) This is the Mathematica Program that we ran to evaluate the Variance Reduction. After comparing Monte Carlo Simulation with Variance Reduction by Antithetic Variates, We made a table of data.

7

100

80

60

40

20

0.2

0.4

0.6

0.8

1.0

Figure 7: S254 =47.084, σ = 0.7

k(1,000) V (40) = 1.74309 V (45) = 4.4789 V (50) = 7.94642

k(5,000) V (40) = 1.66685 V (45) = 4.305103 V (50) = 7.67707

k(50,000) V (40) = 1.69966 V (45) = 4.1778 V (50) = 7.62341

k (100,000) V (40) = 1.70674 V (45) = 4.099835 V (50) = 7.53716

BS Model BS(40) = 1.71179 BS(45) = 4.11716 BS(50) = 7.49477

While the table for the data of Monte Carlo Simulation we get under the same condition is: MCE(1,000) ˆ C(40) = 1.67263 ˆ C(45) = 4.6343 ˆ C(50) = 8.1409 MCE (5,000) ˆ C(40) = 1.72486 ˆ C(45) = 4.34609 ˆ C(50) = 7.92241 MCE (50,000) ˆ C(40) = 1.76289 ˆ C(45) = 4.3014 ˆ C(50) = 7.7889 MCE (100,000) ˆ C(40) = 1.7496 ˆ C(45) = 4.2097 ˆ C(50) = 7.6864

By the chart of Monte Carlo Estimate and the Variance Reduction by Antithetic Variates, we can tell quite easily that the error gets smaller as the sample size increases. Here we calculated the relative error by the equation M onteCarloEstimate − V arianceReductionbyAntitheticV ariates (13) V arianceReductionbyAntitheticV ariates The results are error(1,000) E(40) = 0.0402 E(45) = 0.0347 E(50) = 0.0245 error (50,000) E(40) = 0.0372 E(45) = 0.0296 E(50) = 0.0217 error (100,000) E(40) = 0.0251 E(45) = 0.0268 E(50) = 0.0198

8

70

60

50

40

0.2

0.4

0.6

0.8

1.0

Figure 8: S254 =54.7495, σ = 0.7

0.7

Random Walk ds = µSdt + σS β dw t (14)

Random Walk

this is the random walk, where ds is S(tk +∆t)−S(tk ), so this can be written as √ S(tk + ∆t) − S(tk ) = µS(tk )∆t + σS(tk )β εk ∆t (15) We randomly picked two different positive numbers, one greater than 1 and one less than one. β as 0.5 and 2 and we made mathematica ran twice while keep every other variables constant. Here are some graphs that we plot using Mathematica. Figure 19 is a graph for β = 0.5 In addition, we plot the graphs when β = 2 as well. During Mathematica’s computation, we got some values that are overflowing. Within the values that are available, we plotted some corresponding graphs. Figure 20 is a graph for the random walk when β = 2. By using 1000 as a sample size, we used the same program to compute ˆM (10), CM (20), CM (40) and CM (160) when β = 0.5, where CM (10) is ˆ ˆ ˆ ˆ C ˆM (20) is the Monte Carlo the Monte Carlo estimate for sample size 1000, C ˆ simulation for sample size 2000 ans so forth. We get CM (10) = 3.80332, ˆM (20) = 4.75116, CM (40) = 6.09321 and CM (160) = 7.63962 . Now we ˆ ˆ C are using the equation log2 ( ˆ ˆ CM (10) − CM (20) ) ˆ ˆ CM (40) − CM (160) (16)

9

50

40

30

20

10

0.2

0.4

0.6

0.8

1.0

Figure 9: S254 =7.95396, σ = 0.7 to see the error. After plugging the numbers, we get this value of the equation to be -0.67621. This value is very close to 1 as desired. 2 By the same procedure, we used Mathematica to generate the value for ˆ ˆ ˆ ˆ CM (10), CM (20), CM (40) and CM (160) when β = 2. There are some overflow in the result which gives us trouble for going further and also tells us that when β > 1, there will be overflow.

10

40

35

30

0.2

0.4

0.6

0.8

1.0

Figure 10: S254 =22.7805, σ = 0.7

Array of S for monte carlo estimate1.png

Figure 11: We first array for S in order to make space of memory for the ˆ data we will get later. Here, P(k) is theC(k)

11

Mathematica program for Monte Carlo Estimate1.png

Figure 12: This is the Matehamatica program we made and data we got for the different S(254), after 10 times of computing, named as M[k]=S[255, k] ˆ in mathematica. Again, here P(k) is the C(k). So, according to the data we got, there are 10 different P[k+1], which in other word, 10 different Monte ˆ Carlo estimate C(k + 1) after 10 times of computing

12

S(0)=40.png

Figure 13: This result is the Monte Carlo estimate corresponding to our current stock prices of S(0) = s = 40 after using 254 steps (the number of N) and M ∼ 10, 000 for each trajectories for each Monte Carlo estimate. =

13

S(0)=45.png

Figure 14: This result is the Monte Carlo estimate corresponding to our current stock prices of S(0) = s = 45 after using 254 steps (the number of N) and M ∼ 10, 000 for each trajectories for each Monte Carlo estimate. =

14

S(0)=50.png

Figure 15: This result is the Monte Carlo estimate corresponding to our current stock prices of S(0) = s = 50 after using 254 steps (the number of N) and M ∼ 10, 000 for each trajectories for each Monte Carlo estimate. =

Figure 16: This the program we wrote with Mathematica for the Variance Reduction by Antithetic Variattes. We have S for the positive ε and Y for the negative ε and we get the average of these value by adding them up and divide by 2.

15

Figure 17: This the program we wrote with Mathematica for different values of β.

46

45

44

43

42

41

0.2

0.4

0.6

0.8

1.0

Figure 18: This is the graph when β = 0.5

16

46

45

44

43

42

41

0.2

0.4

0.6

0.8

1.0

Figure 19: This is how the graph looks like when β = 0.5

7

6

5

4

0.2

0.4

0.6

0.8

1.0

Figure 20: This is the graph when β = 2

17

20

15

10

5

0.2

0.4

0.6

0.8

1.0

Figure 21: This is how the graph looks when β = 2

18