Extending Leveraged ETFs Back In Time¶

plt.rcParams["figure.figsize"] = (16, 9)    # Matplotlib likes this in a separate cell

# Source is in percent.
fedfunds = read_fred('DFF').rename('FEDFUNDS') / 100
#tbill = read_fred('DTB3').rename('TBILL') / 100
#libor1d  = read_fred('USDONTD156N').rename('LIBOR-1d') / 100
#libor1w  = read_fred('USD1WKD156N').rename('LIBOR-1w') / 100
#libor1m  = read_fred('USD1MTD156N').rename('LIBOR-1m') / 100     
#libor12m = read_fred('USD12MD156N').rename('LIBOR-12m') / 100

BORROW = fedfunds  

def leverage(prices, factor=2, expense=0.01, borrow_rate=BORROW, ann_periods=TRADING_DAYS, 
             factor_scale=1.0, factor_incr=0, borrow_scale=1.0, expense_incr=0):
    """:Return: a Series giving the daily leveraged value of `prices` at a given leverage `factor`.
    
    This is basically the per-period change in prices minus the expense ratio and borrow rate.
    The expense ratio and borrow rate are divided (evenly, arithmetically) by `ann_periods`.
    
    :param float factor: The leverage factor by which daily returns are multiplied.
    :param float expense: Net expense ratio per `ann_periods` as a fraction. Deducted proportionally from each period.  
      Example: 0.0095 for a 0.95% annual expense ratio.
    :param Series borrow_rate: The (annualized) interest rate used to finance short-term borrowing for leverage.  
      Deducted from daily returns.  Typically the daily Federal Funds Rate or LIBOR. 
    :param int ann_periods: The number of periods over which rates are given.  E.g., 252 for daily periods in a typical trading year.
    
    The `_scale` and `_incr` parameters are adjustments to the corresponding parameters found through curve-fitting.
    """
    
    # Align borrow rates with prices
    prices = pd.Series(prices, dtype=float)
    if prices.isna().any():
        raise ValueError('NaN in prices')
    if isinstance(borrow_rate, pd.Series):
        if borrow_rate.isna().any():
            raise ValueError('NaN in borrow_rate')
        if prices.index[0] < borrow_rate.index[0]:
            raise ValueError(f'Prices start {prices.index[0]}, before borrow_rate {borrow_rate.index[0]}')
        borrow_rate = borrow_rate.reindex(index=prices.index, method='ffill') 
    name = f'{prices.name or ""}:{round(factor, 3)}X'
    
    # Curve-fitting adjustments 
    borrow_rate *= borrow_scale
    expense += expense_incr  # This functions as an additive constant for the whole equation, since exp isn't scaled ("C" above)
    
    change = prices.pct_change() * factor * factor_scale + 1      # Period-to-period changes as ratios
    # Changes less expenses and borrowing costs, evenly distributed among periods
    net_change = change - (expense + borrow_rate * (factor + factor_incr - 1)) / ann_periods
    net_change.iat[0] = prices.iat[0]   # Start leveraged series at same value so it's easily comparable
    leveraged = net_change.cumprod()
    return leveraged.rename(name)


# All possible leverage() parameters for curve fitting and their ranges
ALL_LEV_PARAMS = {'factor_scale': (0, 3), 'factor_incr': (-2, 2), 'borrow_scale': (-5, 5), 'expense_incr': (-2, 2)}
# Actually used leverage() parameters, determined in the 'Model Selection' section
LEV_PARAMS = {p: ALL_LEV_PARAMS[p] for p in ('expense_incr',)}

def norm(prices):
    """:Return: prices normalized to start at 1.0."""
    return prices / prices.iloc[0]


def cat(*dfs, dropna=True):
    """:Return: the column-wise concatenation of a sequence of Series or DataFrames.
    
    :param bool dropna: If True, remove rows with any NaN from the result.
    """
    result = pd.concat(dfs, axis=1)
    if dropna:
        result = result.dropna()
    return result


def align(*prices, dropna=True, norm=False):
    """:Return: The `prices` Series with only the dates in common to all of them, as a sequence.
    
    :param bool norm: If True, normalize each series of prices to start at 1.0.
    """
    aligned = cat(*prices, dropna=dropna)
    if norm:
        aligned = globals()['norm'](aligned)   # Calls norm() function, since we shadowed the name. Bit naughty.
    return tuple(col for _, col in aligned.iteritems())


