File 1: option_montecarlo_simulation.m % option_montecarlo_simulation % Author: Ryan Dutton. Project: Bachelor Thesis % The main routine that calculates the payoff from the sample % asset paths, and computes the price of the call and put options. % It also prints out the price of call and put options using the % Black-Scholes formula with the same parameters. S0 = 100 % Initial value of stock mu = .05 % Drift rate of the return, should be the same % as the risk-free interest rate sigma = .12 % Standard deviation of the periodic return TPeriod = 2 % Time period in year num_path = 2000 % Number of asset paths for the simulation step_length = 1/255 % Time length of one step inside a period % We are using 255 steps int_rate= .05 % Risk-free interest rate strike = 125 % Strike price [sample, time] = simulate_asset_path(sigma,TPeriod,num_path,S0,mu,step_length); figure('Color','white'); plot(time,sample,'linewidth',2); axis([time(1) time(end) min(sample(:)) max(sample(:))]); xlabel('Time');ylabel('Stock Price'); title(sprintf('No. of Asset Path Simulated: %d', size(sample,2))); positive_payoff = (sample(511,:) - strike) > 0; MCE_call_payoff = (sample(511,:) - strike).*positive_payoff; negative_payoff = (strike - sample(511,:)) > 0; MCE_put_payoff = (strike - sample(511,:)).*negative_payoff; MCE_call = mean(MCE_call_payoff) * exp(-int_rate * TPeriod) MCE_put = mean(MCE_put_payoff) * exp(-int_rate * TPeriod) [bs_call, bs_put] = blsprice (S0,strike,int_rate,TPeriod,sigma) --------------------------------------------------------------------- File 2: simulate_asset_path.m function [sample, time] = simulate_asset_path(sigma,TPeriod,num_path,S0,mu,step_length) % Author: Ryan Dutton. Project: Bachelor Thesis % Generates the asset paths from simulation, and returns them in % matrix Sample. Example: % [Sample, time] = simulate_asset_path(.12,5,50,100,.05,1/255); num_step = ceil(TPeriod/step_length); RandomVectors = sigma*sqrt(step_length)*cumsum(randn(num_step, num_path)); Drifts = repmat((mu - sigma^2/2) *step_length * (1:num_step)',1,num_path); sample = [repmat(S0,1,num_path); S0 * exp( Drifts + RandomVectors)]; if nargout > 1 time = [0;step_length * (1:num_step)']; end