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()