def telltale(reference, *dfs, **layout_kws):
    """Plot the growth of several dataframes or series `dfs` relative to a `reference` series.
    
    https://www.bogleheads.org/wiki/Telltale_chart
    """
    tell = norm(cat(reference, *dfs))
    if tell.columns.nunique() < len(tell.columns):
        raise ValueError('Column names must be unique.')
    tell = tell.apply(lambda c: c / tell.iloc[:, 0])   # straight division doesn't work for some reason
    fig = tell.plot(title=f'Telltale Chart: {", ".join(tell.columns)}')
    if pd.options.plotting.backend == "plotly":
        fig.update_layout(margin=dict(t=50), **layout_kws).show()
               
      
def plotret(*prices, title=None):
    """Nice Plotly cumulative returns plot."""
    return norm(cat(*prices)).sub(1).plot(title=title).update_layout(
        yaxis=dict(tickformat=".0%"), 
        margin=dict(t=50), 
        legend_title_text='',
        yaxis_title='Cumulative Return',
        width=950,
        height=450,
    )
    
    
def color_leverage(factor, alpha=1.0, max_factor=3):
    """:Return: a plotly color string for a given leverage `factor`, more green for more long, more red for more short."""
    intensity = int(abs(factor) / max_factor * 255)
    if factor >= 0:
        return f'rgba(0,{intensity},0,{alpha})'
    else:
        return f'rgba({intensity},0,0,{alpha})'
    
    
def rasterize(figure, raster=False, width=1100, height=600, filename=None):
    """Maybe render a plotly `figure` as a static image, to save space and time."""
    if raster:
        if filename:
            figure.write_image(filename, width=width, height=height, scale=1, engine="kaleido")
            return Image(url=f'{filename}?cache_bust={np.random.randint(100000)}')
        else:
            return Image(figure.to_image(format="png", width=width, height=height, scale=1, engine="kaleido"))
    else:
        return figure    

    
def splice(old, new):
    """Splice together `new` prices with `old` prices before them, adjusted so new prices don't change."""
    if old.index[-1] < new.index[0]:
        raise ValueError(f'Last old index {old.index[-1]} and first new index {new.index[0]} must overlap')
    if old.index[0] > new.index[0]:
        warnings.warn(f'Old has no data older than new; old starts {old.index[0]}, new starts {new.index[0]}')
        return new
    first = old.index.get_loc(new.index[0], method='ffill')  # Find previous value if no exact match
    ratio = old.iloc[first] / new.iloc[0]
    return pd.concat((old.iloc[:first] / ratio, new), verify_integrity=True).rename(new.name)

def rmse(a, b):
    return np.sqrt(np.mean((a - b) ** 2))


def mae(a, b):
    return np.mean(np.abs(a - b))


def rel_rmse(a, b):
    """Relative RMSE between 1.0 and b / a, both aligned and normalized to start at 1.0."""
    a, b = align(a, b, norm=True)
    return rmse(1.0, b / a)
    

def rel_mae(a, b):
    """Relative MAE between 1.0 and b / a, both aligned and normalized to start at 1.0."""
    a, b = align(a, b, norm=True)
    return mae(1.0, b / a)
    
    
def return_rmse(a, b):
    """RMSE between simple periodic returns of price series `a` and `b`."""
    a, b = align(a.pct_change(), b.pct_change())
    return rmse(a, b)


def return_mae(a, b):
    """RMSE between simple periodic returns of price series `a` and `b`."""
    a, b = align(a.pct_change(), b.pct_change())
    return mae(a, b)


# This is, as you can imagine, sensitive to the most recent price.
def cumret_diff(a, b):
    """Absolute difference between the cumulative return of `a` and `b`."""
    a, b = align(a, b)
    return abs(cumret(a) - cumret(b))


def errstats(reference, leveraged, ann_periods=TRADING_DAYS):
    """:Return: dict of error metrics between expected `reference` price series and actual `leveraged` series
    at the dates (indices) they have in common."""
    reference, leveraged = align(reference, leveraged)
    return {
        'RMSE': rel_rmse(reference, leveraged), 
        'MAE': rel_mae(reference, leveraged),
        'RETRMSE': return_rmse(reference, leveraged),
        'RETMAE': return_mae(reference, leveraged),
        'CAGR': annret(leveraged, ann_periods) - annret(reference, ann_periods), 
        'VOL': annvol(leveraged, ann_periods) - annvol(reference, ann_periods), 
        'P99': norm(leveraged).div(norm(reference)).sub(1).abs().quantile(.99),
    }


def roundvals(d, digits=4):
    """Round the values of dict `d` to `digits` digits."""
    return {k: round(v, digits) for k, v in d.items()}

