Parsiad Azimzadeh

I have recently taken on the role of maintaining the GNU Octave financial package. I recently implemented stochastic differential equation simulation, which I discuss in this post.

If you don't already have a copy of GNU Octave and the financial package, see this post for installation instructions.

Consider the following SDE: $$dX_t = F(t, X_t) dt + G(t, X_t) dW_t$$ where $W$ is a standard Brownian motion. With the Octave financial package, you can now simulate sample paths of this SDE. Let's take a moment to motivate why you might want to do this...

A simple example: pricing a European call

A classic problem in finance is that of pricing a European option, for which one needs to compute an expectation involving the random variable $X_T$ at some fixed time $T>0$, referred to as the expiry time. For example, the Black-Scholes call pricing problem is that of computing $$E\left[e^{-rT} (X_T-K)^+\right]$$ where $X$ follows a geometric Brownian motion (GBM): $$dX_t = (r-\delta) X_t dt + \sigma X_t dW_t.$$ Here, $r$ is the risk-free rate, $\delta$ is the dividend rate, and $\sigma$ is the volatility.

Approximating such an expectation using a sample mean is referred to as Monte Carlo integration (a.k.a. Monte Carlo simulation). Though the Black-Scholes pricing problem happens to be one in which a closed-form solution is known, as an expository example, let's perform Monte Carlo integration to approximate it using an SDE simulation:

% Test parameters
X_0 = 100.; K = 100.; r = 0.04; delta = 0.01; sigma = 0.2; T = 1.;
Simulations = 1e6; Timesteps = 10;

SDE = gbm (r - delta, sigma, "StartState", X_0);
[Paths, ~, ~] = simByEuler (SDE, 1, "DeltaTime", T, "NTRIALS", Simulations, "NSTEPS", Timesteps, "Antithetic", true);

% Monte Carlo price
CallPrice_MC = exp (-r * T) * mean (max (Paths(end, 1, :) - K, 0.));

% Compare with the exact answer (Black-Scholes formula): 9.3197

The gbm function is used to generate an object describing geometric Brownian motion (GBM). Under the hood, it invokes the sde constructor, which is capable of constructing more general SDEs.

Timing

The GNU Octave financial implementation uses broadcasting to speed up computation. Here is a speed comparison of the above with the MATLAB Financial Toolbox, under a varying number of timesteps:

Timesteps MATLAB Financial Toolbox (secs.) GNU Octave financial package (secs.)
16 0.543231 0.048691
32 1.053423 0.064110
64 2.167072 0.097092
128 4.191894 0.162552
256 8.361655 0.294098
512 16.609718 0.568558
1024 32.839757 1.136864

While both implementations scale more-or-less linearly, the GNU Octave financial package implementation greatly outperforms its MATLAB counterpart.

A multi-dimensional example: pricing a Basket call

Consider two assets, $X^1$ and $X^2$, and the basket call option pricing problem $$E\left[e^{-rT}\left(\max\left(X^1_T,X^2_T\right)-K\right)^+\right]$$ where $$dX_t^i = (r-\delta^i) X^i_t dt + \sigma^i X^i_t dW^i_t \text{ for } i = 1,2$$ and the correlation between the Wiener processes is $dW^1_t dW^2_t = \rho dt$. Sample code for this example is below:

% Test parameters
X1_0 = 40.; X2_0 = 40.; K = 40.; r = 0.05; delta1 = 0.; delta2 = 0.; sigma1 = 0.5; sigma2 = 0.5; T = 0.25; rho = 0.3;
Simulations = 1e5; Timesteps = 10;

SDE = gbm ([r-delta1 0.; 0. r-delta2], [sigma1 0.; 0. sigma2], "StartState", [X1_0; X2_0], "Correlation", [1 rho; rho 1]);
[Paths, ~, ~] = simulate (SDE, 1, "DeltaTime", T, "NTRIALS", Simulations, "NSTEPS", Timesteps, "Antithetic", true);

Max_Xi_T = max (Paths(end, :, :));
BasketCallPrice_MC = exp (-r * T) * mean (max (Max_Xi_T - K, 0.));

% Compare with the exact answer (Stulz 1982): 6.8477