portfolio optimization example - strange results.

[andrewczgithub]

Hi @mmckerns !!

Thank you so much for your help. I am trying to use mystic for a basic portfolio optimzation problem but i am getting strange results. observe the accurate example below and mystic's example.

def sharpe_ratio(x):
    return np.mean(x) / np.std(x)

#change to more accurate sharp ration - method of moments sharp ratio.

def zscore(x, window):
    r = x.rolling(window=window)
    m = r.mean().shift(1)
    s = r.std(ddof=0).shift(1)
    return (x-m)/s

class Portfolio:

    def __init__(self, returns, regularized=False, c=1e-6):
        self.returns = returns
        self.regularized = regularized
        self.c = c

    def fitness_function(self, weights):
        weighted_returns = calc_weighted_returns(self.returns, weights)
        neg_sharpe_ratio = -1 * sharpe_ratio(weighted_returns)
        penalty = self.c * np.linalg.norm(weights, ord=1)
        if self.regularized:
            return neg_sharpe_ratio
        else:
            return neg_sharpe_ratio + penalty
    #need to change to mystic
    def optimize_portfolio(self, display=False):
        num_strategies = self.returns.shape[1]
        initial_guess = [0.5] * num_strategies
        optimization_results = cma.fmin(
                self.fitness_function, initial_guess, 0.25, 
                options={'bounds':[0,1], 'verb_disp':display, 'tolfun':0.0001}
        )
        return optimization_results[0]

#mystic test

from mystic.solvers import fmin

def sharpe_ratio(x):
    return np.mean(x) / np.std(x)

#change to more accurate sharp ration - method of moments sharp ratio.

def zscore(x, window):
    r = x.rolling(window=window)
    m = r.mean().shift(1)
    s = r.std(ddof=0).shift(1)
    return (x-m)/s

class Portfolio:

    def __init__(self, returns, regularized=False, c=1e-6):
        self.returns = returns
        self.regularized = regularized
        self.c = c

    def fitness_function(self, weights):
        weighted_returns = calc_weighted_returns(self.returns, weights)
        neg_sharpe_ratio = -1 * sharpe_ratio(weighted_returns)
        penalty = self.c * np.linalg.norm(weights, ord=1)
        if self.regularized:
            return neg_sharpe_ratio
        else:
            return neg_sharpe_ratio + penalty
    #need to change to mystic
    def optimize_portfolio(self, display=False):
        num_strategies = self.returns.shape[1]
        initial_guess = [0.5] * num_strategies
        optimization_results = fmin(
                self.fitness_function, initial_guess
               # options={'bounds':[0,1], 'verb_disp':display, 'tolfun':0.0001}
        )
        return optimization_results[0]

import numpy as np
from astropy.stats import median_absolute_deviation

def modified_z_score(ys, window):
    ys_arr = np.array(ys)
    ys_arrr = np.roll(ys_arr, window)
    median_y = np.median(ys_arrr)
    median_absolute_deviation_y = median_absolute_deviation(ys_arrr)
    z_score = 0.6745 * (ys_arrr - median_y) / median_absolute_deviation_y
    return z_score

strategies = pd.DataFrame(index=etf_returns.index)
strategy_returns = pd.DataFrame(index=etf_returns.index)

spy = etfs['SPY']
spy_returns = spy['Open'].pct_change().shift(-2)
strategies['SPY_5'] = modified_z_score(spy['Close'], 5) * -1
strategies['SPY_50'] = modified_z_score(spy['Close'], 50)
strategies['SPY_100'] = modified_z_score(spy['Close'], 100)

tlt = etfs['TLT']
tlt_returns = tlt['Open'].pct_change().shift(-2)
strategies['TLT_5'] = modified_z_score(tlt['Close'], 5) * -1
strategies['TLT_50'] = modified_z_score(tlt['Close'], 50)
strategies['TLT_100'] = modified_z_score(tlt['Close'], 100)

