12 July, 2015

Gaussian Form of Geometric Brownian Motion

In numerical analysis, it can in some cases be more efficient to evaluate a function when it is expressed in a certain form.  For example, it is numerically more efficient to compute a second degree polynomial
$$P_{2}(x)=a_{2}x^{2}+a_{1}x+a_{0}$$
 when it is in the form
$$P(x)=x(a_{2}x+a_{1})+a_{0}.$$
That is because evaluating the former requires three multiplication operations and two addition operations, whereas the latter requires only two multiplication operations and two addition operations.  This way of evaluating polynomials numerically is known as Horner's Method, which describes the general formulation for a polynomial of degree $n$ and coefficients $\{a_{j}\}_{0\leq j\leq n}.$  A simple implementation of the algorithm in C can be found here.  In general, evaluation of a $n$ degree polynomial $P(x)$ will require $2n-1$ multiplication operations and $n$ addition operations using the naïve method, whereas it will require just $n$ multiplication operations and $n$ addition operations using Horner's Method.  Although from the point of view of algorithmic analysis both methods are equivalent (they both require $O(n)$ operations to compute $P(x)$), on a practical/implementation level, Horner's Method requires only $2n$ operations whereas the naïve approach requires $3n-1$ operations.  On the other hand, since the savings come from only the multiplication operations, the performance gain (unless each $a_{j}$ and $x$ are all integers) is somewhat tempered by the fact that for most computing hardware floating point multiplication is much faster than floating point addition (the opposite is generally true for integral arithmetic).  Note also that if $x^{2}$ is expressed explicitly as $x*x$, then a modern compiler will very likely optimize the naïve expression into something similar to Horner's Method; however, for high degree polynomials, expanding the powers with the multiplication operator like this is not feasible without producing very ugly, static, and difficult to maintain code - a power function of some sort is likely to be used and it is then significantly less likely the compiler will be able (or if it can, choose) to optimize the resulting code).




Now if the function in question is evaluated repeatedly, as in the case of finite difference methods and especially Monte Carlo simulations, then it is sometimes possible to achieve significant performance gains.  In financial applications, the formula for the exact solution to the SDE
$$dS_{t}=rSdt+\sigma SdW_{t},$$
where $W_{t}$ is a standard Brownian motion process, frequently needs to be evaluated repeatedly in the context Monte Carlo simulations.  The formula is given by
$$S_{t}=S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma W_{t}\right\}.$$
Frequently, certain parameters in the formula are known and can be synthesized as constants before entering into the main control loop, which could be on the order of millions or even billions of iterations.  In the case of a very simple Monte Carlo simulation for the pricing of a European vanilla call option, as discussed in this post, all the parameters except for the samples $W_{t}$ are known ($t$ is fixed in this situation).  Thus we can rewrite $S_{t}$ as
$$S_{t}=Ae^{bN},$$
where
$$A=S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t\right\},\;\;\;\;b=\sigma\sqrt{t},$$
and $N$ is a random variable with a standard normal distribution that is re-sampled at each iteration.  $A$ and $b$ are both computed before entering into the Monte Carlo simulation loop, thereby realizing a significant gain in computational efficiency when the number of iterations (trials) is large.

In other situations, different functional forms of $S_{t}$ can be useful.  A number of optimized numerical routines are available for Gaussian functions, which are functions of the form
$$G(\alpha,\beta,\gamma;x)=\alpha\exp\left\{-\frac{(x-\beta)^{2}}{2\gamma^{2}}\right\}.$$
If we match $S_{t}$ with $G(x)$, then the obvious candidate for identification with the argument $x$ is $\sigma$.  If we complete the square in the argument of the exponential function in $S_{t}$, we get
$$-\frac{t}{2}\left(\sigma^{2}-\frac{2W_{t}}{t}\sigma\right)=-\frac{t}{2}\left(\sigma-\frac{W_{t}}{t}\right)^{2}+\frac{W_{t}^{2}}{2t}.$$
Therefore,
$$S_{t}=S_{0}\exp\left\{rt+\frac{W_{t}^{2}}{2t}\right\}\exp\left\{-\frac{\left(\sigma-\frac{W_{t}}{t}\right)^{2}}{2/t}\right\}.$$
Matching the triple of parameters $(\alpha,\beta,\gamma)$ in $G(\sigma)$, we find
$$S_{t}(\sigma)=G(\sigma)$$
where
$$\alpha=S_{0}\exp\left\{rt+\frac{W_{t}^{2}}{2t}\right\}, \;\;\;\;  \beta=\frac{W_{t}}{t}, \;\;\;\; \gamma=\frac{1}{\sqrt{t}}.$$
If we map $W_{t}\mapsto\sqrt{t}{N}$ where $N$ is a standard normal random variable, then
$$\alpha=S_{0}\exp\left\{rt+\frac{N^{2}}{2}\right\}, \;\;\;\; \beta=\frac{N}{\sqrt{t}}, \;\;\;\; \gamma=\frac{1}{\sqrt{t}}.$$

This is the Gaussian form of $S_{t}$, with independent variable $\sigma$.  This representation is particularly useful for numerical pricing routines when a stochastic/local volatility model is being used.