Efficient Frontier Portfolio Monte Carlo Simulation in Python

'''
Created on 29 Sep 2013
@author: deniz turan (denizstij AT gmail DOT com)
'''

import numpy as np
import QSTK.qstkutil.qsdateutil as du
import QSTK.qstkutil.tsutil as tsu
import QSTK.qstkutil.DataAccess as da

import datetime as dt
import matplotlib.pyplot as plt
import pandas as pd


class PortfolioEfficientFrontierMonteCarloSimulator():
    
    
    def simulate(self, dt_start, endDate, ls_symbols, weights):        
        dt_timeofday = dt.timedelta(hours=16)
        ldt_timestamps = du.getNYSEdays(dt_start, dt_end, dt_timeofday)

        c_dataobj = da.DataAccess('Yahoo')
        ls_keys = ['open', 'high', 'low', 'close', 'volume', 'actual_close']
        ldf_data = c_dataobj.get_data(ldt_timestamps, ls_symbols, ls_keys)
        d_data = dict(zip(ls_keys, ldf_data))        
        closePrices = d_data['close'].values
        
        rows = closePrices.shape[0]
        
        dailyReturns = (closePrices[0:rows, :] / closePrices[0, :])
        dailyPortfolioCum = dailyReturns * weights;
        dailyPortfolioSum = np.sum(dailyPortfolioCum, axis=1)
        dailyPortfolioRet = (dailyPortfolioSum[1:rows] / dailyPortfolioSum[0:rows - 1]) - 1
        dailyPortfolioRet = np.append(0, dailyPortfolioRet) 

        mean = np.mean(dailyPortfolioRet, axis=0)
        vol = np.std(dailyPortfolioRet, axis=0)
        cumReturn = dailyPortfolioSum[-1]
        T = 252;
        riskFree = 0
        sharpRatio = np.sqrt(T) * (mean - riskFree) / vol
        return [mean, vol, sharpRatio,cumReturn]

    def generateShortWeights(self, numWeight):        
        weights = np.zeros(numWeight)
        weights[0:numWeight - 1] = np.random.randint(-101, 101, numWeight - 1)
        weights[numWeight - 1] = 100 - np.sum(weights)
        np.random.shuffle(weights) # eliminate bias on last element
        weights = weights / np.sum(weights)    
        return weights

    def generateNoShortWeights(self, numWeight):        
        weights = np.random.rand(numWeight)
        weights = weights / np.sum(weights)
        return weights

    
    def genareteWeight(self, numWeight, isShortAllowed):
        if isShortAllowed :
            return self.generateShortWeights(numWeight)
        else:
            return self.generateNoShortWeights(numWeight)
        
                
myClass = PortfolioEfficientFrontierMonteCarloSimulator()
dt_start = dt.datetime(2011, 1, 1)
dt_end = dt.datetime(2011, 12, 31)
stockList = ['AAPL', 'GLD', 'GOOG', 'XOM']
# weight=[0.4, 0.4, 0.0, 0.2]
numTrial=10000
res= np.zeros((numTrial,4))
weights= np.zeros((numTrial,4))

for i in range(0, numTrial):
    weights[i] = myClass.genareteWeight(4, True)    
    res[i,:]=myClass.simulate(dt_start, dt_end, stockList, weights[i])
    
#  find index of min/max of mean and vol
minMeanIndex=res[:,0].argmin()
maxMeanIndex=res[:,0].argmax()
minVolIndex=res[:,1].argmin()
maxVolIndex=res[:,1].argmax()

# min and max mean and vol
maxVol=res[maxVolIndex,1]
maxMean=res[maxMeanIndex,0]
minVol=res[minVolIndex,1]
minMean=res[minMeanIndex,0]

# lets plot now 
plt.clf()
plt.scatter(res[:,1],res[:,0],marker="+", linewidths=0.5)
# Plot global mean variance portfolio
plt.scatter(res[minVolIndex,1],res[minVolIndex,0],c='m',marker='x',linewidths=3) 
plt.xlim([minVol*0.8,maxVol*1.2])
plt.ylim([minMean*0.8,maxMean*1.2])
plt.ylabel('Return')
plt.xlabel('Vol')
plt.savefig('efficientFrontier.png', format='png')

# lets print some stats now
print "Global Mean-Variance, mean=%s, vol=%s, weights: %s" %( res[minVolIndex,0],res[minVolIndex,1], weights[minVolIndex,:])
print "minMeanIndex  mean=%s, vol=%s, weights: %s" %( res[minMeanIndex,0],res[minMeanIndex,1], weights[minMeanIndex,:])
print "maxMeanIndex  mean=%s, vol=%s, weights: %s" %( res[maxMeanIndex,0],res[maxMeanIndex,1], weights[maxMeanIndex,:])
print "maxVolIndex mean=%s, vol=%s, weights: %s" %( res[maxVolIndex,0],res[maxVolIndex,1], weights[maxVolIndex,:])

Efficient Frontier Portfolio (no short allowed)