for column in strategies.columns:
    symbol = column[:3]
    if symbol == 'SPY':
        strategy_returns[column] = strategies[column].multiply(spy_returns)
    elif symbol == 'TLT':
        strategy_returns[column] = strategies[column].multiply(tlt_returns)

portfolio = Portfolio(strategy_returns)
results = portfolio.optimize_portfolio()
optimal_weights = pd.DataFrame(np.round(results, 2), index=strategy_returns.columns, columns=[''])
optimal_weights

[andrewczgithub]

The result i get for the cma solver is -

  |  
-- | --
0.23
0.64
0.00
1.00
0.00
0.52

[andrewczgithub]

The result i get for mystic is -

  |  
-- | --
0.13
0.13
0.13
0.13
0.13
0.13

[andrewczgithub]

How do I use mystic to perform a global optimziation solver.

[mmckerns]

@andrewczgithub : fmin is not a global optimizer, so any result you are seeing is going to be due to at least the initial conditions (unless the function you are minimizing is convex). Within the use of fmin, there's a number of settings that can change how the solver progresses along a surface, for example you can change the size of the simplex, or the gradient tolerance, etc -- and if you use the class interface to the fmin solver, you can customize the solver even more, creating custom termination conditions and the like (see NelderMeadSimplexSolver). If you are looking for a global optimizer, you might want to try diffev (differential evolution) or buckshot (a pseudo-global ensemble optimizer).

I also encourage you to use a LoggingMonitor object and one of the monitor viewers to visualize how the solver is progressing -- this should give you a lot of insight into why you ar getting the results you are getting, and how to tweak the solver or solver parameters.

[andrewczgithub]

Hi @mmckerns !! Thank you so much for your help!!

I will give the buckshot solver a try over the weekend.

Would you be able to give some advice abour how i would code the solver in the below class?

#mystic test

from mystic.solvers import fmin

def sharpe_ratio(x):
    return np.mean(x) / np.std(x)

#change to more accurate sharp ration - method of moments sharp ratio.

def zscore(x, window):
    r = x.rolling(window=window)
    m = r.mean().shift(1)
    s = r.std(ddof=0).shift(1)
    return (x-m)/s

class Portfolio:

    def __init__(self, returns, regularized=False, c=1e-6):
        self.returns = returns
        self.regularized = regularized
        self.c = c

    def fitness_function(self, weights):
        weighted_returns = calc_weighted_returns(self.returns, weights)
        neg_sharpe_ratio = -1 * sharpe_ratio(weighted_returns)
        penalty = self.c * np.linalg.norm(weights, ord=1)
        if self.regularized:
            return neg_sharpe_ratio
        else:
            return neg_sharpe_ratio + penalty
    #need to change to mystic
    def optimize_portfolio(self, display=False):
        num_strategies = self.returns.shape[1]
        initial_guess = [0.5] * num_strategies
        optimization_results = fmin(
                self.fitness_function, initial_guess
               # options={'bounds':[0,1], 'verb_disp':display, 'tolfun':0.0001}
        )
        return optimization_results[0]

[mmckerns]

The documentation for fmin is here: https://github.com/uqfoundation/mystic/blob/master/mystic/scipy_optimize.py#L364. Your options above correspond to bounds, disp, and ftol.

[andrewczgithub]

Hi @mmckerns !!

Thank you so much your help!! this is so interesting!! I am trying to use the buckshot algorithm for global optimzation.

This is my attempt below

#mystic test

from mystic.solvers import BuckshotSolver
lb = [0.0, 0.0, 0.0]
ub = [2.0, 2.0, 2.0]

# get monitors and termination condition objects
from mystic.monitors import Monitor
stepmon = Monitor()
from mystic.termination import CandidateRelativeTolerance as CRT

# select the parallel launch configuration
from pyina.launchers import Mpi as Pool

def sharpe_ratio(x):
    return np.mean(x) / np.std(x)

#change to more accurate sharp ration - method of moments sharp ratio.

def zscore(x, window):
    r = x.rolling(window=window)
    m = r.mean().shift(1)
    s = r.std(ddof=0).shift(1)
    return (x-m)/s

