National registry data on income and wealth from Scandinavian countries (esp. Norway) have recently become available (with a lot of security) to some (lucky!) researchers. These data offer a uniquely powerful tool for testing (and improving) our models of consumption and saving behavior over the life cycle.
This notebook is an example of how to construct a life cycle model with the HARK toolkit that makes predictions that can be compared to the raw data statistics=.
For example, some papers have tabulated information about the growth rate of assets at different ages over the life cycle.
The default parameters of the HARK life cycle model have not been optmized to match features of the Norwegian data; a first step in a real "structural" estimation would be to use Norwegian calibrate the inputs to the model (like the profile of income, and the magnitude of income shocks, over the life cycle), and then to find the values of parameters like the time preference rate that allow the model to fit the data best. (See SolvingMicroDSOPs for how this can be done, and search for the corresponding HARK content using our documentation).
# Initial imports and notebook setup, click arrow to show
import HARK.ConsumptionSaving.ConsIndShockModel as cShksModl # The consumption-saving micro model
from HARK.utilities import plot_funcs_der, plot_funcs # Some tools
import pandas as pd
import numpy as np
# ---------------------------------------------------------------------------------
# - Define all of the model parameters for SolvingMicroDSOPs and ConsumerExamples -
# ---------------------------------------------------------------------------------
exp_nest = 3 # Number of times to "exponentially nest" when constructing a_grid
aXtraMin = 0.001 # Minimum end-of-period "assets above minimum" value
aXtraMax = 20 # Maximum end-of-period "assets above minimum" value
aXtraHuge = None # A very large value of assets to add to the grid, not used
aXtraExtra = None # Some other value of assets to add to the grid, not used
aXtraCount = 8 # Number of points in the grid of "assets above minimum"
BoroCnstArt = 0.0 # Artificial borrowing constraint; imposed minimum level of end-of period assets
CubicBool = True # Use cubic spline interpolation when True, linear interpolation when False
vFuncBool = False # Whether to calculate the value function during solution
Rfree = 1.03 # Interest factor on assets
PermShkCount = 7 # Number of points in discrete approximation to permanent income shocks
TranShkCount = 7 # Number of points in discrete approximation to transitory income shocks
UnempPrb = 0.005 # Probability of unemployment while working
UnempPrbRet = 0.000 # Probability of "unemployment" while retired
IncUnemp = 0.0 # Unemployment benefits replacement rate
IncUnempRet = 0.0 # "Unemployment" benefits when retired
final_age = 90 # Age at which the problem ends (die with certainty)
retirement_age = 65 # Age at which the consumer retires
initial_age = 25 # Age at which the consumer enters the model
TT = final_age - initial_age # Total number of periods in the model
retirement_t = retirement_age - initial_age - 1
CRRA_start = 4.0 # Initial guess of the coefficient of relative risk aversion during estimation (rho)
DiscFacAdj_start = 0.99 # Initial guess of the adjustment to the discount factor during estimation (beth)
DiscFacAdj_bound = [0.0001,15.0] # Bounds for beth; if violated, objective function returns "penalty value"
CRRA_bound = [0.0001,15.0] # Bounds for rho; if violated, objective function returns "penalty value"
# Expected growth rates of permanent income over the lifecycle, starting from age 25
PermGroFac = [ 1.025, 1.025, 1.025, 1.025, 1.025, 1.025, 1.025, 1.025,
1.025, 1.025, 1.025, 1.025, 1.025, 1.025, 1.025, 1.025,
1.025, 1.025, 1.025, 1.025, 1.025, 1.025, 1.025, 1.025,
1.025, 1.01 , 1.01 , 1.01 , 1.01 , 1.01 , 1.01 , 1.01 ,
1.01 , 1.01 , 1.01 , 1.01 , 1.01 , 1.01 , 1.01 , 0.7 , # <-- This represents retirement
1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. , 1. ]
# Age-varying discount factors over the lifecycle, lifted from Cagetti (2003)
DiscFac_timevary = [1.064914 , 1.057997 , 1.051422 , 1.045179 , 1.039259 ,
1.033653 , 1.028352 , 1.023348 , 1.018632 , 1.014198 ,
1.010037 , 1.006143 , 1.002509 , 0.9991282, 0.9959943,
0.9931012, 0.9904431, 0.9880143, 0.9858095, 0.9838233,
0.9820506, 0.9804866, 0.9791264, 0.9779656, 0.9769995,
0.9762239, 0.9756346, 0.9752274, 0.9749984, 0.9749437,
0.9750595, 0.9753422, 0.9757881, 0.9763936, 0.9771553,
0.9780698, 0.9791338, 0.9803439, 0.981697 , 0.8287214,
0.9902111, 0.9902111, 0.9902111, 0.9902111, 0.9902111,
0.9902111, 0.9902111, 0.9902111, 0.9902111, 0.9902111,
0.9902111, 0.9902111, 0.9902111, 0.9902111, 0.9902111,
0.9902111, 0.9902111, 0.9902111, 0.9902111, 0.9902111,
0.9902111, 0.9902111, 0.9902111, 0.9902111, 0.9902111]
# Survival probabilities over the lifecycle, starting from age 25
LivPrb = [ 1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. , # <-- automatic survival to age 65
0.98438596, 0.98438596, 0.98438596, 0.98438596, 0.98438596,
0.97567062, 0.97567062, 0.97567062, 0.97567062, 0.97567062,
0.96207901, 0.96207901, 0.96207901, 0.96207901, 0.96207901,
0.93721595, 0.93721595, 0.93721595, 0.93721595, 0.93721595,
0.63095734, 0.63095734, 0.63095734, 0.63095734, 0.63095734]
# Standard deviations of permanent income shocks by age, starting from age 25
PermShkStd = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0, 0.0, 0.0, # <-- no permanent income shocks after retirement
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
# Standard deviations of transitory income shocks by age, starting from age 25
TranShkStd = [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1,
0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0, 0.0, 0.0, # <-- no transitory income shocs after retirement
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
# Age groups for the estimation: calculate average wealth-to-permanent income ratio
# for consumers within each of these age groups, compare actual to simulated data
empirical_cohort_age_groups = [[ 26,27,28,29,30 ],
[ 31,32,33,34,35 ],
[ 36,37,38,39,40 ],
[ 41,42,43,44,45 ],
[ 46,47,48,49,50 ],
[ 51,52,53,54,55 ],
[ 56,57,58,59,60 ]]
initial_wealth_income_ratio_vals = [0.17, 0.5, 0.83] # Three point discrete distribution of initial w
initial_wealth_income_ratio_probs = [0.33333, 0.33333, 0.33334] # Equiprobable discrete distribution of initial w
num_agents = 10000 # Number of agents to simulate
bootstrap_size = 50 # Number of re-estimations to do during bootstrap
seed = 31382 # Just an integer to seed the estimation
# Dictionary that can be passed to ConsumerType to instantiate
init_consumer_objects = {"CRRA":CRRA_start,
"Rfree":Rfree,
"PermGroFac":PermGroFac,
"BoroCnstArt":BoroCnstArt,
"PermShkStd":PermShkStd,
"PermShkCount":PermShkCount,
"TranShkStd":TranShkStd,
"TranShkCount":TranShkCount,
"T_cycle":TT,
"UnempPrb":UnempPrb,
"UnempPrbRet":UnempPrbRet,
"T_retire":retirement_t,
"T_age":TT+1,
"IncUnemp":IncUnemp,
"IncUnempRet":IncUnempRet,
"aXtraMin":aXtraMin,
"aXtraMax":aXtraMax,
"aXtraCount":aXtraCount,
"aXtraExtra":[aXtraExtra,aXtraHuge],
"aXtraNestFac":exp_nest,
"LivPrb":LivPrb,
"DiscFac":DiscFac_timevary,
'AgentCount':num_agents,
'seed':seed,
'tax_rte':0.0,
'vFuncBool':vFuncBool,
'CubicBool':CubicBool
}
# Set up default values for CRRA, DiscFac, and simulation variables in the dictionary
init_consumer_objects["CRRA"]= 2.00 # Default coefficient of relative risk aversion (rho)
init_consumer_objects["DiscFac"]= 0.97 # Default intertemporal discount factor (beta)
init_consumer_objects["PermGroFacAgg"]= 1.0 # Aggregate permanent income growth factor
init_consumer_objects["aNrmInitMean"]= -10.0 # Mean of log initial assets
init_consumer_objects["aNrmInitStd"]= 1.0 # Standard deviation of log initial assets
init_consumer_objects["pLvlInitMean"]= 0.0 # Mean of log initial permanent income
init_consumer_objects["pLvlInitStd"]= 0.0 # Standard deviation of log initial permanent income
# Make an instance of a lifecycle consumer to be used for estimation
LifeCyclePop = cShksModl.IndShockConsumerType(**init_consumer_objects)
# Solve and simulate the model (ignore the "warning" message)
LifeCyclePop.solve() # Obtain consumption rules by age
LifeCyclePop.unpack('cFunc') # Expose the consumption rules
# Which variables do we want to track
LifeCyclePop.track_vars = ['aNrm','aLvl','pLvl','mNrm','cNrm','TranShk']
LifeCyclePop.T_sim = 120 # Nobody lives to be older than 145 years (=25+120)
LifeCyclePop.initialize_sim() # Construct the age-25 distribution of income and assets
LifeCyclePop.simulate() # Simulate a population behaving according to this model
{'aNrm': array([[0.15127186, 0.15125337, 0.15137639, ..., 0.15125232, 0.15129283, 0.15126489], [0.15720283, 0.2640323 , 0.27109329, ..., 0.2484313 , 0.24846159, 0.22613419], [0.29023075, 0.27589679, 0.45803117, ..., 0.34902784, 0.29628105, 0.21172888], ..., [0. , 0.3464155 , 0.49892497, ..., 0.33524317, 0.15485401, 0.40136673], [0.15125389, 0.41011742, 0.7044481 , ..., 0.49148864, 0.10016124, 0.39917186], [0.27813171, 0.36985598, 0.61672583, ..., 0.46655336, 0.05614685, 0.39320286]]), 'aLvl': array([[0.18085553, 0.13184619, 0.1425343 , ..., 0.14869039, 0.15431001, 0.16007211], [0.18476294, 0.21671094, 0.30517834, ..., 0.21288725, 0.22090092, 0.20859615], [0.34791542, 0.23963374, 0.48550224, ..., 0.31650528, 0.2686692 , 0.23350374], ..., [0. , 0.55354082, 0.60636705, ..., 0.3272484 , 0.1127017 , 0.41864968], [0.16006047, 0.64423066, 0.74629729, ..., 0.48933561, 0.07289667, 0.44060239], [0.32520793, 0.64194659, 0.72191836, ..., 0.45664162, 0.0408633 , 0.42666246]]), 'pLvl': array([[1.19556632, 0.87169091, 0.94158876, ..., 0.98306182, 1.01994264, 1.05822383], [1.17531561, 0.82077437, 1.12573181, ..., 0.85692606, 0.88907473, 0.9224441 ], [1.1987545 , 0.868563 , 1.05997643, ..., 0.90681958, 0.90680522, 1.1028431 ], ..., [1.50144644, 1.59791009, 1.21534715, ..., 0.97615233, 0.72779323, 1.04306025], [1.05822383, 1.5708444 , 1.05940707, ..., 0.99561938, 0.72779323, 1.10379122], [1.16925872, 1.73566637, 1.17056611, ..., 0.9787554 , 0.72779323, 1.08509501]]), 'mNrm': array([[1.00006893, 1.00002717, 1.00030507, ..., 1.0000248 , 1.00011634, 1.00005319], [1.0131983 , 1.20305693, 1.21380651, ..., 1.1787878 , 1.17883567, 1.14264056], [1.24214623, 1.22089454, 1.46881532, ..., 1.32519867, 1.25097791, 1.1180577 ], ..., [1.01134323, 1.31772662, 1.5172135 , ..., 1.30621764, 1.23465048, 1.3942102 ], [1.00002834, 1.40055728, 1.76180299, ..., 1.51081641, 1.15949963, 1.39072821], [1.22439054, 1.34621133, 1.65674534, ..., 1.47885994, 1.10316608, 1.38213531]]), 'cNrm': array([[0.84879708, 0.8487738 , 0.84892869, ..., 0.84877248, 0.8488235 , 0.84878831], [0.85599547, 0.93902463, 0.94271322, ..., 0.9303565 , 0.93037408, 0.91650637], [0.95191548, 0.94499775, 1.01078415, ..., 0.97617083, 0.95469687, 0.90632881], ..., [1.01134323, 0.97131112, 1.01828852, ..., 0.97097446, 1.07979647, 0.99284347], [0.84877445, 0.99043985, 1.05735489, ..., 1.01932777, 1.05933838, 0.99155636], [0.94625882, 0.97635535, 1.04001951, ..., 1.01230658, 1.04701923, 0.98893245]]), 'TranShk': array([[1. , 1. , 1. , ..., 1. , 1. , 1. ], [0.85470368, 1.0376015 , 1.08339327, ..., 1.00006632, 1.00006632, 0.96390423], [1.08339327, 0.96390423, 1.1722675 , ..., 1.08339327, 1.00006632, 0.92323938], ..., [1. , 0.85470368, 1.08339327, ..., 1.08339327, 1. , 1.00006632], [1. , 1.0376015 , 1.1722675 , ..., 1.1722675 , 1. , 1.00006632], [1.08339327, 0.96390423, 1.00006632, ..., 0.96390423, 1. , 0.96390423]])}
# Plot the consumption functions during working life
print('Consumption as a function of market resources while working:')
mMin = min([LifeCyclePop.solution[t].mNrmMin for t in range(LifeCyclePop.T_cycle)])
plot_funcs(LifeCyclePop.cFunc[:LifeCyclePop.T_retire],mMin,5)
Consumption as a function of market resources while working:
# Define the saving rate function
def savRteFunc(SomeType, m, t):
"""
Parameters:
----------
SomeType:
Agent type that has been solved and simulated.
m:
normalized market resources of agent
t:
age of agent (from starting in the workforce)
Returns:
--------
savRte: float
"""
inc = (SomeType.Rfree -1.)*(m-1.)+1. # Normalized by permanent labor income
cns = SomeType.solution[t].cFunc(m) # Consumption (normalized)
sav = inc - cns # Flow of saving this period
savRte = sav / inc # Saving Rate
return savRte
# Create a matrix gathering useful data:
# 't_now', 'aNrm_hist', 'cNrm_hist', employment-status in date t and date t-1,
# aLvlGro_hist, Saving rate
w, h = 1, LifeCyclePop.T_cycle
giant_list = [[0 for x in range(w)] for y in range(h)]
savRte_list = []
import warnings
warnings.filterwarnings("ignore") # Suppress some disturbing but harmless warnings
for t in range(1,LifeCyclePop.T_cycle+1):
#aLvlGro[0] = 0 # set the first growth rate to 0, since there is no data for period 0
aLvlGroNow = np.log((LifeCyclePop.history['aNrm'][t] *LifeCyclePop.history['pLvl'][t])/ \
LifeCyclePop.history['aNrm'][t-1] *LifeCyclePop.history['pLvl'][t-1]) # (10000,)
# Call the saving rate function defined above
savRte = savRteFunc(LifeCyclePop, LifeCyclePop.history['mNrm'][t] , t)
savRte_list.append(savRte) # Add this period's saving rate to the list
# Create elements of matrix list
matrix_list = [0 for number in range(7)]
matrix_list[0] = t
matrix_list[1] = LifeCyclePop.history['aNrm'][t]
matrix_list[2] = LifeCyclePop.history['cNrm'][t]
matrix_list[3] = LifeCyclePop.history['TranShk'][t]
matrix_list[4] = LifeCyclePop.history['TranShk'][t-1]
matrix_list[5] = aLvlGroNow
matrix_list[6] = savRte
giant_list[t-1] = matrix_list
# Construct the level of assets A from a*p where a is the ratio to permanent income p
# Remember 41 is "years after entering workforce" (=age 25); 66 is the year right after retirement
LifeCyclePop.history['aLvl'] = LifeCyclePop.history['aNrm']*LifeCyclePop.history['pLvl']
aGro41=LifeCyclePop.history['aLvl'][41]/LifeCyclePop.history['aLvl'][40]
aGro41NoU=aGro41[aGro41[:]>0.2] # Throw out extreme outliers; don't want growth rates relative to 0 income!
# Plot the (truncated) distribution of growth rates of wealth between age 65 and 66 (=25 + 41)
from matplotlib import pyplot as plt
n, bins, patches = plt.hist(aGro41NoU,50,density=True)
Perhaps more interesting than the distribution of asset growth rates over the life cycle is the distribution of the level of assets, or the ratio of assets to permanent income.
Construct a plot similar to the one above for the disributions of $\texttt{aNrm}$ and $\texttt{aLev}$ in the period just BEFORE retirement (44 periods from the start).
# PROBLEM: soln here
# (rename all-caps problem in line above to all-caps solution)
In this model, each consumer experiences a set of draws of permanent income shocks over their lifetime. Some will be lucky and draw a mostly positive series of shocks (and unlucky people experience negative shocks).
This problem asks you to examine the consequences of these shocks for the lifetime pattern of saving.
The first step is to recalibrate the model so that there is no difference in initial assets, then reconstruct the initial conditions and simulate the model:
# PROBLEM: soln here
# (rename all-caps problem in line above to all-caps solution)
Now we are interested in comparing the people who were "lucky" vs those who were "unlucky"
The easiest way to measure this is by the cumulated level of noncapital (labor) income they have experienced over their working life.
For consumer in period 41 (age 66), calculate this object, then plot it against the $\texttt{aNrm}$ ratio at age 66.
# PROBLEM: soln here
# (rename all-caps problem in line above to all-caps solution)
You can have luck in transitory income shocks or in permanent income shocks. Their consequences are quite different. With a permanent shock, you expect your (noncapital) income to change forever, and (according to Friedman (1957)) you should adjust your consumption nearly one-for-one. With a transitory shock, you expect your income to return to its "permanent" level so you do not consume. So if you get a positive transitory shock, you will mostly save it.
The existence of transitory shocks therefore means that people who have on average experienced positive transitory shocks over their lifetimes should have higher saving rates. That would bias the relationship between lifetime income and the $\texttt{aNrm}$ ratio upward.
To see how important this might be, redo the same exercise as before, but using the level of (noncapital) permanent income (rather than overall income including transitory and permanent) over the lifetime. Comment on the result
# PROBLEM: soln here
# (rename all-caps problem in line above to all-caps solution)
The Haig-Simons definition of "saving" is basically the amount by which your wealth changes from one period to the next. This definition includes the consequences of any capital gains (or losses) for your wealth.
In recent work, Faegering, Holm, Natvik, and Moll have proposed that instead households largely ignore the consequences of capital gains and losses. That is, their consumption is largely unchanged by asset price movements.
Specifically, they define "active saving" as the difference between income and consumption neglecting any contriubutions from "buy and hold" assets like houses or stocks. The "active saving rate" is the quantity of active saving divided by the level of income. They find that the "active saving rate" is remarkably stable over the range from roughly the 20th percentile to the 95th percentile of the wealth distribution (see the figures below from their paper).
The basic model considered above does not allow for capital gains or losses, so it can be used to calculate directly the saving behavior of people who do not anticipate capital gains and losses. So, the saving rate computed by the $\texttt{savRte}$ function above should correspond to their "active saving rate."
Your problem: For the entire population simulated above, calculate what the model predicts about the saving rate they measure. You will do this by grouping the population into vigntile bins, and calculating the average active saving rate for all the households in each vigntile, and then plotting the wealth vigntiles against their saving rates.
# PROBLEM: soln here
# (rename all-caps problem in line above to all-caps solution)
We are interested in how income growth over the lifetime of the agent affects their saving rate and asset ratio $a=A/P$.
cumulative_income_first_half = np.sum(LifeCyclePop.history['pLvl'][0:20,:]*LifeCyclePop.history['TranShk'][0:20,:],0)
cumulative_income_second_half = np.sum(LifeCyclePop.history['pLvl'][20:40,:]*LifeCyclePop.history['TranShk'][20:40,:],0)
lifetime_growth = cumulative_income_second_half/cumulative_income_first_half
t=39
vigntiles = pd.qcut(lifetime_growth,20,labels=False)
savRte = savRteFunc(LifeCyclePop, LifeCyclePop.history['mNrm'][t] , t)
savRtgueseByVigtile = np.zeros(20)
assetsByVigtile = np.zeros(20)
assetsNrmByVigtile = np.zeros(20)
savRteByVigtile = np.zeros(20)
for i in range(20):
savRteByVigtile[i] = np.mean(savRte[vigntiles==i])
assetsByVigtile[i] = np.mean(LifeCyclePop.history['aLvl'][t][vigntiles==i])
assetsNrmByVigtile[i] = np.mean(LifeCyclePop.history['aNrm'][t][vigntiles==i])
plt.plot(np.array(range(20)), savRteByVigtile)
plt.title("Saving Rate at age 65, by Vigntile of Lifetime Income Growth")
plt.xlabel("Vigntile of Lifetime Income Growth")
plt.ylabel("Savings Rate")
plt.figure()
plt.plot(np.array(range(20)), assetsByVigtile)
plt.title("Assets at age 65, by Vigntile of Lifetime Income Growth")
plt.xlabel("Vigntile of Lifetime Income Growth")
plt.ylabel("Assets")
plt.figure()
plt.plot(np.array(range(20)), assetsNrmByVigtile)
plt.title("Normalized Assets at age 65, by Vigntile of Lifetime Income Growth")
plt.xlabel("Vigntile of Lifetime Income Growth")
plt.ylabel("Normalized Assets")