diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml
new file mode 100644
index 0000000..f537855
--- /dev/null
+++ b/.github/workflows/build_and_test.yml
@@ -0,0 +1,44 @@
+name: Build Package and Test Source Code [Python 3.7, 3.8, 3.9]
+
+on: [push, pull_request]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: [3.7, 3.8, 3.9]
+
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v2
+        with:
+          persist-credentials: false
+
+      - name: Setup Miniconda using Python ${{ matrix.python-version }}
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+          activate-environment: behresp-dev
+          environment-file: environment.yml
+          python-version: ${{ matrix.python-version }}
+          auto-activate-base: false
+
+      - name: Build
+        shell: bash -l {0}
+        run: |
+          pip install -e .
+          pip install pytest-cov
+          pip install pytest-pycodestyle
+      - name: Test
+        shell: bash -l {0}
+        working-directory: ./
+        run: |
+          pytest -m 'not requires_pufcsv and not pre_release and not local' --pycodestyle --cov=./ --cov-report=xml
+      - name: Upload coverage to Codecov
+        uses: codecov/codecov-action@v1
+        with:
+          file: ./coverage.xml
+          flags: unittests
+          name: codecov-umbrella
+          fail_ci_if_error: true
+          verbose: true
diff --git a/.travis.yml b/.travis.yml
index 194c1ed..ec62fb5 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,9 +1,9 @@
 dist: xenial
 language: python
 python:
-  - "3.6"
   - "3.7"
   - "3.8"
+  - "3.9"
 
 install:
   # Install conda
@@ -15,7 +15,7 @@ install:
   - conda create -n behresp-dev python=$TRAVIS_PYTHON_VERSION;
   - source activate behresp-dev
   - conda env update -f environment.yml
-  - pip install pytest-pycodestyle 
+  - pip install pytest-pycodestyle
   - pip install pyyaml
   - pip install coverage
   - pip install codecov
diff --git a/behresp/behavior.py b/behresp/behavior.py
index aa82028..4d900b1 100644
--- a/behresp/behavior.py
+++ b/behresp/behavior.py
@@ -10,20 +10,23 @@
 import taxcalc as tc
 
 
