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

Issue:
Add constructor architecture
Beginning work on the `construct` method of `AgentType`. As of this moment, a basic prototype works for income distributions. Still a lot more to develop and test.

- [x] Tests for new functionality/models or Tests to reproduce the bug-fix in code.
- [x] Updated documentation of features that add new functionality.
- [x] Update CHANGELOG.md with major/minor changes.



Discussion / hints:
As noted in the issue where I proposed this, there are two significant design points to work out:

1) How to handle situations where it's natural to construct *multiple* objects in the same function. The most obvious case is the first example I implemented: `IncShkDstn`, `PermShkDstn`, and `TranShkDstn`. I put in an ugly workaround, but something else should be done before this is merged.

2) Whether and how we want constructors to indicate whether they've made a time-varying or time-invariant object. I'm beginning to come around to the idea that these concepts (to the extent they're in HARK now and will be until several more overhaul steps are completed) should *not* be intended to be changed at the instance level. Rather, they're a fixed fact of the `AgentType`. There's a very specific reason they were defined at the instance level, but that reason no longer exists.
Can someone tell me how to run the pre-commit routine locally? I can run black, but I don't know how to make lint work.
`$ pre-commit run --all-files`
more info in the developer guide here: https://docs.econ-ark.org/guides/contributing.html#pre-commit-hooks
```
pip install pre-commit # install 
pre-commit install # start-up
git add * # stage your changes
pre-commit # run the hooks
```

> `$ pre-commit run --all-files`

you probably don't want to run all files, only your changes 
A bunch of tests needed syntax changes, but everything passes now. The other things that I'm moving to constructors will be much easier from here on out, and should be completed on Monday.

For this PR, I'm going to set the goal at getting everything in constructor functions that should be (plus some alternates that can be swapped in). After this is merged, in a separate branch/PR, I'm going to tackle undoing/simplifying the custom init methods for each AgentType subclass, and get them all just using AgentType.__init__ ... That also means that the default dictionary should be a class attribute.
We should be doing super().init everywhere 
Even where `super().__init__` is used, there can be problems when `update()` is called by a parent class as part of init before the child class is "ready" for it, leading to unexpected results. In the not-too-distant future, I hope all our AgentType subclasses inherit directly from AgentType, and don't have any special init stuff.
If tests pass on this commit, the next commit will change *exactly one line* of code, and revise tests to fit. As is, the "construct pLvlGrid by simulation" method draws on AgentCount, which is fine *as long as AgentCount is actually large*. The number of eventual simulated agents shouldn't be used in the *solution* of the problem, even extremely indirectly. So I'm just going to have it use a *large number* of agents no matter what. This is also why the hardcoded seeds don't matter in this function/method.
The remaining tasks on this PR are pretty small, almost optional:

- make some example MrkvArray constructors
- make terminal solution constructors for types that are missing it
- consider moving one method in ConsAggShockType to a constructor (but maybe not)
- make a method for AgentType that reports what the user *needs* to make a valid type
- a few "alternate example" constructor functions, to demonstrate swappability

I'm going to try to finish this by Friday, as it has fallen behind schedule.
Also: Is this ready for use?  Ryan Edwards is working on making a new Aiyagari REMARK using HARK and I want him to use the latest and greatest tools ...


> This also allows for transitory shocks, I assume?
> 
> Someday I want a tool that can handle the case that I think captures almost all of what can be said about annual income processes: It's hard to do better than a process with a fully permanent shock and a transitory shock with an MA(1) term (as time aggregation would imply). I take it that your setup here allows for purely transitory terms but probably not for MA(1) transitory terms?

This PR is just a restructuring. I've made *very* small changes to the constructor methods when converting them to functions, but they're almost trivial adjustments. This isn't about economics or new capabilities at all. So all the same income processes are in there that we had before.

It *might* be possible to do an MA(1) transitory process using `MarkovConsumerType`, but I'd have to work out whether there's a better way than having N^2 discrete states, where N is the number of transitory shock values.
Meant to say before posting: If it's not, we can add something later that can handle this income process. Let me think about it.
The simplest way to do it in HARK (read: the way that requires the least new code) is to make a subclass of `MarkovConsumerType` with an adjusted simulation procedure. When solving the model, the Markov state can be "today's transitory shock", which simply gives the forward-looking agent a little bit of information about next period's *transitory income* (weighted average of today's shock and tomorrow's shock).

When *simulating* the model, you just need a mechanism that enforces that (the index of) today's transitory shock is the same as (the index of) today's discrete state.
Still behind schedule, as there were more terminal solvers than I remembered. Still to do the following small tasks:

- make some example MrkvArray constructors
- a few "alternate example" constructor functions, to demonstrate swappability
Once the current tests finish, I'll be pushing a commit that only has documentation / changelog updates. I've added a small number of "alternate" constructors that can be substituted in, but there's an arbitrarily large number of other ones that *could* be written as well.

Following yesterday's discussion with Seb and Alan, I added an additional features to construct: the *option* to `force` it to continue running even if it encounters exceptions. In the default (`force=False`), any exception *other than* simply not finding an argument will raise an exception as normal; reaching the end of the method without having completed its work *is* an exception. Under `force=True`, it records each exception to `_constructor_errors` and then ignores it, terminating only when it runs out of work to do, as usual. This option allows the user to (e.g.) run `construct(force=True)` and then simply examine `ThisType._constructor_errors` (and do `ThisType.describe_constructors()`) to get a quick read of what happened when they tried to build all of their constructed inputs.

This PR is ready for others to look at. No external behavior has changed. The only changes in this PR are how complex objects get created.
This all looks very good to me.

The one thing which I think could be improved is the confusing way in which constructed objects are stored, both (!) as an entry in `parameters` and as a new attribute on the model object. Neither of these seem quite right to me.

What about a `self.constructs` dictionary?

Fix patch:
diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md
index c065c4d89..ebf09a221 100644
--- a/Documentation/CHANGELOG.md
+++ b/Documentation/CHANGELOG.md
@@ -12,6 +12,8 @@ For more information on HARK, see [our Github organization](https://github.com/e
 
 Release Date: TBA
 
+Note: Due to major changes on this release, you may need to adjust how AgentTypes are instantiated in your projects using HARK. If you are manually constructing "complicated" objects like MrkvArray, they should be assigned to your instances *after* initialization, not passed as part of the parameter dictionary. See also the new constructor methodology for how to pass parameters for such constructed inputs.
+
 ### Major Changes
 
 - Replace object-oriented solvers with single function versions. [#1394](https://github.com/econ-ark/HARK/pull/1394)
@@ -29,9 +31,10 @@ Release Date: TBA
 - Moves computation of "stable points" from inside of ConsIndShock solver to a post-solution method. [#1349](https://github.com/econ-ark/HARK/pull/1349)
 - Corrects calculation of "human wealth" under risky returns, providing correct limiting linear consumption function. [#1403](https://github.com/econ-ark/HARK/pull/1403)
 - Removed 'parameters' from new block definitions; these are now 'calibrations' provided separately.
-- Create functions for well-known and repeated calculations in single-function solvers. [#1395](https://github.com/econ-ark/HARK/pull/1395)
 - Create functions for well-known and repeated calculations in single-function solvers. [1395](https://github.com/econ-ark/HARK/pull/1395)
 - Re-work WealthPortfolioSolver to use approximate EGM method [#1404](https://github.com/econ-ark/HARK/pull/1404)
+- Default parameter dictionaries for AgentType subclasses have been "flattened": all parameters appear in one place for each model, rather than inheriting from parent models' dictionaries. The only exception is submodels *within* a file when only 1 or 2 parameters are added or changed. (#1425)[https://github.com/econ-ark/HARK/pull/1425]
+
 
 ### 0.14.1
 
diff --git a/HARK/Calibration/Income/IncomeProcesses.py b/HARK/Calibration/Income/IncomeProcesses.py
index dfe067e4b..82fd021e1 100644
--- a/HARK/Calibration/Income/IncomeProcesses.py
+++ b/HARK/Calibration/Income/IncomeProcesses.py
@@ -540,7 +540,68 @@ def make_AR1_style_pLvlNextFunc(T_cycle, pLvlInitMean, PermGroFac, PrstIncCorr):
     return pLvlNextFunc
 
 
-def construct_pLvlGrid_by_simulation(
+###############################################################################
+
+
+def make_basic_pLvlPctiles(
+    pLvlPctiles_count,
+    pLvlPctiles_bound=[0.001, 0.999],
+    pLvlPctiles_tail_count=0,
+    pLvlPctiles_tail_order=np.e,
+):
+    """
+    Make a relatively basic specification for pLvlPctiles by choosing the number
+    of uniformly spaced nodes in the "body", the percentile boundaries for the
+    body, the number of nodes in each tail, and the order/factor by which the
+    tail percentiles approach 0 and 1 respectively.
+
+    Parameters
+    ----------
+    pLvlPctile_count : int
+        Number of nodes in the "body" of the percentile set.
+    pLvlPctile_bound : [float,float], optional
+        Percentile bounds for the "body" of the set. The default is [0.0, 1.0].
+    pLvlPctile_tail_count : int, optional
+       Number of nodes in each extant tail of the set. The default is 0.
+    pLvlPctile_tail_order : float, optional
+        Factor by which tail percentiles shrink toward 0 and 1. The default is np.e.
+
+    Returns
+    -------
+    pLvlPctiles : np.array
+        Array of percentiles of pLvl, usually used to construct pLvlGrid using
+        the function below.
+    """
+    bound = pLvlPctiles_bound
+    fac = 1.0 / pLvlPctiles_tail_order
+    body = np.linspace(bound[0], bound[1], num=pLvlPctiles_count)
+
+    if bound[0] > 0.0:
+        lower = []
+        val = bound[0]
+        for i in range(pLvlPctiles_tail_count):
+            val *= fac
+            lower.append(val)
+        lower.reverse()
+        lower = np.array(lower)
+    else:
+        lower = np.array([])
+
+    if bound[1] < 1.0:
+        upper = []
+        val = 1.0 - bound[1]
+        for i in range(pLvlPctiles_tail_count):
+            val *= fac
+            upper.append(val)
+        upper = 1.0 - np.array(upper)
+    else:
+        upper = np.array([])
+
+    pLvlPctiles = np.concatenate((lower, body, upper))
+    return pLvlPctiles
+
+
+def make_pLvlGrid_by_simulation(
     cycles,
     T_cycle,
     PermShkDstn,
diff --git a/HARK/ConsumptionSaving/ConsAggShockModel.py b/HARK/ConsumptionSaving/ConsAggShockModel.py
index 3b1bb7183..13d799050 100644
--- a/HARK/ConsumptionSaving/ConsAggShockModel.py
+++ b/HARK/ConsumptionSaving/ConsAggShockModel.py
@@ -11,6 +11,11 @@
 import scipy.stats as stats
 
 from HARK import AgentType, Market
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+)
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
     IndShockConsumerType,
@@ -45,7 +50,7 @@
     CRRAutilityP_inv,
     CRRAutilityPP,
 )
-from HARK.utilities import make_grid_exp_mult
+from HARK.utilities import make_grid_exp_mult, make_assets_grid
 
 __all__ = [
     "AggShockConsumerType",
@@ -101,20 +106,73 @@ def make_aggshock_solution_terminal(CRRA):
 aggshock_constructor_dict = indshk_constructor_dict.copy()
 aggshock_constructor_dict["solution_terminal"] = make_aggshock_solution_terminal
 
-# Interest factor is endogenous in agg shocks model
-del init_agg_shocks["Rfree"]
-del init_agg_shocks["CubicBool"]  # Not supported yet for agg shocks model
-del init_agg_shocks["vFuncBool"]  # Not supported yet for agg shocks model
-init_agg_shocks["PermGroFac"] = [1.0]
-# Grid of capital-to-labor-ratios (factors)
+# Make a dictionary of constructors for the aggregate income shocks model
+aggshock_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "solution_terminal": make_aggshock_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 24,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Choose a grid of capital-to-labor-ratios (factors relative to SS)
+# TODO: Make a constructor for this
 MgridBase = np.array(
     [0.1, 0.3, 0.6, 0.8, 0.9, 0.98, 1.0, 1.02, 1.1, 1.2, 1.6, 2.0, 3.0]
 )
-init_agg_shocks["MgridBase"] = MgridBase
-init_agg_shocks["aXtraCount"] = 24
-init_agg_shocks["aNrmInitStd"] = 0.0
-init_agg_shocks["LivPrb"] = [0.98]
-init_agg_shocks["constructors"] = aggshock_constructor_dict
+
+# Make a dictionary to specify an aggregate income shocks consumer type
+init_agg_shocks = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": aggshock_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.00],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "MgridBase": MgridBase,
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 0.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_agg_shocks.update(default_IncShkDstn_params)
+init_agg_shocks.update(default_aXtraGrid_params)
 
 
 class AggShockConsumerType(IndShockConsumerType):
@@ -435,17 +493,11 @@ def make_euler_error_func(self, mMax=100, approx_inc_dstn=True):
         raise NotImplementedError()
 
 
+###############################################################################
+
+
 # This example makes a high risk, low growth state and a low risk, high growth state
 MrkvArray = np.array([[0.90, 0.10], [0.04, 0.96]])
-PermShkAggStd = [
-    0.012,
-    0.006,
-]  # Standard deviation of log aggregate permanent shocks by state
-TranShkAggStd = [
-    0.006,
-    0.003,
-]  # Standard deviation of log aggregate transitory shocks by state
-PermGroFacAgg = [0.98, 1.02]  # Aggregate permanent income growth factor
 
 # Make a dictionary to specify a Markov aggregate shocks consumer
 init_agg_mrkv_shocks = init_agg_shocks.copy()
@@ -621,6 +673,9 @@ def getMrkvNow(self):  # This function exists to be overwritten in StickyE model
         return self.shocks["Mrkv"] * np.ones(self.AgentCount, dtype=int)
 
 
+##############################################################################
+
+# Make a dictionary for Krusell-Smith agents
 init_KS_agents = {
     "T_cycle": 1,
     "DiscFac": 0.99,
diff --git a/HARK/ConsumptionSaving/ConsBequestModel.py b/HARK/ConsumptionSaving/ConsBequestModel.py
index d68d0e7d0..6919535e0 100644
--- a/HARK/ConsumptionSaving/ConsBequestModel.py
+++ b/HARK/ConsumptionSaving/ConsBequestModel.py
@@ -15,18 +15,27 @@
 import numpy as np
 
 from HARK import NullFunc
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+)
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
     IndShockConsumerType,
-    init_idiosyncratic_shocks,
     make_basic_CRRA_solution_terminal,
 )
+from HARK.Calibration.Assets.AssetProcesses import (
+    make_lognormal_RiskyDstn,
+    combine_IncShkDstn_and_RiskyDstn,
+    calc_ShareLimit_for_CRRA,
+)
 from HARK.ConsumptionSaving.ConsPortfolioModel import (
     PortfolioConsumerType,
     PortfolioSolution,
-    init_portfolio,
     make_portfolio_solution_terminal,
 )
+from HARK.ConsumptionSaving.ConsRiskyAssetModel import make_simple_ShareGrid
 from HARK.distribution import expected
 from HARK.interpolation import (
     BilinearInterp,
@@ -41,6 +50,7 @@
     ValueFuncCRRA,
 )
 from HARK.rewards import UtilityFuncCRRA, UtilityFuncStoneGeary
+from HARK.utilities import make_assets_grid
 
 
 def make_bequest_solution_terminal(
@@ -177,6 +187,90 @@ def make_warmglow_portfolio_solution_terminal(
 ###############################################################################
 
 
+# Make a dictionary of constructors for the warm glow bequest model
+warmglow_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "solution_terminal": make_bequest_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Make a dictionary to specify awarm glow bequest consumer type
+init_warm_glow = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": warmglow_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion on consumption
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "BeqCRRA": 2.0,  # Coefficient of relative risk aversion for bequest motive
+    "BeqFac": 40.0,  # Scaling factor for bequest motive
+    "BeqShift": 0.0,  # Stone-Geary shifter term for bequest motive
+    "BeqCRRATerm": 2.0,  # Coefficient of relative risk aversion for bequest motive, terminal period only
+    "BeqFacTerm": 40.0,  # Scaling factor for bequest motive, terminal period only
+    "BeqShiftTerm": 0.0,  # Stone-Geary shifter term for bequest motive, terminal period only
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_warm_glow.update(default_IncShkDstn_params)
+init_warm_glow.update(default_aXtraGrid_params)
+
+# Make a dictionary with bequest motives turned off
+init_accidental_bequest = init_warm_glow.copy()
+init_accidental_bequest["BeqFac"] = 0.0
+init_accidental_bequest["BeqShift"] = 0.0
+init_accidental_bequest["BeqFacTerm"] = 0.0
+init_accidental_bequest["BeqShiftTerm"] = 0.0
+
+# Make a dictionary that has *only* a terminal period bequest
+init_warm_glow_terminal_only = init_warm_glow.copy()
+init_warm_glow_terminal_only["BeqFac"] = 0.0
+init_warm_glow_terminal_only["BeqShift"] = 0.0
+
+
 class BequestWarmGlowConsumerType(IndShockConsumerType):
     time_inv_ = IndShockConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"]
 
@@ -188,18 +282,6 @@ def __init__(self, **kwds):
         self.solve_one_period = solve_one_period_ConsWarmBequest
 
 
-class BequestWarmGlowPortfolioType(PortfolioConsumerType):
-    time_inv_ = PortfolioConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"]
-
-    def __init__(self, **kwds):
-        params = init_portfolio_bequest.copy()
-        params.update(kwds)
-
-        self.IndepDstnBool = True
-        super().__init__(**params)
-        self.solve_one_period = solve_one_period_ConsPortfolioWarmGlow
-
-
 ###############################################################################
 
 
@@ -968,29 +1050,104 @@ def calc_EndOfPrd_v(S, a, z):
     return solution_now
 
 
-init_accidental_bequest = init_idiosyncratic_shocks.copy()
-init_accidental_bequest["BeqCRRA"] = init_idiosyncratic_shocks["CRRA"]
-init_accidental_bequest["BeqFac"] = 0.0
-init_accidental_bequest["BeqShift"] = 0.0
-init_accidental_bequest["BeqCRRATerm"] = init_idiosyncratic_shocks["CRRA"]
-init_accidental_bequest["BeqFacTerm"] = 0.0
-init_accidental_bequest["BeqShiftTerm"] = 0.0
-temp_dict = init_accidental_bequest["constructors"].copy()
-temp_dict["solution_terminal"] = make_bequest_solution_terminal
-init_accidental_bequest["constructors"] = temp_dict
-
-init_warm_glow_terminal_only = init_accidental_bequest.copy()
-init_warm_glow_terminal_only["BeqCRRATerm"] = init_idiosyncratic_shocks["CRRA"]
-init_warm_glow_terminal_only["BeqFacTerm"] = 40.0  # kid lives 40yr after bequest
-init_warm_glow_terminal_only["BeqShiftTerm"] = 0.0
-
-init_warm_glow = init_warm_glow_terminal_only.copy()
-init_warm_glow["BeqCRRA"] = init_idiosyncratic_shocks["CRRA"]
-init_warm_glow["BeqFac"] = 40.0
-init_warm_glow["BeqShift"] = 0.0
-
-init_portfolio_bequest = init_accidental_bequest.copy()
-init_portfolio_bequest.update(init_portfolio)
-temp_dict_two = init_portfolio["constructors"].copy()
-temp_dict_two["solution_terminal"] = make_warmglow_portfolio_solution_terminal
-init_portfolio_bequest["constructors"] = temp_dict_two
+###############################################################################
+
+
+# Make a dictionary of constructors for the portfolio choice consumer type
+portfolio_bequest_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "RiskyDstn": make_lognormal_RiskyDstn,
+    "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
+    "ShareLimit": calc_ShareLimit_for_CRRA,
+    "ShareGrid": make_simple_ShareGrid,
+    "solution_terminal": make_warmglow_portfolio_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 100,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 1,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 200,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make RiskyDstn with make_lognormal_RiskyDstn (and uniform ShareGrid)
+default_RiskyDstn_and_ShareGrid_params = {
+    "RiskyAvg": 1.08,  # Mean return factor of risky asset
+    "RiskyStd": 0.20,  # Stdev of log returns on risky asset
+    "RiskyCount": 5,  # Number of integration nodes to use in approximation of risky returns
+    "ShareCount": 25,  # Number of discrete points in the risky share approximation
+}
+
+# Make a dictionary to specify a risky asset consumer type
+init_portfolio_bequest = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": portfolio_bequest_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 5.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Return factor on risk free asset
+    "DiscFac": 0.90,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "BeqCRRA": 2.0,  # Coefficient of relative risk aversion for bequest motive
+    "BeqFac": 40.0,  # Scaling factor for bequest motive
+    "BeqShift": 0.0,  # Stone-Geary shifter term for bequest motive
+    "BeqCRRATerm": 2.0,  # Coefficient of relative risk aversion for bequest motive, terminal period only
+    "BeqFacTerm": 40.0,  # Scaling factor for bequest motive, terminal period only
+    "BeqShiftTerm": 0.0,  # Stone-Geary shifter term for bequest motive, terminal period only
+    "DiscreteShareBool": False,  # Whether risky asset share is restricted to discrete values
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    "AdjustPrb": 1.0,  # Probability that the agent can update their risky portfolio share each period
+    "sim_common_Rrisky": True,  # Whether risky returns have a shared/common value across agents
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_portfolio_bequest.update(default_IncShkDstn_params)
+init_portfolio_bequest.update(default_aXtraGrid_params)
+init_portfolio_bequest.update(default_RiskyDstn_and_ShareGrid_params)
+
+
+class BequestWarmGlowPortfolioType(PortfolioConsumerType):
+    time_inv_ = PortfolioConsumerType.time_inv_ + ["BeqCRRA", "BeqShift", "BeqFac"]
+
+    def __init__(self, **kwds):
+        params = init_portfolio_bequest.copy()
+        params.update(kwds)
+
+        self.IndepDstnBool = True
+        super().__init__(**params)
+        self.solve_one_period = solve_one_period_ConsPortfolioWarmGlow
diff --git a/HARK/ConsumptionSaving/ConsGenIncProcessModel.py b/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
index 67042e6f8..47c361036 100644
--- a/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
+++ b/HARK/ConsumptionSaving/ConsGenIncProcessModel.py
@@ -9,17 +9,19 @@
 
 from HARK import AgentType, NullFunc
 from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
     pLvlFuncAR1,
     make_trivial_pLvlNextFunc,
     make_explicit_perminc_pLvlNextFunc,
     make_AR1_style_pLvlNextFunc,
-    construct_pLvlGrid_by_simulation,
+    make_pLvlGrid_by_simulation,
+    make_basic_pLvlPctiles,
 )
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
     IndShockConsumerType,
-    init_idiosyncratic_shocks,
-    indshk_constructor_dict,
 )
 from HARK.distribution import Lognormal, expected
 from HARK.interpolation import (
@@ -46,6 +48,7 @@
     CRRAutilityPP,
     UtilityFuncCRRA,
 )
+from HARK.utilities import make_assets_grid
 
 __all__ = [
     "pLvlFuncAR1",
@@ -533,10 +536,6 @@ def calc_vPP_next(S, a, p):
 
 ###############################################################################
 
-# -----------------------------------------------------------------------------
-# ----- Define additional parameters for the persistent shocks model ----------
-# -----------------------------------------------------------------------------
-
 pLvlPctiles = np.concatenate(
     (
         [0.001, 0.005, 0.01, 0.03],
@@ -545,20 +544,90 @@ def calc_vPP_next(S, a, p):
     )
 )
 
-geninc_constructor_dict = indshk_constructor_dict.copy()
-geninc_constructor_dict["pLvlNextFunc"] = make_trivial_pLvlNextFunc
-geninc_constructor_dict["solution_terminal"] = make_2D_CRRA_solution_terminal
-geninc_constructor_dict["pLvlGrid"] = construct_pLvlGrid_by_simulation
-
-# Make a dictionary for the "explicit permanent income" idiosyncratic shocks model
-init_general_inc = init_idiosyncratic_shocks.copy()
-init_general_inc["pLvlPctiles"] = pLvlPctiles
-init_general_inc["pLvlInitStd"] = 0.4  # This *must* be nonzero
-# long run permanent income growth doesn't work yet
-init_general_inc["PermGroFac"] = [1.0]
-init_general_inc["aXtraMax"] = 30
-init_general_inc["aXtraExtra"] = np.array([0.005, 0.01])
-init_general_inc["constructors"] = geninc_constructor_dict
+# Make a constructor dictionary for the general income process consumer type
+geninc_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "pLvlPctiles": make_basic_pLvlPctiles,
+    "pLvlGrid": make_pLvlGrid_by_simulation,
+    "pLvlNextFunc": make_trivial_pLvlNextFunc,
+    "solution_terminal": make_2D_CRRA_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 30,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": [0.005, 0.01],  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make pLvlGrid using make_basic_pLvlPctiles
+default_pLvlPctiles_params = {
+    "pLvlPctiles_count": 19,  # Number of points in the "body" of the grid
+    "pLvlPctiles_bound": [0.05, 0.95],  # Percentile bounds of the "body"
+    "pLvlPctiles_tail_count": 4,  # Number of points in each tail of the grid
+    "pLvlPctiles_tail_order": np.e,  # Scaling factor for points in each tail
+}
+
+# Default parameters to make pLvlGrid using make_trivial_pLvlNextFunc
+default_pLvlGrid_params = {
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.4,  # Standard deviation of log initial permanent income *MUST BE POSITIVE*
+    # "pLvlPctiles": pLvlPctiles,  # Percentiles of permanent income to use for the grid
+    "pLvlExtra": None,  # Additional permanent income points to automatically add to the grid, optional
+}
+
+# Make a dictionary to specify a general income process consumer type
+init_general_inc = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": geninc_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_general_inc.update(default_IncShkDstn_params)
+init_general_inc.update(default_aXtraGrid_params)
+init_general_inc.update(default_pLvlPctiles_params)
+init_general_inc.update(default_pLvlGrid_params)
 
 
 class GenIncProcessConsumerType(IndShockConsumerType):
@@ -661,7 +730,7 @@ def update_pLvlGrid(self):
         -------
         None
         """
-        self.construct("pLvlGrid")
+        self.construct("pLvlPctiles", "pLvlGrid")
         self.add_to_time_vary("pLvlGrid")
 
     def sim_birth(self, which_agents):
@@ -778,10 +847,19 @@ def get_poststates(self):
 
 ###############################################################################
 
+# Make a dictionary for the "explicit permanent income" consumer type; see parent dictionary above.
 explicit_constructor_dict = geninc_constructor_dict.copy()
 explicit_constructor_dict["pLvlNextFunc"] = make_explicit_perminc_pLvlNextFunc
 init_explicit_perm_inc = init_general_inc.copy()
 init_explicit_perm_inc["constructors"] = explicit_constructor_dict
+init_explicit_perm_inc["PermGroFac"] = [1.0]
+
+# NB: Permanent income growth was not in the default dictionary for GenIncProcessConsumerType
+# because its pLvlNextFunc constructor was *trivial*: no permanent income dynamics at all!
+# For the "explicit permanent income" model, this parameter is added back into the dictionary.
+# However, note that if this model is used in an *infinite horizon* setting, it will work
+# best if the product of PermGroFac (across all periods) is 1. If it is far from 1, then the
+# pLvlGrid that is constructed by the default method might not be appropriate.
 
 
 class IndShockExplicitPermIncConsumerType(GenIncProcessConsumerType):
@@ -800,13 +878,12 @@ class IndShockExplicitPermIncConsumerType(GenIncProcessConsumerType):
 
 ###############################################################################
 
-PrstIncCorr = 0.98  # Serial correlation coefficient for permanent income
-
-# Make a dictionary for the "persistent idiosyncratic shocks" model
+# Make a dictionary for the "persistent idiosyncratic shocks" consumer type; see parent dictionary above.
 persistent_constructor_dict = geninc_constructor_dict.copy()
 persistent_constructor_dict["pLvlNextFunc"] = make_AR1_style_pLvlNextFunc
 init_persistent_shocks = init_explicit_perm_inc.copy()
-init_persistent_shocks["PrstIncCorr"] = PrstIncCorr
+init_persistent_shocks["PrstIncCorr"] = 0.98
+# Serial correlation coefficient for permanent income, which is used by make_AR1_style_pLvlNextFunc
 init_persistent_shocks["constructors"] = persistent_constructor_dict
 
 
diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py
index e4f79367e..6fdc15bf9 100644
--- a/HARK/ConsumptionSaving/ConsIndShockModel.py
+++ b/HARK/ConsumptionSaving/ConsIndShockModel.py
@@ -56,7 +56,7 @@
     UtilityFuncCRRA,
 )
 from HARK.utilities import (
-    construct_assets_grid,
+    make_assets_grid,
     gen_tran_matrix_1D,
     gen_tran_matrix_2D,
     jump_to_grid_1D,
@@ -1040,38 +1040,37 @@ def make_basic_CRRA_solution_terminal(CRRA):
 # == Classes for representing types of consumer agents (and things they do) ==
 # ============================================================================
 
+# Make a dictionary of constructors (very simply for perfect foresight model)
 perf_foresight_constructors = {
     "solution_terminal": make_basic_CRRA_solution_terminal,
 }
 
 # Make a dictionary to specify a perfect foresight consumer type
 init_perfect_foresight = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
     "cycles": 1,  # Finite, non-cyclic model
-    "CRRA": 2.0,  # Coefficient of relative risk aversion,
-    "Rfree": 1.03,  # Interest factor on assets
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": perf_foresight_constructors,  # See dictionary above
+    # PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Interest factor on retained assets
     "DiscFac": 0.96,  # Intertemporal discount factor
-    "LivPrb": [0.98],  # Survival probability
+    "LivPrb": [0.98],  # Survival probability after each period
     "PermGroFac": [1.01],  # Permanent income growth factor
     "BoroCnstArt": None,  # Artificial borrowing constraint
-    # Maximum number of grid points to allow in cFunc (should be large)
-    "MaxKinks": 400,
-    # Number of agents of this type (only matters for simulation)
-    "AgentCount": 10000,
-    # Mean of log initial assets (only matters for simulation)
-    "aNrmInitMean": 0.0,
-    # Standard deviation of log initial assets (only for simulation)
-    "aNrmInitStd": 1.0,
-    # Mean of log initial permanent income (only matters for simulation)
-    "pLvlInitMean": 0.0,
-    # Standard deviation of log initial permanent income (only matters for simulation)
-    "pLvlInitStd": 0.0,
-    # Aggregate permanent income growth factor: portion of PermGroFac attributable to aggregate productivity growth (only matters for simulation)
-    "PermGroFacAgg": 1.0,
+    "MaxKinks": 400,  # Maximum number of grid points to allow in cFunc
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
     "T_age": None,  # Age after which simulated agents are automatically killed
-    "T_cycle": 1,  # Number of periods in the cycle for this agent type
-    "PerfMITShk": False,
-    # Do Perfect Foresight MIT Shock: Forces Newborns to follow solution path of the agent he/she replaced when True
-    "constructors": perf_foresight_constructors,
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
 }
 
 
@@ -1820,49 +1819,70 @@ def calc_stable_points(self):
 
 ###############################################################################
 
+# Make a dictionary of constructors for the idiosyncratic income shocks model
 indshk_constructor_dict = {
     "IncShkDstn": construct_lognormal_income_process_unemployment,
     "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
     "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
-    "aXtraGrid": construct_assets_grid,
+    "aXtraGrid": make_assets_grid,
     "solution_terminal": make_basic_CRRA_solution_terminal,
 }
 
-# Make a dictionary to specify an idiosyncratic income shocks consumer
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Make a dictionary to specify an idiosyncratic income shocks consumer type
 init_idiosyncratic_shocks = {
-    **init_perfect_foresight,
-    **{  # assets above grid parameters
-        "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
-        "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
-        # Exponential nesting factor when constructing "assets above minimum" grid
-        "aXtraNestFac": 3,
-        "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
-        "aXtraExtra": None,  # Some other values to add to the grid, not used
-        # Income process variables
-        # Standard deviation of log permanent income shocks
-        "PermShkStd": [0.1],
-        "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
-        # Standard deviation of log transitory income shocks
-        "TranShkStd": [0.1],
-        "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
-        "UnempPrb": 0.05,  # Probability of unemployment while working
-        "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
-        "IncUnemp": 0.3,  # Unemployment benefits replacement rate
-        "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
-        # Artificial borrowing constraint; imposed minimum level of end-of period assets
-        "BoroCnstArt": 0.0,
-        "tax_rate": 0.0,  # Flat income tax rate
-        "T_retire": 0,  # Period of retirement (0 --> no retirement)
-        "vFuncBool": False,  # Whether to calculate the value function during solution
-        # Use cubic spline interpolation when True, linear interpolation when False
-        "CubicBool": False,
-        "neutral_measure": False,
-        # Use permanent income neutral measure (see Harmenberg 2021) during simulations when True.
-        # Whether Newborns have transitory shock. The default is False.
-        "NewbornTransShk": False,
-        "constructors": indshk_constructor_dict,
-    },
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": indshk_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
 }
+init_idiosyncratic_shocks.update(default_IncShkDstn_params)
+init_idiosyncratic_shocks.update(default_aXtraGrid_params)
 
 
 class IndShockConsumerType(PerfForesightConsumerType):
@@ -3197,22 +3217,17 @@ def check_conditions(self, verbose=None):
 
 ###############################################################################
 
+# Specify default parameters that differ in "kinked R" model compared to base IndShockConsumerType
+kinked_R_different_params = {
+    "Rboro": 1.20,  # Interest factor on assets when borrowing, a < 0
+    "Rsave": 1.02,  # Interest factor on assets when saving, a > 0
+    "BoroCnstArt": None,  # Kinked R only matters if borrowing is allowed
+}
 
 # Make a dictionary to specify a "kinked R" idiosyncratic shock consumer
-init_kinked_R = dict(
-    init_idiosyncratic_shocks,
-    **{
-        "Rboro": 1.20,  # Interest factor on assets when borrowing, a < 0
-        "Rsave": 1.02,  # Interest factor on assets when saving, a > 0
-        "BoroCnstArt": None,
-        # kinked R is a bit silly if borrowing not allowed
-        "CubicBool": True,
-        # kinked R is now compatible with linear cFunc and cubic cFunc
-        "aXtraCount": 48,
-        # ...so need lots of extra gridpoints to make up for it
-    },
-)
-del init_kinked_R["Rfree"]  # get rid of constant interest factor
+init_kinked_R = init_idiosyncratic_shocks.copy()  # See base dictionary above
+init_kinked_R.update(kinked_R_different_params)  # Update with some parameters
+del init_kinked_R["Rfree"]  # Get rid of constant interest factor
 
 
 class KinkedRconsumerType(IndShockConsumerType):
diff --git a/HARK/ConsumptionSaving/ConsLaborModel.py b/HARK/ConsumptionSaving/ConsLaborModel.py
index 5291d298b..f3cc33a33 100644
--- a/HARK/ConsumptionSaving/ConsLaborModel.py
+++ b/HARK/ConsumptionSaving/ConsLaborModel.py
@@ -14,12 +14,14 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-
-from HARK.Calibration.Income.IncomeProcesses import get_TranShkGrid_from_TranShkDstn
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+    get_TranShkGrid_from_TranShkDstn,
+)
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     IndShockConsumerType,
-    init_idiosyncratic_shocks,
-    indshk_constructor_dict,
 )
 from HARK.interpolation import (
     BilinearInterp,
@@ -32,6 +34,7 @@
 )
 from HARK.metric import MetricObject
 from HARK.rewards import CRRAutilityP, CRRAutilityP_inv
+from HARK.utilities import make_assets_grid
 
 
 class ConsumerLaborSolution(MetricObject):
@@ -453,6 +456,83 @@ def uPinv(X):
 ###############################################################################
 
 
+# Make a dictionary of constructors for the intensive margin labor model
+labor_int_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "LbrCost": make_log_polynomial_LbrCost,
+    "TranShkGrid": get_TranShkGrid_from_TranShkDstn,
+    "solution_terminal": make_labor_intmarg_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 16,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 15,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.0,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 80.0,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 200,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Default parameter to make LbrCost using make_log_polynomial_LbrCost
+defualt_LbrCost_params = {
+    "LbrCostCoeffs": [
+        -1.0
+    ]  # Polynomial coefficients (for age) on log labor utility cost
+}
+
+# Make a dictionary to specify an intensive margin labor supply choice consumer type
+init_labor_intensive = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": labor_int_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "WageRte": [1.0],  # Wage rate paid on labor income
+    "BoroCnstArt": None,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_labor_intensive.update(default_IncShkDstn_params)
+init_labor_intensive.update(default_aXtraGrid_params)
+init_labor_intensive.update(defualt_LbrCost_params)
+
+
 class LaborIntMargConsumerType(IndShockConsumerType):
     """
     A class representing agents who make a decision each period about how much
@@ -676,24 +756,7 @@ def plot_LbrFunc(self, t, bMin=None, bMax=None, ShkSet=None):
         plt.show()
 
 
-# Make a default dictionary for the intensive margin labor supply model
-labor_int_constructor_dict = indshk_constructor_dict.copy()
-labor_int_constructor_dict["solution_terminal"] = make_labor_intmarg_solution_terminal
-labor_int_constructor_dict["LbrCost"] = make_log_polynomial_LbrCost
-labor_int_constructor_dict["TranShkGrid"] = get_TranShkGrid_from_TranShkDstn
-
-init_labor_intensive = copy(init_idiosyncratic_shocks)
-init_labor_intensive["constructors"] = labor_int_constructor_dict
-init_labor_intensive["LbrCostCoeffs"] = [-1.0]
-init_labor_intensive["WageRte"] = [1.0]
-init_labor_intensive["IncUnemp"] = 0.0
-init_labor_intensive["TranShkCount"] = 15
-# Crank up permanent shock count - Number of points in discrete approximation to transitory income shocks
-init_labor_intensive["PermShkCount"] = 16  # Crank up permanent shock count
-init_labor_intensive["aXtraCount"] = 200
-# May be important to have a larger number of gridpoints (than 48 initially)
-init_labor_intensive["aXtraMax"] = 80.0
-init_labor_intensive["BoroCnstArt"] = None
+###############################################################################
 
 # Make a dictionary for intensive margin labor supply model with finite lifecycle
 init_labor_lifecycle = init_labor_intensive.copy()
diff --git a/HARK/ConsumptionSaving/ConsMarkovModel.py b/HARK/ConsumptionSaving/ConsMarkovModel.py
index 94a4cc526..7ea8b2246 100644
--- a/HARK/ConsumptionSaving/ConsMarkovModel.py
+++ b/HARK/ConsumptionSaving/ConsMarkovModel.py
@@ -8,13 +8,16 @@
 import numpy as np
 
 from HARK import AgentType, NullFunc
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+)
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
     IndShockConsumerType,
     PerfForesightConsumerType,
     make_basic_CRRA_solution_terminal,
-    indshk_constructor_dict,
-    init_idiosyncratic_shocks,
 )
 from HARK.distribution import MarkovProcess, Uniform, expected
 from HARK.interpolation import (
@@ -35,6 +38,7 @@
     CRRAutilityP_invP,
     CRRAutilityPP,
 )
+from HARK.utilities import make_assets_grid
 
 __all__ = ["MarkovConsumerType"]
 
@@ -651,11 +655,78 @@ def calc_vPPnext(S, a, R):
 ####################################################################################################
 ####################################################################################################
 
-markov_constructor_dict = indshk_constructor_dict.copy()
-markov_constructor_dict["solution_terminal"] = make_markov_solution_terminal
-
-init_indshk_markov = init_idiosyncratic_shocks.copy()
-init_indshk_markov["constructors"] = markov_constructor_dict
+# Make a dictionary of constructors for the markov consumption-saving model
+markov_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "MrkvArray": make_simple_binary_markov,
+    "solution_terminal": make_markov_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make MrkvArray using make_simple_binary_markov
+default_MrkvArray_params = {
+    "Mrkv_p11": [0.9],  # Probability of remaining in binary state 1
+    "Mrkv_p22": [0.4],  # Probability of remaining in binary state 2
+}
+
+# Make a dictionary to specify an idiosyncratic income shocks consumer type
+init_indshk_markov = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": markov_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": np.array([1.03, 1.03]),  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [np.array([0.98, 0.98])],  # Survival probability after each period
+    "PermGroFac": [np.array([0.99, 1.03])],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_indshk_markov.update(default_IncShkDstn_params)
+init_indshk_markov.update(default_aXtraGrid_params)
+init_indshk_markov.update(default_MrkvArray_params)
 
 
 class MarkovConsumerType(IndShockConsumerType):
@@ -668,7 +739,7 @@ class MarkovConsumerType(IndShockConsumerType):
 
     time_vary_ = IndShockConsumerType.time_vary_ + ["MrkvArray"]
 
-    # Is "Mrkv" a shock or a state?
+    # Mrkv is both a shock and a state
     shock_vars_ = IndShockConsumerType.shock_vars_ + ["Mrkv"]
     state_vars = IndShockConsumerType.state_vars + ["Mrkv"]
 
@@ -682,6 +753,22 @@ def __init__(self, **kwds):
         if not hasattr(self, "global_markov"):
             self.global_markov = False
 
+    def update(self):
+        """
+        Update the Markov array, the income process, the assets grid, and
+        the terminal solution.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        """
+        self.construct("MrkvArray")
+        super().update()
+
     def check_markov_inputs(self):
         """
         Many parameters used by MarkovConsumerType are arrays.  Make sure those arrays are the
@@ -700,7 +787,7 @@ def check_markov_inputs(self):
         # Check that arrays are the right shape
         if not isinstance(self.Rfree, np.ndarray) or self.Rfree.shape != (StateCount,):
             raise ValueError(
-                "Rfree not the right shape, it should an array of Rfree of all the states."
+                "Rfree not the right shape, it should be an array of Rfree of all the states."
             )
 
         # Check that arrays in lists are the right shape
@@ -756,6 +843,7 @@ def pre_solve(self):
         """
         AgentType.pre_solve(self)
         self.check_markov_inputs()
+        self.update_solution_terminal()
 
     def initialize_sim(self):
         self.shocks["Mrkv"] = np.zeros(self.AgentCount, dtype=int)
diff --git a/HARK/ConsumptionSaving/ConsMedModel.py b/HARK/ConsumptionSaving/ConsMedModel.py
index 7bb8f641d..aeed4efe0 100644
--- a/HARK/ConsumptionSaving/ConsMedModel.py
+++ b/HARK/ConsumptionSaving/ConsMedModel.py
@@ -8,11 +8,17 @@
 from scipy.optimize import brentq
 
 from HARK import AgentType
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+    make_AR1_style_pLvlNextFunc,
+    make_pLvlGrid_by_simulation,
+    make_basic_pLvlPctiles,
+)
 from HARK.ConsumptionSaving.ConsGenIncProcessModel import (
     PersistentShockConsumerType,
     VariableLowerBoundFunc2D,
-    init_persistent_shocks,
-    persistent_constructor_dict,
 )
 from HARK.ConsumptionSaving.ConsIndShockModel import ConsumerSolution
 from HARK.distribution import Lognormal, add_discrete_outcome_constant_mean, expected
@@ -40,7 +46,7 @@
     CRRAutilityPP,
     UtilityFuncCRRA,
 )
-from HARK.utilities import NullFunc, make_grid_exp_mult
+from HARK.utilities import NullFunc, make_grid_exp_mult, make_assets_grid
 
 __all__ = [
     "MedShockPolicyFunc",
@@ -706,32 +712,111 @@ def make_MedShock_solution_terminal(
 
 ###############################################################################
 
-# -----------------------------------------------------------------------------
-# ----- Define additional parameters for the medical shocks model -------------
-# -----------------------------------------------------------------------------
-CRRA = 2.0
-CRRAmed = 1.5 * CRRA  # Coefficient of relative risk aversion for medical care
-MedShkAvg = [0.001]  # Average of medical need shocks
-MedShkStd = [5.0]  # Standard deviation of (log) medical need shocks
-MedShkCount = 5  # Number of medical shock points in "body"
-MedShkCountTail = 15  # Number of medical shock points in "tail" (upper only)
-MedPrice = [1.5]  # Relative price of a unit of medical care
-
-medshock_constructor_dict = persistent_constructor_dict.copy()
-medshock_constructor_dict["MedShkDstn"] = make_lognormal_MedShkDstn
-medshock_constructor_dict["solution_terminal"] = make_MedShock_solution_terminal
-
-# Make a dictionary for the "medical shocks" model
-init_medical_shocks = init_persistent_shocks.copy()
-init_medical_shocks["CRRAmed"] = CRRAmed
-init_medical_shocks["MedShkAvg"] = MedShkAvg
-init_medical_shocks["MedShkStd"] = MedShkStd
-init_medical_shocks["MedShkCount"] = MedShkCount
-init_medical_shocks["MedShkCountTail"] = MedShkCountTail
-init_medical_shocks["MedPrice"] = MedPrice
-init_medical_shocks["aXtraCount"] = 32
-init_medical_shocks["pLvlExtra"] = [0.0001]
-init_medical_shocks["constructors"] = medshock_constructor_dict
+# Make a constructor dictionary for the general income process consumer type
+medshock_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "pLvlPctiles": make_basic_pLvlPctiles,
+    "pLvlGrid": make_pLvlGrid_by_simulation,
+    "pLvlNextFunc": make_AR1_style_pLvlNextFunc,
+    "MedShkDstn": make_lognormal_MedShkDstn,
+    "solution_terminal": make_MedShock_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 30,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 32,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": [0.005, 0.01],  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make pLvlGrid using make_basic_pLvlPctiles
+default_pLvlPctiles_params = {
+    "pLvlPctiles_count": 19,  # Number of points in the "body" of the grid
+    "pLvlPctiles_bound": [0.05, 0.95],  # Percentile bounds of the "body"
+    "pLvlPctiles_tail_count": 4,  # Number of points in each tail of the grid
+    "pLvlPctiles_tail_order": np.e,  # Scaling factor for points in each tail
+}
+
+# Default parameters to make pLvlGrid using make_trivial_pLvlNextFunc
+default_pLvlGrid_params = {
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.4,  # Standard deviation of log initial permanent income *MUST BE POSITIVE*
+    # "pLvlPctiles": pLvlPctiles,  # Percentiles of permanent income to use for the grid
+    "pLvlExtra": [
+        0.0001
+    ],  # Additional permanent income points to automatically add to the grid, optional
+}
+
+# Default parameters to make MedShkDstn using make_lognormal_MedShkDstn
+default_MedShkDstn_params = {
+    "MedShkAvg": [0.001],  # Average of medical need shocks
+    "MedShkStd": [5.0],  # Standard deviation of (log) medical need shocks
+    "MedShkCount": 5,  # Number of medical shock points in "body"
+    "MedShkCountTail": 15,  # Number of medical shock points in "tail" (upper only)
+    "MedPrice": [1.5],  # Relative price of a unit of medical care
+}
+
+# Default parameters to make pLvlNextFunc using make_AR1_style_pLvlNextFunc
+default_pLvlNextFunc_params = {
+    "PermGroFac": [1.0],  # Permanent income growth factor
+    "PrstIncCorr": 0.98,  # Correlation coefficient on (log) persistent income
+}
+
+# Make a dictionary to specify a medical shocks consumer type
+init_medical_shocks = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": medshock_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion on consumption
+    "CRRAmed": 3.0,  # Coefficient of relative risk aversion on medical care
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_medical_shocks.update(default_IncShkDstn_params)
+init_medical_shocks.update(default_aXtraGrid_params)
+init_medical_shocks.update(default_pLvlPctiles_params)
+init_medical_shocks.update(default_pLvlGrid_params)
+init_medical_shocks.update(default_MedShkDstn_params)
+init_medical_shocks.update(default_pLvlNextFunc_params)
 
 
 class MedShockConsumerType(PersistentShockConsumerType):
diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py
index a7b2cb36d..588972343 100644
--- a/HARK/ConsumptionSaving/ConsPortfolioModel.py
+++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py
@@ -12,10 +12,19 @@
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     IndShockConsumerType,
 )
+from HARK.Calibration.Assets.AssetProcesses import (
+    make_lognormal_RiskyDstn,
+    combine_IncShkDstn_and_RiskyDstn,
+    calc_ShareLimit_for_CRRA,
+)
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+)
 from HARK.ConsumptionSaving.ConsRiskyAssetModel import (
     RiskyAssetConsumerType,
-    init_risky_asset,
-    risky_constructor_dict,
+    make_simple_ShareGrid,
 )
 from HARK.distribution import expected
 from HARK.interpolation import (
@@ -30,6 +39,7 @@
 )
 from HARK.metric import MetricObject
 from HARK.rewards import UtilityFuncCRRA
+from HARK.utilities import make_assets_grid
 
 
 # Define a class to represent the single period solution of the portfolio choice problem
@@ -147,112 +157,6 @@ def __init__(
         self.AdjPrb = AdjPrb
 
 
-class PortfolioConsumerType(RiskyAssetConsumerType):
-    """
-    A consumer type with a portfolio choice. This agent type has log-normal return
-    factors. Their problem is defined by a coefficient of relative risk aversion,
-    intertemporal discount factor, risk-free interest factor, and time sequences of
-    permanent income growth rate, survival probability, and permanent and transitory
-    income shock standard deviations (in logs).  The agent may also invest in a risky
-    asset, which has a higher average return than the risk-free asset.  He *might*
-    have age-varying beliefs about the risky-return; if he does, then "true" values
-    of the risky asset's return distribution must also be specified.
-    """
-
-    time_inv_ = deepcopy(RiskyAssetConsumerType.time_inv_)
-    time_inv_ = time_inv_ + ["AdjustPrb", "DiscreteShareBool"]
-
-    def __init__(self, verbose=False, quiet=False, **kwds):
-        params = init_portfolio.copy()
-        params.update(kwds)
-        kwds = params
-        self.PortfolioBool = True
-
-        # Initialize a basic consumer type
-        RiskyAssetConsumerType.__init__(self, verbose=verbose, quiet=quiet, **kwds)
-
-        # Set the solver for the portfolio model, and update various constructed attributes
-        self.solve_one_period = solve_one_period_ConsPortfolio
-
-    def initialize_sim(self):
-        """
-        Initialize the state of simulation attributes.  Simply calls the same method
-        for IndShockConsumerType, then sets the type of AdjustNow to bool.
-
-        Parameters
-        ----------
-        None
-
-        Returns
-        -------
-        None
-        """
-        # these need to be set because "post states",
-        # but are a control variable and shock, respectively
-        self.controls["Share"] = np.zeros(self.AgentCount)
-        RiskyAssetConsumerType.initialize_sim(self)
-
-    def sim_birth(self, which_agents):
-        """
-        Create new agents to replace ones who have recently died; takes draws of
-        initial aNrm and pLvl, as in ConsIndShockModel, then sets Share and Adjust
-        to zero as initial values.
-        Parameters
-        ----------
-        which_agents : np.array
-            Boolean array of size AgentCount indicating which agents should be "born".
-
-        Returns
-        -------
-        None
-        """
-        IndShockConsumerType.sim_birth(self, which_agents)
-
-        self.controls["Share"][which_agents] = 0
-        # here a shock is being used as a 'post state'
-        self.shocks["Adjust"][which_agents] = False
-
-    def get_controls(self):
-        """
-        Calculates consumption cNrmNow and risky portfolio share ShareNow using
-        the policy functions in the attribute solution.  These are stored as attributes.
-
-        Parameters
-        ----------
-        None
-
-        Returns
-        -------
-        None
-        """
-        cNrmNow = np.zeros(self.AgentCount) + np.nan
-        ShareNow = np.zeros(self.AgentCount) + np.nan
-
-        # Loop over each period of the cycle, getting controls separately depending on "age"
-        for t in range(self.T_cycle):
-            these = t == self.t_cycle
-
-            # Get controls for agents who *can* adjust their portfolio share
-            those = np.logical_and(these, self.shocks["Adjust"])
-            cNrmNow[those] = self.solution[t].cFuncAdj(self.state_now["mNrm"][those])
-            ShareNow[those] = self.solution[t].ShareFuncAdj(
-                self.state_now["mNrm"][those]
-            )
-
-            # Get Controls for agents who *can't* adjust their portfolio share
-            those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
-            cNrmNow[those] = self.solution[t].cFuncFxd(
-                self.state_now["mNrm"][those], ShareNow[those]
-            )
-            ShareNow[those] = self.solution[t].ShareFuncFxd(
-                self.state_now["mNrm"][those], ShareNow[those]
-            )
-
-        # Store controls as attributes of self
-        self.controls["cNrm"] = cNrmNow
-        self.controls["Share"] = ShareNow
-
-
 ###############################################################################
 
 
@@ -1090,24 +994,194 @@ def solve_one_period_ConsPortfolio(
     return solution_now
 
 
-# Make a dictionary to specify a portfolio choice consumer type
-init_portfolio = init_risky_asset.copy()
-portfolio_constructor_dict = risky_constructor_dict.copy()
-portfolio_constructor_dict["solution_terminal"] = make_portfolio_solution_terminal
-init_portfolio["constructors"] = portfolio_constructor_dict
-init_portfolio["AdjustPrb"] = 1.0
-# Flag for whether to optimize risky share on a discrete grid only
-init_portfolio["DiscreteShareBool"] = False
-
-# Adjust some of the existing parameters in the dictionary
-init_portfolio["RiskyAvg"] = 1.08  # Average return of the risky asset
-init_portfolio["RiskyStd"] = 0.20  # Standard deviation of (log) risky returns
-init_portfolio["aXtraMax"] = 100  # Make the grid of assets go much higher...
-init_portfolio["aXtraCount"] = 200  # ...and include many more gridpoints...
-# ...which aren't so clustered at the bottom
-init_portfolio["aXtraNestFac"] = 1
-# Artificial borrowing constraint must be turned on
-init_portfolio["BoroCnstArt"] = 0.0
-# Results are more interesting with higher risk aversion
-init_portfolio["CRRA"] = 5.0
-init_portfolio["DiscFac"] = 0.90  # And also lower patience
+###############################################################################
+
+# Make a dictionary of constructors for the portfolio choice consumer type
+portfolio_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "RiskyDstn": make_lognormal_RiskyDstn,
+    "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
+    "ShareLimit": calc_ShareLimit_for_CRRA,
+    "ShareGrid": make_simple_ShareGrid,
+    "solution_terminal": make_portfolio_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 100,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 1,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 200,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make RiskyDstn with make_lognormal_RiskyDstn (and uniform ShareGrid)
+default_RiskyDstn_and_ShareGrid_params = {
+    "RiskyAvg": 1.08,  # Mean return factor of risky asset
+    "RiskyStd": 0.20,  # Stdev of log returns on risky asset
+    "RiskyCount": 5,  # Number of integration nodes to use in approximation of risky returns
+    "ShareCount": 25,  # Number of discrete points in the risky share approximation
+}
+
+# Make a dictionary to specify a risky asset consumer type
+init_portfolio = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": portfolio_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 5.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Return factor on risk free asset
+    "DiscFac": 0.90,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "DiscreteShareBool": False,  # Whether risky asset share is restricted to discrete values
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    "AdjustPrb": 1.0,  # Probability that the agent can update their risky portfolio share each period
+    "sim_common_Rrisky": True,  # Whether risky returns have a shared/common value across agents
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_portfolio.update(default_IncShkDstn_params)
+init_portfolio.update(default_aXtraGrid_params)
+init_portfolio.update(default_RiskyDstn_and_ShareGrid_params)
+
+
+class PortfolioConsumerType(RiskyAssetConsumerType):
+    """
+    A consumer type with a portfolio choice. This agent type has log-normal return
+    factors. Their problem is defined by a coefficient of relative risk aversion,
+    intertemporal discount factor, risk-free interest factor, and time sequences of
+    permanent income growth rate, survival probability, and permanent and transitory
+    income shock standard deviations (in logs).  The agent may also invest in a risky
+    asset, which has a higher average return than the risk-free asset.  He *might*
+    have age-varying beliefs about the risky-return; if he does, then "true" values
+    of the risky asset's return distribution must also be specified.
+    """
+
+    time_inv_ = deepcopy(RiskyAssetConsumerType.time_inv_)
+    time_inv_ = time_inv_ + ["AdjustPrb", "DiscreteShareBool"]
+
+    def __init__(self, verbose=False, quiet=False, **kwds):
+        params = init_portfolio.copy()
+        params.update(kwds)
+        kwds = params
+        self.PortfolioBool = True
+
+        # Initialize a basic consumer type
+        RiskyAssetConsumerType.__init__(self, verbose=verbose, quiet=quiet, **kwds)
+
+        # Set the solver for the portfolio model, and update various constructed attributes
+        self.solve_one_period = solve_one_period_ConsPortfolio
+
+    def initialize_sim(self):
+        """
+        Initialize the state of simulation attributes.  Simply calls the same method
+        for IndShockConsumerType, then sets the type of AdjustNow to bool.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        """
+        # these need to be set because "post states",
+        # but are a control variable and shock, respectively
+        self.controls["Share"] = np.zeros(self.AgentCount)
+        RiskyAssetConsumerType.initialize_sim(self)
+
+    def sim_birth(self, which_agents):
+        """
+        Create new agents to replace ones who have recently died; takes draws of
+        initial aNrm and pLvl, as in ConsIndShockModel, then sets Share and Adjust
+        to zero as initial values.
+        Parameters
+        ----------
+        which_agents : np.array
+            Boolean array of size AgentCount indicating which agents should be "born".
+
+        Returns
+        -------
+        None
+        """
+        IndShockConsumerType.sim_birth(self, which_agents)
+
+        self.controls["Share"][which_agents] = 0
+        # here a shock is being used as a 'post state'
+        self.shocks["Adjust"][which_agents] = False
+
+    def get_controls(self):
+        """
+        Calculates consumption cNrmNow and risky portfolio share ShareNow using
+        the policy functions in the attribute solution.  These are stored as attributes.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        """
+        cNrmNow = np.zeros(self.AgentCount) + np.nan
+        ShareNow = np.zeros(self.AgentCount) + np.nan
+
+        # Loop over each period of the cycle, getting controls separately depending on "age"
+        for t in range(self.T_cycle):
+            these = t == self.t_cycle
+
+            # Get controls for agents who *can* adjust their portfolio share
+            those = np.logical_and(these, self.shocks["Adjust"])
+            cNrmNow[those] = self.solution[t].cFuncAdj(self.state_now["mNrm"][those])
+            ShareNow[those] = self.solution[t].ShareFuncAdj(
+                self.state_now["mNrm"][those]
+            )
+
+            # Get Controls for agents who *can't* adjust their portfolio share
+            those = np.logical_and(these, np.logical_not(self.shocks["Adjust"]))
+            cNrmNow[those] = self.solution[t].cFuncFxd(
+                self.state_now["mNrm"][those], ShareNow[those]
+            )
+            ShareNow[those] = self.solution[t].ShareFuncFxd(
+                self.state_now["mNrm"][those], ShareNow[those]
+            )
+
+        # Store controls as attributes of self
+        self.controls["cNrm"] = cNrmNow
+        self.controls["Share"] = ShareNow
+
+
+###############################################################################
diff --git a/HARK/ConsumptionSaving/ConsPrefShockModel.py b/HARK/ConsumptionSaving/ConsPrefShockModel.py
index ec7eb0d51..f99290301 100644
--- a/HARK/ConsumptionSaving/ConsPrefShockModel.py
+++ b/HARK/ConsumptionSaving/ConsPrefShockModel.py
@@ -14,8 +14,13 @@
     ConsumerSolution,
     IndShockConsumerType,
     KinkedRconsumerType,
-    init_idiosyncratic_shocks,
-    init_kinked_R,
+    make_assets_grid,
+    make_basic_CRRA_solution_terminal,
+)
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
 )
 from HARK.distribution import MeanOneLogNormal, expected
 from HARK.interpolation import (
@@ -87,28 +92,92 @@ def make_lognormal_PrefShkDstn(
 
 ###############################################################################
 
-pref_shock_params = {
+# Make a dictionary of constructors for the preference shock model
+prefshk_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "PrefShkDstn": make_lognormal_PrefShkDstn,
+    "solution_terminal": make_basic_CRRA_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using construct_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make PrefShkDstn using make_lognormal_PrefShkDstn
+default_PrefShkDstn_params = {
     "PrefShkCount": 12,  # Number of points in discrete approximation to preference shock dist
     "PrefShk_tail_N": 4,  # Number of "tail points" on each end of pref shock dist
     "PrefShkStd": [0.30],  # Standard deviation of utility shocks
-    "aXtraCount": 48,
-    "CubicBool": False,  # pref shocks currently only compatible with linear cFunc
 }
 
-# Make a dictionary to specify a preference shock consumer
-init_preference_shocks = dict(
-    init_idiosyncratic_shocks,
-    **pref_shock_params,
-)
-init_preference_shocks["constructors"]["PrefShkDstn"] = make_lognormal_PrefShkDstn
+# Make a dictionary to specify an preference shocks consumer type
+init_preference_shocks = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": prefshk_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_preference_shocks.update(default_IncShkDstn_params)
+init_preference_shocks.update(default_aXtraGrid_params)
+init_preference_shocks.update(default_PrefShkDstn_params)
+
+
+# Specify default parameters that differ in "kinky preference" model compared to base PrefShockConsumerType
+kinky_pref_different_params = {
+    "Rboro": 1.20,  # Interest factor on assets when borrowing, a < 0
+    "Rsave": 1.02,  # Interest factor on assets when saving, a > 0
+    "BoroCnstArt": None,  # Kinked R only matters if borrowing is allowed
+}
 
 # Make a dictionary to specify a "kinky preference" consumer
-init_kinky_pref = dict(
-    init_kinked_R,
-    **pref_shock_params,
-)
-init_kinky_pref["BoroCnstArt"] = None
-init_kinky_pref["constructors"]["PrefShkDstn"] = make_lognormal_PrefShkDstn
+init_kinky_pref = init_preference_shocks.copy()  # Start with base above
+init_kinky_pref.update(kinky_pref_different_params)  # Change just a few params
+del init_kinky_pref["Rfree"]
 
 
 class PrefShockConsumerType(IndShockConsumerType):
diff --git a/HARK/ConsumptionSaving/ConsRepAgentModel.py b/HARK/ConsumptionSaving/ConsRepAgentModel.py
index c365b7c6e..217faa4f6 100644
--- a/HARK/ConsumptionSaving/ConsRepAgentModel.py
+++ b/HARK/ConsumptionSaving/ConsRepAgentModel.py
@@ -6,15 +6,20 @@
 """
 
 import numpy as np
-
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+)
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
     IndShockConsumerType,
-    init_idiosyncratic_shocks,
+    make_basic_CRRA_solution_terminal,
 )
 from HARK.ConsumptionSaving.ConsMarkovModel import MarkovConsumerType
 from HARK.distribution import MarkovProcess
 from HARK.interpolation import LinearInterp, MargValueFuncCRRA
+from HARK.utilities import make_assets_grid
 
 __all__ = ["RepAgentConsumerType", "RepAgentMarkovConsumerType"]
 
@@ -34,14 +39,13 @@ def solve_ConsRepAgent(
     CRRA : float
         Coefficient of relative risk aversion.
     IncShkDstn : distribution.Distribution
-        A discrete
-        approximation to the income process between the period being solved
-        and the one immediately following (in solution_next). Order:
+        A discrete approximation to the income process between the period being
+        solved and the one immediately following (in solution_next). Order:
         permanent shocks, transitory shocks.
     CapShare : float
         Capital's share of income in Cobb-Douglas production function.
     DeprFac : float
-        Depreciation rate of capital.
+        Depreciation rate for capital.
     PermGroFac : float
         Expected permanent income growth factor at the end of this period.
     aXtraGrid : np.array
@@ -134,10 +138,9 @@ def solve_ConsRepAgentMarkov(
     CRRA : float
         Coefficient of relative risk aversion.
     IncShkDstn : [distribution.Distribution]
-        A list of discrete
-        approximations to the income process between the period being solved
-        and the one immediately following (in solution_next). Order: event
-        probabilities, permanent shocks, transitory shocks.
+        A list of discrete approximations to the income process between the
+        period being solved and the one immediately following (in solution_next).
+        Order: event probabilities, permanent shocks, transitory shocks.
     CapShare : float
         Capital's share of income in Cobb-Douglas production function.
     DeprFac : float
@@ -224,6 +227,76 @@ def solve_ConsRepAgentMarkov(
     return solution_now
 
 
+###############################################################################
+
+# Make a dictionary of constructors for the representative agent model
+repagent_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "solution_terminal": make_basic_CRRA_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.00,  # Probability of unemployment while working
+    "IncUnemp": 0.0,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Make a dictionary to specify a representative agent consumer type
+init_rep_agent = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 0,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": repagent_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Interest factor on retained assets
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [1.0],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "DeprFac": 0.05,  # Depreciation rate for capital
+    "CapShare": 0.36,  # Capital's share in Cobb-Douglas production function
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 1,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_rep_agent.update(default_IncShkDstn_params)
+init_rep_agent.update(default_aXtraGrid_params)
+
+
 class RepAgentConsumerType(IndShockConsumerType):
     """
     A class for representing representative agents with inelastic labor supply.
@@ -286,6 +359,15 @@ def get_states(self):
         self.mNrmNow = self.Rfree * self.kNrmNow + self.wRte * self.shocks["TranShk"]
 
 
+###############################################################################
+
+# Define the default dictionary for a markov representative agent type
+init_markov_rep_agent = init_rep_agent.copy()
+init_markov_rep_agent["PermGroFac"] = [[0.97, 1.03]]
+init_markov_rep_agent["MrkvArray"] = np.array([[0.99, 0.01], [0.01, 0.99]])
+init_markov_rep_agent["Mrkv"] = 0
+
+
 class RepAgentMarkovConsumerType(RepAgentConsumerType):
     """
     A class for representing representative agents with inelastic labor supply
@@ -380,18 +462,3 @@ def get_controls(self):
         t = self.t_cycle[0]
         i = self.shocks["Mrkv"]
         self.controls["cNrm"] = self.solution[t].cFunc[i](self.mNrmNow)
-
-
-# Define the default dictionary for a representative agent type
-init_rep_agent = init_idiosyncratic_shocks.copy()
-init_rep_agent["cycles"] = 0
-init_rep_agent["DeprFac"] = 0.05
-init_rep_agent["CapShare"] = 0.36
-init_rep_agent["UnempPrb"] = 0.0
-init_rep_agent["LivPrb"] = [1.0]
-
-# Define the default dictionary for a markov representative agent type
-init_markov_rep_agent = init_rep_agent.copy()
-init_markov_rep_agent["PermGroFac"] = [[0.97, 1.03]]
-init_markov_rep_agent["MrkvArray"] = np.array([[0.99, 0.01], [0.01, 0.99]])
-init_markov_rep_agent["Mrkv"] = 0
diff --git a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py
index d75b4361e..cff779063 100644
--- a/HARK/ConsumptionSaving/ConsRiskyAssetModel.py
+++ b/HARK/ConsumptionSaving/ConsRiskyAssetModel.py
@@ -1,9 +1,7 @@
 """
 This file contains a class that adds a risky asset with a log-normal return
-factor to IndShockConsumerType.
-This class is not a fully specified model and therefore has no solution or
-simulation methods. It is meant as a container of methods for dealing with
-risky assets that will be useful to models what will inherit from it.
+factor to IndShockConsumerType. It is meant as a container of methods for dealing
+with risky assets that will be useful to models what will inherit from it.
 """
 
 import numpy as np
@@ -11,15 +9,19 @@
 from HARK import NullFunc
 from HARK.ConsumptionSaving.ConsIndShockModel import (
     ConsumerSolution,
-    IndShockConsumerType,  # PortfolioConsumerType inherits from it;
-    init_idiosyncratic_shocks,  # Baseline dictionary to build on
-    indshk_constructor_dict,
+    IndShockConsumerType,
+    make_basic_CRRA_solution_terminal,
 )
 from HARK.Calibration.Assets.AssetProcesses import (
     make_lognormal_RiskyDstn,
     combine_IncShkDstn_and_RiskyDstn,
     calc_ShareLimit_for_CRRA,
 )
+from HARK.Calibration.Income.IncomeProcesses import (
+    construct_lognormal_income_process_unemployment,
+    get_PermShkDstn_from_IncShkDstn,
+    get_TranShkDstn_from_IncShkDstn,
+)
 from HARK.distribution import (
     Bernoulli,
     expected,
@@ -37,6 +39,7 @@
     ValueFuncCRRA,
 )
 from HARK.rewards import UtilityFuncCRRA
+from HARK.utilities import make_assets_grid
 
 ###############################################################################
 
@@ -58,6 +61,92 @@ def make_simple_ShareGrid(ShareCount):
     return ShareGrid
 
 
+###############################################################################
+
+# Make a dictionary of constructors for the risky asset model
+risky_constructor_dict = {
+    "IncShkDstn": construct_lognormal_income_process_unemployment,
+    "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
+    "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
+    "aXtraGrid": make_assets_grid,
+    "RiskyDstn": make_lognormal_RiskyDstn,
+    "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
+    "ShareLimit": calc_ShareLimit_for_CRRA,
+    "ShareGrid": make_simple_ShareGrid,
+    "solution_terminal": make_basic_CRRA_solution_terminal,
+}
+
+# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
+default_IncShkDstn_params = {
+    "PermShkStd": [0.1],  # Standard deviation of log permanent income shocks
+    "PermShkCount": 7,  # Number of points in discrete approximation to permanent income shocks
+    "TranShkStd": [0.1],  # Standard deviation of log transitory income shocks
+    "TranShkCount": 7,  # Number of points in discrete approximation to transitory income shocks
+    "UnempPrb": 0.05,  # Probability of unemployment while working
+    "IncUnemp": 0.3,  # Unemployment benefits replacement rate while working
+    "T_retire": 0,  # Period of retirement (0 --> no retirement)
+    "UnempPrbRet": 0.005,  # Probability of "unemployment" while retired
+    "IncUnempRet": 0.0,  # "Unemployment" benefits when retired
+}
+
+# Default parameters to make aXtraGrid using make_assets_grid
+default_aXtraGrid_params = {
+    "aXtraMin": 0.001,  # Minimum end-of-period "assets above minimum" value
+    "aXtraMax": 20,  # Maximum end-of-period "assets above minimum" value
+    "aXtraNestFac": 3,  # Exponential nesting factor for aXtraGrid
+    "aXtraCount": 48,  # Number of points in the grid of "assets above minimum"
+    "aXtraExtra": None,  # Additional other values to add in grid (optional)
+}
+
+# Default parameters to make RiskyDstn with make_lognormal_RiskyDstn (and uniform ShareGrid)
+default_RiskyDstn_and_ShareGrid_params = {
+    "RiskyAvg": 1.080370891,  # Mean return factor of risky asset
+    "RiskyStd": 0.177196585,  # Stdev of log returns on risky asset
+    "RiskyCount": 5,  # Number of integration nodes to use in approximation of risky returns
+    "ShareCount": 25,  # Number of discrete points in the risky share approximation
+}
+# Risky return factor moments are based on SP500 real returns from Shiller's
+# "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
+
+# Make a dictionary to specify a risky asset consumer type
+init_risky_asset = {
+    # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "cycles": 1,  # Finite, non-cyclic model
+    "T_cycle": 1,  # Number of periods in the cycle for this agent type
+    "constructors": risky_constructor_dict,  # See dictionary above
+    # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
+    "CRRA": 2.0,  # Coefficient of relative risk aversion
+    "Rfree": 1.03,  # Return factor on risk free asset (not used by this type)
+    "DiscFac": 0.96,  # Intertemporal discount factor
+    "LivPrb": [0.98],  # Survival probability after each period
+    "PermGroFac": [1.01],  # Permanent income growth factor
+    "BoroCnstArt": 0.0,  # Artificial borrowing constraint
+    "vFuncBool": False,  # Whether to calculate the value function during solution
+    "CubicBool": False,  # Whether to use cubic spline interpolation when True
+    # (Uses linear spline interpolation for cFunc when False)
+    "AdjustPrb": 1.0,  # Probability that the agent can update their risky portfolio share each period
+    # TODO: This is not used in this file and should be moved to ConsPortfolioModel.py
+    "sim_common_Rrisky": True,  # Whether risky returns have a shared/common value across agents
+    # PARAMETERS REQUIRED TO SIMULATE THE MODEL
+    "AgentCount": 10000,  # Number of agents of this type
+    "T_age": None,  # Age after which simulated agents are automatically killed
+    "aNrmInitMean": 0.0,  # Mean of log initial assets
+    "aNrmInitStd": 1.0,  # Standard deviation of log initial assets
+    "pLvlInitMean": 0.0,  # Mean of log initial permanent income
+    "pLvlInitStd": 0.0,  # Standard deviation of log initial permanent income
+    "PermGroFacAgg": 1.0,  # Aggregate permanent income growth factor
+    # (The portion of PermGroFac attributable to aggregate productivity growth)
+    "NewbornTransShk": False,  # Whether Newborns have transitory shock
+    # ADDITIONAL OPTIONAL PARAMETERS
+    "PerfMITShk": False,  # Do Perfect Foresight MIT Shock
+    # (Forces Newborns to follow solution path of the agent they replaced if True)
+    "neutral_measure": False,  # Whether to use permanent income neutral measure (see Harmenberg 2021)
+}
+init_risky_asset.update(default_IncShkDstn_params)
+init_risky_asset.update(default_aXtraGrid_params)
+init_risky_asset.update(default_RiskyDstn_and_ShareGrid_params)
+
+
 class IndShockRiskyAssetConsumerType(IndShockConsumerType):
     """
     A consumer type that has access to a risky asset for his savings. The
@@ -186,7 +275,7 @@ def update_AdjustDstn(self):
     def update_ShareLimit(self):
         """
         Creates the attribute ShareLimit, representing the limiting lower bound of
-        risky portfolio share as mNrm goes to infinity.
+        risky portfolio share as mNrm goes to infinity, if it is allowed to adjust.
 
         Parameters
         ----------
@@ -341,6 +430,15 @@ def get_shocks(self):
 RiskyAssetConsumerType = IndShockRiskyAssetConsumerType
 
 
+###############################################################################
+
+# Make a dictionary to specify a consumer type with a fixed risky asset share
+init_risky_share_fixed = init_risky_asset.copy()
+init_risky_share_fixed["RiskyShareFixed"] = [
+    0.0
+]  # Fixed share of assets in the risky asset
+
+
 class FixedPortfolioShareRiskyAssetConsumerType(IndShockRiskyAssetConsumerType):
     time_vary_ = IndShockRiskyAssetConsumerType.time_vary_ + ["RiskyShareFixed"]
 
@@ -1789,33 +1887,3 @@ def calc_vPPnext(S, a):
 # Initial parameter sets
 
 # Base risky asset dictionary
-
-risky_constructor_dict = indshk_constructor_dict.copy()
-risky_constructor_dict["RiskyDstn"] = make_lognormal_RiskyDstn
-risky_constructor_dict["ShockDstn"] = combine_IncShkDstn_and_RiskyDstn
-risky_constructor_dict["ShareLimit"] = calc_ShareLimit_for_CRRA
-risky_constructor_dict["ShareGrid"] = make_simple_ShareGrid
-
-risky_asset_params = {
-    # Risky return factor moments. Based on SP500 real returns from Shiller's
-    # "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
-    "RiskyAvg": 1.080370891,
-    "RiskyStd": 0.177196585,
-    "ShareCount": 25,  # Number of discrete points in the risky share approximation
-    # Number of integration nodes to use in approximation of risky returns
-    "RiskyCount": 5,
-    # Probability that the agent can adjust their portfolio each period
-    "AdjustPrb": 1.0,
-    # When simulating the model, should all agents get the same risky return in
-    # a given period?
-    "sim_common_Rrisky": True,
-    "constructors": risky_constructor_dict,
-}
-
-
-# Make a dictionary to specify a risky asset consumer type
-init_risky_asset = init_idiosyncratic_shocks.copy()
-init_risky_asset.update(risky_asset_params)
-
-init_risky_share_fixed = init_risky_asset.copy()
-init_risky_share_fixed["RiskyShareFixed"] = [0.0]
diff --git a/HARK/ConsumptionSaving/ConsRiskyContribModel.py b/HARK/ConsumptionSaving/ConsRiskyContribModel.py
index 8414cb50a..aceda63ab 100644
--- a/HARK/ConsumptionSaving/ConsRiskyContribModel.py
+++ b/HARK/ConsumptionSaving/ConsRiskyContribModel.py
@@ -43,7 +43,7 @@
 from HARK.ConsumptionSaving.ConsRiskyAssetModel import (
     RiskyAssetConsumerType,
     init_risky_asset,
-    risky_asset_params,
+    risky_constructor_dict,
 )
 from HARK.distribution import calc_expectation
 from HARK.interpolation import BilinearInterp  # 2D interpolator
@@ -2025,6 +2025,21 @@ def solveRiskyContrib(
     # Grid for finding the optimal rebalancing flow
     "dCount": 20,
 }
+risky_asset_params = {
+    # Risky return factor moments. Based on SP500 real returns from Shiller's
+    # "chapter 26" data, which can be found at https://www.econ.yale.edu/~shiller/data.htm
+    "RiskyAvg": 1.080370891,
+    "RiskyStd": 0.177196585,
+    "ShareCount": 25,  # Number of discrete points in the risky share approximation
+    # Number of integration nodes to use in approximation of risky returns
+    "RiskyCount": 5,
+    # Probability that the agent can adjust their portfolio each period
+    "AdjustPrb": 1.0,
+    # When simulating the model, should all agents get the same risky return in
+    # a given period?
+    "sim_common_Rrisky": True,
+    "constructors": risky_constructor_dict,
+}
 
 # Infinite horizon version
 init_risky_contrib = init_risky_asset.copy()
diff --git a/HARK/ConsumptionSaving/TractableBufferStockModel.py b/HARK/ConsumptionSaving/TractableBufferStockModel.py
index 94846dc70..08f5c6bb8 100644
--- a/HARK/ConsumptionSaving/TractableBufferStockModel.py
+++ b/HARK/ConsumptionSaving/TractableBufferStockModel.py
@@ -299,6 +299,19 @@ def add_to_stable_arm_points(
     return solution_now
 
 
+###############################################################################
+
+# Define a dictionary for the tractable buffer stock model
+init_tractable = {
+    "cycles": 0,  # infinite horizon
+    "UnempPrb": 0.00625,  # Probability of becoming permanently unemployed
+    "DiscFac": 0.975,  # Intertemporal discount factor
+    "Rfree": 1.01,  # Risk-free interest factor on assets
+    "PermGroFac": 1.0025,  # Permanent income growth factor (uncompensated)
+    "CRRA": 1.0,  # Coefficient of relative risk aversion
+}
+
+
 class TractableConsumerType(AgentType):
     """
     Parameters
@@ -716,13 +729,3 @@ def get_poststates(self):
         """
         self.state_now["aLvl"] = self.state_now["mLvl"] - self.controls["cLvlNow"]
         return None
-
-
-init_tractable = {
-    "cycles": 0,  # infinite horizon
-    "UnempPrb": 0.00625,  # Probability of becoming unemployed
-    "DiscFac": 0.975,  # Intertemporal discount factor
-    "Rfree": 1.01,  # Risk-free interest factor on assets
-    "PermGroFac": 1.0025,  # Permanent income growth factor (uncompensated)
-    "CRRA": 1.0,  # Coefficient of relative risk aversion
-}
diff --git a/HARK/utilities.py b/HARK/utilities.py
index 0424fd96f..7a633c7cb 100644
--- a/HARK/utilities.py
+++ b/HARK/utilities.py
@@ -117,7 +117,7 @@ def distance(self, other):
 # =======================================================
 
 
-def construct_assets_grid(aXtraMin, aXtraMax, aXtraCount, aXtraExtra, aXtraNestFac):
+def make_assets_grid(aXtraMin, aXtraMax, aXtraCount, aXtraExtra, aXtraNestFac):
     """
     Constructs the base grid of post-decision states, representing end-of-period
     assets above the absolute minimum. Has three modes, depending on aXtraNestFac:
diff --git a/examples/ConsumptionSaving/example_ConsGenIncProcessModel.ipynb b/examples/ConsumptionSaving/example_ConsGenIncProcessModel.ipynb
index d80dc70e0..02fc41c42 100644
--- a/examples/ConsumptionSaving/example_ConsGenIncProcessModel.ipynb
+++ b/examples/ConsumptionSaving/example_ConsGenIncProcessModel.ipynb
@@ -15,8 +15,8 @@
     "    IndShockConsumerType,\n",
     "    PersistentShockConsumerType,\n",
     "    init_explicit_perm_inc,\n",
-    "    indshk_constructor_dict,\n",
     ")\n",
+    "from HARK.ConsumptionSaving.ConsIndShockModel import indshk_constructor_dict\n",
     "\n",
     "\n",
     "def mystr(number):\n",
@@ -32,27 +32,6 @@
     "do_simulation = True"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Display information about the pLvlGrid used in these examples\n",
-    "print(\n",
-    "    \"The infinite horizon examples presented here use a grid of persistent income levels (pLvlGrid)\"\n",
-    ")\n",
-    "print(\n",
-    "    \"based on percentiles of the long run distribution of pLvl for the given parameters. These percentiles\"\n",
-    ")\n",
-    "print(\n",
-    "    \"are specified in the attribute pLvlPctiles. Here, the lowest percentile is \"\n",
-    "    + str(init_explicit_perm_inc[\"pLvlPctiles\"][0] * 100)\n",
-    "    + \" and the highest\"\n",
-    ")\n",
-    "print(\"percentile is \" + str(init_explicit_perm_inc[\"pLvlPctiles\"][-1] * 100) + \".\\n\")"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -73,6 +52,27 @@
     ")"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Display information about the pLvlGrid used in these examples\n",
+    "print(\n",
+    "    \"The infinite horizon examples presented here use a grid of persistent income levels (pLvlGrid)\"\n",
+    ")\n",
+    "print(\n",
+    "    \"based on percentiles of the long run distribution of pLvl for the given parameters. These percentiles\"\n",
+    ")\n",
+    "print(\n",
+    "    \"are specified in the attribute pLvlPctiles. Here, the lowest percentile is \"\n",
+    "    + mystr(ExplicitExample.pLvlPctiles[0] * 100)\n",
+    "    + \" and the highest\"\n",
+    ")\n",
+    "print(\"percentile is \" + mystr(ExplicitExample.pLvlPctiles[-1] * 100) + \".\\n\")"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
diff --git a/examples/ConsumptionSaving/example_ConsMarkovModel.ipynb b/examples/ConsumptionSaving/example_ConsMarkovModel.ipynb
index 963c656db..60cd050a6 100644
--- a/examples/ConsumptionSaving/example_ConsMarkovModel.ipynb
+++ b/examples/ConsumptionSaving/example_ConsMarkovModel.ipynb
@@ -184,14 +184,15 @@
     "        employed_income_dist,\n",
     "        unemployed_income_dist,\n",
     "    ]\n",
-    "]"
+    "]\n",
+    "SerialUnemploymentExample.assign_parameters(MrkvArray=[MrkvArray])"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Note that $\\texttt{MarkovConsumerType}$ currently has no method to automatically construct a valid IncomeDstn - $\\texttt{IncomeDstn}$ is manually constructed in each case. Writing a method to supersede $\\texttt{IndShockConsumerType.update_income_process}$ for the “Markov model” would be a welcome contribution!"
+    "Note that $\\texttt{MarkovConsumerType}$ currently has no method to automatically construct a valid IncomeDstn - $\\texttt{IncomeDstn}$ is manually constructed in each case. Writing a method to supersede $\\texttt{IndShockConsumerType.update\\_income_process}$ for the “Markov model” would be a welcome contribution!"
    ]
   },
   {
@@ -201,13 +202,9 @@
    "outputs": [],
    "source": [
     "# Interest factor, permanent growth rates, and survival probabilities are constant arrays\n",
-    "SerialUnemploymentExample.assign_parameters(\n",
-    "    Rfree=np.array(4 * [SerialUnemploymentExample.Rfree])\n",
-    ")\n",
-    "SerialUnemploymentExample.PermGroFac = [\n",
-    "    np.array(4 * SerialUnemploymentExample.PermGroFac)\n",
-    "]\n",
-    "SerialUnemploymentExample.LivPrb = [SerialUnemploymentExample.LivPrb * np.ones(4)]"
+    "SerialUnemploymentExample.assign_parameters(Rfree=np.array(1.03 * np.ones(4)))\n",
+    "SerialUnemploymentExample.PermGroFac = [1.0 * np.ones(4)]\n",
+    "SerialUnemploymentExample.LivPrb = [0.98 * np.ones(4)]"
    ]
   },
   {
@@ -321,6 +318,7 @@
     "init_unemployment_immunity = copy(init_indshk_markov)\n",
     "init_unemployment_immunity[\"MrkvArray\"] = [MrkvArray]\n",
     "ImmunityExample = MarkovConsumerType(**init_unemployment_immunity)\n",
+    "ImmunityExample.assign_parameters(MrkvArray=[MrkvArray])\n",
     "ImmunityExample.assign_parameters(\n",
     "    Rfree=np.array(np.array(StateCount * [1.03])),  # Interest factor same in all states\n",
     "    PermGroFac=[\n",
@@ -415,6 +413,7 @@
     "init_serial_growth = copy(init_indshk_markov)\n",
     "init_serial_growth[\"MrkvArray\"] = [MrkvArray]\n",
     "SerialGroExample = MarkovConsumerType(**init_serial_growth)\n",
+    "SerialGroExample.assign_parameters(MrkvArray=[MrkvArray])\n",
     "SerialGroExample.assign_parameters(\n",
     "    Rfree=np.array(\n",
     "        np.array(StateCount * [1.03])\n",


Test patch:
diff --git a/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py b/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
index b9c915aed..b740c1ff8 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsGenIncProcessModel.py
@@ -44,13 +44,6 @@
     "CubicBool": False,  # Use cubic spline interpolation when True, linear interpolation when False
     "vFuncBool": True,  # Whether to calculate the value function during solution
     # More parameters specific to "Explicit Permanent income" shock model
-    "pLvlPctiles": np.concatenate(
-        (
-            [0.001, 0.005, 0.01, 0.03],
-            np.linspace(0.05, 0.95, num=19),
-            [0.97, 0.99, 0.995, 0.999],
-        )
-    ),
     "PermGroFac": [
         1.0
     ],  # Permanent income growth factor - long run permanent income growth doesn't work yet
@@ -65,7 +58,7 @@ def setUp(self):
     def test_solution(self):
         pLvlGrid = self.agent.pLvlGrid[0]
         self.assertAlmostEqual(
-            self.agent.pLvlGrid[0][0], 0.283792, places=HARK_PRECISION
+            self.agent.pLvlGrid[0][0], 0.28063, places=HARK_PRECISION
         )
 
         self.assertAlmostEqual(self.agent.solution[0].mLvlMin(pLvlGrid[0]), 0.0)
@@ -93,7 +86,7 @@ def test_solution(self):
 
         self.assertAlmostEqual(
             self.agent.solution[0].cFunc(10, pLvlGrid[1]).tolist(),
-            5.29136,
+            5.27723,
             places=HARK_PRECISION,
         )
 
@@ -103,7 +96,7 @@ def test_value(self):
         self.assertTrue(self.agent.vFuncBool)
         self.assertAlmostEqual(
             self.agent.solution[0].vFunc(10, pLvlGrid[3]),
-            -0.36557,
+            -0.36683,
             places=HARK_PRECISION,
         )
 
diff --git a/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py b/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py
index 895b8d88c..c5f531e07 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsMarkovModel.py
@@ -58,13 +58,17 @@ def setUp(self):
 
         init_serial_unemployment = copy(init_indshk_markov)
         init_serial_unemployment["MrkvArray"] = [MrkvArray]
-        init_serial_unemployment["UnempPrb"] = (
-            0.0  # to make income distribution when employed
-        )
+        init_serial_unemployment["UnempPrb"] = 0.0
+        # to make income distribution when employed
         init_serial_unemployment["global_markov"] = False
+        init_serial_unemployment["Rfree"] = np.array([1.03, 1.03, 1.03, 1.03])
+        init_serial_unemployment["LivPrb"] = [np.array([0.98, 0.98, 0.98, 0.98])]
+        init_serial_unemployment["PermGroFac"] = [np.array([1.01, 1.01, 1.01, 1.01])]
+
         self.model = MarkovConsumerType(**init_serial_unemployment)
         self.model.cycles = 0
         self.model.vFuncBool = False  # for easy toggling here
+        self.model.MrkvArray = [MrkvArray]
 
         # Replace the default (lognormal) income distribution with a custom one
         employed_income_dist = DiscreteDistributionLabeled(
@@ -88,6 +92,7 @@ def setUp(self):
 
     def test_check_markov_inputs(self):
         # check Rfree
+        self.model.Rfree = 1.03
         self.assertRaises(ValueError, self.model.check_markov_inputs)
         # fix Rfree
         self.model.Rfree = np.array(4 * [self.model.Rfree])
@@ -98,24 +103,20 @@ def test_check_markov_inputs(self):
         # then fix it back
         self.model.MrkvArray = self.MrkvArray
         # check LivPrb
+        self.model.LivPrb = [0.98]
         self.assertRaises(ValueError, self.model.check_markov_inputs)
         # fix LivPrb
         self.model.LivPrb = [np.array(4 * self.model.LivPrb)]
         # check PermGroFac
+        self.model.PermGroFac = [1.01]
         self.assertRaises(ValueError, self.model.check_markov_inputs)
         # fix PermGroFac
         self.model.PermGroFac = [np.array(4 * self.model.PermGroFac)]
 
     def test_solve(self):
-        self.model.Rfree = np.array(4 * [self.model.Rfree])
-        self.model.LivPrb = [np.array(4 * self.model.LivPrb)]
-        self.model.PermGroFac = [np.array(4 * self.model.PermGroFac)]
         self.model.solve()
 
     def test_simulation(self):
-        self.model.Rfree = np.array(4 * [self.model.Rfree])
-        self.model.LivPrb = [np.array(4 * self.model.LivPrb)]
-        self.model.PermGroFac = [np.array(4 * self.model.PermGroFac)]
         self.model.solve()
         self.model.T_sim = 120
         self.model.MrkvPrbsInit = [0.25, 0.25, 0.25, 0.25]
diff --git a/HARK/ConsumptionSaving/tests/test_ConsMedModel.py b/HARK/ConsumptionSaving/tests/test_ConsMedModel.py
index 1f55a0af0..64311c847 100644
--- a/HARK/ConsumptionSaving/tests/test_ConsMedModel.py
+++ b/HARK/ConsumptionSaving/tests/test_ConsMedModel.py
@@ -16,7 +16,7 @@ def test_solution(self):
         pLvl = 2.0
         Shk = 1.5
         self.assertAlmostEqual(
-            cFunc(mLvl, pLvl, Shk).tolist(), 4.16230, places=HARK_PRECISION
+            cFunc(mLvl, pLvl, Shk).tolist(), 4.16194, places=HARK_PRECISION
         )
         self.assertAlmostEqual(
             MedFunc(mLvl, pLvl, Shk).tolist(), 2.45432, places=HARK_PRECISION
diff --git a/HARK/ConsumptionSaving/tests/test_KinkedRconsumerType.py b/HARK/ConsumptionSaving/tests/test_KinkedRconsumerType.py
index c53d0e19a..bde126727 100644
--- a/HARK/ConsumptionSaving/tests/test_KinkedRconsumerType.py
+++ b/HARK/ConsumptionSaving/tests/test_KinkedRconsumerType.py
@@ -17,20 +17,20 @@ def test_liquidity_constraint(self):
         KinkyExample.solve()
 
         self.assertAlmostEqual(
-            KinkyExample.solution[0].cFunc(1).tolist(), 0.96233, places=HARK_PRECISION
+            KinkyExample.solution[0].cFunc(1).tolist(), 0.96161, places=HARK_PRECISION
         )
 
         self.assertAlmostEqual(
-            KinkyExample.solution[0].cFunc(4).tolist(), 1.34377, places=HARK_PRECISION
+            KinkyExample.solution[0].cFunc(4).tolist(), 1.34274, places=HARK_PRECISION
         )
 
         KinkyExample.BoroCnstArt = -0.2
         KinkyExample.solve()
 
         self.assertAlmostEqual(
-            KinkyExample.solution[0].cFunc(1).tolist(), 0.93469, places=HARK_PRECISION
+            KinkyExample.solution[0].cFunc(1).tolist(), 0.93444, places=HARK_PRECISION
         )
 
         self.assertAlmostEqual(
-            KinkyExample.solution[0].cFunc(4).tolist(), 1.34014, places=HARK_PRECISION
+            KinkyExample.solution[0].cFunc(4).tolist(), 1.33927, places=HARK_PRECISION
         )
diff --git a/HARK/ConsumptionSaving/tests/test_modelInits.py b/HARK/ConsumptionSaving/tests/test_modelInits.py
index 616e789ec..d4a55050a 100644
--- a/HARK/ConsumptionSaving/tests/test_modelInits.py
+++ b/HARK/ConsumptionSaving/tests/test_modelInits.py
@@ -4,7 +4,6 @@
 
 # Bring in modules we need
 import unittest
-from copy import copy
 
 import numpy as np
 
@@ -12,7 +11,6 @@
     IndShockConsumerType,
     KinkedRconsumerType,
     PerfForesightConsumerType,
-    init_idiosyncratic_shocks,
     init_kinked_R,
     init_lifecycle,
 )
@@ -88,7 +86,7 @@ def test_MarkovConsumerType(self):
             )
 
             # Make a consumer with serially correlated unemployment, subject to boom and bust cycles
-            init_serial_unemployment = copy(init_idiosyncratic_shocks)
+            init_serial_unemployment = {}
             init_serial_unemployment["MrkvArray"] = [MrkvArray]
             init_serial_unemployment["UnempPrb"] = (
                 0.0  # to make income distribution when employed
diff --git a/HARK/ConsumptionSaving/tests/test_modelcomparisons.py b/HARK/ConsumptionSaving/tests/test_modelcomparisons.py
index 27da5e63c..0b4669e16 100644
--- a/HARK/ConsumptionSaving/tests/test_modelcomparisons.py
+++ b/HARK/ConsumptionSaving/tests/test_modelcomparisons.py
@@ -108,14 +108,6 @@ def setUp(self):
         TBSType.solve()
 
         # Set up and solve Markov
-        MrkvArray = [
-            np.array(
-                [
-                    [1.0 - base_primitives["UnempPrb"], base_primitives["UnempPrb"]],
-                    [0.0, 1.0],
-                ]
-            )
-        ]
         Markov_primitives = {
             "CRRA": base_primitives["CRRA"],
             "Rfree": np.array(2 * [base_primitives["Rfree"]]),
@@ -128,6 +120,8 @@ def setUp(self):
                     ]
                 )
             ],
+            "Mrkv_p11": [1.0 - base_primitives["UnempPrb"]],
+            "Mrkv_p22": [1.0],
             "BoroCnstArt": None,
             "PermShkStd": [0.0],
             "PermShkCount": 1,
@@ -154,7 +148,6 @@ def setUp(self):
             "tax_rate": 0.0,
             "vFuncBool": False,
             "CubicBool": True,
-            "MrkvArray": MrkvArray,
             "T_cycle": 1,
         }
 
diff --git a/HARK/tests/test_HARKutilities.py b/HARK/tests/test_HARKutilities.py
index 0a4445973..172f18413 100644
--- a/HARK/tests/test_HARKutilities.py
+++ b/HARK/tests/test_HARKutilities.py
@@ -15,7 +15,7 @@
     CRRAutilityPPP,
     CRRAutilityPPPP,
 )
-from HARK.utilities import construct_assets_grid, make_figs
+from HARK.utilities import make_assets_grid, make_figs
 
 
 class testsForHARKutilities(unittest.TestCase):
@@ -76,7 +76,7 @@ def test_asset_grid(self):
             "aXtraNestFac": -1,
         }
 
-        aXtraGrid = construct_assets_grid(**params)
+        aXtraGrid = make_assets_grid(**params)
         test = np.unique(np.diff(aXtraGrid).round(decimals=3))
         self.assertEqual(test.size, 1)