Challenge introduction

Brownian motions play a fundamental role in modeling and simulating various financial processes. Among these models, the geometric Brownian model is widely employed to describe how stock prices evolve over time. This challenge offers a valuable opportunity to gain familiarity with the properties of this model and to develop an efficient implementation. Such knowledge can have practical applications across different business domains, including derivatives pricing and stress testing.

This challenge comes from my own experience in the field of quantitative finance. I recommend attempting to solve the challenge independently before reviewing the provided solution.

Problem statement

You have been provided with the following input parameters:

• A vector of floating-point values of annualized asset returns, with a length of $\small N$:
$\small M = (\mu_0, \mu_1, ..., \mu_n)$
• A vector of floating-point values of annualized asset volatilities, also of length $\small N$:
$\small \Sigma = (\sigma_0, \sigma_1, ..., \sigma_n)$
• A square matrix of floating-point values of asset correlations, with dimensions $\small NxN$.
$\small P = \begin{bmatrix} \rho_{0, 0} & \rho_{0, 1} & ... & \rho_{0, n} \\ \rho_{1, 0} & \rho_{1, 1} & ... & \rho_{1, n} \\ ... & ... & ... & ... \\ \rho_{n, 0} & \rho_{n, 1} & ... & \rho_{n, n} \end{bmatrix}$

Your objective is to develop a class that computes $\small N$ geometric Brownian motions. Each of these motions should exhibit returns specified by the corresponding values of mu and sigma. Additionally, the correlations between these geometric Brownian motions should align with the correlation matrix, $\small P$. The class should offer flexibility to users by allowing them to configure the number of steps in the simulation and the count of independent simulations for each asset.

You are restricted from utilizing loops or iteration methods like "apply." In other words, all calculations must be fully implemented using vectorization with NumPy. The returns within vector $\small M$ are computed using discrete compounding. Therefore, it is necessary to make the requisite adjustments to transform them into the continuous equivalent required for further calculations.

Solution

The solution to this challenge begins by creating an auxiliary class that implements a basic Wiener process. This process serves as the foundational element upon which the geometric Brownian motions are constructed.

The Wiener process is simulated as a regular random walk, achieved by accumulating white noise perturbations sampled from a standard Gaussian distribution.

Geometric Brownian motion is then built on top of this standard Brownian motion using the following formula:

$\displaystyle S_t = S_0 e^{(\mu - \frac{\sigma^2}{2})t + \sigma W_t}$

where $\small W_t$ represents a Wiener process. The key to vectorizing the entire process lies in conducting Cholesky decomposition on the original correlation matrix, transforming it into a lower triangular matrix. This hierarchical approach ensures that the stochastic component of each asset simulation (the Brownian motion) is generated in a manner where each stochastic component depends only on the assets simulated thus far.

class WienerProcess:
    """
    This class implements an arbitrary number of independent Wiener processes.
    The output is a 3D numpy ndarray with the following axes:
        · n_steps: number of steps for each brownian motion.
        · n_processes: number of independent processes to simulate.
        · n_assets: number of independent simulations with n_steps and n_processes.
    """

    def __init__(self, periods: int = 252) -> None:
        self.periods = periods

    def calculate(self, n_steps: int, n_processes: int, 
                  n_assets: int = 1) -> Tuple[np.ndarray, np.ndarray]:
        shape = (n_steps, n_processes, n_assets)
        u_t = np.random.normal(loc=0, scale=1 / np.sqrt(self.periods), size=shape)
        w_t = np.concatenate((np.zeros(shape=(1, n_processes, n_assets)), 
                              u_t.cumsum(axis=0)), axis=0)
        return u_t, w_t


class CorrelatedGeometricBrownianMotion:
    """
    This class implements an arbitrary number of correlated geometric brownian motions.
    The output is a 3D numpy ndarray with the following axes:
        · n_steps: number of steps for each brownian motion.
        · n_processes: number of independent processes to simulate.
        · n_assets: number of correlated simulations with n_steps and n_processes.
    """

    def __init__(self, mu: np.ndarray, sigma: np.ndarray, rho: np.ndarray, s_0: np.ndarray, 
                 periods: int = 252) -> None:
        self.periods = periods
        self.mu, self.sigma, self.rho = np.log(1 + mu), sigma, rho
        self.s_0 = s_0

    def calculate(self, n_steps: int, n_processes: int) -> np.ndarray:
        shape = (n_steps, n_processes, len(self.mu))
        _, w_t = WienerProcess(periods=self.periods).calculate(*shape)
        rho_decomposed = np.linalg.cholesky(self.rho)
        wcorr_t = np.tensordot(w_t, rho_decomposed, axes=(-1, -1))
        t = np.arange(0, n_steps + 1)[:, np.newaxis] / self.periods
        drift = (self.mu - .5 * self.sigma ** 2) * t
        s_t = self.s_0 * np.exp(drift[:, np.newaxis] + (self.sigma * wcorr_t))
        return s_t