In this tutorial we will implement the Monte Carlo Method to simulate a stock portfolio over time.

To be clear, we will be assuming daily returns are distributed by a Multivariate Normal Distribution \(R_t \sim MVN(\mu,\sum)\)

**Cholesky Decomposition** is used to determine Lower Triangular Matrix \(L \in LL’=\sum\)

\(R_t= \mu+LZ_t\)

\(Z_t \sim N(0,I)\)

Where \(Z_t\) are the samples from a normal distribution (Ι represents the Identity matrix).

Firstly let’s import our financial data from yahoo and compute the mean returns and covariance matrix. The weights of the portfolio are assigned at random and normalised to ensure the summation equals 1.

import pandas as pd import numpy as np import matplotlib.pyplot as plt import datetime as dt from pandas_datareader import data as pdr # import data def get_data(stocks, start, end): stockData = pdr.get_data_yahoo(stocks, start, end) stockData = stockData['Close'] returns = stockData.pct_change() meanReturns = returns.mean() covMatrix = returns.cov() return meanReturns, covMatrix stockList = ['CBA', 'BHP', 'TLS', 'NAB', 'WBC', 'STO'] stocks = [stock + '.AX' for stock in stockList] endDate = dt.datetime.now() startDate = endDate - dt.timedelta(days=300) meanReturns, covMatrix = get_data(stocks, startDate, endDate) weights = np.random.random(len(meanReturns)) weights /= np.sum(weights)

The monte carlo method as described above is implemented below. It is important to remember that as the number of simulations increases, the computation required to complete the method increases, you are aiming for a convergence to the ‘exact’ solution – within reason!

# Monte Carlo Method mc_sims = 400 # number of simulations T = 100 #timeframe in days meanM = np.full(shape=(T, len(weights)), fill_value=meanReturns) meanM = meanM.T portfolio_sims = np.full(shape=(T, mc_sims), fill_value=0.0) initialPortfolio = 10000 for m in range(0, mc_sims): Z = np.random.normal(size=(T, len(weights)))#uncorrelated RV's L = np.linalg.cholesky(covMatrix) #Cholesky decomposition to Lower Triangular Matrix dailyReturns = meanM + np.inner(L, Z) #Correlated daily returns for individual stocks portfolio_sims[:,m] = np.cumprod(np.inner(weights, dailyReturns.T)+1)*initialPortfolio plt.plot(portfolio_sims) plt.ylabel('Portfolio Value ($)') plt.xlabel('Days') plt.title('MC simulation of a stock portfolio') plt.show()