<jupyter_start><jupyter_text>Treasury Bond Fund SimulatorThis is based off of [longinvest's post on bogleheads][1]. He implemented it all in a spreadsheet (which is linked in the thread).The goal is to calculate returns of a simulated bond fund given a bunch of interest rates. By having returns (instead of rates) we can perform backtesting with historical data.Keep in mind that this a simulation so it has some simplifying things that will make it slightly different from the real world. This is not a comprehensive list:* It assumes no transaction costs. This is mostly truish for Treasuries, though.* It assumes an extremely liquid marketplace where there is no margin when selling. This is mostly truish for Treasuries, as well.* Because bonds are only rolled once a year, [the duration of the fund shifts a bit][2].* Real funds [shift their durations all the time][3], even if only by small amounts. Our simulated shifts and their actual shifts are unlikely to line up exactly.* All interest is reinvested rather than being paid out. This actually works in our favor, since total returns are what we care about for backtesting anyway.First we need to import some libraries....[1]: https://www.bogleheads.org/forum/viewtopic.php?f=10&t=179425[2]: https://www.bogleheads.org/forum/viewtopic.php?f=10&t=179425&start=100p2998433[3]: https://www.bogleheads.org/forum/viewtopic.php?f=10&t=179425&start=100p2998031<jupyter_code>import numpy
from collections import deque
import pandas
import math<jupyter_output><empty_output><jupyter_text>Simulating the Bond FundSimulating the bond fund is conceptually straightforward.You have a ladder of bonds, one for each year. Something like this:    deque([Maturity: 1 | Yield: 5.00% | Face Value: $50.00,       Maturity: 2 | Yield: 5.00% | Face Value: $52.50,       Maturity: 3 | Yield: 5.00% | Face Value: $55.12,       Maturity: 4 | Yield: 5.00% | Face Value: $57.88,       Maturity: 5 | Yield: 5.00% | Face Value: $60.78,       Maturity: 6 | Yield: 5.00% | Face Value: $63.81,       Maturity: 7 | Yield: 5.00% | Face Value: $67.00,       Maturity: 8 | Yield: 5.00% | Face Value: $70.36,       Maturity: 9 | Yield: 5.00% | Face Value: $73.87])           Every year three things will happen:1. All of the bonds will pay out their cash coupon. This is based on their yield and their face value.1. When a bond gets "too young" (I'll come back to this) we sell it. The exact price will also be explained later. Every year you will sell one bond of the youngest maturity.1. Now you've got a pile of cash and one fewer bond. Use the cash to buy a new bond of the longest maturity. Youngest maturity & oldest maturityWhen you create the bond fund, you can select the youngest maturity and the oldest maturity. Say that you want fund where the oldest bond has a 10-year maturity and the youngest bond has a 2-year maturity. As a shorthand, we'll call that a 10-2 fund. Every year a 2-year bond becomes a 1-year bond and will be sold and replaced with a brand new 10-year bond.<jupyter_code>def iterate_fund(ladder, yield_curve, max_maturity):
    reduce_maturity(ladder)
    
    payments = get_payments(ladder)

    sold_bond = ladder.popleft()
    payments += sold_bond.value(yield_curve)

    new_bond = Bond(payments, yield_curve[max_maturity-1], max_maturity)
    ladder.append(new_bond)
    
    # This happens *after* we sell the shortest bond and buy a new long one
    # (at least, that's what longinvest does...)
    nav = get_nav(ladder, yield_curve)

    return (ladder, payments, nav)

def get_nav(ladder, rates):
    return sum((b.value(rates) for b in ladder))

def get_payments(ladder):
    return sum((b.gen_payment() for b in ladder))

def reduce_maturity(ladder):
    for b in ladder:
        b.maturity -= 1
    return ladder<jupyter_output><empty_output><jupyter_text>Bond MechanicsA bond is just three things: a yield, a face value, and a maturity. If you called up your broker you would say, "I want to buy $100 of the 10-year Treasury that is yielding 3.2%." The maturity is 10-years; the face value is $100; and the yield is 3.2%.There are only two things you can do with a bond. Receive your paymentEvery year the bond will generate a payment -- a "coupon" in bond-speak. This is simply the **yield × face value**. Going back to the previous example, with a face value of $100 and a yield of 3.2%, every year you would get a payment of $3.20. (Not very impressive, admittedly.) Check the current value of the bondBonds are designed to be held until their maturity. At that point you'll receive a payment for the face value. In our example, that would mean after holding the bond for 10 years you would get your full $100 back.But what if you wanted to sell the bond **before** maturity? That's (usually) possible but the exact price will depend on current rates. Say we want to sell our bond after 9 years. In essence, we have a 1-year bond that yields 3.2%. What if the current going yield for 1-year bonds was 2.5%? Then our bond will be worth a little more. If the current going yield for 1-year bonds is 4.2% then our bond will be worth a little less.Here's how it actually gets calculated:* take the current maturity remaining on the bond* take the current yield on bonds of that maturity* take the bond face valueThen mix all of them into present value calculation: **pv(current yield, current maturity, face value)**. (The pv function is found in every spreadsheet and many calculators.)If the current yields are 2.5% and you have a face value of $100 then the present value is $97.56. This calculates the "exact" value of the bond. In the real world, someone would probably offer you slightly less, maybe $97.25, because they are offering you liquidity and taking some risk. For Treasuries that is a very small number, so ignoring it is probably good enough. From rates to returnsAt the end of the day, checking the current value of the bonds we hold is what we're trying to achieve. By adding up the value of all the bonds we hold we can figure out the Net Asset Value (NAV) of our fund. And then we compare that NAV over time. This is what we wanted: to be able to calculate the returns of a (simulated) bond fund.<jupyter_code>class Bond:
    def __init__(self, face_value, yield_pct, maturity):
        self.face_value = face_value
        self.yield_pct = yield_pct
        self.maturity = maturity
        
    def __repr__(self):
        return ('Maturity: %d | Yield: %.2f%% | Face Value: $%.2f' % (self.maturity, self.yield_pct * 100, self.face_value))
        
    def gen_payment(self):
        return self.face_value * self.yield_pct
    
    def value(self, rates):
        value = numpy.pv(rates[self.maturity - 1], self.maturity, self.gen_payment(), self.face_value)
        return -value<jupyter_output><empty_output><jupyter_text>Bootstrapping the LadderOur bond ladder is straightforward enough. Sell the youngest bond and buy another one of the old bonds, using whatever cash we currently have available.But how do you get the ladder **started**? Where do those first bonds come from?Here's where things get a little bit unavoidably hacky. In the real world, you could slowly build up a ladder over time. For instance, buy 1/10th of the ladder every year for a decade. That takes, well, a decade. Which means there's an entire decade in our simulation with no results. We can shortcut that at the cost of a slight loss of accuracy for those first few years.If we're building a 10-2 ladder then we have 9 bonds (we don't have a bond with 1-year maturity, hence only 9 bonds). We bootstrap the ladder by buying all 9 instantly. That means they will all have the same yield -- whatever the current yield is.<jupyter_code>def bootstrap(yield_curve, max_bonds, min_maturity):
    bond_yield = yield_curve[max_bonds - 1]
    ladder = deque()
    starting_face_value = 50 # chosen arbitrarily (to match longinvest)

    for i, j in zip(range(max_bonds), range(min_maturity, max_bonds+1)):
        face_value = pow(1 + bond_yield, i) * starting_face_value
        b = Bond(face_value, bond_yield, j)
        ladder.append(b)
    return ladder
bootstrap([.0532]*10, 10, 2)<jupyter_output><empty_output><jupyter_text>Why do we have a different face value for each one? Why not just $50 for each? That's how we ensure that each rung of the ladder has equivalent value. If they each had a face value of $50 then the longer-term bonds would actually be worth less than the younger-term bonds. [For more explanation, refer to this post by longinvest][1].[1]: https://www.bogleheads.org/forum/viewtopic.php?p=3142416p3142416  RatesNow that we understand how the ladder works and how to bootstrap it, we need a source of rates in order to drive the engine.We have a number of sources of rate data.* Shiller provides 10 year yields on Treasuries, going back to 1871* Shiller provides 1 year interest rates, going back to 1871* [FRED provides 1-, 2-, 3-, 5-, 7-, 20-, and 30-year rates][1]. The data begins in the 1954-1977 range. When available, we prefer the FRED data over Shiller data.* NBER provides [historical data (1942-1962)][2] for 20-year rates* NBER provides [yields on long-term bonds (1919-1944)][3]So we will start by importing those. (I've spliced them all into a single CSV file to make importing things simpler; [you can get the CSV from github][4])[1]: https://fred.stlouisfed.org/categories/115[2]: https://fred.stlouisfed.org/series/M13058USM156NNBR[3]: https://fred.stlouisfed.org/series/M1333AUSM156NNBR[4]: https://github.com/hoostus/prime-harvesting/blob/master/bond_rates.csv<jupyter_code>HISTORICAL_RATES = pandas.read_csv('bond_rates.csv', index_col=0)
HISTORICAL_RATES.head()<jupyter_output><empty_output><jupyter_text>Rate interpolationFor a given year, we will have **some** rate data. At the very least we will have the 1-year and 10-year rates; the data on those go back the further thanks to Shiller.However, we may *also* have other rate data from FRED.But we need to have rate data for every year on the yield curve. That is: 1-, 2-, 3-, 4-, 5-, 6-, 7-, 9-, and 10-year rates. When we don't have the data available we will perform linear interpolation from data we *do* have to fill in the gaps.So if we only have the 1- and 10-year data then we need to do a linear interpolation for the other 8 years. If we have 1-, 3-, and 10-year data then we do linear interpolation between the 1- and 3-year data to fill in the 2-year data. And we'll do linear interpolation between the 3- and 10-year data for the rest. Missing data & flat yield curveWe can only do linear interpolation if we have two sets of data. We *always* have data for 1-year and 10-year rates, so we can fill in that part of the yield curve. But before the FRED data series of the mid-20th century we don't have any rates beyond 10 years. How do we fill that part of the rate yield curve when we don't have anything to interpolate?We don't interpolate: we just create a flat yield curve. That is, the 11-year rate is the same as the 10-year rate. And the 20-year rate is *also* the same as the 10-year rate. And the 30-year rate is **also** the same as the 10-year rate.This is obviously far from idea. So you should take with a large grain of salt any results before 1954. longinvest has [some comments on the longer terms starting approximately here][2]. Potential problems with linear interpolationThis linear interpolation is not perfect: it assumes that the yield curve is linear and that may not be the case. In particular, [look at this post from Fryxell][1] where he notes that before the 1920s the yield curve may have looked very different from what it did today.Still, trying to handle that is beyond the scope of this simulation. The more historical data (like those extra FRED data points) that we have, the less of this linear interpolation we need to do. That makes our post-1954 numbers better than the earlier numbers.[1]: https://www.bogleheads.org/forum/viewtopic.php?f=10&t=179425&start=100p2973643[2]: https://www.bogleheads.org/forum/viewtopic.php?f=10&t=179425&start=100p3013350<jupyter_code>def splice_data(raw_rates, series):
    # Start by loading the data we get from Shiller.
    # This will always exist.
    series.iloc[0] = raw_rates['1 year']
    series.iloc[9] = raw_rates['10 year']
    
    # Try to load any FRED rates.
    series.iloc[1] = raw_rates['GS2']
    series.iloc[2] = raw_rates['GS3']
    series.iloc[4] = raw_rates['GS5']
    series.iloc[6] = raw_rates['GS7']
    series.iloc[19] = raw_rates['GS20']
    series.iloc[29] = raw_rates['GS30']
    
    def safe_add(series_index, rate_index):
        # Don't overwrite any data we already have.
        if math.isnan(series.iloc[series_index]):
            series.iloc[series_index] = raw_rates[rate_index]

    # These are in order of preference. This is try to use M13058 before
    # trying to use M1333.
    safe_add(19, 'M13058')
    safe_add(19, 'M1333')
    
    # See the note below under "Going Beyond 30 Years" about how longinvest got these numbers
    safe_add(19, 'longinvest 20')
    safe_add(29, 'longinvest 30')

