Fix Sabelhaus-Song variances

Documentation/CHANGELOG.md

@@ -28,7 +28,7 @@ interpolation for problems with CRRA utility. See [#888](https://github.com/econ
 * Adds a module for extracting initial distributions of permanent income (`pLvl`) and normalized assets (`aNrm`) from the SCF [#932](https://github.com/econ-ark/HARK/pull/932).
 * Fix the return fields of `dcegm/calcCrossPoints`[#909](https://github.com/econ-ark/HARK/pull/909).
 * Corrects location of constructor documentation to class string for Sphinx rendering [#908](https://github.com/econ-ark/HARK/pull/908)
-* Adds a module with tools for parsing and using various income calibrations from the literature. It includes the option of using life-cycle profiles of income shock variances from [Sabelhaus and Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358). See [#921](https://github.com/econ-ark/HARK/pull/921), [#941](https://github.com/econ-ark/HARK/pull/941).
+* Adds a module with tools for parsing and using various income calibrations from the literature. It includes the option of using life-cycle profiles of income shock variances from [Sabelhaus and Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358). See [#921](https://github.com/econ-ark/HARK/pull/921), [#941](https://github.com/econ-ark/HARK/pull/941), [#980](https://github.com/econ-ark/HARK/pull/980).
 * remove "Now" from model variable names [#936](https://github.com/econ-ark/HARK/pull/936)
 * remove Model.__call__; use Model init in Market and AgentType init to standardize on parameters dictionary [#947](https://github.com/econ-ark/HARK/issues/947)
 * Moves state MrkvNow to shocks['Mrkv'] in AggShockMarkov and KrusellSmith models [#935](https://github.com/econ-ark/HARK/pull/935)

HARK/Calibration/Income/IncomeTools.py

@@ -463,13 +463,13 @@ def sabelhaus_song_var_profile(age_min=27, age_max=54, cohort=None, smooth=True)
     # Construct variances
     # They use 1926 as the base year for cohort effects.
     ages = np.arange(age_min, age_max + 1)
-    tran_std = tran_dummy_interp(ages) + (cohort - 1926) * beta_eps
-    perm_std = perm_dummy_interp(ages) + (cohort - 1926) * beta_eta
+    tran_var = tran_dummy_interp(ages) + (cohort - 1926) * beta_eps
+    perm_var = perm_dummy_interp(ages) + (cohort - 1926) * beta_eta
 
     profiles = {
         "Age": list(ages),
-        "TranShkStd": list(tran_std),
-        "PermShkStd": list(perm_std),
+        "TranShkStd": list(np.sqrt(tran_var)),
+        "PermShkStd": list(np.sqrt(perm_var)),
     }
 
     return profiles

HARK/Calibration/Income/tests/test_IncomeTools.py

@@ -423,13 +423,13 @@ def test_paper_results(self):
 
         self.assertTrue(
             np.allclose(
-                self.Fig6Coh1940Tran, np.array(stds1940["TranShkStd"]), atol=1e-03
+                self.Fig6Coh1940Tran, np.array(stds1940["TranShkStd"])**2, atol=1e-03
             )
         )
 
         self.assertTrue(
             np.allclose(
-                self.Fig6Coh1940Perm, np.array(stds1940["PermShkStd"]), atol=1e-03
+                self.Fig6Coh1940Perm, np.array(stds1940["PermShkStd"])**2, atol=1e-03
             )
         )
 
@@ -440,13 +440,13 @@ def test_paper_results(self):
 
         self.assertTrue(
             np.allclose(
-                self.Fig6Coh1965Tran, np.array(stds1965["TranShkStd"]), atol=1e-03
+                self.Fig6Coh1965Tran, np.array(stds1965["TranShkStd"])**2, atol=1e-03
             )
         )
 
         self.assertTrue(
             np.allclose(
-                self.Fig6Coh1965Perm, np.array(stds1965["PermShkStd"]), atol=1e-03
+                self.Fig6Coh1965Perm, np.array(stds1965["PermShkStd"])**2, atol=1e-03
             )
         )
 
@@ -457,11 +457,11 @@ def test_aggregate_results(self):
         )
 
         self.assertTrue(
-            np.allclose(self.AggTran, np.array(stds_agg["TranShkStd"]), atol=1e-03)
+            np.allclose(self.AggTran, np.array(stds_agg["TranShkStd"])**2, atol=1e-03)
         )
 
         self.assertTrue(
-            np.allclose(self.AggPerm, np.array(stds_agg["PermShkStd"]), atol=1e-03)
+            np.allclose(self.AggPerm, np.array(stds_agg["PermShkStd"])**2, atol=1e-03)
         )
 
     def test_smoothing(self):
@@ -484,34 +484,34 @@ def test_smoothing(self):
         # 1940
         self.assertTrue(
             np.allclose(
-                np.array(smooth1940["TranShkStd"]), self.Fig6Coh1940Tran, rtol=rtol
+                np.array(smooth1940["TranShkStd"])**2, self.Fig6Coh1940Tran, rtol=rtol
             )
         )
 
         self.assertTrue(
             np.allclose(
-                np.array(smooth1940["PermShkStd"]), self.Fig6Coh1940Perm, atol=rtol
+                np.array(smooth1940["PermShkStd"])**2, self.Fig6Coh1940Perm, atol=rtol
             )
         )
 
         # 1965
         self.assertTrue(
             np.allclose(
-                np.array(smooth1965["TranShkStd"]), self.Fig6Coh1965Tran, rtol=rtol
+                np.array(smooth1965["TranShkStd"])**2, self.Fig6Coh1965Tran, rtol=rtol
             )
         )
 
         self.assertTrue(
             np.allclose(
-                np.array(smooth1965["PermShkStd"]), self.Fig6Coh1965Perm, atol=rtol
+                np.array(smooth1965["PermShkStd"])**2, self.Fig6Coh1965Perm, atol=rtol
             )
         )
 
         # Aggregate
         self.assertTrue(
-            np.allclose(np.array(smoothAgg["TranShkStd"]), self.AggTran, rtol=rtol)
+            np.allclose(np.array(smoothAgg["TranShkStd"])**2, self.AggTran, rtol=rtol)
         )
 
         self.assertTrue(
-            np.allclose(np.array(smoothAgg["PermShkStd"]), self.AggPerm, atol=rtol)
+            np.allclose(np.array(smoothAgg["PermShkStd"])**2, self.AggPerm, atol=rtol)
         )

HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py

@@ -35,18 +35,18 @@ def test_ConsIndShockSolverBasic(self):
         # test the solution_terminal
         self.assertAlmostEqual(LifecycleExample.solution[-1].cFunc(2).tolist(), 2)
 
-        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.89066194)
-        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.89144313)
-        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.89210133)
+        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.79429538)
+        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.79391692)
+        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.79253095)
 
         self.assertAlmostEqual(
-            LifecycleExample.solution[0].cFunc(1).tolist(), 0.8928547282397321
+            LifecycleExample.solution[0].cFunc(1).tolist(), 0.7506184692092213
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[1].cFunc(1).tolist(), 0.8930303445748624
+            LifecycleExample.solution[1].cFunc(1).tolist(), 0.7586358637239385
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[2].cFunc(1).tolist(), 0.8933075371183773
+            LifecycleExample.solution[2].cFunc(1).tolist(), 0.7681247572911291
         )
 
         solver = ConsIndShockSolverBasic(
@@ -66,22 +66,22 @@ def test_ConsIndShockSolverBasic(self):
         solver.prepare_to_solve()
 
         self.assertAlmostEqual(solver.DiscFacEff, 0.9586233599999999)
-        self.assertAlmostEqual(solver.PermShkMinNext, 0.9037554719886154)
+        self.assertAlmostEqual(solver.PermShkMinNext, 0.6554858756904397)
         self.assertAlmostEqual(solver.cFuncNowCnst(4).tolist(), 4.0)
-        self.assertAlmostEqual(solver.prepare_to_calc_EndOfPrdvP()[0], -0.2732742703949109)
-        self.assertAlmostEqual(solver.prepare_to_calc_EndOfPrdvP()[-1], 19.72572572960506)
+        self.assertAlmostEqual(solver.prepare_to_calc_EndOfPrdvP()[0], -0.19792871012285213)
+        self.assertAlmostEqual(solver.prepare_to_calc_EndOfPrdvP()[-1], 19.801071289877118)
 
         EndOfPrdvP = solver.calc_EndOfPrdvP()
 
-        self.assertAlmostEqual(EndOfPrdvP[0], 6710.672670733023)
-        self.assertAlmostEqual(EndOfPrdvP[-1], 0.14122987153089447)
+        self.assertAlmostEqual(EndOfPrdvP[0], 6657.839372100613)
+        self.assertAlmostEqual(EndOfPrdvP[-1], 0.2606075215645896)
 
         solution = solver.make_basic_solution(
             EndOfPrdvP, solver.aNrmNow, solver.make_linear_cFunc
         )
         solver.add_MPC_and_human_wealth(solution)
 
-        self.assertAlmostEqual(solution.cFunc(4).tolist(), 1.484118342351686)
+        self.assertAlmostEqual(solution.cFunc(4).tolist(), 1.0028005137373956)
 
     def test_simulated_values(self):
         self.agent.initialize_sim()

HARK/ConsumptionSaving/tests/test_IndShockConsumerTypeFast.py

@@ -39,18 +39,18 @@ def test_ConsIndShockSolverBasic(self):
         # test the solution_terminal
         self.assertAlmostEqual(LifecycleExample.solution[-1].cFunc(2).tolist(), 2)
 
-        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.89066194)
-        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.89144313)
-        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.89210133)
+        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.79429538)
+        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.79391692)
+        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.79253095)
 
         self.assertAlmostEqual(
-            LifecycleExample.solution[0].cFunc(1).tolist(), 0.8928547282397321
+            LifecycleExample.solution[0].cFunc(1).tolist(), 0.7506184692092213
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[1].cFunc(1).tolist(), 0.8930303445748624
+            LifecycleExample.solution[1].cFunc(1).tolist(), 0.7586358637239385
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[2].cFunc(1).tolist(), 0.8933075371183773
+            LifecycleExample.solution[2].cFunc(1).tolist(), 0.7681247572911291
         )
 
     def test_simulated_values(self):

examples/Calibration/Sabelhaus_Song_var_profiles.py

@@ -15,6 +15,7 @@
 
 import matplotlib.pyplot as plt
 from HARK.Calibration.Income.IncomeTools import sabelhaus_song_var_profile
+import numpy as np
 
 # Set up ages and cohorts at which we will get the variances
 age_min = 27
@@ -36,7 +37,7 @@
     coh_label = "aggregate" if cohorts[i] is None else cohorts[i]
     plt.plot(
         variances[i]["Age"],
-        variances[i]["TranShkStd"],
+        np.power(variances[i]["TranShkStd"],2),
         label="Tran. {} cohort".format(coh_label),
     )
 
@@ -49,7 +50,7 @@
     coh_label = "aggregate" if cohorts[i] is None else cohorts[i]
     plt.plot(
         variances[i]["Age"],
-        variances[i]["PermShkStd"],
+        np.power(variances[i]["PermShkStd"],2),
         label="Perm. {} cohort".format(coh_label),
     )