class Portfolio:

    def __init__(self, returns, regularized=False, c=1e-6):
        self.returns = returns
        self.regularized = regularized
        self.c = c

    def fitness_function(self, weights):
        weighted_returns = calc_weighted_returns(self.returns, weights)
        neg_sharpe_ratio = -1 * sharpe_ratio(weighted_returns)
        penalty = self.c * np.linalg.norm(weights, ord=1)
        if self.regularized:
            return neg_sharpe_ratio
        else:
            return neg_sharpe_ratio + penalty
    #need to change to mystic
    def optimize_portfolio(self, display=False):
        num_strategies = self.returns.shape[1]
        initial_guess = [0.5] * num_strategies
        solver = BuckshotSolver(len(lb), npts)
        solver.SetMapper(Pool(NNODES).map)
        solver.SetGenerationMonitor(stepmon)
        solver.SetTermination(CRT())
        optimization_results = solver.Solve(

                self.fitness_function, initial_guess
               # options={'bounds':[0,1], 'verb_disp':display, 'tolfun':0.0001}
        )
        return solver.Solution()

[andrewczgithub]

Hi @mmckerns ,

Just a question, should in not be using the below the get the global optimzation?

https://github.com/uqfoundation/mystic/blob/master/examples3/test_searcher.py

Best, Andrew

[mmckerns]

Buckshot (or the one-liner, buckshot) is probably a good choice. It's a pseudo-global solver. You will want to increase the number of solvers in the ensemble. The searcher you are referring to is typically used as an improved replacement for parameter sweeps. The searcher codes are used to find all minima and maxima, so you can efficiently and accurately learn/interpolate/estimate the unknown function surface.

[andrewczgithub]

Hi @mmckerns !!

Thank you so much for your help!!

This is my first attempt at buckshot.

#mystic test

from mystic.solvers import BuckshotSolver
from mystic.termination import NormalizedChangeOverGeneration as NCOG

lb = [0.0, 0.0, 0.0]
ub = [2.0, 2.0, 2.0]
npts = 20

# get monitors and termination condition objects
from mystic.monitors import Monitor
stepmon = Monitor()
from mystic.termination import CandidateRelativeTolerance as CRT

# select the parallel launch configuration
from pyina.launchers import Mpi as Pool

def sharpe_ratio(x):
    return np.mean(x) / np.std(x)

#change to more accurate sharp ration - method of moments sharp ratio.

def zscore(x, window):
    r = x.rolling(window=window)
    m = r.mean().shift(1)
    s = r.std(ddof=0).shift(1)
    return (x-m)/s

class Portfolio:

    def __init__(self, returns, regularized=False, c=1e-6):
        self.returns = returns
        self.regularized = regularized
        self.c = c

    def fitness_function(self, weights):
        weighted_returns = calc_weighted_returns(self.returns, weights)
        neg_sharpe_ratio = -1 * sharpe_ratio(weighted_returns)
        penalty = self.c * np.linalg.norm(weights, ord=1)
        if self.regularized:
            return neg_sharpe_ratio
        else:
            return neg_sharpe_ratio + penalty
    #need to change to mystic
    def optimize_portfolio(self, display=False):
        num_strategies = self.returns.shape[1]
        initial_guess = [0.5] * num_strategies
        solver = BuckshotSolver(len(lb), npts)
        #solver.SetMapper(Pool(NNODES).map)
        #solver.SetGenerationMonitor(stepmon)
        solver.SetTermination(CRT())
        optimization_results = solver.Solve(

                self.fitness_function, initial_guess, 0.25,
            bounds=[0,1], ftol=0.005, disp =True

        )
        return solver.Solution()

[andrewczgithub]

import numpy as np
from astropy.stats import median_absolute_deviation

def modified_z_score(ys, window):
    ys_arr = np.array(ys)
    ys_arrr = np.roll(ys_arr, window)
    median_y = np.median(ys_arrr)
    median_absolute_deviation_y = median_absolute_deviation(ys_arrr)
    z_score = 0.6745 * (ys_arrr - median_y) / median_absolute_deviation_y
    return z_score