def build_yield_curve(raw_rates, yield_curve_size=30):
    s = pandas.Series(math.nan, index=numpy.arange(yield_curve_size))

    # We use NaN to indicate "the data needs to be interpolated"
    # We have a few different data series that we splice together.
    splice_data(raw_rates, s)
    
    def left_number(series, index):
        """ Find the index of first number to the left """
        if not math.isnan(series.iloc[index]):
            return index
        else:
            return left_number(series, index-1)
        
    def right_number(series, index):
        """ Find the index of the first number to the right """
        if not math.isnan(series.iloc[index]):
            return index
        else:
            return right_number(series, index+1)
                
    # now fill in the gaps with linear interpolation.
    for i in range(yield_curve_size):
        if math.isnan(s.iloc[i]):
            # First, try to find any existing data on the left and right.
            # We might not find any, for instance when we look beyond 10-years
            # before we have FRED data.

            try:
                left = left_number(s, i)
            except IndexError:
                left = None
                
            try:
                right = right_number(s, i)
            except IndexError:
                right = None

            if (left is None) and (right is None):
                raise IndexError("Couldn't find any rate data to fill out the yield curve.")

            if left is None:
                # If we can't find any data to the left then we can't do any linear interpolation
                # So just fill from the right
                s.iloc[i] = s.iloc[right]
            elif right is None:
                # If we can't find any data to the right then fill from the left
                # Both of these will result in a flat yield curve, which isn't ideal
                s.iloc[i] = s.iloc[left]
            else:
                # We can actually do linear interpolation
                steps = right - left
                rate = s.iloc[left] + ((s.iloc[right] - s.iloc[left]) * (i - left) / steps)
                s.iloc[i] = rate

    return s.tolist()
