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

Fix patch:
diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md
index bff8479b8..db2c01d82 100644
--- a/Documentation/CHANGELOG.md
+++ b/Documentation/CHANGELOG.md
@@ -23,7 +23,7 @@ interpolation for problems with CRRA utility. See [#888](https://github.com/econ
 * Adds a SSA life tables and methods to extract survival probabilities from them [#986](https://github.com/econ-ark/HARK/pull/906).
 * 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 for producing 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).
 #### Minor Changes
 
 * Move AgentType constructor parameters docs to class docstring so it is rendered by Sphinx.
diff --git a/HARK/Calibration/Income/README.md b/HARK/Calibration/Income/README.md
new file mode 100644
index 000000000..e83abe6cd
--- /dev/null
+++ b/HARK/Calibration/Income/README.md
@@ -0,0 +1,9 @@
+# Tools for calibrating income profiles
+
+This folder contains various tools that are used for replicating calibrations/estimates of income processes from the literature.
+
+## `SabelhausSongProfiles.py`
+
+This file contains tools for producing life-cycle profiles of the volatilities of transitory and permanent shocks to income using the results from [Sabelhaus & Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358).
+
+The estimates used to reproduce the results from the paper were generously provided by [John Sabelhaus](https://www.sites.google.com/view/johnsabelhaus).
diff --git a/HARK/Calibration/Income/SabelhausSongProfiles.py b/HARK/Calibration/Income/SabelhausSongProfiles.py
new file mode 100644
index 000000000..a9b145dce
--- /dev/null
+++ b/HARK/Calibration/Income/SabelhausSongProfiles.py
@@ -0,0 +1,303 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 14 15:49:43 2021
+
+@author: Mateo
+
+This scripts contains representations of the results in:
+
+[1] Sabelhaus, J., & Song, J. (2010). The great moderation in micro labor
+    earnings. Journal of Monetary Economics, 57(4), 391-403.
+
+It provides functions to produce life-cycle profiles of the variances of 
+income shocks from the raw results from the paper.
+
+The raw estimates were generously shared by John Sabelhaus.
+
+"""
+
+import numpy as np
+from warnings import warn
+from HARK.interpolation import LinearInterp
+
+
+# The raw results shared by John Sabelhaus contain the following two
+# sets of estimates (with and without cohor trends), which we will
+# use for constructing the age profiles.
+
+# The first specification contains a cohort trend. The variance of
+# (transitory or permanent) shocks to income of a person born in year
+# "cohort" and who is now age "age" is
+# age_dummy(age) + beta * (cohort - 1926)
+# Where we have dummies for ages 27 to 54
+Sabelhaus_Song_cohort_trend = {
+    "Ages": np.arange(27, 55),
+    "AgeDummiesPrm": np.array(
+        [
+            0.0837941,
+            0.0706855,
+            0.0638561,
+            0.0603879,
+            0.0554693,
+            0.0532388,
+            0.0515262,
+            0.0486079,
+            0.0455297,
+            0.0456573,
+            0.0417433,
+            0.0420146,
+            0.0391508,
+            0.0395776,
+            0.0369826,
+            0.0387158,
+            0.0365356,
+            0.036701,
+            0.0364236,
+            0.0358601,
+            0.0348528,
+            0.0362901,
+            0.0373366,
+            0.0372724,
+            0.0401297,
+            0.0415868,
+            0.0434772,
+            0.046668,
+        ]
+    ),
+    "AgeDummiesTrn": np.array(
+        [
+            0.1412842,
+            0.1477754,
+            0.1510265,
+            0.1512203,
+            0.1516837,
+            0.151412,
+            0.1489388,
+            0.148521,
+            0.1470632,
+            0.143514,
+            0.1411806,
+            0.1378733,
+            0.135245,
+            0.1318365,
+            0.1299689,
+            0.1255799,
+            0.1220823,
+            0.1178995,
+            0.1148793,
+            0.1107577,
+            0.1073337,
+            0.102347,
+            0.0962066,
+            0.0918819,
+            0.0887777,
+            0.0835057,
+            0.0766663,
+            0.0698848,
+        ]
+    ),
+    "CohortCoefPrm": -0.0005966,
+    "CohortCoefTrn": -0.0017764 / 2,
+}
+
+# The second specification contains no cohort trend. The variance of
+# (transitory or permanent) shocks to income of a person born in year
+# "cohort" and who is now age "age" is: age_dummy(age)
+# Where we have dummies for ages 27 to 54. We use this "aggregate"
+# specification if no cohort is provided.
+Sabelhaus_Song_all_years = {
+    "Ages": np.arange(27, 55),
+    "AgeDummiesPrm": np.array(
+        [
+            0.0599296,
+            0.0474176,
+            0.0411848,
+            0.0383132,
+            0.0339912,
+            0.0323573,
+            0.0312414,
+            0.0289196,
+            0.0264381,
+            0.0271623,
+            0.0238449,
+            0.0247128,
+            0.0224456,
+            0.0234691,
+            0.0214706,
+            0.0238005,
+            0.0222169,
+            0.0229789,
+            0.0232982,
+            0.0233312,
+            0.0229205,
+            0.0249545,
+            0.0265975,
+            0.02713,
+            0.0305839,
+            0.0326376,
+            0.0351246,
+            0.038912,
+        ]
+    ),
+    "AgeDummiesTrn": np.array(
+        [
+            0.1061999,
+            0.1135794,
+            0.1177187,
+            0.1188007,
+            0.1201523,
+            0.1207688,
+            0.1191838,
+            0.1196542,
+            0.1190846,
+            0.1164236,
+            0.1149784,
+            0.1125594,
+            0.1108192,
+            0.1082989,
+            0.1073195,
+            0.1038188,
+            0.1012094,
+            0.0979148,
+            0.0957829,
+            0.0925495,
+            0.0900136,
+            0.0859151,
+            0.0806629,
+            0.0772264,
+            0.0750105,
+            0.0706267,
+            0.0646755,
+            0.0587821,
+        ]
+    ),
+    "CohortCoefPrm": 0,
+    "CohortCoefTrn": 0,
+}
+
+
+def sabelhaus_song_var_profile(age_min=27, age_max=54, cohort=None,
+                               smooth = True):
+    """
+    This is a function to find the life-cycle profiles of the volatilities
+    of transitory and permanent shocks to income using the estimates in 
+    [1] Sabelhaus and Song (2010).
+
+    Parameters
+    ----------
+    age_min : int, optional
+        Minimum age at which to construct volatilities. The default is 27.
+    age_max : int, optional
+        Maximum age at which to construct volatilities. The default is 54.
+    cohort : int, optional
+        Birth year of the hypothetical person for which the volatilities will
+        be constructed. The default is None, and in this case the we will
+        use the specification that does not have cohort trends.
+    smooth: bool, optional
+        Boolean indicating whether to smooth the variance profile estimates
+        using third degree polynomials for the age dummies estimated by
+        Sabelhaus and Song. If False, the original dummies are used.
+        
+    Returns
+    -------
+    profiles : dict
+        Dictionary with entries:
+            - Ages: array of ages for which we found income volatilities in
+                ascending order
+            - TranShkStd: array of standard deviations of transitory income
+                shocks. Position n corresponds to Ages[n].
+            - PermShkStd: array of standard deviations of permanent income
+                shocks. Position n corresponds to Ages[n].
+        
+        Note that TransShkStd[n] and PermShkStd[n] are the volatilities of
+        shocks _experienced_ at age Age[n], (not those expected at Age[n+1]
+        from the perspective of Age[n]).
+        
+        Note that Sabelhaus and Song work in discrete time and with periods
+        that represent one year. Therefore, the outputs must be interpreted
+        at the yearly frequency.
+    """
+
+    assert age_max >= age_min, (
+        "The maximum age can not be lower than the " + "minimum age."
+    )
+
+    # Determine which set of estimates to use based on wether a cohort is
+    # provided or not.
+    if cohort is None:
+
+        spec = Sabelhaus_Song_all_years
+        cohort = 0
+        warn("No cohort was provided. Using aggregate specification.")
+
+    else:
+
+        spec = Sabelhaus_Song_cohort_trend
+
+    # Extract coefficients
+    beta_eps = spec["CohortCoefTrn"]
+    beta_eta = spec["CohortCoefPrm"]
+    tran_age_dummies = spec["AgeDummiesTrn"]
+    perm_age_dummies = spec["AgeDummiesPrm"]
+
+    # Smooth out dummies using a 3rd degree polynomial if requested
+    if smooth:
+        
+        # Fit polynomials
+        tran_poly = np.poly1d(np.polyfit(spec['Ages'],
+                                         tran_age_dummies, deg = 3))
+        perm_poly = np.poly1d(np.polyfit(spec['Ages'],
+                                         perm_age_dummies, deg = 3))
+        
+        # Replace dummies
+        tran_age_dummies = tran_poly(spec['Ages'])
+        perm_age_dummies = perm_poly(spec['Ages'])
+        
+    # Make interpolators for transitory and permanent dummies. Alter to use
+    # flat extrapolation.
+
+    # We use Sabelhaus and Song (2010) dummies for ages 27-54 and extrapolate
+    # outside of that just using the endpoints.
+    tran_dummy_interp = LinearInterp(
+        np.arange(min(spec['Ages'])-1, max(spec['Ages'])+2),
+        np.concatenate(
+            [[tran_age_dummies[0]], tran_age_dummies, [tran_age_dummies[-1]]]
+        ),
+        lower_extrap=True,
+    )
+
+    perm_dummy_interp = LinearInterp(
+        np.arange(min(spec['Ages'])-1, max(spec['Ages'])+2),
+        np.concatenate(
+            [[perm_age_dummies[0]], perm_age_dummies, [perm_age_dummies[-1]]]
+        ),
+        lower_extrap=True,
+    )
+
+    if age_min < 27 or age_max > 54:
+        warn(
+            "Sabelhaus and Song (2010) provide variance profiles for ages "
+            + "27 to 54. Extrapolating variances using the extreme points."
+        )
+
+    if cohort < 1926 or cohort > 1980:
+        warn(
+            "Sabelhaus and Song (2010) use data from birth cohorts "
+            + "[1926,1980]. Extrapolating variances."
+        )
+
+        cohort = max(min(cohort, 1980), 1926)
+
+    # 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
+
+    profiles = {
+        "Age": list(ages),
+        "TranShkStd": list(tran_std),
+        "PermShkStd": list(perm_std),
+    }
+
+    return profiles
diff --git a/HARK/Calibration/Income/__init__.py b/HARK/Calibration/Income/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/HARK/Calibration/__init__.py b/HARK/Calibration/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/examples/Calibration/Sabelhaus_Song_var_profiles.py b/examples/Calibration/Sabelhaus_Song_var_profiles.py
new file mode 100644
index 000000000..636fc9b2c
--- /dev/null
+++ b/examples/Calibration/Sabelhaus_Song_var_profiles.py
@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 14 16:44:09 2021
+
+@author: Mateo
+
+This short script demonstrates how to use the module for computing
+[1] Sabelhaus & Song (2010) age profiles of income volatility.
+It does so by replicating the results from the original paper (Figure 6 in [1])
+
+[1] Sabelhaus, J., & Song, J. (2010). The great moderation in micro labor
+    earnings. Journal of Monetary Economics, 57(4), 391-403.
+
+"""
+
+import matplotlib.pyplot as plt
+from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile
+
+# Set up ages and cohorts at which we will get the variances
+age_min = 27
+age_max = 54
+cohorts = [1940, 1965, None]
+
+# Find volatility profiles using the module
+variances = [
+    sabelhaus_song_var_profile(age_min=age_min, age_max=age_max, cohort=c)
+    for c in cohorts
+]
+
+# %% Plots
+
+# Plot transitory shock variances
+plt.figure()
+for i in range(len(cohorts)):
+
+    coh_label = "aggregate" if cohorts[i] is None else cohorts[i]
+    plt.plot(
+        variances[i]["Age"],
+        variances[i]["TranShkStd"],
+        label="Tran. {} cohort".format(coh_label),
+    )
+
+plt.legend()
+
+# Plot permanent shock variances
+plt.figure()
+for i in range(len(cohorts)):
+
+    coh_label = "aggregate" if cohorts[i] is None else cohorts[i]
+    plt.plot(
+        variances[i]["Age"],
+        variances[i]["PermShkStd"],
+        label="Perm. {} cohort".format(coh_label),
+    )
+
+plt.legend()


Test patch:
diff --git a/HARK/Calibration/Income/tests/__init__.py b/HARK/Calibration/Income/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/HARK/Calibration/Income/tests/test_SabelhausSongProfiles.py b/HARK/Calibration/Income/tests/test_SabelhausSongProfiles.py
new file mode 100644
index 000000000..9e781d3ca
--- /dev/null
+++ b/HARK/Calibration/Income/tests/test_SabelhausSongProfiles.py
@@ -0,0 +1,320 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 21 15:04:24 2021
+
+@author: Mateo
+"""
+
+# Bring in modules we need
+import unittest
+import numpy as np
+from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile
+
+
+class test_SabelhausSongProfiles(unittest.TestCase):
+    def setUp(self):
+
+        # Write results from Figure 6 in the original paper
+        self.Fig6Coh1940Tran = np.array(
+            [
+                0.1288,
+                0.1353,
+                0.1386,
+                0.1388,
+                0.1392,
+                0.1390,
+                0.1365,
+                0.1361,
+                0.1346,
+                0.1311,
+                0.1287,
+                0.1254,
+                0.1228,
+                0.1194,
+                0.1175,
+                0.1131,
+                0.1096,
+                0.1055,
+                0.1024,
+                0.0983,
+                0.0949,
+                0.0899,
+                0.0838,
+                0.0794,
+                0.0763,
+                0.0711,
+                0.0642,
+                0.0575,
+            ]
+        )
+
+        self.Fig6Coh1940Perm = np.array(
+            [
+                0.0754,
+                0.0623,
+                0.0555,
+                0.0520,
+                0.0471,
+                0.0449,
+                0.0432,
+                0.0403,
+                0.0372,
+                0.0373,
+                0.0334,
+                0.0337,
+                0.0308,
+                0.0312,
+                0.0286,
+                0.0304,
+                0.0282,
+                0.0283,
+                0.0281,
+                0.0275,
+                0.0265,
+                0.0279,
+                0.0290,
+                0.0289,
+                0.0318,
+                0.0332,
+                0.0351,
+                0.0383,
+            ]
+        )
+
+        self.Fig6Coh1965Tran = np.array(
+            [
+                0.1066,
+                0.1131,
+                0.1164,
+                0.1166,
+                0.1170,
+                0.1168,
+                0.1143,
+                0.1139,
+                0.1124,
+                0.1089,
+                0.1065,
+                0.1032,
+                0.1006,
+                0.0972,
+                0.0953,
+                0.0909,
+                0.0874,
+                0.0833,
+                0.0802,
+                0.0761,
+                0.0727,
+                0.0677,
+                0.0616,
+                0.0572,
+                0.0541,
+                0.0489,
+                0.0420,
+                0.0352,
+            ]
+        )
+        self.Fig6Coh1965Perm = np.array(
+            [
+                0.0605,
+                0.0474,
+                0.0406,
+                0.0371,
+                0.0322,
+                0.0300,
+                0.0283,
+                0.0253,
+                0.0223,
+                0.0224,
+                0.0185,
+                0.0187,
+                0.0159,
+                0.0163,
+                0.0137,
+                0.0154,
+                0.0133,
+                0.0134,
+                0.0132,
+                0.0126,
+                0.0116,
+                0.0130,
+                0.0141,
+                0.0140,
+                0.0169,
+                0.0183,
+                0.0202,
+                0.0234,
+            ]
+        )
+
+        # Aggregate result from Sabelhaus' excel file
+        self.AggTran = np.array(
+            [
+                0.1062,
+                0.1136,
+                0.1177,
+                0.1188,
+                0.1202,
+                0.1208,
+                0.1192,
+                0.1197,
+                0.1191,
+                0.1164,
+                0.1150,
+                0.1126,
+                0.1108,
+                0.1083,
+                0.1073,
+                0.1038,
+                0.1012,
+                0.0979,
+                0.0958,
+                0.0925,
+                0.0900,
+                0.0859,
+                0.0807,
+                0.0772,
+                0.0750,
+                0.0706,
+                0.0647,
+                0.0588,
+            ]
+        )
+        self.AggPerm = np.array(
+            [
+                0.0599,
+                0.0474,
+                0.0412,
+                0.0383,
+                0.0340,
+                0.0324,
+                0.0312,
+                0.0289,
+                0.0264,
+                0.0272,
+                0.0238,
+                0.0247,
+                0.0224,
+                0.0235,
+                0.0215,
+                0.0238,
+                0.0222,
+                0.0230,
+                0.0233,
+                0.0233,
+                0.0229,
+                0.0250,
+                0.0266,
+                0.0271,
+                0.0306,
+                0.0326,
+                0.0351,
+                0.0389,
+            ]
+        )
+
+    def test_paper_results(self):
+        # Test own function against the profiles from Figure 6 in the paper
+
+        # 1940 cohort
+        stds1940 = sabelhaus_song_var_profile(
+            age_min=27, age_max=54, cohort=1940, smooth=False
+        )
+
+        self.assertTrue(
+            np.allclose(
+                self.Fig6Coh1940Tran, np.array(stds1940["TranShkStd"]), atol=1e-03
+            )
+        )
+
+        self.assertTrue(
+            np.allclose(
+                self.Fig6Coh1940Perm, np.array(stds1940["PermShkStd"]), atol=1e-03
+            )
+        )
+
+        # 1965 cohort
+        stds1965 = sabelhaus_song_var_profile(
+            age_min=27, age_max=54, cohort=1965, smooth=False
+        )
+
+        self.assertTrue(
+            np.allclose(
+                self.Fig6Coh1965Tran, np.array(stds1965["TranShkStd"]), atol=1e-03
+            )
+        )
+
+        self.assertTrue(
+            np.allclose(
+                self.Fig6Coh1965Perm, np.array(stds1965["PermShkStd"]), atol=1e-03
+            )
+        )
+
+    def test_aggregate_results(self):
+        # Tests own function against the aggregate profiles provided by Sabelhaus
+        stds_agg = sabelhaus_song_var_profile(
+            age_min=27, age_max=54, cohort=None, smooth=False
+        )
+
+        self.assertTrue(
+            np.allclose(self.AggTran, np.array(stds_agg["TranShkStd"]), atol=1e-03)
+        )
+
+        self.assertTrue(
+            np.allclose(self.AggPerm, np.array(stds_agg["PermShkStd"]), atol=1e-03)
+        )
+
+    def test_smoothing(self):
+        # Tests the smoothing approximations that we make. The test is wether
+        # the smoothed profiles are "close" to the exact estimates from Sabelhaus
+
+        # Ensure smoothed versions are within +-10% of the exact estimates
+        rtol = 1e-1
+
+        smooth1940 = sabelhaus_song_var_profile(
+            age_min=27, age_max=54, cohort=1940, smooth=True
+        )
+        smooth1965 = sabelhaus_song_var_profile(
+            age_min=27, age_max=54, cohort=1965, smooth=True
+        )
+        smoothAgg = sabelhaus_song_var_profile(
+            age_min=27, age_max=54, cohort=None, smooth=True
+        )
+
+        # 1940
+        self.assertTrue(
+            np.allclose(
+                np.array(smooth1940["TranShkStd"]), self.Fig6Coh1940Tran, rtol= rtol
+            )
+        )
+
+        self.assertTrue(
+            np.allclose(
+                np.array(smooth1940["PermShkStd"]), self.Fig6Coh1940Perm, atol= rtol
+            )
+        )
+        
+        # 1965
+        self.assertTrue(
+            np.allclose(
+                np.array(smooth1965["TranShkStd"]), self.Fig6Coh1965Tran, rtol= rtol
+            )
+        )
+
+        self.assertTrue(
+            np.allclose(
+                np.array(smooth1965["PermShkStd"]), self.Fig6Coh1965Perm, atol= rtol
+            )
+        )
+        
+        # Aggregate
+        self.assertTrue(
+            np.allclose(
+                np.array(smoothAgg["TranShkStd"]), self.AggTran, rtol= rtol
+            )
+        )
+
+        self.assertTrue(
+            np.allclose(
+                np.array(smoothAgg["PermShkStd"]), self.AggPerm, atol= rtol
+            )
+        )