Geometric Brownian Motion

 Share 

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 $ \mu $ is the drift rate, $ \sigma $ the volatility of $ S $ and $ dW^P(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) = r B(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 $ d W^Q(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