Brownian simulation of correlated assets

  |   Source

When using Monte Carlo methods to price options dependent on a basket of underlying assets (multidimensional stochastic simulations), the correlations between assets should be considered. Here I will show an example of how this can be simulated using pandas.

The Geometric Brownian Motion model can be described by the stochastic differential equation:

$$dS_t = \mu S_t dt + \sigma S_t dW_t$$

In practice the solution to this equation is normally implemented as:

$$S_{t+∆t}=S_t \exp(\mu−\frac{\sigma^2}{2}) ∆t+\sigma \sqrt{\Delta t} W_t$$

So we just need to be careful how to generate the $W_t$ correctly in order to preserve the correlations between the assets.

Download and prepare the data

First we download some data from Yahoo:

In [7]:
from pandas.io.data import DataReader
from pandas import Panel, DataFrame

symbols = ['AAPL',    # Apple Inc.
           'GLD',     # SPDR Gold Trust ETF
           'SPX',     # S&P 500 Index
           'MCD']     # McDonald's Corporation
data = dict((symbol, DataReader(symbol, "yahoo", pause=1)) for symbol in symbols)
panel = Panel(data).swapaxes('items', 'minor')
closing = panel['Close'].dropna()
closing.head()
Out[7]:
AAPL GLD MCD SPX
Date
2010-01-04 214.01 109.80 62.78 1132.99
2010-01-05 214.38 109.70 62.30 1136.52
2010-01-06 210.97 111.51 61.45 1137.14
2010-01-07 210.58 110.82 61.90 1141.69
2010-01-08 211.98 111.37 61.84 1144.98

5 rows × 4 columns

Now we can calculate the log returns:

In [8]:
rets = log(closing / closing.shift(1)).dropna()
rets.head()
Out[8]:
AAPL GLD MCD SPX
Date
2010-01-05 0.001727 -0.000911 -0.007675 0.003111
2010-01-06 -0.016034 0.016365 -0.013738 0.000545
2010-01-07 -0.001850 -0.006207 0.007296 0.003993
2010-01-08 0.006626 0.004951 -0.000970 0.002878
2010-01-11 -0.008861 0.013202 0.007732 0.001745

5 rows × 4 columns

The correlation matrix has information about the historical correlations between stocks in the group. We work under the assumption that this quantity is conserved, so the generated stocks will need to satisfy this condition:

In [3]:
corr_matrix = rets.corr()
corr_matrix
Out[3]:
AAPL GLD MCD SPX
AAPL 1.000000 0.121700 0.328659 0.548297
GLD 0.121700 1.000000 0.042045 0.103287
MCD 0.328659 0.042045 1.000000 0.603910
SPX 0.548297 0.103287 0.603910 1.000000

4 rows × 4 columns

So the most correlated assets are MCD (McDonald's Corporation) and the SPX (S&P 500 Index). Pandas has a nice utility to plot the correlations:

In [6]:
from pandas.tools.plotting import scatter_matrix
import seaborn as sns

scatter_matrix(rets, figsize=(8,8));

Simulation

The simulation procedure for generating random variables will go like this:

  1. Calculate the Cholesky Decomposition matrix, this step will return an upper triangular matrix $L^T$.
  2. Generate random vector $X ∼ N(0,I)$.
  3. Obtain a correlated random vector $Z = XL^T$.

The Cholesky decomposition of the correlation matrix corr_matrix is impemented in scipy:

In [20]:
from scipy.linalg import cholesky

upper_cholesky = cholesky(corr_matrix, lower=False)
upper_cholesky
Out[20]:
array([[ 1.        ,  0.12169957,  0.3286587 ,  0.54829741],
       [ 0.        ,  0.99256698,  0.00206288,  0.03683347],
       [ 0.        ,  0.        ,  0.94444651,  0.44854975],
       [ 0.        ,  0.        ,  0.        ,  0.70485202]])

We set up the parameters for the simulation:

In [21]:
import numpy as np 
from pandas import bdate_range   # business days

n_days = 21
dates = bdate_range(start=closing.ix[-1].name, periods=n_days)
n_assets = len(symbols)
n_sims = 1000
dt = 1./252
mu = rets.mean().values
sigma = rets.std().values*sqrt(252)
np.random.seed(1234)            # init random number generator for reproducibility

Now we generate the correlated random values $X$:

In [22]:
rand_values = np.random.standard_normal(size = (n_days * n_sims, n_assets)) #
corr_values = rand_values.dot(upper_cholesky)*sigma
corr_values
Out[22]:
array([[ 0.13087785, -0.21044768,  0.22093827,  0.10789383],
       [-0.20004681,  0.14835154,  0.08464768, -0.07208435],
       [ 0.00435756, -0.41614387,  0.15946412,  0.19329902],
       ..., 
       [-0.20106783, -0.35762663, -0.185891  , -0.47797347],
       [ 0.12112218, -0.2046531 , -0.08907654, -0.11178027],
       [ 0.17662738,  0.17000015, -0.08653295,  0.06223719]])

With everything set up we can start iterating through the time interval. The results for each specific time are saved along the third axis of a pandas Panel.

In [23]:
prices = Panel(items=range(n_sims), minor_axis=symbols, major_axis=dates)
prices.ix[:, 0, :] = closing.ix[-1].values.repeat(1000).reshape(4,1000).T # set initial values

for i in range(1,n_days):
    prices.ix[:, i, :] = prices.ix[:, i-1,:] * (exp((mu-0.5*sigma**2)*dt +  sqrt(dt)*corr_values[i::n_days])).T    

prices.ix[123, :, :].head()   # show random path
Out[23]:
AAPL GLD SPX MCD
2014-01-10 532.940000 120.260000 95.800000 1842.370000
2014-01-13 545.450785 121.020824 96.720283 1862.052264
2014-01-14 532.016843 119.696718 96.785642 1845.683695
2014-01-15 542.829893 121.037367 96.631135 1857.364592
2014-01-16 550.590683 121.234086 96.694849 1841.950857

5 rows × 4 columns

And thats all! Now it is time to check our results. First a plot of all random paths for AAPL (Apple Inc.).

In [24]:
prices.ix[::10, :, 'AAPL'].plot(title='AAPL', legend=False);

We can take a look at the statistics for the last day:

In [26]:
prices.ix[:, -1, :].T.describe()
Out[26]:
AAPL GLD SPX MCD
count 1000.000000 1000.000000 1000.000000 1000.000000
mean 532.191090 120.200223 95.792754 1841.582015
std 42.769326 6.415497 4.047702 88.367691
min 418.318399 102.907420 81.802915 1538.856420
25% 502.136623 115.641946 93.065585 1780.649973
50% 529.483890 120.160960 95.678727 1837.775000
75% 560.921198 124.402412 98.176784 1900.917199
max 654.817826 140.270510 109.348568 2118.691896

8 rows × 4 columns

Finally a quick check of the correlations matrix:

In [29]:
simulated_prices = prices.ix[:, -1, :].T
predicted_rets = log(predicted_prices / predicted_prices.shift(1)).dropna()
corr_matrix = predicted_rets.corr()
corr_matrix
Out[29]:
AAPL GLD SPX MCD
AAPL 1.000000 0.155397 0.382274 0.603628
GLD 0.155397 1.000000 0.059997 0.178133
SPX 0.382274 0.059997 1.000000 0.630910
MCD 0.603628 0.178133 0.630910 1.000000

4 rows × 4 columns

Comments powered by Disqus