# Monte Carlo European Option Pricing

I've written code below that simulates GBM paths for determining the price of a given European call option and put option. The stock is priced at 150 USD, strike price at 155 USD, risk-free rate was assumed to be 0.02, expected return was equal to 0.05, volatility at 0.1 and it's one year to maturity.

function [call, put] = monte_carlo_price(S_init, K, T, r, mu, sigma, n)
% Computes European call and put options using Monte Carlo simulations
% 'S_init' is the current underlying stock price
% 'K' is the strike price
% 'T' is years to maturity
% 'r' is the risk-free rate
% 'mu' is the expected return
% 'sigma' is the volatility
%----------------------------------------------------------------------

Paths = n;
Steps = T*365;
MC_scenes = zeros(Steps+1, Paths);
dT = T/Steps;

% First asset price is the initial price
MC_scenes(1,:) = S_init;

% Generate paths
for iPath = 1:Paths
for iStep = 1:Steps
MC_scenes(iStep+1, iPath) = MC_scenes(iStep, iPath) * exp((mu - 0.5*sigma^2)*dT + sigma * sqrt(dT)*normrnd(0,1));
end
end

% Calculate put and call option payoffs
putPayoff = max(K-MC_scenes(end,:),0);
callPayoff = max(MC_scenes(end,:)-K,0);

% Discount prices back to present day
put = mean(putPayoff)* exp(-r*T) ;
call = mean(callPayoff) * exp(-r*T);

end

My problem is that I'm trying to come up with an efficient way to get my simulation results to within 10 cents of the Black-Scholes output (5.10 USD for the call and 7.04 USD for the put). I tried increasing my 'Paths' (number of scenarios) through iterating and I got what I needed, but it took way too long to run. Is there way I can optimize setting the 'Paths' and even 'Steps' get to within the Black-Scholes outputs??
Thanks.

• In the line MC_scenes(iStep+1, iPath) = paths(iStep, iPath) * exp((mu - 0.5*sigma^2)*dT + sigma * sqrt(dT)*normrnd(0,1)); - shouldn't the first term on the right hand side be MC_scenes(iStep, iPath) instead? – LocalVolatility Mar 25 '17 at 21:58
• MC_scenes(iStep, :) holds the initial stock price 'S_init' where iStep = 1. That is why my loops does MC_scenes(iSteps + 1, iPath). – Marcus L Mar 25 '17 at 22:07
• My point is that you never defined the variable paths and I think you meant to write MC_scenes instead..? – LocalVolatility Mar 25 '17 at 22:08
• Ah sorry about that, 'Paths' is the number of scenarios I will run (I initially set 10,000 scenarios) – Marcus L Mar 25 '17 at 22:46
• I am not referring to Paths (capital) but paths (lower case) on the right hand side of the line I referenced. – LocalVolatility Mar 25 '17 at 22:48

A few suggestions:

1. As your underlying follows a geometric Brownian motion and you are solely interested in pricing European options, there is no need to simulate intermediate steps. Since your solution is exact, you can directly sample $S_T$ as $$S_T = S_0 \exp \left\{ \left( r - \frac{1}{2} \sigma^2 \right) T + \sigma \sqrt{T} Z \right\},$$ where $Z \sim \mathcal{N}(0, 1)$.

2. Try replacing the remaining for-loop with matrix operations. Vectorized code tends to be significantly faster in MATLAB.

3. Regarding stopping criterion. In general you don't know the closed-form solution for the option that you are pricing via the Monte Carlo simulation. Otherwise you wouldn't do it in the first place, except for as an exercise. So using the distance to the true analytical solution as a stopping criterion is off the table. Instead, I suggest to compute the Monte Carlo standard error $$s_n = \sqrt{\frac{1}{n - 1} \sum_{i = 1}^n \left( X_i - \bar{X}_n \right)^2},$$ where $\left\{ X_i \right\}_{i = 1}^n$ are the sequence of sample prices for your $n$ paths. You can then construct a confidence interval for your Monte Carlo estimate for some probability and stop your simulation once it is sufficiently small.

To see why your distance-based stopping criterion is problematic, consider a very small number of sample paths that happen to be such that the corresponding average discounted payoff is very close to the analytical price. This is by pure chance and the standard error of your estimate will be very large.

4. Note that you should be simulation the asset under the risk-neutral probability measure. Here, the drift of the asset is $r$ and not $\mu$ as in your code. The latter is irrelevant.

5. In order to further improve the convergence, you could use antithetics. I.e. for a sample of $S_T$ generated using the random normal variate $Z$, you generate a second sample of $S_T$ using $-Z$.

• Hey LocalVolatility, thanks for your quick responds. 1) I am directly sampling 'S' using GBM whereby I keep generating random numbers for 'Z' (or am I not doing this right?) 2) Yes, I am aware of the advantage of matrix operations and I will try to implement that over my for loops instead – Marcus L Mar 25 '17 at 22:14
• 3) Is it not possible to set up a stopping criteria to be that the code stops when the difference between the monte carlo put(call) option pricing is within +/- 0.10 of the Black-Scholes solution (so they are off by a dime); I' did that initially and keep increasing my number of simulation paths, but the code took forever to complete and is not very viable – Marcus L Mar 25 '17 at 22:14
• Please read my answer a bit more carefully: 1) My point is that you don't need to simulate any intermediate steps but can directly go from $t = 0$ to $t = T$ in one step as the solution is exact and you are only interested in the terminal value $S_T$ anyways. This will significantly speed up your code. 3) You can do that. However, as I wrote, this is unrealistic in all interesting cases. If you had the analytical solution, you wouldn't run the MC in the first place. – LocalVolatility Mar 25 '17 at 22:17
• 1) Ah yes, I understand what you mean and that is actually what I did (since maturity is only one year); however I do not think I can modify that aspect since my instructor specifically said that we must implement GBM in this manner........ 3) The only reason I am using Monte Carlo is to compare the performance between it and Black-Scholes equation. Having the difference between results be within 10 cents of one another is one of the criteria to do that (instead of simulating a large amount of scenarios for the monte carlo to converge) – Marcus L Mar 25 '17 at 22:34
• 1) This contradicts your question. You are asking how to optimize the number of steps Steps. I suggested you should use Steps = 1 as choosing higher values doesn't improve accuracy. Now you say, you need to use 365 steps per year. 3) Again, there is not much I can suggest except for vectorizing given that the stopping criteria have already been set by your instructor. Also - please flag homework questions as such in the future. – LocalVolatility Mar 25 '17 at 22:40