Welcome to THETAWIKI. If you like to create or edit a page please make sure to login or register an account. All registered users please make sure to provide a valid email address.

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:

$d S(t) = \mu S(t) dt \ + \ \sigma S(t) \ dW^P(t), \qquad (1)$

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

$d \frac{S(t)}{B(t)}= S(t)d\frac{1}{B(t)} \ + \ \frac{1}{B(t)}dS(t) \ + \ \left \langle dS(t), d\frac{1}{B(t)} \right \rangle$
$= S(t)\left [ -\frac{dB(t)}{B(t)^2} \ + \ \frac{1}{B(t)^3}\left \langle dB(t), dB(t) \right \rangle \right ] \ + \ \frac{1}{B(t)} \left [\mu S(t) d t \ + \ \sigma S(t) d W^P(t) \right ]$
$= S(t)\left [ -\frac{r B(t)dt}{B(t)^2} \right ] \ + \ \frac{1}{B(t)} \left [\mu S(t) d t \ + \ \sigma S(t) d W^P(t) \right ]$
$= \frac{S(t)}{B(t)}\left [ (\mu \ - \ r)dt \ + \ \sigma d W^P(t) \right ], \qquad (2)$

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

$d \frac{S(t)}{B(t)} = \sigma \frac{S(t)}{B(t)} d W^Q(t), \qquad (3)$

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.

$d W^Q(t) = d W^P(t) \ + \ \frac {\mu \ - \ r}{\sigma} dt. \qquad (4)$

Using equations (1) and (4), it is easy to obtain the stock price process under the risk-neutral measure Q:

$d S(t) = \mu S(t) dt \ + \ \sigma S(t) \ d W^P(t)$
$= \mu S(t) dt \ + \ \sigma S(t) \left [ d W^Q(t) \ - \ \frac {\mu - r}{\sigma} dt \right ]$
$= r S(t) dt \ + \ \sigma S(t) \ d W^Q(t), \qquad (5)$

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

1-Dimensional

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

N-Dimensional

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