def history(csvfile):
    df = pd.read_csv(csvfile, index_col='Date', parse_dates=True)
    df = df.dropna()
    df.index = df.index.to_period('M').to_timestamp('M')
    display(df.head(1).append(df.tail(1)))
    annual_price = df.resample('Y').agg({
        'Open': 'first',
        'Adj Close': 'last',
    })
    monthly_price = df[['Open', 'Adj Close']]
    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(14, 4))
    annual_price['Adj Close'].diff().cumsum().plot(ax=axes[0])
    annual_price.loc['2000':'2019', 'Adj Close'].diff().cumsum().plot(ax=axes[1])
    return [
        r.pct_change(axis=1).dropna(axis=1).rename(columns={'Adj Close': 'value'}) 
        for r in (annual_price, monthly_price)
    ]

annual_return, monthly_return = history('./^GSPC.csv')
display(annual_return['value'].mean())
display(annual_return.loc['2000':'2019', 'value'].mean())
plt.show(block=True)

0.07677447259284405

0.05639410820846351

fund_df = pd.DataFrame({
    'SBIバンガードS&P500': [0.000938],
    'eMAXIS Slim米国株式(S&P500)': [0.000968],
    '農林中金<パートナーズ>つみたてNISA米国株式S&P500': [0.00495],
}, index=['信託報酬年率']).T
display(fund_df)

def calculate(return_df, n_div=1, freq='Y'):
    base_dfs = []
    for year_limit, total_limit in (
        (1_200_000, 6_000_000),  # 一般NISA年間上限金額 / 最大金額
        (  400_000, 8_000_000),  # つみたてNISA年間上限金額 / 最大金額
    ):
        df = pd.DataFrame({
            'base': [year_limit / n_div] * len(return_df),
        }, index=return_df.index)
        df = df[df.cumsum() <= total_limit].fillna(0)
        base_dfs.append(df)
    fund_dfs = []
    for fund_name, fee in fund_df.itertuples():
        dfs = []
        for base_df in base_dfs:
            df = pd.DataFrame({'earn': 0}, index=base_df.index)
            for i, index in enumerate(df.index):
                base = base_df.loc[index, 'base']
                earn = df.loc[df.index[i - 1], 'earn']
                ret = return_df.loc[index, 'value']
                df.loc[index, 'earn'] = \
                    (base + earn) * (1 + ret) / (1 + fee / n_div)
            dfs.append(df.astype('int64'))
        df = pd.concat(dfs, keys=['一般', 'つみたて'], axis=1)
        fund_dfs.append(df)
    base_df = pd.concat(
        [df.cumsum() for df in base_dfs],
        keys=fund_dfs[0].columns.get_level_values(0),
        axis=1).astype('int64')
    df = pd.concat(
        [base_df] + fund_dfs,
        keys=['', 'SBIバンガード', 'eMAXIS', '農林中金'], 
        axis=1)
    display(df)
    plot_df = df.copy()
    plot_df.columns = df.columns.to_flat_index()
    ax = plot_df.plot(
        xticks=df.index[::n_div],
        figsize=(12, 6),
        markevery=n_div,
        style=['yo--', 'cs--'])
    ax.xaxis.set_major_formatter(FixedFormatter(df.index.year.unique()))
    plt.show(block=True)
    return df

annual_df = calculate(return_df=annual_return['2000':'2019'])

annual_df[annual_df[('', 'つみたて', 'base')] >= 6000000]

last = annual_df.last('Y')
base, earn = last[''], last['SBIバンガード']
base.columns = earn.columns = ['一般NISA投資利益率', 'つみたてNISA投資利益率']
display(earn.T / base.T)

pd.concat([
    annual_df.last('Y').loc[:, pd.IndexSlice[:, c, 'earn']].diff(axis=1).fillna(0)
    for c in annual_df.columns.get_level_values(1).unique()
], axis=1)

monthly_df = calculate(return_df=monthly_return['2000':'2019'], n_div=12, freq='M')

display(monthly_df[monthly_df[('', 'つみたて', 'base')] == 5999999])

annual_last = annual_df.last('Y')
monthly_last = monthly_df.last('M')
pd.concat([annual_last, monthly_last, annual_last - monthly_last])