def find_params(reference, proxy, factor=2, expense=0.01, borrow_rate=BORROW, ann_periods=TRADING_DAYS, 
                params=LEV_PARAMS, errfunc=rel_rmse):
    """Find `params` that minimize the error between a `reference` series and its leveraged `proxy`.
    
    :param dict params: Maps parameters to `leverage()` to the range to search for optimal values.
    :param func errfunc: Error function that will be minimized; takes two price series and returns a distance metric between them.
    """
    if not params:
        return {}   # Well that was easy
    
    reference, proxy = align(reference, proxy)
    def obj(x):
        param_dict = dict(zip(params.keys(), x))   # param name: value
        return errfunc(reference, leverage(proxy, factor, expense, borrow_rate=borrow_rate, ann_periods=ann_periods, **param_dict))
    
    # Find params x that minimize obj(x)
    x0 = tuple(map(np.mean, params.values()))  # Initial guess = midpoint of bounds
    res = minimize(obj, x0, bounds=list(params.values()))
    best = dict(zip(params.keys(), res.x))   # param name: optimal value
    return best

def plotbest(reference, proxy, factor=2, expense=0.01, borrow_rate=BORROW, ann_periods=TRADING_DAYS, plot=True, errfunc=rel_rmse, params=LEV_PARAMS):
    """Find leverage parameters that minimize error between `reference` and leveraged `proxy`, plot a telltale
    chart, and return the new leveraged series."""
    best = find_params(reference, proxy, factor=factor, expense=expense, borrow_rate=borrow_rate, ann_periods=ann_periods, errfunc=errfunc, params=params)
    print(reference.name + ':' + proxy.name, '\tparams:', ', '.join(f'{k}={v}' for k, v in roundvals(best, 4).items()))

    # Get leveraged series with best params
    leveraged = leverage(proxy, factor=factor, expense=expense, borrow_rate=borrow_rate, ann_periods=ann_periods, **best)
    ref, lev = align(reference, leveraged, norm=True)  # Might be superfluous
    error = errstats(ref, lev, ann_periods=ann_periods)
    print(', '.join(f'{k}: {v}' for k, v in roundvals(error).items()))
    sim = leveraged[:reference.index[0]]
    simret = cumret(sim) if not sim.empty else 0
    print(f'CUMRET: sim {simret:.4f} + actual {cumret(reference):.4f} = {(simret + 1) * (cumret(reference) + 1) - 1:.4f}')
    if plot:
        telltale(ref, lev)
    return leveraged

letf, proxy = yget('UPRO'), yget('^SP500TR')

