Repository: econ-ark/HARK (language: python)

Issue:
`init_lifecycle` should be replaced with a real life cycle calibration
At present, the `init_lifecycle` calibration does not reflect a realistic description of a life cycle; it's just a 10 period finite-horizon model.  We should replace it with something based on the calibration in SolvingMicroDSOPs, but updated to be a more realistic calibration (the SolvingMicroDSOPs calibration dates from 1997, and a lot has been learned since then).   

I didn't go ahead and do this myself right now because:

1. `init_lifecycle` is imported via for a couple of tests, and so I want @sbenthall, who created those tests, and @mnwhite who wrote the original `init_lifecycle` to be involved in the rewrite.  
   * When done, we should also use the new `init_lifecycle` in the `LifeCycleModel` notebook in `HARK/examples`.  
1. I want to do a more serious job of parameterizing the model.  In particular, we need 
   * a better representation of the age profile of transitory and permanent shocks, from Sabelhaus and Song or Panousi et al or some other source
    * To incorporate aggregate growth
    * To construct a plausible calibration of distributions at age of workforce entry of:
        * Permanent income
        * Assets
    * To have a parallel calibration for a version of a `persistent shock` version of the model
        * lp_{t+1} = \gamma_{0} + \gamma_{1} lp_{t} + \psi_{t+1}
        * ly_{t+1} = lp_{t+1} + l\theta_{t+1}
        * where l indicates log, \psi and \theta are IID, and typical calibrations would have \gamma_{1} \approx 0.97

This would be a nice mini-project for an econ PhD student.
```
$ cd HARK ; grep -rs --include="*.py" 'init_life' * | grep import 
HARK/tests/test_simulate.py:from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType, init_lifecycle
HARK/ConsumptionSaving/tests/test_modelInits.py:from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType, init_idiosyncratic_shocks, init_lifecycle
```