-def response(calc_1, calc_2, elasticities, dump=False):
+def response(calc_1, calc_2, elasticities, dump=False, inplace=False):
     """
     Implements TaxBrain "Partial Equilibrium Simulation" dynamic analysis
     returning results as a tuple of Pandas DataFrame objects (df1, df2) where:
-    df1 is extracted from a baseline-policy calc_1 copy, and
-    df2 is extracted from a reform-policy calc_2 copy that incorporates the
+    df1 is extracted from a baseline-policy calc_1, and
+    df2 is extracted from a reform-policy calc_2 that incorporates the
         behavioral responses given by the nature of the baseline-to-reform
         change in policy and elasticities in the specified behavior dictionary.
 
-    Note: this function internally modifies a copy of calc_2 records to account
-      for behavioral responses that arise from the policy reform that involves
-      moving from calc1 policy to calc2 policy.  Neither calc_1 nor calc_2 need
-      to have had calc_all() executed before calling the response function.
-      And neither calc_1 nor calc_2 are affected by this response function.
+    Note: By default, this function internally modifies a copy of calc_2
+      records to account for behavioral responses that arise from the
+      policy reform that involves moving from calc1 policy to calc2
+      policy.  Neither calc_1 nor calc_2 need to have had calc_all()
+      executed before calling the response function. By default, neither
+      calc_1 nor calc_2 are affected by this response function. To
+      perform in-place calculations that affect calc_1 and calc_2, set
+      inplace equal to True.
 
     The elasticities argument is a dictionary containing the assumed response
     elasticities.  Omitting an elasticity key:value pair in the dictionary
@@ -93,8 +96,12 @@ def response(calc_1, calc_2, elasticities, dump=False):
     # pylint: disable=too-many-locals,too-many-statements,too-many-branches
 
     # Check function argument types and elasticity values
-    calc1 = copy.deepcopy(calc_1)
-    calc2 = copy.deepcopy(calc_2)
+    if inplace:
+        calc1 = calc_1
+        calc2 = calc_2
+    else:
+        calc1 = copy.deepcopy(calc_1)
+        calc2 = copy.deepcopy(calc_2)
     assert isinstance(calc1, tc.Calculator)
     assert isinstance(calc2, tc.Calculator)
     assert isinstance(elasticities, dict)
@@ -220,10 +227,14 @@ def _mtr12(calc__1, calc__2, mtr_of='e00200p', tax_type='combined'):
         df1['mtr_combined'] = wage_mtr1 * 100
     else:
         df1 = calc1.dataframe(tc.DIST_VARIABLES)
-    del calc1
+    if not inplace:
+        del calc1
     # Add behavioral-response changes to income sources
-    calc2_behv = copy.deepcopy(calc2)
-    del calc2
+    if inplace:
+        calc2_behv = calc2
+    else:
+        calc2_behv = copy.deepcopy(calc2)
+        del calc2
     if not zero_sub_and_inc:
         calc2_behv = _update_ordinary_income(si_chg, calc2_behv)
     calc2_behv = _update_cap_gain_income(ltcg_chg, calc2_behv)
@@ -237,7 +248,8 @@ def _mtr12(calc__1, calc__2, mtr_of='e00200p', tax_type='combined'):
         df2['mtr_combined'] = wage_mtr2 * 100
     else:
         df2 = calc2_behv.dataframe(tc.DIST_VARIABLES)
-    del calc2_behv
+    if not inplace:
+        del calc2_behv
     # Return the two dataframes
     return (df1, df2)
 
diff --git a/behresp/tests/test_behavior.py b/behresp/tests/test_behavior.py
index d55b5bc..254dfd7 100644
--- a/behresp/tests/test_behavior.py
+++ b/behresp/tests/test_behavior.py
@@ -13,7 +13,8 @@
 from behresp import response, quantity_response, labor_response
 
 
-def test_default_response_function(cps_subsample):
+@pytest.mark.parametrize("inplace", [False, True])
+def test_default_response_function(cps_subsample, inplace):
     """
     Test that default behavior parameters produce static results.
     """
@@ -37,10 +38,32 @@ def test_default_response_function(cps_subsample):
     calc2s.calc_all()
     df2s = calc2s.dataframe(['iitax', 's006'])
     itax2s = round((df2s['iitax'] * df2s['s006']).sum() * 1e-9, 3)
+
+    # Keep track of some of the variables that will be modifed
+    # in response if inplace is True.
+    df_before_1 = calc1.dataframe(tc.DIST_VARIABLES)
+    df_before_2 = calc2d.dataframe(tc.DIST_VARIABLES)
+
     # ... calculate aggregate inctax using zero response elasticities
-    _, df2d = response(calc1, calc2d, elasticities={}, dump=True)
+    _, df2d = response(calc1, calc2d, elasticities={}, dump=True,
+                       inplace=inplace)
     itax2d = round((df2d['iitax'] * df2d['s006']).sum() * 1e-9, 3)
     assert np.allclose(itax2d, itax2s)
+
+    # Grab the same variables after.
+    df_after_1 = calc1.dataframe(tc.DIST_VARIABLES)
+    df_after_2 = calc2d.dataframe(tc.DIST_VARIABLES)
+
+    # If inplace is True, the variables should have been modified.
+    # If inplace is False, response only modified a copy of calc_1
+    # and calc_2.
+    if inplace:
+        assert not df_before_1.equals(df_after_1)
+        assert not df_before_2.equals(df_after_2)
+    else:
+        assert df_before_1.equals(df_after_1)
+        assert df_before_2.equals(df_after_2)
+
     # ... clean up
     del calc1
     del calc2s
@@ -76,9 +99,9 @@ def test_nondefault_response_function(be_inc, cps_subsample):
     del df1
     del df2
     if be_inc == 0.0:
-        assert np.allclose([itax1, itax2], [1355.556, 1304.984])
+        assert np.allclose([itax1, itax2], [1354.7, 1304.166])
     elif be_inc == -0.1:
-        assert np.allclose([itax1, itax2], [1355.556, 1303.898])
+        assert np.allclose([itax1, itax2], [1354.7, 1303.08])
 
 
 def test_alternative_behavior_parameters(cps_subsample):
@@ -107,7 +130,7 @@ def test_alternative_behavior_parameters(cps_subsample):
     itax2 = round((df2['iitax'] * df2['s006']).sum() * 1e-9, 3)
     del df1
     del df2
-    assert np.allclose([itax1, itax2], [1355.556, 1302.09])
+    assert np.allclose([itax1, itax2], [1354.7, 1301.281])
 
 
 def test_quantity_response():
diff --git a/conda.recipe/meta.yaml b/conda.recipe/meta.yaml
index 5ab9da4..30b8b1b 100755
--- a/conda.recipe/meta.yaml
+++ b/conda.recipe/meta.yaml
@@ -5,12 +5,12 @@ package:
 requirements:
   build:
     - python
-    - "taxcalc>=3.0.0"
+    - "taxcalc>=3.1.0"
     - pytest
 
   run:
     - python
-    - "taxcalc>=3.0.0"
+    - "taxcalc>=3.1.0"
     - pytest
 
 test:
diff --git a/environment.yml b/environment.yml
index c1470d9..26c2686 100644
--- a/environment.yml
+++ b/environment.yml
@@ -1,8 +1,7 @@
 name: behresp-dev
 channels:
-- PSLmodels
 - conda-forge
 dependencies:
 - python
-- "taxcalc>=3.0.0"
+- "taxcalc>=3.1.0"
 - pytest