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
```