[andrewczgithub]

strategies = pd.DataFrame(index=etf_returns.index)
strategy_returns = pd.DataFrame(index=etf_returns.index)

spy = etfs['SPY']
spy_returns = spy['Open'].pct_change().shift(-2)
strategies['SPY_5'] = modified_z_score(spy['Close'], 5) * -1
strategies['SPY_50'] = modified_z_score(spy['Close'], 50)
strategies['SPY_100'] = modified_z_score(spy['Close'], 100)

tlt = etfs['TLT']
tlt_returns = tlt['Open'].pct_change().shift(-2)
strategies['TLT_5'] = modified_z_score(tlt['Close'], 5) * -1
strategies['TLT_50'] = modified_z_score(tlt['Close'], 50)
strategies['TLT_100'] = modified_z_score(tlt['Close'], 100)

for column in strategies.columns:
    symbol = column[:3]
    if symbol == 'SPY':
        strategy_returns[column] = strategies[column].multiply(spy_returns)
    elif symbol == 'TLT':
        strategy_returns[column] = strategies[column].multiply(tlt_returns)

ortfolio = Portfolio(strategy_returns)
results = portfolio.optimize_portfolio()
optimal_weights = pd.DataFrame(np.round(results, 2), index=strategy_returns.columns, columns=[''])
optimal_weights

[andrewczgithub]

i get the below error

TypeError Traceback (most recent call last)

in 1 portfolio = Portfolio(strategy_returns) ----> 2 results = portfolio.optimize_portfolio() 3 optimal_weights = pd.DataFrame(np.round(results, 2), index=strategy_returns.columns, columns=['']) 4 optimal_weights in optimize_portfolio(self, display) 56 57 self.fitness_function, initial_guess, 0.25, ---> 58 bounds=[0,1], ftol=0.005, disp =True 59 60 ~/mystic/mystic/abstract_solver.py in Solve(self, cost, termination, ExtraArgs, **kwds) 930 931 # register: cost, termination, ExtraArgs --> 932 cost = self._bootstrap_objective(cost, ExtraArgs) 933 if termination is not None: self.SetTermination(termination) 934 #XXX: self.Step(cost, termination, ExtraArgs, **settings) ? ~/mystic/mystic/abstract_solver.py in _bootstrap_objective(self, cost, ExtraArgs) 734 return _cost 735 # 'wrap' the 'new' cost function with _decorate --> 736 self.SetObjective(cost, ExtraArgs) 737 return self._decorate_objective(*self._cost[1:]) 738 ~/mystic/mystic/abstract_solver.py in SetObjective(self, cost, ExtraArgs) 589 args = () if args is None else args 590 # quick validation check (so doesn't screw up internals) --> 591 if not isvalid(cost, [0]*self.nDim, *args): 592 try: name = cost.__name__ 593 except AttributeError: # raise new error for non-callables TypeError: isvalid() argument after * must be an iterable, not float

[andrewczgithub]

The reference for the above strategy is below -

http://quantfiction.com/2018/12/20/position-sizing-for-practitioners-part-3-a-portfolio-approach/

[andrewczgithub]

I am not sure how to code the solver to be able to get the global minimum.

Best, Andrew

[mmckerns]

@andrewczgithub: if you check the documentation for Solve, the second argument is the termination condition, not the initial value. mystic.solvers.buckshot and Solve from a mystic.solvers.BuckshotSolver instance have different interfaces. You might want to check the examples directory in the GitHub repo for examples of how to use either of the two interfaces to the buckshot solver.

[andrewczgithub]

hey @mmckerns !! Ah ok Though which one am i supposed to use? What is the one liiner solver?

[mmckerns]

the one-liner is a function call wrapped around the solver class, where the function call mimics the spicy.optimize interface for solvers.

[mmckerns]

@andrewczgithub: I'm closing this issue. Please reopen if there's more to do here.