Back to Community
Markowitz Portfolio Construction

This algorithm performs a standard mean-variance optimization over a group of large cap US stocks. The algorithm constructs an efficient frontier of allocations and allows the user to choose an allocation based on risk preference.

Ryan

Clone Algorithm
457
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

# This algorithm performs a Markowitz-style mean-variance optimization with a 
# universe of the 14 highest market capitalization stocks in the S&P 500. The
# algorithm constructs an efficient frontier and then selects a weighting from
# the efficient frontier. Uses ideas from the Global Minimum Variance Portfolio 
# algorithm posted on Quantopian.

# This is the lwindow that governs the pulling of historical data and
# portfolio rebalancing. 
window = 100
refresh_rate = 10

# Compute the expected return of the portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
        return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = 0.1*abs(mean-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean, var = compute_mean_var(W, R, C)
    utility = (mean - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(min(R), max(R), num=20): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data
def assets_meanvar(daily_returns):    
    
    # Calculate expected returns
    expreturns = array([])
    (rows, cols) = daily_returns.shape
    for r in range(rows):
        expreturns = append(expreturns, mean(daily_returns[r]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    
    return expreturns, covars
    
def initialize(context):
  
  # Set day
  context.day = 0
    
  # Set risk-free rate
 
  # 14 largest stocks by market cap in S&P in 2010
  context.securities = [sid(24),sid(26578),sid(3149), sid(5061), sid(3766),sid(23112),sid(4151),sid(5938),sid(5923),sid(6653),sid(700), sid(1900), sid(8347), sid(8229)]
  # Set Max and Min positions in security
  context.max_notional = 1000000.1
  context.min_notional = -1000000.0
  # Set commission

def handle_data(context, data):

  # Get 40 days of prices for each security
  all_prices = get_past_prices(data)
    
  # Circuit breaker in case transform returns none
  if all_prices is None:
        return
  # Circuit breaker, only calculate every 20 days
  if context.day%refresh_rate is not 0:
        context.day = context.day+1
        return
  
  daily_returns = np.zeros((len(context.securities),window))
    
  # Calculate daily returns into daily_returns
  security_index = 0;
  for security in context.securities:
      if data.has_key(security):
          for day in range(0,window):
              day_of = all_prices[security][day]
              day_before = all_prices[security][day-1]
              daily_returns[security_index][day] = (day_of-day_before)/day_before
          security_index = security_index + 1  
    
  expreturns, covars = assets_meanvar(daily_returns)
  R = expreturns
  C = covars
  rf = 0.015
  frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf) 
  
  f_w = array(frontier_weights)          
  (row_1, col_1) = f_w.shape
  log.info(row_1)
  log.info(col_1)          
  
  wts = frontier_weights[5]
  new_weights = wts            
  #log.info(new_weights)          
            
  #set leverage to 1
  leverage = sum(abs(new_weights))
  portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
    
  #reweight portfolio 
  security_index = 0
  for security in context.securities:
        current_position = context.portfolio.positions[security].amount
        new_position = (portfolio_value*new_weights[security_index])/all_prices[security][window-1]
        order(security,new_position-current_position)
        security_index = security_index+1   
  context.day = context.day+1

@batch_transform(refresh_period=refresh_rate, window_length=window)  
def get_past_prices(data):  
    prices = data['price']    
    return prices
This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.
Disclaimer

The material on this website is provided for informational purposes only and does not constitute an offer to sell, a solicitation to buy, or a recommendation or endorsement for any security or strategy, nor does it constitute an offer to provide investment advisory services by Quantopian. In addition, the material offers no opinion with respect to the suitability of any security or specific investment. No information contained herein should be regarded as a suggestion to engage in or refrain from any investment-related course of action as none of Quantopian nor any of its affiliates is undertaking to provide investment advice, act as an adviser to any plan or entity subject to the Employee Retirement Income Security Act of 1974, as amended, individual retirement account or individual retirement annuity, or give advice in a fiduciary capacity with respect to the materials presented herein. If you are an individual retirement or other investor, contact your financial advisor or other fiduciary unrelated to Quantopian about whether any given investment idea, strategy, product or service described herein may be appropriate for your circumstances. All investments involve risk, including loss of principal. Quantopian makes no guarantees as to the accuracy or completeness of the views expressed in the website. The views are subject to change, and may have become unreliable for various reasons, including changes in market conditions or economic circumstances.

20 responses

http://en.wikipedia.org/wiki/Modern_portfolio_theory

Modern portfolio theory (MPT) is a theory of finance that attempts to maximize portfolio expected return for a given amount of portfolio risk, or equivalently minimize risk for a given level of expected return, by carefully choosing the proportions of various assets. MPT is a mathematical formulation of the concept of diversification in investing, with the aim of selecting a collection of investment assets that has collectively lower risk than any individual asset. This is possible, intuitively speaking, because different types of assets often change in value in opposite ways. For example, to the extent prices in the stock market move differently from prices in the bond market, a collection of both types of assets can in theory face lower overall risk than either individually. But diversification lowers risk even if assets' returns are not negatively correlated—indeed, even if they are positively correlated.

MPT assumes that investors are risk averse, meaning that given two portfolios that offer the same expected return, investors will prefer the less risky one. Thus, an investor will take on increased risk only if compensated by higher expected returns. Conversely, an investor who wants higher expected returns must accept more risk. The exact trade-off will be the same for all investors, but different investors will evaluate the trade-off differently based on individual risk aversion characteristics. The implication is that a rational investor will not invest in a portfolio if a second portfolio exists with a more favorable risk-expected return profile – i.e., if for that level of risk an alternative portfolio exists that has better expected returns.

See the attached backtest for a few updates. The algo is designed so that anyone can input any group of stocks he or she likes to be put into the allocation engine. Within the initialize method, the field risk_tolerance is added to the context data frame. One can choose a risk/return profile on a scale of 1 to 20 (1 being the lowest risk and 20 being the highest risk). Higher risk portfolios will have a higher expected return.

Ryan

Clone Algorithm
89
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

# This algorithm performs a Markowitz-style mean-variance optimization with a 
# universe of the 14 highest market capitalization stocks in the S&P 500. The
# algorithm constructs an efficient frontier and then selects a weighting from
# the efficient frontier. Uses ideas from the Global Minimum Variance Portfolio 
# algorithm posted on Quantopian. Ideas also adopted from 
# http://www.quantandfinancial.com/2013/07/mean-variance-portfolio-optimization.html

# window gives the number of data points that are extracted from historical data 
# to estimate the expected returns and covariance matrix. 
window = 100
refresh_rate = 10

# Compute the expected return of the portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
    return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean_1, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = 0.1*abs(mean_1-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean_1, var = compute_mean_var(W, R, C)
    utility = (mean_1 - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(min(R), max(R), num=50): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data
def assets_meanvar(daily_returns):    
    
    # Calculate expected returns
    expreturns = array([])
    (rows, cols) = daily_returns.shape
    for r in range(rows):
        expreturns = append(expreturns, mean(daily_returns[r]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    
    return expreturns, covars
    
def initialize(context):
  
  # Set day
  context.day = 0
    
  # Enter the risk tolerance on a scale of 1 to 20.
  # 1 is the lowest risk tolerance (lowest risk portfolio)
  # 20 is the highest risk tolerance (highest risk portfolio)
  context.risk_tolerance = 10
     
  # Select equities for the portfolio
  # The default selection of equities is the fourteen highest market cap stocks in the S&P 500
  context.securities = [sid(24),sid(26578),sid(3149), sid(5061), sid(3766),sid(23112),sid(4151),sid(5938),sid(5923),sid(6653),sid(700), sid(1900), sid(8347), sid(8229)]
  # Set Max and Min positions in security
  context.max_notional = 1000000.1
  context.min_notional = -1000000.0

def handle_data(context, data):

  # Get 40 days of prices for each security
  all_prices = get_past_prices(data)
    
  # Circuit breaker in case transform returns none
  if all_prices is None:
        return
  # Circuit breaker, only calculate every 20 days
  if context.day%refresh_rate is not 0:
        context.day = context.day+1
        return
  
  daily_returns = np.zeros((len(context.securities),window))
    
  # Calculate daily returns into daily_returns
  security_index = 0;
  for security in context.securities:
      if data.has_key(security):
          for day in range(0,window):
              day_of = all_prices[security][day]
              day_before = all_prices[security][day-1]
              daily_returns[security_index][day] = (day_of-day_before)/day_before
          security_index = security_index + 1  
    
  expreturns, covars = assets_meanvar(daily_returns)
  R = expreturns
  C = covars
  rf = 0.015
  frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf) 
  
  new_weights = frontier_weights[context.risk_tolerance]            
            
  #set leverage to 1
  leverage = sum(abs(new_weights))
  portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
    
  #reweight portfolio 
  security_index = 0
  for security in context.securities:
        current_position = context.portfolio.positions[security].amount
        new_position = (portfolio_value*new_weights[security_index])/all_prices[security][window-1]
        order(security,new_position-current_position)
        security_index = security_index+1   
  context.day = context.day+1

@batch_transform(refresh_period=refresh_rate, window_length=window)  
def get_past_prices(data):  
    prices = data['price']    
    return prices
This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.

Please see the attached improvement of the algorithm. Now includes history() instead of batch_transform() and re-balances every month. Comments welcome.

Clone Algorithm
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
from pytz import timezone  
from datetime import datetime, timedelta  
from zipline.utils.tradingcalendar import get_early_closes
from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

# Compute expected return on portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
    return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean_1, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = 50**abs(mean_1-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean_1, var = compute_mean_var(W, R, C)
    utility = (mean_1 - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(min(R), max(R), num=50): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data

def assets_meanvar(daily_returns, context):    
    
    # Calculate expected returns
    expreturns = array([])
    daily_returns = daily_returns.transpose()
    
    for i in range(0, len(context.securities)):
        expreturns = append(expreturns, mean(daily_returns[i,:]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    expreturns = np.array(expreturns)
    
    return expreturns, covars


def initialize(context):
    # Set day
    context.day = 0
    
    # Enter the risk tolerance on a scale of 1 to 20.
    # 1 is the lowest risk tolerance (lowest risk portfolio)
    # 20 is the highest risk tolerance (highest risk portfolio)
    context.risk_tolerance = 10
     
    # Select equities for the portfolio
    # The default selection of equities is the fourteen highest market cap stocks in the S&P 500
    context.securities = [sid(24),sid(26578),sid(3149), sid(5061), sid(3766),sid(23112),sid(4151),sid(5938),sid(5923),sid(6653),sid(700), sid(1900), sid(8347), sid(8229)]
    # Set Max and Min positions in security
    context.max_notional = 1000000.1
    context.min_notional = -1000000.0
    context.previous_month = 0


def handle_data(context, data):    
    
    loc_dt = get_datetime().astimezone(timezone('US/Eastern'))
    date = get_datetime().date()
    
    if loc_dt.month != context.previous_month and context.previous_month != 0:
        all_prices = history(200, '1d', 'price')
        daily_returns = all_prices.pct_change().dropna()
                
        dr = np.array(daily_returns)
        
        (rr,cc) = dr.shape
        
        expreturns, covars = assets_meanvar(dr, context)
        R = expreturns
        C = covars
        rf = 0.015
        expreturns = np.array(expreturns)
        
        frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf)
        
        f_w = array(frontier_weights)          
        (row_1, col_1) = f_w.shape         
  
        # Choose an allocation along the efficient frontier
        wts = frontier_weights[10]
        new_weights = wts  
            
        # Set leverage to 1
        leverage = sum(abs(new_weights))
        portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
    
        # Reweight portfolio 
        i = 0
        for sec in context.securities:
            order_target_percent(sec, wts[i])
            log.info(get_datetime())
            i=i+1
        
    context.previous_month = loc_dt.month
    
    
This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.

I was playing around with this and made a few changes and some notes. Interesting.

Clone Algorithm
18
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
from pytz import timezone
from datetime import datetime, timedelta  
from zipline.utils.tradingcalendar import get_early_closes
from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

def initialize(context):
    # INPUTS:
    # 1) Enter the risk tolerance on a scale of 1 to 20.
    #    1 is the lowest risk tolerance (lowest risk portfolio)
    #    20 is the highest risk tolerance (highest risk portfolio)
    # 2) Stock listing (manually or from set_universe)
    #    NOTE: To meet the 50 second requirement the number of 
    #    securities must be not more than 120 or about 1.5 % of the 
    #    universe.
    # 3) Do or do not use margin
    #    1 is no margin
    #    2 is double cash
    # 4) Frequency to rebalance

    # Configuration parameters
    context.risk_tolerance = 15
    set_universe(universe.DollarVolumeUniverse(99.5, 100))
    context.allowableMargin        = 2.0  # 
    context.requiredMargin         = 0.0  # In Dollars
    context.usedMargin             = 0.0  # In Dollars
    context.lastPortfolioUpdate    = datetime(1989,1,30,9,15,0,tzinfo=timezone('US/Eastern')) # randomly chosen init value
    context.trailingStopPortfolioValue  = 18000 # Based on starting cash of $20,000
    context.disableMarginPortfolioValue = 19000 # Based on starting cash of $20,000
    context.enableMarginPortfolioValue  = 22000 # Based on starting cash of $20,000
    context.rebalanceFrequency = 1        # Number of weeks
    
def handle_data(context, data):    
    # Morning Margin Check (UTC timezone) - LONG ONLY
    if get_datetime().hour == 14 and get_datetime().minute == 35:
        context.requiredMargin = marginRequirements(context.portfolio) * 3.0
        if context.portfolio.cash < 0.:
            context.usedMargin = abs(context.portfolio.cash)
        else:
            context.usedMargin = 0.
        if context.requiredMargin < context.usedMargin:
                log.warn('MARGIN REQUIREMENTS EXCEEDED. ' +\
                         'Used Margin = ' + str(context.usedMargin) +\
                         ' Allowable Margin = ' + str(context.requiredMargin))
        # Liquidate if total value falls 10% or more (disable margin use after 5% loss)
        if 0.9*(context.portfolio.positions_value+context.portfolio.cash) > context.trailingStopPortfolioValue:
            context.trailingStopPortfolioValue  = 0.90*(context.portfolio.positions_value+context.portfolio.cash)
            context.disableMarginPortfolioValue = 0.95*(context.portfolio.positions_value+context.portfolio.cash)
        if (context.portfolio.positions_value+context.portfolio.cash) < context.trailingStopPortfolioValue:
            log.warn('*** L I Q U I D A T E ***')
            liquidate(context.portfolio)
            context.trailingStopPortfolioValue  = 0.90*(context.portfolio.positions_value+context.portfolio.cash)
        if (context.portfolio.positions_value+context.portfolio.cash) < context.disableMarginPortfolioValue:
            log.info('*** MARGIN USE DISABLED ***')
            context.allowableMargin = 1.
            context.enableMarginPortfolioValue = 1.10*(context.portfolio.positions_value+context.portfolio.cash)
        elif (context.portfolio.positions_value+context.portfolio.cash) > context.enableMarginPortfolioValue:
            log.info('*** MARGIN USE ENABLED ***')
            context.allowableMargin = 2.
    
    # End of Day
    if get_datetime().hour == 20 and get_datetime().minute == 55:
        for stock in data.keys():
            closeAnyOpenOrders(stock)
    
    #if loc_dt.month != context.previous_month:
    if (get_datetime() - context.lastPortfolioUpdate) >= timedelta(weeks=context.rebalanceFrequency):
        context.lastPortfolioUpdate = get_datetime()
        log.debug('Number of secruities to be considered: ' + str(len(data.keys())))
        all_prices = history(200, '1d', 'price')
        daily_returns = all_prices.pct_change().dropna()
                
        dr = np.array(daily_returns)
        (rr,cc) = dr.shape
        
        expreturns, covars = assets_meanvar(dr, data.keys())
        R = expreturns
        C = covars
        rf = 0.015
        expreturns = np.array(expreturns)
        
        frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf,context)
        
        f_w = array(frontier_weights)          
        (row_1, col_1) = f_w.shape         
  
        # Choose an allocation along the efficient frontier
        wts = frontier_weights[context.risk_tolerance]
        new_weights = wts  
            
        # Set leverage to 1
        leverage = sum(abs(new_weights))
        portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
        record(PV=portfolio_value)
        record(Cash=context.portfolio.cash)
    
        # Reweight portfolio 
        i = 0
        for sec in data.keys():
            if wts[i] < 0.01:
                wts[i] = 0.0
            if wts[i] < 0.01 and context.portfolio.positions[sec].amount == 0:
                i = i+1
                continue
            order_target_percent(sec, wts[i],None,None)
            log.info('Adjusting ' + str(sec) + ' to ' + str(wts[i]*100.0) + '%')
            i=i+1

#################################################################################
#
# HELPER FUNCTIONS ( This algorithm's math functions)
#
#################################################################################                     
# Compute expected return on portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
    return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean_1, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = 50**abs(mean_1-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean_1, var = compute_mean_var(W, R, C)
    utility = (mean_1 - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf,context):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf,context):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(min(R), max(R), num=50): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data

def assets_meanvar(daily_returns, stockList):    
    
    # Calculate expected returns
    expreturns = array([])
    daily_returns = daily_returns.transpose()
    
    for i in range(0, len(stockList)):
        expreturns = append(expreturns, mean(daily_returns[i,:]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    expreturns = np.array(expreturns)
    
    return expreturns, covars
            
#################################################################################
#
# HELPER FUNCTIONS ( Portfolio / Stock Management)
#
#################################################################################         
def closeAnyOpenOrders(stock):
    orders = get_open_orders(stock)
    if orders:
        for order in orders:
             message = 'Canceling order for {amount} shares in {stock}'  
             message = message.format(amount=order.amount, stock=stock)  
             #log.debug(message)
             cancel_order(order)
def hasOrders(stock):
    hasOrders = False
    orders = get_open_orders(stock)
    if orders:
        for order in orders:
             message = 'Open order for {amount} shares in {stock}'  
             message = message.format(amount=order.amount, stock=stock)  
             log.info(message)
             hasOrders = True
    return hasOrders   
def liquidate(portfolio):
    for stock in portfolio.positions:
        closeAnyOpenOrders(stock)
        if portfolio.positions[stock].amount > 0:
         order(stock,-portfolio.positions[stock].amount,None,None)
         log.info('Sold stake in ' + str(stock) +
                  ' @ $' + str(portfolio.positions[stock].last_sale_price) +
                  ' Cost Basis: ' + str(portfolio.positions[stock].cost_basis) )
def marginRequirements(portfolio):  
    req = 0  
    for stock in portfolio.positions:  
        amount = portfolio.positions[stock].amount  
        last_price = portfolio.positions[stock].last_sale_price  
        if amount > 0:  
            req += .25 * amount * last_price  
        elif amount < 0:  
            if last_price < 5:  
                req += max(2.5 * amount, abs(amount * last_price))  
            else:  
                req += max(5 * amount, abs(0.3 * amount * last_price))  
    return req
   
    
This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.

It seems to work much better when given a good list of stocks/ETFs than using the set_universe.

From what I gather you are implementing an efficient frontier/MPT algo. I don't think that the optimization process is necessary, there's closed form solutions to most MPT problems. Here is a reference. If nothing else, it would speed up the calculations.

Brian, thank you for your improvements to the algorithm. I noticed you generated 50 different points along the efficient frontier

for r in linspace(min(R), max(R), num=50):  

while the documentation tells the user to select a risk profile using one of 20 points along the efficient frontier:

# 1) Enter the risk tolerance on a scale of 1 to 20.  
#    1 is the lowest risk tolerance (lowest risk portfolio)  
#    20 is the highest risk tolerance (highest risk portfolio)  

I modified the code to generate only twenty points along the efficient frontier.

David, yes, in most formulations, the portfolio allocation problem of MPT does have a closed form solution. However, I would like to allow for maximum flexibility when Quantopian community members use the algorithm. The objective function could include a term for tracking error, for example, which may necessitate numerical methods for the problem.

In the fitness function for the optimization, the penalty function tends to dominate the optimization. Estimating returns can often be far more difficult than estimating volatilities and correlation, so it may be more effective to more heavily weight the covariance term in relation to the estimated return term.

def fitness(W, R, C, r):  
    # For given level of return r, find weights which minimizes portfolio variance.  
    mean_1, var = compute_mean_var(W, R, C)  
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint  
    # Here, r is the 'target' return  
    penalty = 50**abs(mean_1-r)  
    return var + penalty  

Instead of using exponentiation (**), we can use multiplication.

def fitness(W, R, C, r):  
    # For given level of return r, find weights which minimizes portfolio variance.  
    mean_1, var = compute_mean_var(W, R, C)  
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint  
    # Here, r is the 'target' return  
    penalty = (1/100)*abs(mean_1-r)  
    return var + penalty  
Clone Algorithm
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
from pytz import timezone
from datetime import datetime, timedelta  
from zipline.utils.tradingcalendar import get_early_closes
from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

def initialize(context):
    # INPUTS:
    # 1) Enter the risk tolerance on a scale of 1 to 20.
    #    1 is the lowest risk tolerance (lowest risk portfolio)
    #    20 is the highest risk tolerance (highest risk portfolio)
    # 2) Stock listing (manually or from set_universe)
    #    NOTE: To meet the 50 second requirement the number of 
    #    securities must be not more than 120 or about 1.5 % of the 
    #    universe.
    # 3) Do or do not use margin
    #    1 is no margin
    #    2 is double cash
    # 4) Frequency to rebalance

    # Configuration parameters
    context.risk_tolerance = 15
    set_universe(universe.DollarVolumeUniverse(99.0, 100))
    context.allowableMargin        = 2.0  # 
    context.requiredMargin         = 0.0  # In Dollars
    context.usedMargin             = 0.0  # In Dollars
    context.lastPortfolioUpdate    = datetime(1989,1,30,9,15,0,tzinfo=timezone('US/Eastern')) # randomly chosen init value
    context.trailingStopPortfolioValue  = 18000 # Based on starting cash of $20,000
    context.disableMarginPortfolioValue = 19000 # Based on starting cash of $20,000
    context.enableMarginPortfolioValue  = 22000 # Based on starting cash of $20,000
    context.rebalanceFrequency = 1        # Number of weeks
    
def handle_data(context, data):    
    # Morning Margin Check (UTC timezone) - LONG ONLY
    if get_datetime().hour == 14 and get_datetime().minute == 35:
        context.requiredMargin = marginRequirements(context.portfolio) * 3.0
        if context.portfolio.cash < 0.:
            context.usedMargin = abs(context.portfolio.cash)
        else:
            context.usedMargin = 0.
        if context.requiredMargin < context.usedMargin:
                log.warn('MARGIN REQUIREMENTS EXCEEDED. ' +\
                         'Used Margin = ' + str(context.usedMargin) +\
                         ' Allowable Margin = ' + str(context.requiredMargin))
        # Liquidate if total value falls 10% or more (disable margin use after 5% loss)
        if 0.9*(context.portfolio.positions_value+context.portfolio.cash) > context.trailingStopPortfolioValue:
            context.trailingStopPortfolioValue  = 0.90*(context.portfolio.positions_value+context.portfolio.cash)
            context.disableMarginPortfolioValue = 0.95*(context.portfolio.positions_value+context.portfolio.cash)
        if (context.portfolio.positions_value+context.portfolio.cash) < context.trailingStopPortfolioValue:
            log.warn('*** L I Q U I D A T E ***')
            liquidate(context.portfolio)
            context.trailingStopPortfolioValue  = 0.90*(context.portfolio.positions_value+context.portfolio.cash)
        if (context.portfolio.positions_value+context.portfolio.cash) < context.disableMarginPortfolioValue:
            log.info('*** MARGIN USE DISABLED ***')
            context.allowableMargin = 1.
            context.enableMarginPortfolioValue = 1.10*(context.portfolio.positions_value+context.portfolio.cash)
        elif (context.portfolio.positions_value+context.portfolio.cash) > context.enableMarginPortfolioValue:
            log.info('*** MARGIN USE ENABLED ***')
            context.allowableMargin = 2.
    
    # End of Day
    if get_datetime().hour == 20 and get_datetime().minute == 55:
        for stock in data.keys():
            closeAnyOpenOrders(stock)
    
    #if loc_dt.month != context.previous_month:
    if (get_datetime() - context.lastPortfolioUpdate) >= timedelta(weeks=context.rebalanceFrequency):
        context.lastPortfolioUpdate = get_datetime()
        log.debug('Number of secruities to be considered: ' + str(len(data.keys())))
        all_prices = history(250, '1d', 'price')
        daily_returns = all_prices.pct_change().dropna()
                
        dr = np.array(daily_returns)
        (rr,cc) = dr.shape
        
        expreturns, covars = assets_meanvar(dr, data.keys())
        R = expreturns
        C = covars
        rf = 0.015
        expreturns = np.array(expreturns)
        
        frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf,context)
        
        f_w = array(frontier_weights)          
        (row_1, col_1) = f_w.shape         
  
        # Choose an allocation along the efficient frontier
        wts = frontier_weights[context.risk_tolerance]
        new_weights = wts  
            
        # Set leverage to 1
        leverage = sum(abs(new_weights))
        portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
        record(PV=portfolio_value)
        record(Cash=context.portfolio.cash)
    
        # Reweight portfolio 
        i = 0
        for sec in data.keys():
            if wts[i] < 0.01:
                wts[i] = 0.0
            if wts[i] < 0.01 and context.portfolio.positions[sec].amount == 0:
                i = i+1
                continue
            order_target_percent(sec, wts[i],None,None)
            log.info('Adjusting ' + str(sec) + ' to ' + str(wts[i]*100.0) + '%')
            i=i+1

#################################################################################
#
# HELPER FUNCTIONS ( This algorithm's math functions)
#
#################################################################################                     

# Compute expected return on portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
    return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean_1, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = (1/100)*abs(mean_1-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean_1, var = compute_mean_var(W, R, C)
    utility = (mean_1 - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf,context):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf,context):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(max(min(R), rf), max(R), num=20): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)s
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data

def assets_meanvar(daily_returns, stockList):    
    
    # Calculate expected returns
    expreturns = array([])
    daily_returns = daily_returns.transpose()
    
    for i in range(0, len(stockList)):
        expreturns = append(expreturns, mean(daily_returns[i,:]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    expreturns = np.array(expreturns)
    
    return expreturns, covars
            
#################################################################################
#
# HELPER FUNCTIONS ( Portfolio / Stock Management)
#
#################################################################################         

def closeAnyOpenOrders(stock):
    orders = get_open_orders(stock)
    if orders:
        for order in orders:
             message = 'Canceling order for {amount} shares in {stock}'  
             message = message.format(amount=order.amount, stock=stock)  
             #log.debug(message)
             cancel_order(order)
def hasOrders(stock):
    hasOrders = False
    orders = get_open_orders(stock)
    if orders:
        for order in orders:
             message = 'Open order for {amount} shares in {stock}'  
             message = message.format(amount=order.amount, stock=stock)  
             log.info(message)
             hasOrders = True
    return hasOrders   
def liquidate(portfolio):
    for stock in portfolio.positions:
        closeAnyOpenOrders(stock)
        if portfolio.positions[stock].amount > 0:
         order(stock,-portfolio.positions[stock].amount,None,None)
         log.info('Sold stake in ' + str(stock) +
                  ' @ $' + str(portfolio.positions[stock].last_sale_price) +
                  ' Cost Basis: ' + str(portfolio.positions[stock].cost_basis) )
def marginRequirements(portfolio):  
    req = 0  
    for stock in portfolio.positions:  
        amount = portfolio.positions[stock].amount  
        last_price = portfolio.positions[stock].last_sale_price  
        if amount > 0:  
            req += .25 * amount * last_price  
        elif amount < 0:  
            if last_price < 5:  
                req += max(2.5 * amount, abs(amount * last_price))  
            else:  
                req += max(5 * amount, abs(0.3 * amount * last_price))  
    return req
   
    
This backtest was created using an older version of the backtester. Please re-run this backtest to see results using the latest backtester. Learn more about the recent changes.
There was a runtime error.

That is an error. I was playing around. I've changed it to 10. That provides enough resolution in my opinion.

Here is original Ryan Davis version on sector etf products from state street plus TLT

Clone Algorithm
52
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
# Markowitz Portfolio Construction
# https://www.quantopian.com/posts/markowitz-portfolio-construction

from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

# This algorithm performs a Markowitz-style mean-variance optimization with a 
# universe of sector spyder etf products from state street plus TLT. The
# algorithm constructs an efficient frontier and then selects a weighting from
# the efficient frontier. Uses ideas from the Global Minimum Variance Portfolio 
# algorithm posted on Quantopian.

# This is the lwindow that governs the pulling of historical data and
# portfolio rebalancing. 
window = 100
refresh_rate = 10

def initialize(context):  
  # Set day
  context.day = 0    
    
  # Set risk-free rate
 
  # based on the sector spyder etf products from state street plus TLT
  context.securities =symbols( 'tlt','xlf', 'xle', 'xlu', 'xlk', 'xlb', 'xlp', 'xly','xli', 'xlv')
  # Set Max and Min positions in security
  context.max_notional = 1000000.1
  context.min_notional = -1000000.0
  # Set commission
    
@batch_transform(refresh_period=refresh_rate, window_length=window)  
def get_past_prices(data):  
    prices = data['price']    
    return prices    

# Compute the expected return of the portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
        return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = 0.1*abs(mean-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean, var = compute_mean_var(W, R, C)
    utility = (mean - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(min(R), max(R), num=20): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1. }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data
def assets_meanvar(daily_returns):    
    
    # Calculate expected returns
    expreturns = array([])
    (rows, cols) = daily_returns.shape
    for r in range(rows):
        expreturns = append(expreturns, mean(daily_returns[r]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    
    return expreturns, covars 

def handle_data(context, data):

  # Get 40 days of prices for each security
  all_prices = get_past_prices(data)
    
  # Circuit breaker in case transform returns none
  if all_prices is None:
        return
  # Circuit breaker, only calculate every 20 days
  if context.day%refresh_rate is not 0:
        context.day = context.day+1
        return
  
  daily_returns = np.zeros((len(context.securities),window))
    
  # Calculate daily returns into daily_returns
  security_index = 0;
  for security in context.securities:
      if data.has_key(security):
          for day in range(0,window):
              day_of = all_prices[security][day]
              day_before = all_prices[security][day-1]
              daily_returns[security_index][day] = (day_of-day_before)/day_before
          security_index = security_index + 1  
    
  expreturns, covars = assets_meanvar(daily_returns)
  R = expreturns
  C = covars
  rf = 0.015
  frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf) 
  
  f_w = array(frontier_weights)          
  (row_1, col_1) = f_w.shape
  log.info(row_1)
  log.info(col_1)          
  
  wts = frontier_weights[5]
  new_weights = wts            
  #log.info(new_weights)          
            
  #set leverage to 1
  leverage = sum(abs(new_weights))
  portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
    
  #reweight portfolio 
  security_index = 0
  for security in context.securities:
        current_position = context.portfolio.positions[security].amount
        new_position = (portfolio_value*new_weights[security_index])/all_prices[security][window-1]
        order(security,new_position-current_position)
        security_index = security_index+1   
  context.day = context.day+1


There was a runtime error.

Brian Vetere version on sector etf products from state street plus TLT

Clone Algorithm
130
Loading...
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
Returns 1 Month 3 Month 6 Month 12 Month
Alpha 1 Month 3 Month 6 Month 12 Month
Beta 1 Month 3 Month 6 Month 12 Month
Sharpe 1 Month 3 Month 6 Month 12 Month
Sortino 1 Month 3 Month 6 Month 12 Month
Volatility 1 Month 3 Month 6 Month 12 Month
Max Drawdown 1 Month 3 Month 6 Month 12 Month
# Markowitz Portfolio Construction Reply ???


from pytz import timezone
from datetime import datetime, timedelta  
from zipline.utils.tradingcalendar import get_early_closes
from numpy import matrix, array, zeros, empty, sqrt, ones, dot, append, mean, cov, transpose, linspace
import numpy as np
import scipy.optimize
import math
import random

def initialize(context):
    # INPUTS:
    # 1) Enter the risk tolerance on a scale of 1 to 20.
    #    1 is the lowest risk tolerance (lowest risk portfolio)
    #    20 is the highest risk tolerance (highest risk portfolio)
    # 2) Stock listing (manually or from set_universe)
    #    NOTE: To meet the 50 second requirement the number of 
    #    securities must be not more than 120 or about 1.5 % of the 
    #    universe.
    # 3) Do or do not use margin
    #    1 is no margin
    #    2 is double cash
    # 4) Frequency to rebalance

    # Configuration parameters
    context.risk_tolerance = 15
    #set_universe(universe.DollarVolumeUniverse(99.0, 100))
    context.securities =symbols( 'tlt','xlf', 'xle', 'xlu', 'xlk', 'xlb', 'xlp', 'xly','xli', 'xlv')
    context.allowableMargin        = 2.0  # 
    context.requiredMargin         = 0.0  # In Dollars
    context.usedMargin             = 0.0  # In Dollars
    context.lastPortfolioUpdate    = datetime(1989,1,30,9,15,0,tzinfo=timezone('US/Eastern')) # randomly chosen init value
    context.trailingStopPortfolioValue  = 18000 # Based on starting cash of $20,000
    context.disableMarginPortfolioValue = 19000 # Based on starting cash of $20,000
    context.enableMarginPortfolioValue  = 22000 # Based on starting cash of $20,000
    context.rebalanceFrequency = 1        # Number of weeks
    
def handle_data(context, data):    
    # Morning Margin Check (UTC timezone) - LONG ONLY
    if get_datetime().hour == 14 and get_datetime().minute == 35:
        context.requiredMargin = marginRequirements(context.portfolio) * 3.0
        if context.portfolio.cash < 0.:
            context.usedMargin = abs(context.portfolio.cash)
        else:
            context.usedMargin = 0.
        if context.requiredMargin < context.usedMargin:
                log.warn('MARGIN REQUIREMENTS EXCEEDED. ' +\
                         'Used Margin = ' + str(context.usedMargin) +\
                         ' Allowable Margin = ' + str(context.requiredMargin))
        # Liquidate if total value falls 10% or more (disable margin use after 5% loss)
        if 0.9*(context.portfolio.positions_value+context.portfolio.cash) > context.trailingStopPortfolioValue:
            context.trailingStopPortfolioValue  = 0.90*(context.portfolio.positions_value+context.portfolio.cash)
            context.disableMarginPortfolioValue = 0.95*(context.portfolio.positions_value+context.portfolio.cash)
        if (context.portfolio.positions_value+context.portfolio.cash) < context.trailingStopPortfolioValue:
            log.warn('*** L I Q U I D A T E ***')
            liquidate(context.portfolio)
            context.trailingStopPortfolioValue  = 0.90*(context.portfolio.positions_value+context.portfolio.cash)
        if (context.portfolio.positions_value+context.portfolio.cash) < context.disableMarginPortfolioValue:
            log.info('*** MARGIN USE DISABLED ***')
            context.allowableMargin = 1.
            context.enableMarginPortfolioValue = 1.10*(context.portfolio.positions_value+context.portfolio.cash)
        elif (context.portfolio.positions_value+context.portfolio.cash) > context.enableMarginPortfolioValue:
            log.info('*** MARGIN USE ENABLED ***')
            context.allowableMargin = 2.
    
    # End of Day
    if get_datetime().hour == 20 and get_datetime().minute == 55:
        for stock in data.keys():
            closeAnyOpenOrders(stock)
    
    #if loc_dt.month != context.previous_month:
###    if (get_datetime() - context.lastPortfolioUpdate) >= timedelta(weeks=context.rebalanceFrequency):
###        context.lastPortfolioUpdate = get_datetime()
###        log.debug('Number of secruities to be considered: ' + str(len(data.keys())))
    if get_datetime().hour == 14 and get_datetime().minute == 35:
        all_prices = history(250, '1d', 'price')
        daily_returns = all_prices.pct_change().dropna()
                
        dr = np.array(daily_returns)
        (rr,cc) = dr.shape
        
        expreturns, covars = assets_meanvar(dr, data.keys())
        R = expreturns
        C = covars
        rf = 0.015
        expreturns = np.array(expreturns)
        
        frontier_mean, frontier_var, frontier_weights = solve_frontier(R, C, rf,context)
        
        f_w = array(frontier_weights)          
        (row_1, col_1) = f_w.shape         
  
        # Choose an allocation along the efficient frontier
        wts = frontier_weights[context.risk_tolerance]
        new_weights = wts  
            
        # Set leverage to 1
        leverage = sum(abs(new_weights))
        portfolio_value = (context.portfolio.positions_value + context.portfolio.cash)/leverage
        record(PV=portfolio_value)
        record(Cash=context.portfolio.cash)
    
        # Reweight portfolio 
        i = 0
        for sec in data.keys():
            if wts[i] < 0.01:
                wts[i] = 0.0
            if wts[i] < 0.01 and context.portfolio.positions[sec].amount == 0:
                i = i+1
                continue
            order_target_percent(sec, wts[i],None,None)
            log.info('Adjusting ' + str(sec) + ' to ' + str(wts[i]*100.0) + '%')
            i=i+1

#################################################################################
#
# HELPER FUNCTIONS ( This algorithm's math functions)
#
#################################################################################                     

# Compute expected return on portfolio.
def compute_mean(W,R):
    return sum(R*W)

# Compute the variance of the portfolio.
def compute_var(W,C):
    return dot(dot(W, C), W)

# Combination of the two functions above - mean and variance of returns calculation. 
def compute_mean_var(W, R, C):
    return compute_mean(W, R), compute_var(W, C)

def fitness(W, R, C, r):
    # For given level of return r, find weights which minimizes portfolio variance.
    mean_1, var = compute_mean_var(W, R, C)
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint
    # Here, r is the 'target' return
    penalty = (1/100)*abs(mean_1-r)
    return var + penalty
    
# Given risk-free rate, asset returns, and covariances, this function
# calculates the weights of the tangency portfolio with regard to Sharpe
# ratio maximization.

def fitness_sharpe(W, R, C, rf):
    mean_1, var = compute_mean_var(W, R, C)
    utility = (mean_1 - rf)/sqrt(var)
    return 1/utility

# Solves for the optimal portfolio weights using the Sharpe ratio 
# maximization objective function.

def solve_weights(R, C, rf,context):
    n = len(R)
    W = ones([n])/n # Start optimization with equal weights
    b_ = [(0.,1.) for i in range(n)] # Bounds for decision variables
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) 
    # Constraints - weights must sum to 1
    optimized = scipy.optimize.minimize(fitness, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
    if not optimized.success:
        raise BaseException(optimized.message)
    return optimized.x    

# Solve for the efficient frontier using the variance + penalty minimization
# function fitness. 

def solve_frontier(R, C, rf,context):
    frontier_mean, frontier_var, frontier_weights = [], [], []
    n = len(R)      # Number of assets in the portfolio
    for r in linspace(max(min(R), rf), max(R), num=20): # Iterate through the range of returns on Y axis
        W = ones([n])/n # Set initial guess for weights
        b_ = [(0,1) for i in range(n)] # Set bounds on weights
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) # Set constraints
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:
        #    raise BaseException(optimized.message)s
        # Add point to the efficient frontier
        frontier_mean.append(r)
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights
        frontier_weights.append(optimized.x)
    return array(frontier_mean), array(frontier_var), frontier_weights
        
# Weights - array of asset weights (derived from market capitalizations)
# Expreturns - expected returns based on historical data
# Covars - covariance matrix of asset returns based on historical data

def assets_meanvar(daily_returns, stockList):    
    
    # Calculate expected returns
    expreturns = array([])
    daily_returns = daily_returns.transpose()
    
    for i in range(0, len(stockList)):
        expreturns = append(expreturns, mean(daily_returns[i,:]))
    
    # Compute covariance matrix
    covars = cov(daily_returns)
    # Annualize expected returns and covariances
    # Assumes 255 trading days per year    
    expreturns = (1+expreturns)**255-1
    covars = covars * 255
    expreturns = np.array(expreturns)
    
    return expreturns, covars
            
#################################################################################
#
# HELPER FUNCTIONS ( Portfolio / Stock Management)
#
#################################################################################         

def closeAnyOpenOrders(stock):
    orders = get_open_orders(stock)
    if orders:
        for order in orders:
             message = 'Canceling order for {amount} shares in {stock}'  
             message = message.format(amount=order.amount, stock=stock)  
             #log.debug(message)
             cancel_order(order)
def hasOrders(stock):
    hasOrders = False
    orders = get_open_orders(stock)
    if orders:
        for order in orders:
             message = 'Open order for {amount} shares in {stock}'  
             message = message.format(amount=order.amount, stock=stock)  
             log.info(message)
             hasOrders = True
    return hasOrders   
def liquidate(portfolio):
    for stock in portfolio.positions:
        closeAnyOpenOrders(stock)
        if portfolio.positions[stock].amount > 0:
         order(stock,-portfolio.positions[stock].amount,None,None)
         log.info('Sold stake in ' + str(stock) +
                  ' @ $' + str(portfolio.positions[stock].last_sale_price) +
                  ' Cost Basis: ' + str(portfolio.positions[stock].cost_basis) )
def marginRequirements(portfolio):  
    req = 0  
    for stock in portfolio.positions:  
        amount = portfolio.positions[stock].amount  
        last_price = portfolio.positions[stock].last_sale_price  
        if amount > 0:  
            req += .25 * amount * last_price  
        elif amount < 0:  
            if last_price < 5:  
                req += max(2.5 * amount, abs(amount * last_price))  
            else:  
                req += max(5 * amount, abs(0.3 * amount * last_price))  
    return req
   
    
There was a runtime error.

Just want to point out that last algorithm trades daily instead of weekly like all others.

Turns out there's a pretty fundamental problem with the code; the weights are the same for each point on the Frontier, just do a log.info(frontier_weights). I'm trying to find the problem, it's somewhere in the scipy.optimize function...

Edit: looks like r remains at 0.015 (rf) throughout the minimization iterations. WIP

Edit 2: I'm kind of baffled, it seems scipy.optimize.minimize doesn't pick up the new iteration values of r here:

def solve_frontier(R, C, rf, context):  
    frontier_mean, frontier_var, frontier_weights, = [], [], []  
    n = len(R)      # Number of assets in the portfolio  
    for r in linspace(max(min(R), rf), max(R), num=20): # Iterate through the range of returns on Y axis  
        W = ones([n])/n # Set initial guess for weights  
        b_ = [(0,1) for i in range(n)] # Set bounds on weights  
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin}) # Set constraints  
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', bounds=b_, constraints=c_) #PROBLEM: the r does not increase!  
        #if not optimized.success:  
            #raise BaseException(optimized.message)s  
        # Add point to the efficient frontier  
        frontier_mean.append(r) #OK  
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights  
        frontier_weights.append(optimized.x)  
    return array(frontier_mean), array(frontier_var), frontier_weights  

Anyone experienced with S.O.M. care to input?

Joel, see:
MAD Portfolio an alternative to Markowitz?

Scipy Optimise Minimise is a pig and the documentation is terrible. I am tempted to abandon it and use the critical line algo.

Edit: looks like r remains at 0.015 (rf) throughout the minimization iterations. WIP

I'm not sure I understand. r is merely an iterator to loop through the range of returns. It has nothing to do with rf (which is the risk free rate ). The risk free rate will remain the same throughout the test.

Correction to mean-variance solver: replaced "fitness" with "fitness_sharpe". For my own purposes I discarded margin provisions: put back in if you need to.

def solve_weights(R, C, rf,context):  
    n = len(R)  
    W = ones([n])/n # Start optimization with equal weights  
    b_ = [(lower_bound,upper_bound) for i in range(n)] # Bounds for decision variables  
    c_ = ({'type':'eq', 'fun': lambda W: sum(W)-1.0 })  
    # Constraints - weights must sum to 1  
    optimized = scipy.optimize.minimize(fitness_sharpe, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)  
    if not optimized.success:  
        raise BaseException(optimized.message)  
    return optimized.x  

Changed the penalty in "fitness" to make more sense. It is now a much blunter instrument and actually works to change the weights as planned:

def fitness(W, R, C, r):  
    # For given level of return r, find weights which minimizes portfolio variance.  
    mean_1, var = compute_mean_var(W, R, C)  
    # Penalty for not meeting stated portfolio return effectively serves as optimization constraint  
    # Here, r is the 'target' return  
    penalty = 100*abs(mean_1-r)  
    #print("mean_1",mean_1,"var",var,"penalty",penalty, "var + penalty",var+penalty)  
    return var + penalty  

I'll come back with more details but I am still finding surprising differences between the scipy.optimize.minimum approach and a monte carlo approach using 25,000 randomly generated portfolio weightings.

I take my words back on scipy.optimise - it isn't such a pig after all. But I intuitively prefer the broader, blunter, quasi brute force approach on the monte carlo method. Perhaps I could overcome that bias if I studied method='SLSQP' line by line.....but there again I can't help thinking a blunderbus approach in general helps to avoid local minima. Not that there are any in this sort of optimisation; but in general.....

But of course for very large portfolios I guess the MC approach is too cumbersome and time consuming.

I notice that from first to last you have all commented out the following code:

#if not optimized.success:  
            #   raise BaseException(optimized.message)  

In this context:

def solve_frontier(R, C, rf,context):  
    frontier_mean, frontier_var, frontier_weights = [], [], []  
    n = len(R)      # Number of assets in the portfolio  
    for r in linspace(max(min(R), rf), max(R), num=20): # Iterate through the range of returns on Y axis  
        W = ones([n])/n # Set initial guess for weights  
        b_ = [(0,1) for i in range(n)] # Set bounds on weights  
        c_ = ({'type':'eq', 'fun': lambda W: sum(W)-context.allowableMargin }) # Set constraints  
        optimized = scipy.optimize.minimize(fitness, W, (R, C, r), method='SLSQP', constraints=c_, bounds=b_)  
        #if not optimized.success:  
        #    raise BaseException(optimized.message)s  
        # Add point to the efficient frontier  
        frontier_mean.append(r)  
        frontier_var.append(compute_var(optimized.x, C))   # Min-variance based on optimized weights  
        frontier_weights.append(optimized.x)  
    return array(frontier_mean), array(frontier_var), frontier_weights  

Can anyone tell me what errors you were encountering and why you could not fix them?

def fitness_sharpe(W, R, C, rf):  
    mean, var = compute_mean_var(W, R, C)  
    utility = (mean - rf)/sqrt(var)  
    return 1/utility  

Likely to lead you into this problem:
scipy.optimize.minimize SLSQP leads to out of bounds solution scipy isues