diff --git a/taxcalc/tests/test_tbi.py b/taxcalc/tests/test_tbi.py
index 7e8c309a8..8d8c9dab2 100644
--- a/taxcalc/tests/test_tbi.py
+++ b/taxcalc/tests/test_tbi.py
@@ -1,6 +1,7 @@
 """
 Test functions in taxcalc/tbi directory using both puf.csv and cps.csv input.
 """
+from __future__ import print_function
 import numpy as np
 import pandas as pd
 import pytest
@@ -262,3 +263,117 @@ def test_reform_warnings_errors():
     msg_dict = reform_warnings_errors(bad2_mods)
     assert len(msg_dict['warnings']) == 0
     assert len(msg_dict['errors']) > 0
+
+
+@pytest.mark.pre_release
+@pytest.mark.tbi_vs_std_behavior
+@pytest.mark.requires_pufcsv
+def test_behavioral_response(puf_subsample):
+    """
+    Test that behavioral-response results are the same
+    when generated from standard Tax-Calculator calls and
+    when generated from tbi.run_nth_year_tax_calc_model() calls
+    """
+    # specify reform and assumptions
+    reform_json = """
+    {"policy": {
+        "_II_rt5": {"2020": [0.25]},
+        "_II_rt6": {"2020": [0.25]},
+        "_II_rt7": {"2020": [0.25]},
+        "_PT_rt5": {"2020": [0.25]},
+        "_PT_rt6": {"2020": [0.25]},
+        "_PT_rt7": {"2020": [0.25]},
+        "_II_em": {"2020": [1000]}
+    }}
+    """
+    assump_json = """
+    {"behavior": {"_BE_sub": {"2013": [0.25]}},
+     "growdiff_baseline": {},
+     "growdiff_response": {},
+     "consumption": {}
+    }
+    """
+    params = Calculator.read_json_param_objects(reform_json, assump_json)
+    # specify keyword arguments used in tbi function call
+    kwargs = {
+        'start_year': 2019,
+        'year_n': 0,
+        'use_puf_not_cps': True,
+        'use_full_sample': False,
+        'user_mods': {
+            'policy': params['policy'],
+            'behavior': params['behavior'],
+            'growdiff_baseline': params['growdiff_baseline'],
+            'growdiff_response': params['growdiff_response'],
+            'consumption': params['consumption']
+        },
+        'return_dict': False
+    }
+    # generate aggregate results two ways: using tbi and standard calls
+    num_years = 9
+    std_res = dict()
+    tbi_res = dict()
+    for using_tbi in [True, False]:
+        for year in range(0, num_years):
+            cyr = year + kwargs['start_year']
+            if using_tbi:
+                kwargs['year_n'] = year
+                tables = run_nth_year_tax_calc_model(**kwargs)
+                tbi_res[cyr] = dict()
+                for tbl in ['aggr_1', 'aggr_2', 'aggr_d']:
+                    tbi_res[cyr][tbl] = tables[tbl]
+            else:
+                rec = Records(data=puf_subsample)
+                pol = Policy()
+                calc1 = Calculator(policy=pol, records=rec)
+                pol.implement_reform(params['policy'])
+                assert not pol.reform_errors
+                beh = Behavior()
+                beh.update_behavior(params['behavior'])
+                calc2 = Calculator(policy=pol, records=rec, behavior=beh)
+                assert calc2.behavior_has_response()
+                calc1.advance_to_year(cyr)
+                calc2.advance_to_year(cyr)
+                calc2 = Behavior.response(calc1, calc2)
+                std_res[cyr] = dict()
+                for tbl in ['aggr_1', 'aggr_2', 'aggr_d']:
+                    if tbl.endswith('_1'):
+                        itax = calc1.weighted_total('iitax')
+                        ptax = calc1.weighted_total('payrolltax')
+                        ctax = calc1.weighted_total('combined')
+                    elif tbl.endswith('_2'):
+                        itax = calc2.weighted_total('iitax')
+                        ptax = calc2.weighted_total('payrolltax')
+                        ctax = calc2.weighted_total('combined')
+                    elif tbl.endswith('_d'):
+                        itax = (calc2.weighted_total('iitax') -
+                                calc1.weighted_total('iitax'))
+                        ptax = (calc2.weighted_total('payrolltax') -
+                                calc1.weighted_total('payrolltax'))
+                        ctax = (calc2.weighted_total('combined') -
+                                calc1.weighted_total('combined'))
+                    cols = ['0_{}'.format(year)]
+                    rows = ['ind_tax', 'payroll_tax', 'combined_tax']
+                    datalist = [itax, ptax, ctax]
+                    std_res[cyr][tbl] = pd.DataFrame(data=datalist,
+                                                     index=rows,
+                                                     columns=cols)
+    # compare the two sets of results
+    # NOTE that the tbi results have been "fuzzed" for PUF privacy reasons,
+    #      so there is no expectation that the results should be identical.
+    no_diffs = True
+    reltol = 2e-3  # std and tbi differ if more than 0.2 percent different
+    for year in range(0, num_years):
+        cyr = year + kwargs['start_year']
+        col = '0_{}'.format(year)
+        for tbl in ['aggr_1', 'aggr_2', 'aggr_d']:
+            tbi = tbi_res[cyr][tbl][col]
+            std = std_res[cyr][tbl][col]
+            if not np.allclose(tbi, std, atol=0.0, rtol=reltol):
+                no_diffs = False
+                print('**** DIFF for year {} (year_n={}):'.format(cyr, year))
+                print('TBI RESULTS:')
+                print(tbi)
+                print('STD RESULTS:')
+                print(std)
+    assert no_diffs