Hey everyone!

Attached below is my code. It has worked well in the Quant book but when put into my algorithm and during backtesting, it keeps telling me the data coming from the ‘_black_litterman_’ function ‘ must be 1 dimensional’, when I'm assured it is. Can someone help me understand my problem and why it keeps failing during backtesting?

NOTE : The problematic line seems to be line 92 : starting ‘mu_bl = …’


import pandas as pd
import numpy as np
from scipy.optimize import minimize
from sklearn.linear_model import ElasticNet
from sklearn.model_selection import GridSearchCV
from numpy.linalg import inv
from AlgorithmImports import *
from datetime import datetime
from datetime import timedelta
from io import StringIO
import matplotlib as mpl
from matplotlib import cm
import matplotlib.pyplot as plt
from arch import arch_model

class ghostyfund(QCAlgorithm):

def Initialize(self):


self.UniverseSettings.Resolution = Resolution.Daily
self.AddUniverse(self.CoarseSelectionFunction, self.FineSelectionFunction)
self.UniverseSettings.Leverage = 3

self._changes = None

def implied_returns(self, delta, sigma, w):
delta = Risk aversion coefficient (scalar) (excess return on the market / variance of the market)
w = portfolio weights (N x 1) as Series
sigma = N x N covariance matrix as df

Returns a (N x 1) vector of implied expected returns by reverse engineering the portfolio weights
ir = delta * sigma.values.dot(w).squeeze()
return ir

def proportional_prior(self, sigma, tau, p, index):
sigma = N x N covariance matrix as df
tau = a scalar
p = a K x N df linking Q & Assets

Returns a (P x P) df, a matrix representing prior uncertainties
helit_omega = p.values.dot(tau * sigma.values).dot(p.values.T)
return pd.DataFrame(np.diag(np.diag(helit_omega)), index=index, columns=index)

def _black_litterman_(self, w_prior, sigma_prior, p, p_not_values, q, omega=None, delta=2.5, tau=0.2):
Computes the black litterman expected returns
W.prior must be an N x 1 vector of weights, series
Sigma.prior is an N x N covariance matrix, a df
P must be an K x N linking Q and the assets, a df
Q must be a K x 1 vector of views, a series
Omega must be a K x K matrix, a df

if omega == None:
its proportional to the variance of the prior

Delta, tau are scalars
index = sigma_prior.index

def as_column_vec(x):
if (x.ndim == 2):
return x
return np.expand_dims(x, axis=1)

if omega is None:
omega = self.proportional_prior(sigma_prior, tau, p_not_values, index)

# no of assets
N = w_prior.shape[0]
# no of views
K = q.shape[0]
# 1) Reverse engineer weights to get pi
pi = self.implied_returns(delta, sigma_prior, w_prior)
# 2) Adjust (scale) sigma by the uncertainty scale factor
sigma_prior_scaled = sigma_prior.values * tau
sigma_prior = sigma_prior.values
# note; .dot == matrix multiplication
mu_bl = pi + sigma_prior_scaled.dot(p.T).dot(inv(p.dot(sigma_prior_scaled).dot(p.T) + omega).dot(q - p.dot(pi)))
# posterior estimate of uncertainty of mu_bl
sigma_bl = sigma_prior + sigma_prior_scaled - sigma_prior_scaled.dot(p.T).dot(inv(p.dot(sigma_prior_scaled).dot(p.T) + omega)).dot(p).dot(sigma_prior_scaled)

return (mu_bl, sigma_bl)

def black_litty_implied_rets(self, tickers):
q = pd.Series(index=tickers, data=[0.03]*int(len(tickers)))
p = pd.DataFrame(columns=tickers, index=range(0, len(tickers)))
w = pd.Series({}, index=tickers).fillna(float(1/len(tickers)))

delta = 2.5
tau = .02

for i in range(0, len(p.columns)):
p.iloc[i, i] = 1