# Fund: (benchmark, leverage factor, expense ratio, issuer, start year for good data (or None to use all))
FUNDS = {
    # Mutual Funds
    'RYNVX': ('^SP500TR', 1.5, .0138, 'Rydex', '2000'),
    'ULPIX': ('^SP500TR', 2, .016, 'ProFunds', '2003'),
    #'RYTPX': ('^SP500TR', -2, .0184, 'Rydex', None),  # Bad data
    #'UOPIX': ('QQQ', 2, .0159, 'ProFunds', None),  # Bad data
    #'RYVNX': ('QQQ', -2, .0187, 'Rydex', None),    # Bad data
    'UAPIX': ('IWM', 2, .0178, 'ProFunds', '2003'),  # Russel 2000
    #'RYIRX': ('IWM', -2, .0191, 'Rydex', None),    # Russel 2000; Bad data
    'UMPIX': ('MDY', 2, .0166, 'ProFunds', '2003'), # S&P MidCap 400
    #'UDPIX': ('DIA', 2, .0172, 'ProFunds', None), # Dow; bad data
    'UTPIX': ('IDU', 1.5, .0173, 'ProFunds', '2004'), # Utilities
    'REPIX': ('IYR', 1.5, .0178, 'ProFunds', '2010'), # Real Estate
    #'SRPIX': ('IYR', -1, .0178, 'ProFunds', None), # Real Estate; Bad Data
    #'RYEUX': ('FEZ', 1.25, .0182, 'Rydex', None),  # EuroSTOXX 50; no good benchmark (data)
    'DXKLX': ('IEF', 2, .0143, 'Direxion', '2013'), # ITT
    'DXKSX': ('IEF', -2, .014, 'Direxion', '2013'), # ITT
    'UNPIX': ('EFA', 2, .0178, 'ProFunds', None),  # MSCI EAFE (large - mid foreign)
    'UUPIX': ('ADRE', 2, .0178, 'ProFunds', '2009'), # Emerging Markets
    
    
    # ETFs
    'SSO': ('^SP500TR', 2, .0091, 'ProShares', '2009'), # S&P 500
    'UPRO': ('^SP500TR', 3, .0093, 'ProShares', None),
    'SPXL': ('^SP500TR', 3, .0101, 'Direxion', '2013'),
    'SH': ('^SP500TR', -1, .009, 'ProShares', '2009'),
    'SDS': ('^SP500TR', -2, .0091, 'ProShares', '2009'),
    'SPXS': ('^SP500TR', -3, .0107, 'Direxion', '2013'),
    
    'QLD': ('QQQ', 2, .0095, 'ProShares', '2009'),     # NASDAQ 100
    'TQQQ': ('QQQ', 3, .0095, 'ProShares', None),
    'PSQ': ('QQQ', -1, .0095, 'ProShares', '2009'),
    'QID': ('QQQ', -2, .0095, 'ProShares', '2009'),
    'SQQQ': ('QQQ', -3, .0095, 'ProShares', None),     # Maybe 2013?
    
    'MVV': ('MDY', 2, .0095, 'ProShares', '2010'),     # MidCap 400
    'MYY': ('MDY', -1, .0095, 'ProShares', '2010'),
    'MZZ': ('MDY', -2, .0095, 'ProShares', '2010'),
    
    'UWM': ('IWM', 2, .0095, 'ProShares', '2010'),     # Russel 2000
    'TNA': ('IWM', 3, .0112, 'Direxion', '2013'),
    'RWM': ('IWM', -1, .0095, 'ProShares', '2010'),
    'TWM': ('IWM', -2, .0095, 'ProShares', '2009'),
    'TZA': ('IWM', -3, .0107, 'Direxion', '2013'),
    
    'URE': ('IYR', 2, .0095, 'ProShares', '2010'),     # Real Estate
    'REK': ('IYR', -1, .0095, 'ProShares', '2011'),
    #'SRS': ('IYR', -2, .0095, 'ProShares', '2010'),     # Poor fit
    
    'UPW': ('IDU', 2, .0095, 'ProShares', '2009'),     # Utilities
    'SDP': ('IDU', -2, .0095, 'ProShares', '2009'),
    
    'EFO': ('EFA', 2, .0095, 'ProShares', None),       # MSCI EAFE
    'EFZ': ('EFA', -1, .0095, 'ProShares', '2010'),    
    'EFU': ('EFA', -2, .0095, 'ProShares', '2010'),
    
    'EET': ('EEM', 2, .0095, 'ProShares', None),       # Emerging Markets
    'EUM': ('EEM', -1, .0095, 'ProShares', '2009'),
    'EEV': ('EEM', -2, .0095, 'ProShares', '2011'),
    
    'UST': ('IEF', 2, .0095, 'ProShares', '2012'),     # 7-10 Yr Treasury
    'TYD': ('IEF', 3, .0109, 'Direxion', '2010'),
    'TBX': ('IEF', -1, .0095, 'ProShares', None),
    'PST': ('IEF', -2, .0095, 'ProShares', '2011'),    
    'TYO': ('IEF', -3, .0108, 'Direxion', '2011'),
    
    'UBT': ('TLT', 2, .0095, 'ProShares', None),       # 20+ Yr Treasury
    'TMF': ('TLT', 3, .0105, 'Direxion', '2011'),
    'TBF': ('TLT', -1, .0094, 'ProShares', '2011'),
    'TBT': ('TLT', -2, .0092, 'ProShares', '2011'),    
    'TMV': ('TLT', -3, .0104, 'Direxion', '2011'),

    'UGL': ('GLD', 2, .0095, 'ProShares', None),       # Gold
    'GLL': ('GLD', -2, .0132, 'ProShares', None),
}

extras = ('VUSTX',)   # We'll use these later

def is_mutual_fund(ticker):
    return len(ticker) == 5 and ticker.endswith('X')

assert False, "The note-buck stops here"

def oos_error(funds, prices, param_ranges):
    """Find the best leverage parameters for the last half of each fund in `funds`, then use them
    to leverage the first half of each fund and compute the error."""
    tells = {}
    params = {}
    err = {}
    for name, (proxy, factor, exp, _, start) in funds.items():
        # Cut out initial bad data from leveraged fund, align with proxy
        proxy = prices[proxy].dropna()
        lfund, _ = align(prices.loc[start:, name], proxy)  # We do *not* want to modify the proxy
        # Find best params for last half of data
        assert len(lfund) > 500, "That's not enough data!"
        mid = len(lfund) // 2
        params[name] = find_params(lfund.iloc[mid:], proxy, factor, exp, params=param_ranges)
        # Use params to leverage all data
        lev = leverage(proxy, factor, exp, **params[name])
        # Compute error on first half (out of sample)
        err[name] = errstats(lfund.iloc[:mid], lev)
        # Plot the whole thing
        lfund, lev = align(prices[name], lev, norm=True)   # Plot whole series including bad initial data
        tells[name] = lev.div(lfund).rename(name)
        
    return pd.DataFrame.from_dict(err, orient='index'), pd.DataFrame.from_dict(params, orient='index'), pd.concat(tells, axis=1, verify_integrity=True)

def powerset(iterable):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))

# Show error stats for best model / fit params
err, params, tells = oos_error(FUNDS, prices, LEV_PARAMS)