Discussion / hints:
Is this a duplicate of #586 ?
It's an elaboration of #586 
I will start working on this. I will start out simple. Here is my plan (we'll see how far I get this winter).

- [ ] Start with the calibration from SolvingMicroDSOPs. Add that one to HARK and make sure it can be called by the examples in an easy enough way.
- [ ] Add the age-varying shock variances.
- [ ] Allow for different income profiles for different education levels.
- [ ] Get a reasonable calibration of assets at "birth". Could be pulled from the SCF, which I've worked with.

An initial question for @llorracc and @sbenthall would be where to put this? I feel like it would not work as well as a dictionary at the end of some model file. Ideally I see it as a function `getLifeCycleCalibration` with different options, e.g., educational level, do you want age-varying variances?, what age-range do you want to span?, etc.

Any initial thoughts on location and format would help.


As a first step, I am trying to collect all the various serious calibrations that have been mentioned. So far I have
- The one from [Micro DSOPs](https://github.com/econ-ark/REMARK/blob/master/REMARKs/SolvingMicroDSOPs/Calibration/EstimationParameters.py).
- The one from the [Pandemic dashboard](https://github.com/econ-ark/Pandemic/blob/master/Code/Python/Parameters.py).
- The one we pulled from [Cocco, Gomes, Maenhout (2005)](https://github.com/econ-ark/CGMPortfolio/blob/master/Code/Python/Calibration/params.py).

Is there any other one? 

@MridulS, in #586 @llorracc says that the age-varying variances from the Sabelhaus-Song paper were an option in the Portfolio Choice Blog Post, but I can't seem to find them in that repo. They are in the Pandemic file listed above, but I'm trying to make sure I collect as much info as I can.
I know that @mnwhite did a good bit of work to make it possible to incorporate the age-varying variances for the Sabelhaus-Song paper in the portfolio model when we were working on the blog post.  My recollection was that in the end they made less difference than I thought they might.  But in any case, if they are in the Pandemic calibration but not the blog post, it would be because Pandemic inherited from the work Matt did for Portfolio -- so Pandemic should be the source.  

PS.  One aspect of the SolvingMicroDSOPs calibration that is not "serious" is that we set the underlying rate of "aggregate productivity" growth to be 0 percent.  That should probably be more like 1 or 1.5 percent.  It doesn't matter much substantively -- the estimated time preference rate will basically move by exactly enough to offset it -- but that will be good because it should get rid of the almost eerily close to 1.0000 estimate of time time preference rate that comes out of the baseline SolvingMicroDSOPs calibration.  Because 1.0 is a round number, let's make that the default calibration.

> Any initial thoughts on location and format would help.
@llorracc says ConsIndShock and function.

I'm thinking about the volatility of income now. We want age-varying volatility of both transitory and permanent shocks.

You flagged Sabelhaus-Song as a starting point for this. I found a more recent paper by Moffitt on the matter: [Income Volatility and the PSID: Past Research and New Results](https://www.aeaweb.org/articles?id=10.1257/pandp.20181048). A drawback is that the volatility estimates are not at the age level, but grouped for ages [30-39], [40-49], and [50-59]. I wanted to put it up for consideration.

For now I'll continue trying to extract the results from Sabelhaus-Song.
I am now looking at survival probabilities. The [Pandemic paper](https://github.com/econ-ark/Pandemic/blob/c6feed4c4e4a3d7636d02b63c0d263ab801ed4a9/Code/Python/Parameters.py) had different survival probabilities by level of educational attainment. However, both the base rates and adjustments for education are loaded from .txt files that are not in the repo.

Does anyone remember where the numbers came from? @llorracc @MridulS 


> An initial question for @llorracc and @sbenthall would be where to put this?

I agree this is a little tricky.

There is pressure from two different sets of conventions:
 - "default" values are normally packaged with software classes so that they can be ran and tested easily
 - the parameter values here have scientific value; their provenance is important.

In my view, it would be better to think of these kinds of parameter collections as if they were included _data to be loaded_, rather than hard-coded into any Python file.

HARK has some utilities for the inclusion of datasets which follows scientific python conventions here:
https://github.com/econ-ark/HARK/tree/master/HARK/datasets

I'm not sure if we should use that dataset loading logic or something else like it for model configurations.

See also #878
SSA survival probabilities were included in https://github.com/econ-ark/HARK/pull/906.

I am working on similar modules that will allow us to use:
- Life-cycle profiles of income shock variances from Sabelhaus-Song.
- `aNrm` and `pLvl` distributions from the SCF.

PRs coming soon.
Get Chris' opinion on the Moffitt paper. I have a vague memory of him
disagreeing with the econometric approach, but that could be a phantom
memory, or a dream that I've misplaced into reality.

On Mon, Jan 18, 2021 at 3:38 PM Mateo Velásquez-Giraldo <
notifications@github.com> wrote:

> SSA survival probabilities were included in #906
> <https://github.com/econ-ark/HARK/pull/906>.
>
> I am working on similar modules that will allow us to use:
>
>    - Life-cycle profiles of income shock variances from Sabelhaus-Song.
>    - aNrm and pLvl distributions from the SCF.
>
> PRs coming soon.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/econ-ark/HARK/issues/592#issuecomment-762460607>, or
> unsubscribe
> <https://github.com/notifications/unsubscribe-auth/ADKRAFKEJ4K2PCEIUH7B3PDS2SL3TANCNFSM4LRJZSXQ>
> .
>

Thank you!

We discussed it. A cumbersome aspect of it is that it does not estimate volatilities for every year, but for age brackets `(30,40]`, `(40, 50]`, etc. iirc. We thought having volatility jumps when you change brackets would cause weird behavior, esp. around risk exposure (portfolio shares).

Also John Sabelhaus kindly provided us the original estimates for his paper, so we're going with that. Progress at #921.
>        * To construct a plausible calibration of distributions at age of workforce entry of:          
>          * Permanent income
>          * Assets
In progress at #922 .



Fix patch:
diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md
index 17999c3f4..a29c1fa66 100644
--- a/Documentation/CHANGELOG.md
+++ b/Documentation/CHANGELOG.md
@@ -30,6 +30,7 @@ interpolation for problems with CRRA utility. See [#888](https://github.com/econ
 * 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).
 * remove "Now" from model variable names [#936](https://github.com/econ-ark/HARK/pull/936)
 * Moves state MrkvNow to shocks['Mrkv'] in AggShockMarkov and KrusellSmith models [#935](https://github.com/econ-ark/HARK/pull/935)
+* Replaces `ConsIndShock`'s `init_lifecycle` with an actual life-cycle calibration [#951](https://github.com/econ-ark/HARK/pull/951).
 
 #### Minor Changes
 
diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index a69abce5e..c7a6f5583 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -49,6 +49,9 @@
 from HARK import _log
 from HARK import set_verbosity_level
 
+from HARK.Calibration.Income.IncomeTools import ParseIncomeSpec, ParseTimeParams, Cagetti_income
+from HARK.datasets.SCF.WealthIncomeDist.SCFDistTools import income_wealth_dists_from_scf
+from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table
 
 __all__ = [
     "ConsumerSolution",
@@ -2825,27 +2828,46 @@ def constructAssetsGrid(parameters):
 
 
 # Make a dictionary to specify a lifecycle consumer with a finite horizon
+
+# Main calibration characteristics
+birth_age = 25
+death_age = 90
+adjust_infl_to = 1992
+# Use income estimates from Cagetti (2003) for High-school graduates
+education = "HS"
+income_calib = Cagetti_income[education]
+
+# Income specification
+income_params = ParseIncomeSpec(
+    age_min=birth_age,
+    age_max=death_age,
+    adjust_infl_to=adjust_infl_to,
+    **income_calib,
+    SabelhausSong=True
+)
+
+# Initial distribution of wealth and permanent income
+dist_params = income_wealth_dists_from_scf(
+    base_year=adjust_infl_to, age=birth_age, education=education, wave=1995
+)
+
+# We need survival probabilities only up to death_age-1, because survival
+# probability at death_age is 1.
+liv_prb = parse_ssa_life_table(
+    female=False, cross_sec=True, year=2004, min_age=birth_age, max_age=death_age - 1
+)
+
+# Parameters related to the number of periods implied by the calibration
+time_params = ParseTimeParams(age_birth=birth_age, age_death=death_age)
+
+# Update all the new parameters
 init_lifecycle = copy(init_idiosyncratic_shocks)
-init_lifecycle["PermGroFac"] = [
-    1.01,
-    1.01,
-    1.01,
-    1.01,
-    1.01,
-    1.02,
-    1.02,
-    1.02,
-    1.02,
-    1.02,
-]
-init_lifecycle["PermShkStd"] = [0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0, 0, 0]
-init_lifecycle["TranShkStd"] = [0.3, 0.2, 0.1, 0.3, 0.2, 0.1, 0.3, 0, 0, 0]
-init_lifecycle["LivPrb"] = [0.99, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]
-init_lifecycle["T_cycle"] = 10
-init_lifecycle["T_retire"] = 7
-init_lifecycle[
-    "T_age"
-] = 11  # Make sure that old people die at terminal age and don't turn into newborns!
+init_lifecycle.update(time_params)
+init_lifecycle.update(dist_params)
+# Note the income specification overrides the pLvlInitMean from the SCF.
+init_lifecycle.update(income_params)
+init_lifecycle.update({"LivPrb": liv_prb})
+
 
 # Make a dictionary to specify an infinite consumer with a four period cycle
 init_cyclical = copy(init_idiosyncratic_shocks)


Test patch:
diff --git a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
index 07c51efbe..3b95eaa8c 100644
--- a/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
+++ b/HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
@@ -33,20 +33,20 @@ def test_ConsIndShockSolverBasic(self):
         LifecycleExample.solve()
 
         # test the solution_terminal
-        self.assertAlmostEqual(LifecycleExample.solution[10].cFunc(2).tolist(), 2)
+        self.assertAlmostEqual(LifecycleExample.solution[-1].cFunc(2).tolist(), 2)
 
-        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.97769632)
-        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.96624445)
-        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.95691449)
+        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[0].cFunc(1).tolist(), 0.87362789
+            LifecycleExample.solution[0].cFunc(1).tolist(), 0.8928547282397321
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[1].cFunc(1).tolist(), 0.9081621
+            LifecycleExample.solution[1].cFunc(1).tolist(), 0.8930303445748624
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[2].cFunc(1).tolist(), 0.9563899
+            LifecycleExample.solution[2].cFunc(1).tolist(), 0.8933075371183773
         )
 
         solver = ConsIndShockSolverBasic(
@@ -65,23 +65,23 @@ def test_ConsIndShockSolverBasic(self):
 
         solver.prepareToSolve()
 
-        self.assertAlmostEqual(solver.DiscFacEff, 0.9503999999999999)
-        self.assertAlmostEqual(solver.PermShkMinNext, 0.850430160026919)
+        self.assertAlmostEqual(solver.DiscFacEff, 0.9586233599999999)
+        self.assertAlmostEqual(solver.PermShkMinNext, 0.9037554719886154)
         self.assertAlmostEqual(solver.cFuncNowCnst(4).tolist(), 4.0)
-        self.assertAlmostEqual(solver.prepareToCalcEndOfPrdvP()[0], -0.2491750859108316)
-        self.assertAlmostEqual(solver.prepareToCalcEndOfPrdvP()[-1], 19.74982491408914)
+        self.assertAlmostEqual(solver.prepareToCalcEndOfPrdvP()[0], -0.2732742703949109)
+        self.assertAlmostEqual(solver.prepareToCalcEndOfPrdvP()[-1], 19.72572572960506)
 
         EndOfPrdvP = solver.calcEndOfPrdvP()
 
-        self.assertAlmostEqual(EndOfPrdvP[0], 6622.251864311334)
-        self.assertAlmostEqual(EndOfPrdvP[-1], 0.026301061207747087)
+        self.assertAlmostEqual(EndOfPrdvP[0], 6710.672670733023)
+        self.assertAlmostEqual(EndOfPrdvP[-1], 0.14122987153089447)
 
         solution = solver.makeBasicSolution(
             EndOfPrdvP, solver.aNrmNow, solver.makeLinearcFunc
         )
         solver.addMPCandHumanWealth(solution)
 
-        self.assertAlmostEqual(solution.cFunc(4).tolist(), 1.7391265696400773)
+        self.assertAlmostEqual(solution.cFunc(4).tolist(), 1.484118342351686)
 
     def test_simulated_values(self):
         self.agent.initializeSim()
diff --git a/HARK/ConsumptionSaving/tests/test_IndShockConsumerTypeFast.py b/HARK/ConsumptionSaving/tests/test_IndShockConsumerTypeFast.py
index e90f7ba24..d9538a7ec 100644
--- a/HARK/ConsumptionSaving/tests/test_IndShockConsumerTypeFast.py
+++ b/HARK/ConsumptionSaving/tests/test_IndShockConsumerTypeFast.py
@@ -37,20 +37,20 @@ def test_ConsIndShockSolverBasic(self):
         LifecycleExample.solve()
 
         # test the solution_terminal
-        self.assertAlmostEqual(LifecycleExample.solution[10].cFunc(2).tolist(), 2)
+        self.assertAlmostEqual(LifecycleExample.solution[-1].cFunc(2).tolist(), 2)
 
-        self.assertAlmostEqual(LifecycleExample.solution[9].cFunc(1), 0.97769632)
-        self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), 0.96624445)
-        self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), 0.95691449)
+        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[0].cFunc(1).tolist(), 0.87362789
+            LifecycleExample.solution[0].cFunc(1).tolist(), 0.8928547282397321
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[1].cFunc(1).tolist(), 0.9081621
+            LifecycleExample.solution[1].cFunc(1).tolist(), 0.8930303445748624
         )
         self.assertAlmostEqual(
-            LifecycleExample.solution[2].cFunc(1).tolist(), 0.9563899
+            LifecycleExample.solution[2].cFunc(1).tolist(), 0.8933075371183773
         )
 
     def test_simulated_values(self):