SDEs and Monte Carlo in Octave Financial
This post is a tutorial on my fist major contribution as maintainer of GNU Octave financial package: a framework to simulate stochastic differential equations (SDEs) of the form
\[dX_t = F(t, X_t) dt + G(t, X_t) dW_t\]where W is a standard n-dimensional Wiener process.
To follow along with the examples in this post youโll need a copy of GNU Octave and the Financial package. See the release announcement post for installation instructions.
A one dimensional example: pricing a European call
A classic problem in finance is that of pricing a European call option. The price of a European call option with strike price K and expiry time T years from now is given by
\[\mathbb{E}\left[ e^{-rT} \max\{X_T - K, 0\} \right].\]In the celebrated Black-Scholes framework, the functions F and G which parameterize the SDE are taken to be F(t, Xt) = (r - ๐ฟ) Xt and G(t, Xt) = ๐ Xt where r, ๐ฟ, and ๐ are the per-annum interest, dividend, and volatility rates of the stock X. In other words,
\[dX_t = \left(r - \delta\right) X_t dt + \sigma X_t dW_t.\]When F and G are linear functions of the state variable (as they are in this case), the SDE is called a Geometric Brownian Motion.
Approximating the above expectation using a sample mean is referred to as Monte Carlo integration or 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 two dimensional example: pricing a European basket call
Consider now the basket call pricing problem
\[\mathbb{E}\left[ e^{-rT} \max\{\max\{X_T^1, X_T^2\} - K, 0\} \right]\]involving two stocks X1 and X2 which follow the SDEs
\[dX_t^i = \left(r^i - \delta^i\right) X_t^i dt + \sigma^i X_t^i dW_t^i \qquad \text{for } i = 1,2.\]To make matters more interesting, we also assume a correlation between the coordinates of the Wiener process:
\[dW_t^1 dW_t^2 = \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