['%.2f' % (s*100) for s in build_yield_curve(HISTORICAL_RATES.iloc[0])]<jupyter_output><empty_output><jupyter_text>Putting it all togetherNow we have all the building blocks. We have a source of rates. We have a way to bootstrap our ladder. We have a way to see how the NAV changes over time.We only have one decision left to make -- what are the youngest & oldest maturities that we care about? Do we want a 10-2 fund? Or 10-4 fund? Or a 3-2 fund? Or how about a 7-4 fund?The maximum you can chose is 30 and the minimum is 1.Do you want a 10-2 fund, or 10-4 fund, or something else? That's actually done by the way you create the bootstrap ladder. This is how you build a 10-4 ladder.<jupyter_code>bootstrap(build_yield_curve(HISTORICAL_RATES.iloc[0]), 10, 4)
def loop(ladder, rates, max_maturity, start_year, end_year):
    df = pandas.DataFrame(columns=['NAV', 'Payments', 'Change'], index=numpy.arange(start_year, end_year + 1))

    for (year, current_rates) in rates:
        (ladder, payments, nav) = iterate_fund(ladder, build_yield_curve(current_rates), max_maturity)
        df.loc[year] = {'NAV' : nav, 'Payments' : payments}

    calculate_returns(df)
    return df

def calculate_returns(df):
    # Longinvest calculates the return based on comparison's to
    # next year's NAV. So I'll do the same. Even though that seems
    # weird to me. Maybe it's because the rates are based on January?
    # Hmmm...that sounds plausible.
    max_row = df.shape[0]

    for i in range(max_row - 1):
        next_nav = df.iloc[i+1]['NAV']
        nav = df.iloc[i]['NAV']
        change = (next_nav - nav) / nav
        df.iloc[i]['Change'] = change
    return df