p = p.fillna(0)

history = self.History(tickers, 60, Resolution.Daily)["close"]
sigma_prior = history.pct_change().unstack(level=0).dropna().cov()
index = sigma_prior.index

pi = self.implied_returns(delta, sigma_prior, w)
omega = self.proportional_prior(sigma_prior, tau, p, index)

bl = self._black_litterman_(w_prior=w, sigma_prior=sigma_prior, p=p.values, p_not_values=p, q=q)

df = pd.DataFrame(bl[0])
df["Ticker"] = tickers
df = df.pivot_table(columns="Ticker")
exp_rets = df.sort_values(by=0, ascending=False, axis=1)
#highest_exp_ret = exp_rets.iloc[0].index[0]

return exp_rets

def monte_carlo_gbm(self, n_years, n_scenarios, mu, sigma, steps_per_year, s_0):
Evolution of Geometric Brownian Motion trajectories, such as for Stock Prices through Monte Carlo
:param n_years: The number of years to generate data for
:param n_paths: The number of scenarios/trajectories
:param mu: Annualized Drift, e.g. Market Return
:param sigma: Annualized Volatility
:param steps_per_year: granularity of the simulation
:param s_0: initial value
:return: a numpy array of n_paths columns and n_years*steps_per_year rows
dt = 1/steps_per_year
n_steps = int(n_years*steps_per_year) + 1
rets_plus_1 = np.random.normal(loc=mu*dt+1, scale=sigma*np.sqrt(dt), size=(n_steps, n_scenarios))
rets_plus_1[0] = 1
prices = s_0*pd.DataFrame(rets_plus_1).cumprod()
return prices

def get_max_drawdown(self, array):
drawdowns = []
for i in range(len(array)):
max_array = max(array[:i+1])
drawdown = max_array - array[i]

return max(drawdowns)

def CoarseSelectionFunction(self, coarse):
sortedByDollarVolume = sorted(coarse, key=lambda x: x.DollarVolume, reverse=True)
return [x.Symbol for x in sortedByDollarVolume]

def FineSelectionFunction(self, fine):
sortedByPeRatio = sorted(fine, key=lambda x: x.ValuationRatios.PERatio)
tickers = [x.Symbol for x in sortedByPeRatio][:5]
return tickers

def OnData(self, data):
if self._changes is None:

tickers = [x.Symbol for x in self._changes.AddedSecurities]

mu = self.black_litty_implied_rets(tickers)
iteration = 0
avg_drawdowns = []

for ticker in tickers:
symbol = self.AddEquity(ticker).Symbol
symbol_hist = self.History(symbol, 252, Resolution.Daily)["close"]
returns = symbol_hist.pct_change().dropna().unstack(level=0)
vol = returns.rolling(window=30, center=False).std().dropna()

am = arch_model(vol, vol='Garch', p=1, o=0, q=1, dist='Normal')
res = am.fit(disp='off')
garch_implied_vol = 0.1 * np.sqrt(res.params['omega'] + res.params['alpha[1]'] * res.resid**2 +
res.conditional_volatility**2 * res.params['beta[1]'])

mc = self.monte_carlo_gbm(n_years=1, n_scenarios=2500,
mu=mu.iloc[:, iteration], sigma=garch_implied_vol.max(), steps_per_year=252, s_0=100)

drawdowns = []
for i in range(0, len(mc.columns)):
drawdown = self.get_max_drawdown(mc[i])

avg_drawd = np.mean(drawdowns)

iteration += 1

df = pd.DataFrame(avg_drawdowns)
df["Ticker"] = tickers
df = df.pivot_table(columns="Ticker")
expected_drawdowns = df.sort_values(by=0, ascending=True, axis=1)

ticker_symbol = expected_drawdowns.columns[0]

self.SetHoldings(ticker_symbol, 1)

self._changes = None

def OnSecuritiesChanged(self, changes):
self._changes = changes