Monte Carlo European Option Pricing

Question

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.

Answer 1

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 \begin{equation} S_T = S_0 \exp \left\{ \left( r - \frac{1}{2} \sigma^2 \right) T + \sigma \sqrt{T} Z \right\}, \end{equation} 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 \begin{equation} s_n = \sqrt{\frac{1}{n - 1} \sum_{i = 1}^n \left( X_i - \bar{X}_n \right)^2}, \end{equation} 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$.