def simulate(max_maturity, min_maturity, rates):
    """ This is just something to save on typing...and make clearer what the bounds on the fund are """
    ladder = bootstrap(build_yield_curve(rates.iloc[0]), max_maturity, min_maturity)
    start_year = int(rates.iloc[0].name)
    end_year = int(rates.iloc[-1].name)
    return loop(ladder, rates.iterrows(), max_maturity, start_year, end_year)<jupyter_output><empty_output><jupyter_text>Simulate All The ThingsNow we can simulate a 10-2 fund.<jupyter_code>simulate(10, 2, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>Or a 10-4 fund.<jupyter_code>simulate(10, 4, HISTORICAL_RATES).head()
simulate(4, 2, HISTORICAL_RATES).head()
simulate(3, 2, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>Holding to maturityOur previous examples all sell prior to bond maturity.If you want to hold until maturity, then specify a minimum duration of 1.<jupyter_code>simulate(10, 1, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>Only holding a single maturityIf you only want to fund to hold a single maturity, specify the same number twice. This constructs fund that only holds bonds of 2-year maturity. That is, every year it sells them and buys new 2-year bonds.<jupyter_code>simulate(2, 2, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>The "naive" approachOther bond return simulations don't hold a ladder, they just sell the bond after a single year. That is, buy a 10-year bond, sell it after 1-year and buy another 10-year bond. That's the same as holding a single maturity (above) but with a 10-year maturity.<jupyter_code>simulate(10, 10, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>Going beyond 30 yearsKeeping in mind the caveats above about results before we have FRED data from the mid-1950s...we can also build longer term bond funds. Warning: Extra interpolationThere are a few gaps in the FRED data. In the FRED 20-year series, there is no data for 1987-1993. (They weren't issued in those years, IIRC.)In the FRED 30-year series, there is no data for 2003-2006.longinvest does a linear interpolation **across years** to fill in the missing data. **I currently do not do this, so my results differ from those in the spreadsheet.**The way it works is:There is no 20-year data for 1987. But there *is* 30-year data. longinvest looks back at 1986 and sees that the spread between 20-year and 30-year in 1986 was -0.19%. Then he looks forward to 1994 and sees that the spread then was -0.10%.  Then he assumes the spread changed linearly over the missing years 1987-1993. That is:* 1987: -0.178%* 1988: -0.167%* 1989: -0.156%* 1990: -0.145%* 1991: -0.133%* 1992: -0.122%* 1993: -0.111%Then he applies that calculated spread to the known 30-year rates.When the 30-year rates are missing, a similar calculation is done based on the years 2002 & 2007 and the known 20-year rates.This is captured in the CSV as two new data series "longinvest 20" and "longinvest 30".<jupyter_code>simulate(30, 11, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>Crazy FundsYou can select any set of maturities you want, since the computer does all the heavy lifting. Want to simulate a fund that is **28-14** or **17-9**? That's easy.<jupyter_code>simulate(28, 14, HISTORICAL_RATES).head()
simulate(17, 9, HISTORICAL_RATES).head()<jupyter_output><empty_output><jupyter_text>Saving to CSVTo do anything useful, you probably want to save the results to a CSV file. Here's a commented-out example of how to do that.<jupyter_code>#simulate(10, 5, HISTORICAL_RATES).to_csv('10-5.csv')
import numpy
from collections import deque
import pandas
import math
import pandas_datareader.data as web
import datetime
import requests
import requests_cache
import xlrd
import tempfile

def get_morningstar(secid):
    url = 'http://mschart.morningstar.com/chartweb/defaultChart?type=getcc&secids=%s&dataid=8225&startdate=1900-01-01&enddate=2016-11-18&currency=&format=1' % secid
    expire_after = datetime.timedelta(days=3)
    session = requests_cache.CachedSession(cache_name='data-cache', backend='sqlite', expire_after=expire_after)

    # TODO: why doesn't this work!?!
    #r = session.get(url)
    r = requests.get(url)
    j = r.json()
    
    # The Morningstar data is pretty deeply nested....
    m = j['data']['r'][0]
    assert m['i'] == secid
    
    actual_data = m['t'][0]['d']
    # convert from strings to real data types
    as_dict = dict([(datetime.datetime.strptime(n['i'], '%Y-%m-%d'), float(n['v'])) for n in m['t'][0]['d']])
    
    # Strip out data?
    # Do we only want start of month, end of month, start of year, end of year, etc?
    s = pandas.Series(as_dict, name=secid)

    return s

barclays_index = get_morningstar('XIUSA000MJ')

# Use only final value for each calendar year
def annual(series):
    return series.groupby(by=lambda x: x.year).last()
# Use only final value for each calendar month
def monthly(series):
    return series.groupby(by=lambda x: datetime.date(x.year, x.month, 1)).last()

def calculate_change_prev(df, column):
    max_row = df.shape[0]
    
    series = pandas.Series()

    for i in range(max_row - 1):
        val = df.iloc[i][column]
        prev_val = df.iloc[i-1][column]
        change = (val - prev_val) / prev_val
        series.loc[df.iloc[i].name] = change
    return series

def calculate_change_next(df, column):
    max_row = df.shape[0]
    
    series = pandas.Series()

    for i in range(max_row - 1):
        val = df.iloc[i][column]
        next_val = df.iloc[i+1][column]
        change = (next_val - val) / val
        series.loc[df.iloc[i].name] = change
    return series

barclays_index = annual(barclays_index)
sim_10_4 = simulate(10, 4, HISTORICAL_RATES)
sim_10_10 = simulate(10, 10, HISTORICAL_RATES)
joined = pandas.concat([sim_10_4, barclays_index], axis=1, join='outer')
s_ind = calculate_change_prev(joined, 'XIUSA000MJ')
s_nav = calculate_change_next(joined, 'NAV')
joined = joined.assign(Change=s_nav, index_change=s_ind)
joined.to_csv('check.csv')
print(joined[["Change", "index_change"]].corr())
joined = pandas.concat([sim_10_10, barclays_index], axis=1, join='outer')
s_ind = calculate_change_prev(joined, 'XIUSA000MJ')
s_nav = calculate_change_next(joined, 'NAV')
joined = joined.assign(Change=s_nav, index_change=s_ind)
print(joined[["Change", "index_change"]].corr())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08980323530667123, 'fraction_numerical': 0.03269452388066072, 'mean_word_length': 3.710494571773221, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 10, 'lorem ipsum': 0, 'www.': 6, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '62609', 'license': 'BlueOak-1.0.0', 'n_tokens_mistral': 7777, 'n_tokens_neox': 6970, 'n_words': 2966, 'path': 'Bond Fund Simulator.ipynb', 'repo_name': 'hoostus/prime-harvesting', 'size': '23429'}