Geometric Brownian Motion
GBM Model overview
The underlying assumption for Black-Scholes Option Pricing is that stock prices follow a Geometric Brownian motion process. The stochastic differential equation(SDE) for the stock prices take the following function form:
where μ is the drift rate, σ the volatility of S and dWP(t) the increment of a standard Wiener process under the physical measure P.
To derive the stock price process under the risk-neutral measure Q, we proceed as follows. Use the bank account B(t) (defined in Discounting, B(t) is here short for B(0,t)) as numeraire, the discounted stock price process takes the following form
where in the above equation, we have substitued the bank account process dB(t) = rB(t)dt for B(t), and the stock price process in equation (1) for dS(t). So far, we have derived the discounted stock price process under the physical measure P. To change the underlying probability measure from the physical measure P to the risk-neutral measure Q, first note that stock price discounted with the numeraire bank account will itself be a Martingale under this numeraire induced probability measure - namely the risk-neutral measure Q. Mathematically, the discounted stock price process in equation (2) re-expressed under the risk-neutral measure Q is
where dWQ(t) is the increment of a wiener process under the risk-neutral measure Q. Compare equation (2) with equation (3), the Girsanov theorm gives the measure change from P to Q, i.e.
Using equations (1) and (4), it is easy to obtain the stock price process under the risk-neutral measure Q:
namely, under the risk-neutral pricing measure Q, stock prices grow at the risk-free rate r.
Having derived the stochastic differential equation for stock prices under the risk-neutral measure Q, we show, in the following, some code examples in Theta Suite for simulating stock price process and for analytical call and put option prices.
GBM Process simulation
Implementation in ThetaML
The following is an ThetaML implementation for the stock price process (5).
model GBM %This model computes future stock prices 'S'; the stock prices 'S' follow %a Geometric Brownian motion process under the risk-neutral measure; %the stock prices 'S' is a process variable, 'S' implicitly incorporate %scenario and time indexes import S0 "Initial stock price" import r "Risk-free interest rate" import sigma "Volatility of stock" export S "GBM stock prices" %initialize the stock prices at 'S0' S = S0 %'loop ... inf' is an infinite loop; this infinite loop computes a stock %price process of an arbitrary length; the lifetime of the infinite loop is %automatically extended to the desired length depending on a specific pricing %application loop inf %the ThetaML command 'theta' passes time by '@dt' units %the ThetaML parameter '@dt' denotes an arbitrary time unit; its specific %value depends on a specific pricing application theta @dt %update the stock prices for the time interval '@dt' S = S * exp( (r - 0.5*sigma^2)*@dt + sigma*sqrt(@dt)*randn() ) end end
The following code examples are Theta.m implementations for the stock price process (5) under the risk-neutral measure Q.
Implementation as a simple Matlab stepping function Theta.m
function state = Theta(dt, state) %This function returns a struct object 'state'; the struct 'state' has subfields %'S' for the stock price, 'sigma' for the stock price volatility, and 'r' for the %risk-free rate. if nargin == 0 %if no function arguments state.S.comment = 'Stock price'; state.sigma.comment = 'Volatility'; state.r.comment = 'Risk-free rate'; else % risk-free rate r = state.r(1); % volatility of stock sigma = state.sigma(1); % update the stock price process state.S = state.S .* exp( (r - 0.5*sigma^2)*dt + sqrt(dt)*sigma*randn(size(state.S))); end end
function state = Theta(dt, state) %This function returns a struct object 'state'; the struct 'state' has subfields %'S' for stock prices, 'rho' for the correlation among the stock prices, %'sigma' for stock price volatilities if nargin == 0 state.S.comment = 'Stock prices'; state.rho.comment = 'Correlation'; state.sigma.comment = 'Volatilities'; else r = 0.05; rho = mean(state.rho); sigma = ones(1,size(state.S,2)) * mean(state.sigma); Z = zeros(length(sigma), length(sigma)); %construct the variance-covariance matrix for i = 1 : length(sigma) for k = 1 : length(sigma) if i == k %matrix diagonal Z(i,i) = sigma(i)^2; else %upper and lower matrix Z(i,k) = rho*sigma(i) * sigma(k); end end end %do the Cholesky decomposition of the variance-covariance matrice 'Z' A = (chol(Z))'; %sample standard normal variates B = randn(size(state.S')); %update the N-dimensional stock price processes for i = 1 : length(sigma) state.S(:,i) = state.S(:,i) .* exp( (r - 0.5*sigma(i)^2)*dt + sqrt(dt)*(A(i,:)*B)'); end end end
Analytic Black-Scholes option pricing in Matlab
The price of European puts and calls can be obtained by calling the Matlab function theta_bls().
function result = theta_bls(S, K, sigma, r, T, call, greek) %This function returns option prices as well as analytic option 'delta' and 'gamma' %respectively for European Call and Put options %function parameters: % S : asset price % K : option strike price % sigma: volatility of asset price % r : risk-free rate % T : option maturity time % call: true for call, otherwise put % greek: a string for 'value' that returns the option value, % 'delta' that returns the option delta, % or 'gamma' that returns the option gamma if T <= 0 if call V = max(S - K, 0); else V = max(K - S, 0); end return end d1 = ( log(S./K) + (r + (sigma.^2)/2)*T ) ./ (sigma * sqrt(T)); d2 = d1 - sigma * sqrt(T); switch call case true switch greek case 'value' result = S.*normcdf(d1,0,1) - K .* exp(-r.*T).*normcdf(d2); case 'delta' result = normcdf(d1,0,1)'; case 'gamma' result = (normpdf(d1,0,1)./(S.*sigma.*sqrt(T)))'; end otherwise switch greek case 'value' result = -S.*normcdf(-d1,0,1) + K .* exp(-r.*T).*normcdf(-d2); case 'delta' result = normcdf(d1,0,1)'-1; case 'gamma' result = (normpdf(d1,0,1)./(S.*sigma.*sqrt(T)))'; end end end