From 337f4fa984e945111e3117433f78104add5456d9 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sat, 19 Dec 2020 16:14:26 -0500 Subject: [PATCH 01/47] Add calibrations from other parts of HARK --- HARK/Calibration/CGMRemark.py | 216 +++++++++++++++++ HARK/Calibration/Pandemic.py | 334 ++++++++++++++++++++++++++ HARK/Calibration/SolvingMicroDSOPs.py | 173 +++++++++++++ 3 files changed, 723 insertions(+) create mode 100644 HARK/Calibration/CGMRemark.py create mode 100644 HARK/Calibration/Pandemic.py create mode 100644 HARK/Calibration/SolvingMicroDSOPs.py diff --git a/HARK/Calibration/CGMRemark.py b/HARK/Calibration/CGMRemark.py new file mode 100644 index 000000000..224454907 --- /dev/null +++ b/HARK/Calibration/CGMRemark.py @@ -0,0 +1,216 @@ +import numpy as np +# %% Preferences + +# Relative risk aversion +CRRA = 10 +# Discount factor +DiscFac = 0.96 + +# Survival probabilities from the author's Fortran code +n = 80 +survprob = np.zeros(n+1) +survprob[1] = 0.99845 +survprob[2] = 0.99839 +survprob[3] = 0.99833 +survprob[4] = 0.9983 +survprob[5] = 0.99827 +survprob[6] = 0.99826 +survprob[7] = 0.99824 +survprob[8] = 0.9982 +survprob[9] = 0.99813 +survprob[10] = 0.99804 +survprob[11] = 0.99795 +survprob[12] = 0.99785 +survprob[13] = 0.99776 +survprob[14] = 0.99766 +survprob[15] = 0.99755 +survprob[16] = 0.99743 +survprob[17] = 0.9973 +survprob[18] = 0.99718 +survprob[19] = 0.99707 +survprob[20] = 0.99696 +survprob[21] = 0.99685 +survprob[22] = 0.99672 +survprob[23] = 0.99656 +survprob[24] = 0.99635 +survprob[25] = 0.9961 +survprob[26] = 0.99579 +survprob[27] = 0.99543 +survprob[28] = 0.99504 +survprob[29] = 0.99463 +survprob[30] = 0.9942 +survprob[31] = 0.9937 +survprob[32] = 0.99311 +survprob[33] = 0.99245 +survprob[34] = 0.99172 +survprob[35] = 0.99091 +survprob[36] = 0.99005 +survprob[37] = 0.98911 +survprob[38] = 0.98803 +survprob[39] = 0.9868 +survprob[40] = 0.98545 +survprob[41] = 0.98409 +survprob[42] = 0.9827 +survprob[43] = 0.98123 +survprob[44] = 0.97961 +survprob[45] = 0.97786 +survprob[46] = 0.97603 +survprob[47] = 0.97414 +survprob[48] = 0.97207 +survprob[49] = 0.9697 +survprob[50] = 0.96699 +survprob[51] = 0.96393 +survprob[52] = 0.96055 +survprob[53] = 0.9569 +survprob[54] = 0.9531 +survprob[55] = 0.94921 +survprob[56] = 0.94508 +survprob[57] = 0.94057 +survprob[58] = 0.9357 +survprob[59] = 0.93031 +survprob[60] = 0.92424 +survprob[61] = 0.91717 +survprob[62] = 0.90922 +survprob[63] = 0.90089 +survprob[64] = 0.89282 +survprob[65] = 0.88503 +survprob[66] = 0.87622 +survprob[67] = 0.86576 +survprob[68] = 0.8544 +survprob[69] = 0.8423 +survprob[70] = 0.82942 +survprob[71] = 0.8154 +survprob[72] = 0.80002 +survprob[73] = 0.78404 +survprob[74] = 0.76842 +survprob[75] = 0.75382 +survprob[76] = 0.73996 +survprob[77] = 0.72464 +survprob[78] = 0.71057 +survprob[79] = 0.6961 +survprob[80] = 0.6809 + +# Fix indexing problem (fortran starts at 1, python at 0) +survprob = np.delete(survprob, [0]) +# Now we have 80 probabilities of death, +# for ages 20 to 99. + +# Labor income + +# They assume its a polinomial of age. Here are the coefficients +a=-2.170042+2.700381 +b1=0.16818 +b2=-0.0323371/10 +b3=0.0019704/100 + +time_params = {'Age_born': 20, 'Age_retire': 65, 'Age_death': 100} +t_start = time_params['Age_born'] +t_ret = time_params['Age_retire'] # We are currently interpreting this as the last period of work +t_end = time_params['Age_death'] + +# They assume retirement income is a fraction of labor income in the +# last working period +repl_fac = 0.68212 + +# Compute average income at each point in (working) life +f = np.arange(t_start, t_ret+1,1) +f = a + b1*f + b2*(f**2) + b3*(f**3) +det_work_inc = np.exp(f) + +# Retirement income +det_ret_inc = repl_fac*det_work_inc[-1]*np.ones(t_end - t_ret) + +# Get a full vector of the deterministic part of income +det_income = np.concatenate((det_work_inc, det_ret_inc)) + +# ln Gamma_t+1 = ln f_t+1 - ln f_t +gr_fac = np.exp(np.diff(np.log(det_income))) + +# Now we have growth factors for T_end-1 periods. + +# Finally define the normalization factor used by CGM, for plots. +# ### IMPORTANT ### +# We adjust this normalization factor for what we believe is a typo in the +# original article. See the REMARK jupyter notebook for details. +norm_factor = det_income * np.exp(0) + +# %% Shocks + +# Transitory and permanent shock variance from the paper +std_tran_shock = np.sqrt(0.0738) +std_perm_shock = np.sqrt(0.0106) + +# Vectorize. (HARK turns off these shocks after T_retirement) +std_tran_vec = np.array([std_tran_shock]*(t_end-t_start)) +std_perm_vec = np.array([std_perm_shock]*(t_end-t_start)) + +# %% Financial instruments + +# Risk-free factor +Rfree = 1.02 + +# Creation of risky asset return distributions + +Mu = 0.06 # Equity premium +Std = 0.157 # standard deviation of rate-of-return shocks + +RiskyAvg = Mu + Rfree +RiskyStd = Std + +# Make a dictionary to specify the rest of params +dict_portfolio = { + # Usual params + 'CRRA': CRRA, + 'Rfree': Rfree, + 'DiscFac': DiscFac, + + # Life cycle + 'T_age' : t_end-t_start+1, # Time of death + 'T_cycle' : t_end-t_start, # Number of non-terminal periods + 'T_retire':t_ret-t_start+1, + 'LivPrb': survprob.tolist(), + 'PermGroFac': gr_fac.tolist(), + 'cycles': 1, + + # Income shocks + 'PermShkStd': std_perm_vec, + 'PermShkCount': 3, + 'TranShkStd': std_tran_vec, + 'TranShkCount': 3, + 'UnempPrb': 0, + 'UnempPrbRet': 0, + 'IncUnemp': 0, + 'IncUnempRet': 0, + 'BoroCnstArt': 0, + 'tax_rate':0.0, + + # Portfolio related params + 'RiskyAvg': RiskyAvg, + 'RiskyStd': RiskyStd, + 'RiskyCount': 3, + 'RiskyShareCount': 30, + + # Grid + 'aXtraMin': 0.001, + 'aXtraMax': 400, + 'aXtraCount': 400, + 'aXtraExtra': [None], + 'aXtraNestFac': 3, + + # General + 'vFuncBool': False, + 'CubicBool': False, + + # Simulation params + 'AgentCount': 10, + 'pLvlInitMean' : np.log(det_income[0]), # Mean of log initial permanent income (only matters for simulation) + 'pLvlInitStd' : std_perm_shock, # Standard deviation of log initial permanent income (only matters for simulation) + 'T_sim': (t_end - t_start+1)*50, + + # Unused params required for simulation + 'PermGroFacAgg': 1, + 'aNrmInitMean': -50.0, # Agents start with 0 assets (this is log-mean) + 'aNrmInitStd' : 0.0 +} + +age_plot_params = [20, 30, 55, 75] diff --git a/HARK/Calibration/Pandemic.py b/HARK/Calibration/Pandemic.py new file mode 100644 index 000000000..84c5f3696 --- /dev/null +++ b/HARK/Calibration/Pandemic.py @@ -0,0 +1,334 @@ +import numpy as np +import matplotlib.pyplot as plt +import os +import csv +from HARK.distribution import Uniform +from importlib import reload +figs_dir = '../../Figures/' + +# Import configurable parameters, and keep them updated through reload +import parameter_config +reload(parameter_config) +from parameter_config import * + +############################################################################### + +# Size of simulations +AgentCountTotal = 1000000 # Total simulated population +T_sim = 13 # Number of quarters to simulate in counterfactuals + +# Basic lifecycle length parameters (don't touch these) +init_age = 24 +working_T = 41*4 # Number of working periods +retired_T = 55*4 # Number of retired periods +T_cycle = working_T + retired_T + +# Define the distribution of the discount factor for each eduation level +DiscFacCount = 7 +DiscFacDstnD = Uniform(DiscFacMeanD-DiscFacSpread, DiscFacMeanD+DiscFacSpread).approx(DiscFacCount) +DiscFacDstnH = Uniform(DiscFacMeanH-DiscFacSpread, DiscFacMeanH+DiscFacSpread).approx(DiscFacCount) +DiscFacDstnC = Uniform(DiscFacMeanC-DiscFacSpread, DiscFacMeanC+DiscFacSpread).approx(DiscFacCount) +DiscFacDstns = [DiscFacDstnD, DiscFacDstnH, DiscFacDstnC] + +# Define permanent income growth rates for each education level (from Cagetti 2003) +PermGroRte_d_ann = [5.2522391e-002, 5.0039782e-002, 4.7586132e-002, 4.5162424e-002, 4.2769638e-002, 4.0408757e-002, 3.8080763e-002, 3.5786635e-002, 3.3527358e-002, 3.1303911e-002, 2.9117277e-002, 2.6968437e-002, 2.4858374e-002, 2.2788068e-002, 2.0758501e-002, 1.8770655e-002, 1.6825511e-002, 1.4924052e-002, 1.3067258e-002, 1.1256112e-002, 9.4915947e-003, 7.7746883e-003, 6.1063742e-003, 4.4876340e-003, 2.9194495e-003, 1.4028022e-003, -6.1326258e-005, -1.4719542e-003, -2.8280999e-003, -4.1287819e-003, -5.3730185e-003, -6.5598280e-003, -7.6882288e-003, -8.7572392e-003, -9.7658777e-003, -1.0713163e-002, -1.1598112e-002, -1.2419745e-002, -1.3177079e-002, -1.3869133e-002, -4.3985368e-001, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003] +PermGroRte_h_ann = [4.1102173e-002, 4.1194381e-002, 4.1117402e-002, 4.0878307e-002, 4.0484168e-002, 3.9942056e-002, 3.9259042e-002, 3.8442198e-002, 3.7498596e-002, 3.6435308e-002, 3.5259403e-002, 3.3977955e-002, 3.2598035e-002, 3.1126713e-002, 2.9571062e-002, 2.7938153e-002, 2.6235058e-002, 2.4468848e-002, 2.2646594e-002, 2.0775369e-002, 1.8862243e-002, 1.6914288e-002, 1.4938576e-002, 1.2942178e-002, 1.0932165e-002, 8.9156095e-003, 6.8995825e-003, 4.8911556e-003, 2.8974003e-003, 9.2538802e-004, -1.0178097e-003, -2.9251214e-003, -4.7894755e-003, -6.6038005e-003, -8.3610250e-003, -1.0054077e-002, -1.1675886e-002, -1.3219380e-002, -1.4677487e-002, -1.6043137e-002, -5.5864350e-001, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002] +PermGroRte_c_ann = [3.9375106e-002, 3.9030288e-002, 3.8601230e-002, 3.8091011e-002, 3.7502710e-002, 3.6839406e-002, 3.6104179e-002, 3.5300107e-002, 3.4430270e-002, 3.3497746e-002, 3.2505614e-002, 3.1456953e-002, 3.0354843e-002, 2.9202363e-002, 2.8002591e-002, 2.6758606e-002, 2.5473489e-002, 2.4150316e-002, 2.2792168e-002, 2.1402124e-002, 1.9983263e-002, 1.8538663e-002, 1.7071404e-002, 1.5584565e-002, 1.4081224e-002, 1.2564462e-002, 1.1037356e-002, 9.5029859e-003, 7.9644308e-003, 6.4247695e-003, 4.8870812e-003, 3.3544449e-003, 1.8299396e-003, 3.1664424e-004, -1.1823620e-003, -2.6640003e-003, -4.1251914e-003, -5.5628564e-003, -6.9739162e-003, -8.3552918e-003, -6.8938447e-001, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004] +PermGroRte_d_ann += 31*[PermGroRte_d_ann[-1]] # Add 31 years of the same permanent income growth rate to the end of the sequence +PermGroRte_h_ann += 31*[PermGroRte_h_ann[-1]] +PermGroRte_c_ann += 31*[PermGroRte_c_ann[-1]] +PermGroRte_d_retire = PermGroRte_d_ann[40] # Store the big shock to permanent income at retirement +PermGroRte_h_retire = PermGroRte_h_ann[40] +PermGroRte_c_retire = PermGroRte_c_ann[40] +PermGroRte_d_ann[40] = PermGroRte_d_ann[39] # Overwrite the "retirement drop" with the adjacent growth rate +PermGroRte_h_ann[40] = PermGroRte_h_ann[39] +PermGroRte_c_ann[40] = PermGroRte_c_ann[39] +PermGroFac_d = [] +PermGroFac_h = [] +PermGroFac_c = [] +for j in range(len(PermGroRte_d_ann)): # Make sequences of quarterly permanent income growth factors from annual permanent income growth rates + PermGroFac_d += 4*[(1 + PermGroRte_d_ann[j])**0.25] + PermGroFac_h += 4*[(1 + PermGroRte_h_ann[j])**0.25] + PermGroFac_c += 4*[(1 + PermGroRte_c_ann[j])**0.25] +PermGroFac_d[working_T-1] = 1 + PermGroRte_d_retire # Put the big shock at retirement back into the sequence +PermGroFac_h[working_T-1] = 1 + PermGroRte_h_retire +PermGroFac_c[working_T-1] = 1 + PermGroRte_c_retire +SkillRot = 0.00125 +for t in range(T_cycle): + PermGroFac_d[t] = PermGroFac_d[t]*np.ones(6) + PermGroFac_d[t][1:3] = PermGroFac_d[t][0] - SkillRot + PermGroFac_d[t][3:6] = PermGroFac_d[t][0] - SkillRot + PermGroFac_h[t] = PermGroFac_h[t]*np.ones(6) + PermGroFac_h[t][1:3] = PermGroFac_h[t][0] - SkillRot + PermGroFac_h[t][3:6] = PermGroFac_h[t][0] - SkillRot + PermGroFac_c[t] = PermGroFac_c[t]*np.ones(6) + PermGroFac_c[t][1:3] = PermGroFac_c[t][0] - SkillRot + PermGroFac_c[t][3:6] = PermGroFac_c[t][0] - SkillRot +PermGroFac_d_small = [PermGroFac_d[t][:2] for t in range(T_cycle)] +PermGroFac_h_small = [PermGroFac_h[t][:2] for t in range(T_cycle)] +PermGroFac_c_small = [PermGroFac_c[t][:2] for t in range(T_cycle)] + +# Define the paths of permanent and transitory shocks (from Sabelhaus and Song) +TranShkStd = (np.concatenate((np.linspace(0.1,0.12,17), 0.12*np.ones(17), np.linspace(0.12,0.075,61), np.linspace(0.074,0.007,68), np.zeros(retired_T+1)))*4)**0.5 +TranShkStd = np.ndarray.tolist(TranShkStd) +PermShkStd = np.concatenate((((0.00011342*(np.linspace(24,64.75,working_T-1)-47)**2 + 0.01)/(11.0/4.0))**0.5,np.zeros(retired_T+1))) +PermShkStd[123:162] = PermShkStd[122] # Don't extrapolate +PermShkStd = np.ndarray.tolist(PermShkStd) + +# Import survival probabilities from SSA data +data_location = os.path.dirname(os.path.abspath(__file__)) +f = open(data_location + '/' + 'USactuarial.txt','r') +actuarial_reader = csv.reader(f,delimiter='\t') +raw_actuarial = list(actuarial_reader) +base_death_probs = [] +for j in range(len(raw_actuarial)): + base_death_probs += [float(raw_actuarial[j][4])] # This effectively assumes that everyone is a white woman +f.close() + +# Import adjustments for education and apply them to the base mortality rates +f = open(data_location + '/' + 'EducMortAdj.txt','r') +adjustment_reader = csv.reader(f,delimiter=' ') +raw_adjustments = list(adjustment_reader) +d_death_probs = [] +h_death_probs = [] +c_death_probs = [] +for j in range(76): + d_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[j][1])] + h_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[j][2])] + c_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[j][3])] +for j in range(76,96): + d_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[75][1])] + h_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[75][2])] + c_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[75][3])] +LivPrb_d = [] +LivPrb_h = [] +LivPrb_c = [] +for j in range(len(d_death_probs)): # Convert annual mortality rates to quarterly survival rates + LivPrb_d += 4*[(1 - d_death_probs[j])**0.25] + LivPrb_h += 4*[(1 - h_death_probs[j])**0.25] + LivPrb_c += 4*[(1 - c_death_probs[j])**0.25] +LivPrb_d_small = [LivPrb_d[t]*np.ones(2) for t in range(T_cycle)] +LivPrb_h_small = [LivPrb_h[t]*np.ones(2) for t in range(T_cycle)] +LivPrb_c_small = [LivPrb_c[t]*np.ones(2) for t in range(T_cycle)] +LivPrb_d = [LivPrb_d[t]*np.ones(6) for t in range(T_cycle)] +LivPrb_h = [LivPrb_h[t]*np.ones(6) for t in range(T_cycle)] +LivPrb_c = [LivPrb_c[t]*np.ones(6) for t in range(T_cycle)] + + +def makeMrkvArray(Urate, Uspell, Dspell, Lspell, Dexit=0): + ''' + Make an age-varying list of Markov transition matrices given the unemployment + rate (in normal times), the average length of an unemployment spell, the average + length of a deep unemployment spell, and the average length of a lockdown. + Parameters + ---------- + Urate: float + Erogodic unemployment rate + Uspell: float + Expected length of unemployment spell + Dspell: float + Expected length of deep unemployment spell + Lspell: float + Expected length of pandemic + Dexit: float + How likely to leave deep unemployment when pandemic ends + 0 => No correlation between exiting deep unemployment and pandemic ending + 1 => Certain to exit deep unemployment when pandemic ends + ''' + U_persist = 1.-1./Uspell + E_persist = 1.-Urate*(1.-U_persist)/(1.-Urate) + D_persist = 1.-1./Dspell + u = U_persist + e = E_persist + d = D_persist + r = 1-1./Lspell + dr = r*(1-Dexit) + + MrkvArray_working = np.array([[e, 1-e, 0.0, 0.0, 0.0, 0.0], + [1-u, u, 0.0, 0.0, 0.0, 0.0], + [0.0, 1-d, d, 0.0, 0.0, 0.0], + [e*(1-r), (1-e)*(1-r), 0.0, e*r, (1-e)*r, 0.0], + [(1-u)*(1-r), u*(1-r), 0.0, (1-u)*r, u*r, 0.0], + [0.0, (1-d)*(1-dr), d*(1-dr), 0.0, (1-d)*dr, d*dr] + ]) + MrkvArray_retired = np.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0, 0.0], + [(1-r), 0.0, 0.0, r, 0.0, 0.0], + [(1-r), 0.0, 0.0, r, 0.0, 0.0], + [(1-r), 0.0, 0.0, r, 0.0, 0.0] + ]) + MrkvArray = (working_T-1)*[MrkvArray_working] + (retired_T+1)*[MrkvArray_retired] + return MrkvArray + + +# Make Markov transition arrays among discrete states in each period of the lifecycle (ACTUAL / SIMULATION) +MrkvArray_real = makeMrkvArray(Urate_normal, Uspell, Dspell_real, Lspell_real) + +# Make Markov transition arrays among discrete states in each period of the lifecycle (PERCEIVED / SIMULATION) +MrkvArray_pcvd = makeMrkvArray(Urate_normal, Uspell, Dspell_pcvd, Lspell_pcvd) + +# Make a two state Markov array ("small") that is only used when generating the initial distribution of states +U_logit_base = np.log(1./(1.-Urate_normal)-1.) +U_persist = 1.-1./Uspell +E_persist = 1.-Urate_normal*(1.-U_persist)/(1.-Urate_normal) +u = U_persist +e = E_persist +MrkvArray_working_small = np.array([[e, 1-e], + [1-u, u] + ]) +MrkvArray_retired_small = np.array([[1., 0.], + [1., 0.] + ]) +MrkvArray_small = (working_T-1)*[MrkvArray_working_small] + (retired_T+1)*[MrkvArray_retired_small] + +# Define a parameter dictionaries for three education levels +init_dropout = {"cycles" : 1, + "T_cycle": T_cycle, + "T_retire": working_T-1, + 'T_sim': T_cycle, + 'T_age': T_cycle+1, + 'AgentCount': 10000, + "PermGroFacAgg": PermGroFacAgg, + "PopGroFac": PopGroFac, + "CRRA": CRRA, + "DiscFac": 0.98, # This will be overwritten at type construction + "Rfree_big" : np.array(6*[1.01]), + "PermGroFac_big": PermGroFac_d, + "LivPrb_big": LivPrb_d, + "uPfac_big" : np.array(3*[1.0] + 3*[uPfac_L]), + "MrkvArray_big" : MrkvArray_pcvd, + "Rfree" : np.array(2*[1.01]), + "PermGroFac": PermGroFac_d_small, + "LivPrb": LivPrb_d_small, + "uPfac" : np.array(2*[1.0]), + "MrkvArray" : MrkvArray_small, # Yes, this is intentional + "MrkvArray_pcvd" : MrkvArray_small, # Yes, this is intentional + "MrkvArray_sim" : MrkvArray_real, + "BoroCnstArt": 0.0, + "PermShkStd": PermShkStd, + "PermShkCount": PermShkCount, + "TranShkStd": TranShkStd, + "TranShkCount": TranShkCount, + "UnempPrb": 0.0, # Unemployment is modeled as a Markov state + "UnempPrbRet": UnempPrbRet, + "IncUnemp": IncUnemp, + "IncUnempRet": IncUnempRet, + "aXtraMin": aXtraMin, + "aXtraMax": aXtraMax, + "aXtraCount": aXtraCount, + "aXtraExtra": aXtraExtra, + "aXtraNestFac": aXtraNestFac, + "CubicBool": False, + "vFuncBool": False, + 'aNrmInitMean': np.log(0.00001), # Initial assets are zero + 'aNrmInitStd': 0.0, + 'pLvlInitMean': pLvlInitMeanD, + 'pLvlInitStd': pLvlInitStd, + "MrkvPrbsInit" : np.array([1-Urate_normal, Urate_normal] + 4*[0.0]), + 'Urate' : Urate_normal, + 'Uspell' : Uspell, + 'L_shared' : L_shared, + 'Lspell_pcvd' : Lspell_pcvd, + 'Lspell_real' : Lspell_real, + 'Dspell_pcvd' : Dspell_pcvd, + 'Dspell_real' : Dspell_real, + 'EducType': 0, + 'UpdatePrb': UpdatePrb, + 'track_vars' : [] + } +if L_shared: + init_dropout['T_lockdown'] = int(Lspell_real) + +adj_highschool = { + "PermGroFac" : PermGroFac_h_small, + "LivPrb" : LivPrb_h_small, + "PermGroFac_big" : PermGroFac_h, + "LivPrb_big" : LivPrb_h, + 'pLvlInitMean' : pLvlInitMeanH, + 'EducType' : 1} +init_highschool = init_dropout.copy() +init_highschool.update(adj_highschool) + +adj_college = { + "PermGroFac" : PermGroFac_c_small, + "LivPrb" : LivPrb_c_small, + "PermGroFac_big" : PermGroFac_c, + "LivPrb_big" : LivPrb_c, + 'pLvlInitMean' : pLvlInitMeanC, + 'EducType' : 2} +init_college = init_dropout.copy() +init_college.update(adj_college) + +# Define a dictionary to represent the baseline scenario +base_dict = { + 'PanShock' : False, + 'StimMax' : 0., + 'StimCut0' : None, + 'StimCut1' : None, + 'BonusUnemp' : 0., + 'BonusDeep' : 0., + 'T_ahead' : 0, + 'UnempD' : U_logit_base, + 'UnempH' : U_logit_base, + 'UnempC' : U_logit_base, + 'UnempP' : 0., + 'UnempA1' : 0., + 'UnempA2' : 0., + 'DeepD' : -np.inf, + 'DeepH' : -np.inf, + 'DeepC' : -np.inf, + 'DeepP' : 0., + 'DeepA1' : 0., + 'DeepA2' : 0., + 'Dspell_pcvd' : Dspell_pcvd, # These five parameters don't do anything in baseline scenario + 'Dspell_real' : Dspell_real, + 'Lspell_pcvd' : Lspell_pcvd, + 'Lspell_real' : Lspell_real, + 'L_shared' : L_shared + } + +# Define a dictionary to mutate baseline for the pandemic +pandemic_changes = { + 'PanShock' : True, + 'UnempD' : UnempD + Unemp0, + 'UnempH' : UnempH + Unemp0, + 'UnempC' : UnempC + Unemp0, + 'UnempP' : UnempP, + 'UnempA1' : UnempA1, + 'UnempA2' : UnempA2, + 'DeepD' : DeepD + Deep0, + 'DeepH' : DeepH + Deep0, + 'DeepC' : DeepC + Deep0, + 'DeepP' : DeepP, + 'DeepA1' : DeepA1, + 'DeepA2' : DeepA2, + } + +# Define a dictionary to mutate baseline for the fiscal stimulus +stimulus_changes = { + 'StimMax' : StimMax, + 'StimCut0' : StimCut0, + 'StimCut1' : StimCut1, + 'BonusUnemp' : BonusUnemp, + 'BonusDeep' : BonusDeep, + 'T_ahead' : T_ahead, + } + +# Define a dictionary to mutate baseline for a deep unemployment pandemic +deep_pandemic_changes = { + 'PanShock' : True, + 'UnempD' : UnempD + DeepPanAdj1, + 'UnempH' : UnempH + DeepPanAdj1, + 'UnempC' : UnempC + DeepPanAdj1, + 'UnempP' : UnempP + DeepPanAdj3, + 'UnempA1' : UnempA1, + 'UnempA2' : UnempA2, + 'DeepD' : DeepD + DeepPanAdj2, + 'DeepH' : DeepH + DeepPanAdj2, + 'DeepC' : DeepC + DeepPanAdj2, + 'DeepP' : DeepP, + 'DeepA1' : DeepA1, + 'DeepA2' : DeepA2, + } + diff --git a/HARK/Calibration/SolvingMicroDSOPs.py b/HARK/Calibration/SolvingMicroDSOPs.py new file mode 100644 index 000000000..c684145fb --- /dev/null +++ b/HARK/Calibration/SolvingMicroDSOPs.py @@ -0,0 +1,173 @@ +''' +Specifies the full set of calibrated values required to estimate the SolvingMicroDSOPs +model. The empirical data is stored in a separate csv file and is loaded in SetupSCFdata. +''' +from __future__ import print_function +import numpy as np +# --------------------------------------------------------------------------------- +# Debugging flags +# --------------------------------------------------------------------------------- + +show_PermGroFacAgg_error = False +# Error Notes: +# This sets a "quick fix" to the error, AttributeError: 'TempConsumerType' object has no attribute 'PermGroFacAgg' +# If you set this flag to "True" you will see the error. A more thorough fix is to +# fix the place where this error was introduced (Set to "True" and it will appear; +# this was almost certainly introduced when the code was extended to be used in the +# GE setting). An even more thorough solution, which moves beyond the scope of +# fixing this error, is adding unit tests to ID when changes to some code will +# break things elsewhere. +# Note: alternatively, decide that the "init_consumer_objects['PermGroFacAgg'] = 1.0" +# line below fixes it properly ('feature not a bug') and remove all this text. + +# --------------------------------------------------------------------------------- +# - 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 = np.array([0.17, 0.5, 0.83]) # Three point discrete distribution of initial w +initial_wealth_income_ratio_probs = np.array([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 + +# ----------------------------------------------------------------------------- +# -- Set up the dictionary "container" for making a basic lifecycle type ------ +# ----------------------------------------------------------------------------- + +# 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_rate':0.0, + 'vFuncBool':vFuncBool, + 'CubicBool':CubicBool + } + +if show_PermGroFacAgg_error: + pass # do nothing +else: + print("***NOTE: using a 'quick fix' for an attribute error. See 'Error Notes' in EstimationParameter.py for further discussion.***") + init_consumer_objects['PermGroFacAgg'] = 1.0 + + +if __name__ == '__main__': + print("Sorry, EstimationParameters doesn't actually do anything on its own.") + print("This module is imported by StructEstimation, providing calibrated ") + print("parameters for the example estimation. Please see that module if you ") + print("want more interesting output.") From 67ade11a76b2a15a82498b10c71d551b89a21761 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sat, 19 Dec 2020 16:14:43 -0500 Subject: [PATCH 02/47] Compute growth factors from poly --- HARK/Calibration/Calibration.py | 59 +++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 HARK/Calibration/Calibration.py diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py new file mode 100644 index 000000000..ffbfa3784 --- /dev/null +++ b/HARK/Calibration/Calibration.py @@ -0,0 +1,59 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat Dec 19 15:08:54 2020 + +@author: Mateo +""" + +import numpy as np + +def AgeLogPolyToGrowthRates(coefs, age_min, age_max): + """ + The deterministic component of permanent income is often expressed as a + log-polynomial of age. In HARK, this part of the income process is + expressed in a sequence of growth factors 'PermGroFac'. + + This function computes growth factors from the coefficients of a + log-polynomial specification + + Parameters + ---------- + coefs : numpy array or list of floats + Coefficients of the income log-polynomial, in ascending degree order + (starting with the constant). + age_min : int + Starting age at which the polynomial applies. + age_max : int + Final age at which the polynomial applies. + + Returns + ------- + GrowthFac : [float] + List of growth factors that replicate the polynomial. + + """ + # Figure out the degree of the polynomial + deg = len(coefs) - 1 + + # Create age matrices + age = np.arange(age_min, age_max + 1).reshape(age_max - age_min +1,1) + age_mat = np.hstack(list(map(lambda n: age**n, range(deg+1)))) + + # Fing the value of the polynomial + lnYDet = np.dot(age_mat, np.array(coefs)) + + # Compute growth factors + GrowthFac = np.exp(np.diff(lnYDet)) + # The last growth factor is nan: we do not know lnYDet(age_max+1) + GrowthFac = np.append(GrowthFac, np.nan) + + return GrowthFac.tolist() + + +a=-2.170042+2.700381 +b1=0.16818 +b2=-0.0323371/10 +b3=0.0019704/100 + +CoefsTest = np.array([a, b1, b2, b3]) +print(AgeLogPolyToGrowthRates(CoefsTest, 20, 65)) \ No newline at end of file From bf13247bd096f58dec086eaf4242e5e6adee0ed2 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 27 Dec 2020 12:20:14 -0500 Subject: [PATCH 03/47] Start adding CGM calibrations --- HARK/Calibration/Calibration.py | 48 ++++++++++++++++++++++++++------- 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index ffbfa3784..fa6e9e744 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -16,6 +16,10 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): This function computes growth factors from the coefficients of a log-polynomial specification + The form of the polynomial is assumed to be + alpha_0 + age/10 * alpha_1 + age^2/100 * alpha_2 + ... + (age/10)^n * alpha_n + Be sure to adjust the coefficients accordingly. + Parameters ---------- coefs : numpy array or list of floats @@ -30,30 +34,54 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): ------- GrowthFac : [float] List of growth factors that replicate the polynomial. - + + Y0 : float + Initial level of income """ # Figure out the degree of the polynomial deg = len(coefs) - 1 # Create age matrices - age = np.arange(age_min, age_max + 1).reshape(age_max - age_min +1,1) - age_mat = np.hstack(list(map(lambda n: age**n, range(deg+1)))) + age_10 = np.arange(age_min, age_max + 1).reshape(age_max - age_min +1,1)/10 + age_mat = np.hstack(list(map(lambda n: age_10**n, range(deg+1)))) # Fing the value of the polynomial lnYDet = np.dot(age_mat, np.array(coefs)) + # Find the starting level + Y0 = np.exp(lnYDet[0]) + # Compute growth factors GrowthFac = np.exp(np.diff(lnYDet)) # The last growth factor is nan: we do not know lnYDet(age_max+1) GrowthFac = np.append(GrowthFac, np.nan) - return GrowthFac.tolist() + return GrowthFac.tolist(), Y0 + +def findProfile(GroFacs, Y0): + + factors = np.array([Y0] + GroFacs[:-1]) + Y = np.cumprod(factors) + + return Y -a=-2.170042+2.700381 -b1=0.16818 -b2=-0.0323371/10 -b3=0.0019704/100 +CGM_income = { + 'NoHS' : {'Poly': [-2.1361, 0.1684/10, -0.0353/10, 0.0023/10], + 'ReplRate': 0.8898}, + + 'HS' : {'Poly': [-2.1700, 0.1682/10, -0.0323/10, 0.0020/10], + 'ReplRate': 0.6821}, + + 'College' : {'Poly': [-4.3148, 0.3194/10, -0.0577/10, 0.0033/10], + 'ReplRate': 0.9389} +} + +PermGroFac, Y0 = AgeLogPolyToGrowthRates(CGM_income['NoHS']['Poly'], 21,65) +Y = findProfile(PermGroFac, Y0) + +import matplotlib.pyplot as plt + +plt.plot(Y) -CoefsTest = np.array([a, b1, b2, b3]) -print(AgeLogPolyToGrowthRates(CoefsTest, 20, 65)) \ No newline at end of file +findProfile(PermGroFac, Y0) \ No newline at end of file From 7fcbbfe4577c8df18b6cc44acff477baea3b388f Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 27 Dec 2020 15:44:07 -0500 Subject: [PATCH 04/47] Fix CGM income calibration --- HARK/Calibration/Calibration.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index fa6e9e744..7cf74eae9 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -67,13 +67,13 @@ def findProfile(GroFacs, Y0): CGM_income = { - 'NoHS' : {'Poly': [-2.1361, 0.1684/10, -0.0353/10, 0.0023/10], + 'NoHS' : {'Poly': [-2.1361 + 2.6275, 0.1684*10, -0.0353*10, 0.0023*10], 'ReplRate': 0.8898}, - 'HS' : {'Poly': [-2.1700, 0.1682/10, -0.0323/10, 0.0020/10], + 'HS' : {'Poly': [-2.1700 + 2.7004, 0.1682*10, -0.0323*10, 0.0020*10], 'ReplRate': 0.6821}, - 'College' : {'Poly': [-4.3148, 0.3194/10, -0.0577/10, 0.0033/10], + 'College' : {'Poly': [-4.3148 + 2.3831, 0.3194*10, -0.0577*10, 0.0033*10], 'ReplRate': 0.9389} } From c893fe9e4eb69a2cbb8d5efc97e9f7a796fffec8 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 27 Dec 2020 16:05:50 -0500 Subject: [PATCH 05/47] Add retirement period and plot CGM --- HARK/Calibration/Calibration.py | 43 +++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 7cf74eae9..ecffce8b4 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -57,6 +57,25 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): GrowthFac = np.append(GrowthFac, np.nan) return GrowthFac.tolist(), Y0 + +def findPermGroFacs(age_min, age_max, age_ret, poly_coefs, repl_rate): + + # First find working age growth rates and starting income + WrkGroFacs, Y0 = AgeLogPolyToGrowthRates(poly, age_min, age_ret) + + # Replace the last item, which must be NaN, with the replacement rate + WrkGroFacs[-1] = repl_rate + + # Now create the retirement phase + n_ret_years = age_max - age_ret + RetGroFacs = [1.0] * (n_ret_years - 1) + [np.nan] + + # Concatenate + GroFacs = WrkGroFacs + RetGroFacs + + return GroFacs, Y0 + + def findProfile(GroFacs, Y0): @@ -77,11 +96,25 @@ def findProfile(GroFacs, Y0): 'ReplRate': 0.9389} } -PermGroFac, Y0 = AgeLogPolyToGrowthRates(CGM_income['NoHS']['Poly'], 21,65) -Y = findProfile(PermGroFac, Y0) - import matplotlib.pyplot as plt -plt.plot(Y) +age_min = 21 +age_max = 100 +age_ret = 65 + +ages = np.arange(age_min, age_max + 1) + +plt.figure() +for spec in CGM_income.items(): + + label = spec[0] + poly = spec[1]['Poly'] + repl = spec[1]['ReplRate'] + + PermGroFac, Y0 = findPermGroFacs(age_min, age_max, age_ret, + poly, repl) + + plt.plot(ages, findProfile(PermGroFac, Y0), label = label) -findProfile(PermGroFac, Y0) \ No newline at end of file +plt.legend() +plt.show() \ No newline at end of file From a677da473ecec5eedc77d854a21a4ffc6c05d3fb Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 28 Dec 2020 18:06:36 -0500 Subject: [PATCH 06/47] Improve parsing and add constant volatilities --- HARK/Calibration/Calibration.py | 103 +++++++++++++++++++++++--------- 1 file changed, 76 insertions(+), 27 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index ecffce8b4..67131fd43 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -58,24 +58,67 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): return GrowthFac.tolist(), Y0 -def findPermGroFacs(age_min, age_max, age_ret, poly_coefs, repl_rate): +def findPermGroFacs(age_min, age_max, age_ret, PolyCoefs, ReplRate): - # First find working age growth rates and starting income - WrkGroFacs, Y0 = AgeLogPolyToGrowthRates(poly, age_min, age_ret) - - # Replace the last item, which must be NaN, with the replacement rate - WrkGroFacs[-1] = repl_rate - - # Now create the retirement phase - n_ret_years = age_max - age_ret - RetGroFacs = [1.0] * (n_ret_years - 1) + [np.nan] - - # Concatenate - GroFacs = WrkGroFacs + RetGroFacs + if age_ret is None: + + GroFacs, Y0 = AgeLogPolyToGrowthRates(PolyCoefs, age_min, age_max) + + else: + + # First find working age growth rates and starting income + WrkGroFacs, Y0 = AgeLogPolyToGrowthRates(PolyCoefs, age_min, age_ret) + + # Replace the last item, which must be NaN, with the replacement rate + WrkGroFacs[-1] = ReplRate + + # Now create the retirement phase + n_ret_years = age_max - age_ret + RetGroFacs = [1.0] * (n_ret_years - 1) + [np.nan] + + # Concatenate + GroFacs = WrkGroFacs + RetGroFacs return GroFacs, Y0 - + +def ParseIncomeSpec(age_min, age_max, + age_ret = None, PolyCoefs = None, ReplRate = None, + PermShkStd = None, TranShkStd = None): + + N_periods = age_max - age_min + 1 + + if age_ret is not None: + N_work_periods = age_ret - age_min + 1 + N_ret_periods = age_max - age_ret + + # Growth factors + if PolyCoefs is not None: + + PermGroFac, P0 = findPermGroFacs(age_min, age_max, age_ret, + PolyCoefs, ReplRate) + + else: + pass + + # Volatilities + if isinstance(PermShkStd, float) and isinstance(TranShkStd, float): + + if age_ret is None: + + PermShkStd = [PermShkStd] * N_periods + TranShkStd = [TranShkStd] * N_periods + + else: + + PermShkStd = [PermShkStd] * N_work_periods + [0.0] * N_ret_periods + TranShkStd = [TranShkStd] * N_work_periods + [0.0] * N_ret_periods + + else: + pass + + return {'PermGroFac': PermGroFac, 'P0': P0, + 'PermShkStd': PermShkStd, 'TranShkStd': TranShkStd} def findProfile(GroFacs, Y0): @@ -86,21 +129,29 @@ def findProfile(GroFacs, Y0): CGM_income = { - 'NoHS' : {'Poly': [-2.1361 + 2.6275, 0.1684*10, -0.0353*10, 0.0023*10], - 'ReplRate': 0.8898}, - - 'HS' : {'Poly': [-2.1700 + 2.7004, 0.1682*10, -0.0323*10, 0.0020*10], - 'ReplRate': 0.6821}, + 'NoHS' : {'PolyCoefs': [-2.1361 + 2.6275, 0.1684*10, -0.0353*10, 0.0023*10], + 'age_ret': 65, + 'ReplRate': 0.8898, + 'PermShkStd': np.sqrt(0.0105), + 'TranShkStd': np.sqrt(0.1056)}, + + 'HS' : {'PolyCoefs': [-2.1700 + 2.7004, 0.1682*10, -0.0323*10, 0.0020*10], + 'age_ret': 65, + 'ReplRate': 0.6821, + 'PermShkStd': np.sqrt(0.0106), + 'TranShkStd': np.sqrt(0.0738)}, - 'College' : {'Poly': [-4.3148 + 2.3831, 0.3194*10, -0.0577*10, 0.0033*10], - 'ReplRate': 0.9389} + 'College' : {'PolyCoefs': [-4.3148 + 2.3831, 0.3194*10, -0.0577*10, 0.0033*10], + 'age_ret': 65, + 'ReplRate': 0.9389, + 'PermShkStd': np.sqrt(0.0169), + 'TranShkStd': np.sqrt(0.0584)} } import matplotlib.pyplot as plt age_min = 21 age_max = 100 -age_ret = 65 ages = np.arange(age_min, age_max + 1) @@ -108,13 +159,11 @@ def findProfile(GroFacs, Y0): for spec in CGM_income.items(): label = spec[0] - poly = spec[1]['Poly'] - repl = spec[1]['ReplRate'] - PermGroFac, Y0 = findPermGroFacs(age_min, age_max, age_ret, - poly, repl) + params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) + MeanY = findProfile(params['PermGroFac'], params['P0']) - plt.plot(ages, findProfile(PermGroFac, Y0), label = label) + plt.plot(ages, MeanY, label = label) plt.legend() plt.show() \ No newline at end of file From 27d955fb1e7dc9c996d05fd08cf44966c7237a78 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sat, 2 Jan 2021 11:41:23 -0500 Subject: [PATCH 07/47] Add base year and Cagetti calibration (base) --- HARK/Calibration/Calibration.py | 62 ++++++++++++++++++++++++++++++--- 1 file changed, 58 insertions(+), 4 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 67131fd43..ed2467842 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -84,7 +84,8 @@ def findPermGroFacs(age_min, age_max, age_ret, PolyCoefs, ReplRate): def ParseIncomeSpec(age_min, age_max, age_ret = None, PolyCoefs = None, ReplRate = None, - PermShkStd = None, TranShkStd = None): + PermShkStd = None, TranShkStd = None, + **unused): N_periods = age_max - age_min + 1 @@ -128,28 +129,59 @@ def findProfile(GroFacs, Y0): return Y +# Processes from Cocco, Gomes, Maenhout (2005). CGM_income = { 'NoHS' : {'PolyCoefs': [-2.1361 + 2.6275, 0.1684*10, -0.0353*10, 0.0023*10], 'age_ret': 65, 'ReplRate': 0.8898, 'PermShkStd': np.sqrt(0.0105), - 'TranShkStd': np.sqrt(0.1056)}, + 'TranShkStd': np.sqrt(0.1056), + 'BaseYear': 1992}, 'HS' : {'PolyCoefs': [-2.1700 + 2.7004, 0.1682*10, -0.0323*10, 0.0020*10], 'age_ret': 65, 'ReplRate': 0.6821, 'PermShkStd': np.sqrt(0.0106), - 'TranShkStd': np.sqrt(0.0738)}, + 'TranShkStd': np.sqrt(0.0738), + 'BaseYear': 1992}, 'College' : {'PolyCoefs': [-4.3148 + 2.3831, 0.3194*10, -0.0577*10, 0.0033*10], 'age_ret': 65, 'ReplRate': 0.9389, 'PermShkStd': np.sqrt(0.0169), - 'TranShkStd': np.sqrt(0.0584)} + 'TranShkStd': np.sqrt(0.0584), + 'BaseYear': 1992} +} + +# Processes from Cagetti (2003). +# - He uses volatilities from Carroll-Samwick (1997) +Cagetti_income = { + 'NoHS' : {'PolyCoefs': [1.2430,0.6941,0.0361,-0.0259,0.0018], + 'age_ret': 65, + 'ReplRate': 0.8344, + 'PermShkStd': np.sqrt(0.0214), # Take 9-12 from CS + 'TranShkStd': np.sqrt(0.0658), # Take 9-12 from CS + 'BaseYear': 1992}, + + 'HS' : {'PolyCoefs': [3.0551,-0.6925,0.4339,-0.0703,0.0035], + 'age_ret': 65, + 'ReplRate': 0.8033, + 'PermShkStd': np.sqrt(0.0277), # Take HS diploma from CS + 'TranShkStd': np.sqrt(0.0431), # Take HS diploma from CS + 'BaseYear': 1992}, + + 'College' : {'PolyCoefs': [2.1684,0.5230,-0.0002,-0.0057,0.0001], + 'age_ret': 65, + 'ReplRate': 0.7215, + 'PermShkStd': np.sqrt(0.0146), # Take College degree from CS + 'TranShkStd': np.sqrt(0.0385), # Take College degree from CS + 'BaseYear': 1992} } import matplotlib.pyplot as plt +# %% CGM calibration + age_min = 21 age_max = 100 @@ -165,5 +197,27 @@ def findProfile(GroFacs, Y0): plt.plot(ages, MeanY, label = label) +plt.title('CGM') +plt.legend() +plt.show() + +# %% Cagetti calibration + +age_min = 25 +age_max = 91 + +ages = np.arange(age_min, age_max + 1) + +plt.figure() +for spec in Cagetti_income.items(): + + label = spec[0] + + params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) + MeanY = findProfile(params['PermGroFac'], params['P0']) + + plt.plot(ages, MeanY, label = label) + +plt.title('Cagetti') plt.legend() plt.show() \ No newline at end of file From 9c007a9c71c5e2bd429810caef9cd137c2ca6d76 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sat, 2 Jan 2021 14:47:36 -0500 Subject: [PATCH 08/47] Post-retirement polynomial in Cagetti --- HARK/Calibration/Calibration.py | 36 +++++++++++++++++++++++++++------ 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index ed2467842..0b3e8d2a3 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -83,7 +83,9 @@ def findPermGroFacs(age_min, age_max, age_ret, PolyCoefs, ReplRate): def ParseIncomeSpec(age_min, age_max, - age_ret = None, PolyCoefs = None, ReplRate = None, + age_ret = None, + PolyCoefs = None, ReplRate = None, + PolyRetir = None, PermShkStd = None, TranShkStd = None, **unused): @@ -96,12 +98,31 @@ def ParseIncomeSpec(age_min, age_max, # Growth factors if PolyCoefs is not None: - PermGroFac, P0 = findPermGroFacs(age_min, age_max, age_ret, - PolyCoefs, ReplRate) + if PolyRetir is None: + + PermGroFac, P0 = findPermGroFacs(age_min, age_max, age_ret, + PolyCoefs, ReplRate) + + else: + + # Working period + PermGroWrk, P0 = findPermGroFacs(age_min, age_ret, None, + PolyCoefs, ReplRate) + PLast = findProfile(PermGroWrk, P0)[-1] + + # Retirement period + PermGroRet, R0 = findPermGroFacs(age_ret+1, age_max, None, + PolyRetir, ReplRate) + + # Input the replacement rate into the Work grow factors + PermGroWrk[-1] = R0/PLast + PermGroFac = PermGroWrk + PermGroRet + else: pass + # Volatilities if isinstance(PermShkStd, float) and isinstance(TranShkStd, float): @@ -156,21 +177,24 @@ def findProfile(GroFacs, Y0): # Processes from Cagetti (2003). # - He uses volatilities from Carroll-Samwick (1997) Cagetti_income = { - 'NoHS' : {'PolyCoefs': [1.2430,0.6941,0.0361,-0.0259,0.0018], + 'NoHS' : {'PolyCoefs': [1.2430, 0.6941, 0.0361, -0.0259, 0.0018], + 'PolyRetir': [3.6872, -0.1034], 'age_ret': 65, 'ReplRate': 0.8344, 'PermShkStd': np.sqrt(0.0214), # Take 9-12 from CS 'TranShkStd': np.sqrt(0.0658), # Take 9-12 from CS 'BaseYear': 1992}, - 'HS' : {'PolyCoefs': [3.0551,-0.6925,0.4339,-0.0703,0.0035], + 'HS' : {'PolyCoefs': [3.0551, -0.6925, 0.4339, -0.0703, 0.0035], + 'PolyRetir': [4.1631, -0.1378], 'age_ret': 65, 'ReplRate': 0.8033, 'PermShkStd': np.sqrt(0.0277), # Take HS diploma from CS 'TranShkStd': np.sqrt(0.0431), # Take HS diploma from CS 'BaseYear': 1992}, - 'College' : {'PolyCoefs': [2.1684,0.5230,-0.0002,-0.0057,0.0001], + 'College' : {'PolyCoefs': [2.1684, 0.5230, -0.0002, -0.0057, 0.0001], + 'PolyRetir': [3.7636, -0.0369], 'age_ret': 65, 'ReplRate': 0.7215, 'PermShkStd': np.sqrt(0.0146), # Take College degree from CS From b014c5dff054eaaf346541099e63ed7c85d7ba12 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sat, 2 Jan 2021 16:13:02 -0500 Subject: [PATCH 09/47] First life-table parser --- HARK/Calibration/Calibration.py | 44 ++++++++++++++++++++++++++++++++- 1 file changed, 43 insertions(+), 1 deletion(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 0b3e8d2a3..4cbc14cc4 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -244,4 +244,46 @@ def findProfile(GroFacs, Y0): plt.title('Cagetti') plt.legend() -plt.show() \ No newline at end of file +plt.show() + +# %% Life probabilities + +import pandas as pd + +def parse_ssa_life_table(filename, sep, sex, min_age, max_age): + + lt = pd.read_csv(filename, sep = sep, header=[0,1,2]) + + # Death probability column depends on sex + if sex == 'female': + death_col = 4 + else: + death_col = 1 + + # Keep only age and death probability + lt = pd.DataFrame({'Age': lt.iloc[:,0], 'DProb': lt.iloc[:,death_col]}) + # And relevant years + lt = lt[lt['Age'] >= min_age] + lt = lt[lt['Age'] <= max_age].sort_values(by = ['Age']) + + # Compute survival probability + LivPrb = 1 - lt['DProb'].to_numpy() + # Make agents die with certainty in the last period + LivPrb[-1] = 0 + + return(list(LivPrb)) + +min_age = 21 +max_age = 100 +ages = np.arange(min_age, max_age + 1) + +plt.figure() +for s in ['male', 'female']: + + LivPrb = parse_ssa_life_table(filename = 'LifeTables/SSA_LifeTable2017.csv', + sep = ',', sex = s, + min_age = min_age, max_age = max_age) + + plt.plot(ages, LivPrb, label = s) + +plt.legend() \ No newline at end of file From e1edb45adc4b4753c97cce85695443b5aba8e140 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 3 Jan 2021 10:25:17 -0500 Subject: [PATCH 10/47] Rm unused replRate from Cagetti --- HARK/Calibration/Calibration.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 4cbc14cc4..f5f5efa6a 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -180,7 +180,6 @@ def findProfile(GroFacs, Y0): 'NoHS' : {'PolyCoefs': [1.2430, 0.6941, 0.0361, -0.0259, 0.0018], 'PolyRetir': [3.6872, -0.1034], 'age_ret': 65, - 'ReplRate': 0.8344, 'PermShkStd': np.sqrt(0.0214), # Take 9-12 from CS 'TranShkStd': np.sqrt(0.0658), # Take 9-12 from CS 'BaseYear': 1992}, @@ -188,7 +187,6 @@ def findProfile(GroFacs, Y0): 'HS' : {'PolyCoefs': [3.0551, -0.6925, 0.4339, -0.0703, 0.0035], 'PolyRetir': [4.1631, -0.1378], 'age_ret': 65, - 'ReplRate': 0.8033, 'PermShkStd': np.sqrt(0.0277), # Take HS diploma from CS 'TranShkStd': np.sqrt(0.0431), # Take HS diploma from CS 'BaseYear': 1992}, @@ -196,7 +194,6 @@ def findProfile(GroFacs, Y0): 'College' : {'PolyCoefs': [2.1684, 0.5230, -0.0002, -0.0057, 0.0001], 'PolyRetir': [3.7636, -0.0369], 'age_ret': 65, - 'ReplRate': 0.7215, 'PermShkStd': np.sqrt(0.0146), # Take College degree from CS 'TranShkStd': np.sqrt(0.0385), # Take College degree from CS 'BaseYear': 1992} From a8c430cb89475d8f4e06b294f0915627dde0d428 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 3 Jan 2021 10:49:42 -0500 Subject: [PATCH 11/47] Port tools from CGM --- HARK/Calibration/CalibTest.ipynb | 55 +++++++++++++++++++++ HARK/Calibration/CalibTest.py | 84 ++++++++++++++++++++++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 HARK/Calibration/CalibTest.ipynb create mode 100644 HARK/Calibration/CalibTest.py diff --git a/HARK/Calibration/CalibTest.ipynb b/HARK/Calibration/CalibTest.ipynb new file mode 100644 index 000000000..25af34153 --- /dev/null +++ b/HARK/Calibration/CalibTest.ipynb @@ -0,0 +1,55 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType" + ] + } + ], + "metadata": { + "jupytext": { + "formats": "ipynb,py:percent" + }, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.9" + }, + "latex_envs": { + "LaTeX_envs_menu_present": true, + "autoclose": false, + "autocomplete": true, + "bibliofile": "biblio.bib", + "cite_by": "apalike", + "current_citInitial": 1, + "eqLabelWithNumbers": true, + "eqNumInitial": 1, + "hotkeys": { + "equation": "Ctrl-E", + "itemize": "Ctrl-I" + }, + "labels_anchors": false, + "latex_user_defs": false, + "report_style_numbering": false, + "user_envs_cfg": false + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/HARK/Calibration/CalibTest.py b/HARK/Calibration/CalibTest.py new file mode 100644 index 000000000..b71516529 --- /dev/null +++ b/HARK/Calibration/CalibTest.py @@ -0,0 +1,84 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.2' +# jupytext_version: 1.2.4 +# kernelspec: +# display_name: Python 3 +# language: python +# name: python3 +# --- + +# %% +from HARK.ConsumptionSaving.ConsIndShockModel import ( +IndShockConsumerType, +init_lifecycle +) + +from HARK.Calibration.Calibration import ( + ParseIncomeSpec, + CGM_income, + Cagetti_income +) + +import matplotlib.pyplot as plt +import pandas as pd + +# %% Alter calibration +birth_age = 21 +death_age = 90 + +income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, + **CGM_income['NoHS']) + + +# %% Create and solve agent +Agent = IndShockConsumerType(**init_lifecycle) +Agent.solve() + +# %% Simulation + +# Setup + +# Number of agents and periods in the simulation. +Agent.AgentCount = 5 +Agent.T_sim = 8 + +# Set up the variables we want to keep track of. +Agent.track_vars = ['aNrmNow','cNrmNow', 'pLvlNow', 't_age','mNrmNow'] + +# Run the simulations +Agent.initializeSim() +Agent.simulate() + +# %% Extract and format simulation results + +raw_data = {'Age': Agent.history['t_age'].flatten()+birth_age - 1, + 'pIncome': Agent.history['pLvlNow'].flatten(), + 'nrmM': Agent.history['mNrmNow'].flatten(), + 'nrmC': Agent.history['cNrmNow'].flatten()} + +Data = pd.DataFrame(raw_data) +Data['Cons'] = Data.nrmC * Data.pIncome +Data['M'] = Data.nrmM * Data.pIncome + +# %% Plots + +# Find the mean of each variable at every age +AgeMeans = Data.groupby(['Age']).mean().reset_index() + +plt.figure() +plt.plot(AgeMeans.Age, AgeMeans.pIncome, + label = 'Permanent Income') +plt.plot(AgeMeans.Age, AgeMeans.M, + label = 'Market resources') +plt.plot(AgeMeans.Age, AgeMeans.Cons, + label = 'Consumption') +plt.legend() +plt.xlabel('Age') +plt.title('Variable Means Conditional on Survival') +plt.grid() From 425744d6ade7d02ced0843f51275afa398534c88 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 3 Jan 2021 11:04:41 -0500 Subject: [PATCH 12/47] Rearrange tools and examples --- HARK/Calibration/Calibration.py | 74 +++++----------------------- HARK/Calibration/__init__.py | 0 examples/Calibration/Calib_plots.py | 75 +++++++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 63 deletions(-) create mode 100644 HARK/Calibration/__init__.py create mode 100644 examples/Calibration/Calib_plots.py diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index f5f5efa6a..ff09e4fbc 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -6,6 +6,15 @@ """ import numpy as np +import pandas as pd + +__all__ = [ + "Cagetti_income", + "CGM_income", + "ParseIncomeSpec", + "findProfile", + "parse_ssa_life_table" +] def AgeLogPolyToGrowthRates(coefs, age_min, age_max): """ @@ -199,53 +208,7 @@ def findProfile(GroFacs, Y0): 'BaseYear': 1992} } -import matplotlib.pyplot as plt - -# %% CGM calibration - -age_min = 21 -age_max = 100 - -ages = np.arange(age_min, age_max + 1) - -plt.figure() -for spec in CGM_income.items(): - - label = spec[0] - - params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) - MeanY = findProfile(params['PermGroFac'], params['P0']) - - plt.plot(ages, MeanY, label = label) - -plt.title('CGM') -plt.legend() -plt.show() - -# %% Cagetti calibration - -age_min = 25 -age_max = 91 - -ages = np.arange(age_min, age_max + 1) - -plt.figure() -for spec in Cagetti_income.items(): - - label = spec[0] - - params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) - MeanY = findProfile(params['PermGroFac'], params['P0']) - - plt.plot(ages, MeanY, label = label) - -plt.title('Cagetti') -plt.legend() -plt.show() - -# %% Life probabilities - -import pandas as pd +# %% Tools for survival probabilities def parse_ssa_life_table(filename, sep, sex, min_age, max_age): @@ -268,19 +231,4 @@ def parse_ssa_life_table(filename, sep, sex, min_age, max_age): # Make agents die with certainty in the last period LivPrb[-1] = 0 - return(list(LivPrb)) - -min_age = 21 -max_age = 100 -ages = np.arange(min_age, max_age + 1) - -plt.figure() -for s in ['male', 'female']: - - LivPrb = parse_ssa_life_table(filename = 'LifeTables/SSA_LifeTable2017.csv', - sep = ',', sex = s, - min_age = min_age, max_age = max_age) - - plt.plot(ages, LivPrb, label = s) - -plt.legend() \ No newline at end of file + return(list(LivPrb)) \ No newline at end of file diff --git a/HARK/Calibration/__init__.py b/HARK/Calibration/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/examples/Calibration/Calib_plots.py b/examples/Calibration/Calib_plots.py new file mode 100644 index 000000000..39854a735 --- /dev/null +++ b/examples/Calibration/Calib_plots.py @@ -0,0 +1,75 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun Jan 3 10:50:02 2021 + +@author: Mateo +""" + +from HARK.Calibration.Calibration import ( + ParseIncomeSpec, + findProfile, + parse_ssa_life_table, + Cagetti_income, + CGM_income +) +import numpy as np +import matplotlib.pyplot as plt + +# %% CGM calibration + +age_min = 21 +age_max = 100 + +ages = np.arange(age_min, age_max + 1) + +plt.figure() +for spec in CGM_income.items(): + + label = spec[0] + + params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) + MeanY = findProfile(params['PermGroFac'], params['P0']) + + plt.plot(ages, MeanY, label = label) + +plt.title('CGM') +plt.legend() +plt.show() + +# %% Cagetti calibration + +age_min = 25 +age_max = 91 + +ages = np.arange(age_min, age_max + 1) + +plt.figure() +for spec in Cagetti_income.items(): + + label = spec[0] + + params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) + MeanY = findProfile(params['PermGroFac'], params['P0']) + + plt.plot(ages, MeanY, label = label) + +plt.title('Cagetti') +plt.legend() +plt.show() + +# %% Life probabilities + +min_age = 21 +max_age = 100 +ages = np.arange(min_age, max_age + 1) + +plt.figure() +for s in ['male', 'female']: + + LivPrb = parse_ssa_life_table(filename = 'LifeTables/SSA_LifeTable2017.csv', + sep = ',', sex = s, + min_age = min_age, max_age = max_age) + + plt.plot(ages, LivPrb, label = s) + +plt.legend() \ No newline at end of file From 8702cb08ca5b8d54220e3ff53f2ff5fbe23e539e Mon Sep 17 00:00:00 2001 From: MateoVG Date: Sun, 3 Jan 2021 11:34:47 -0500 Subject: [PATCH 13/47] Start. Everything is nan --- .../Calibration/CalibTest.ipynb | 0 {HARK => examples}/Calibration/CalibTest.py | 22 ++++++++++++++++--- 2 files changed, 19 insertions(+), 3 deletions(-) rename {HARK => examples}/Calibration/CalibTest.ipynb (100%) rename {HARK => examples}/Calibration/CalibTest.py (75%) diff --git a/HARK/Calibration/CalibTest.ipynb b/examples/Calibration/CalibTest.ipynb similarity index 100% rename from HARK/Calibration/CalibTest.ipynb rename to examples/Calibration/CalibTest.ipynb diff --git a/HARK/Calibration/CalibTest.py b/examples/Calibration/CalibTest.py similarity index 75% rename from HARK/Calibration/CalibTest.py rename to examples/Calibration/CalibTest.py index b71516529..f529a5daf 100644 --- a/HARK/Calibration/CalibTest.py +++ b/examples/Calibration/CalibTest.py @@ -21,12 +21,14 @@ from HARK.Calibration.Calibration import ( ParseIncomeSpec, + parse_ssa_life_table, CGM_income, Cagetti_income ) import matplotlib.pyplot as plt import pandas as pd +from copy import copy # %% Alter calibration birth_age = 21 @@ -35,9 +37,23 @@ income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, **CGM_income['NoHS']) +liv_prb = parse_ssa_life_table(filename = 'LifeTables/SSA_LifeTable2017.csv', + sep = ',', sex = 'male', + min_age = birth_age, max_age = death_age) + +time_params = { + 'T_age': None, # Age at which agents are automatically killed + 'T_cycle': death_age - birth_age + 1, + 'T_retire': 0 +} + +params = copy(init_lifecycle) +params.update(time_params) +params.update(income_params) +params.update({'LivPrb': liv_prb}) # %% Create and solve agent -Agent = IndShockConsumerType(**init_lifecycle) +Agent = IndShockConsumerType(**params) Agent.solve() # %% Simulation @@ -45,8 +61,8 @@ # Setup # Number of agents and periods in the simulation. -Agent.AgentCount = 5 -Agent.T_sim = 8 +Agent.AgentCount = 50 +Agent.T_sim = 200 # Set up the variables we want to keep track of. Agent.track_vars = ['aNrmNow','cNrmNow', 'pLvlNow', 't_age','mNrmNow'] From e9353e496d15f24e80e1f05d73ed90eeb05a531e Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 4 Jan 2021 11:06:03 -0500 Subject: [PATCH 14/47] Return non-terminal periods as hark expects --- HARK/Calibration/Calibration.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index ff09e4fbc..27204a979 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -41,7 +41,7 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): Returns ------- - GrowthFac : [float] + GrowthFac : [float] of length age_max - age_min + 1 List of growth factors that replicate the polynomial. Y0 : float @@ -98,11 +98,13 @@ def ParseIncomeSpec(age_min, age_max, PermShkStd = None, TranShkStd = None, **unused): - N_periods = age_max - age_min + 1 + # There is no income distribution for the last period, as the distributions + # are forward-looking and the last period is terminal. + N_periods = age_max - age_min if age_ret is not None: N_work_periods = age_ret - age_min + 1 - N_ret_periods = age_max - age_ret + N_ret_periods = age_max - age_ret - 1 # Growth factors if PolyCoefs is not None: @@ -127,6 +129,10 @@ def ParseIncomeSpec(age_min, age_max, PermGroWrk[-1] = R0/PLast PermGroFac = PermGroWrk + PermGroRet + # In any case, PermGroFac[-1] will be np.nan, signaling that there is + # no expected growth in the terminal period. Discard it, as HARK expect + # list of growth rates for non-terminal periods + PermGroFac = PermGroFac[:-1] else: pass @@ -224,11 +230,9 @@ def parse_ssa_life_table(filename, sep, sex, min_age, max_age): lt = pd.DataFrame({'Age': lt.iloc[:,0], 'DProb': lt.iloc[:,death_col]}) # And relevant years lt = lt[lt['Age'] >= min_age] - lt = lt[lt['Age'] <= max_age].sort_values(by = ['Age']) + lt = lt[lt['Age'] < max_age].sort_values(by = ['Age']) # Compute survival probability LivPrb = 1 - lt['DProb'].to_numpy() - # Make agents die with certainty in the last period - LivPrb[-1] = 0 return(list(LivPrb)) \ No newline at end of file From 44c2ecd02e5f21ada77fcf54acae480787cdd58a Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 4 Jan 2021 11:06:22 -0500 Subject: [PATCH 15/47] Add initial level of income mean --- HARK/Calibration/Calibration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 27204a979..b80629668 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -154,7 +154,7 @@ def ParseIncomeSpec(age_min, age_max, else: pass - return {'PermGroFac': PermGroFac, 'P0': P0, + return {'PermGroFac': PermGroFac, 'pLvlInitMean': np.log(P0), 'PermShkStd': PermShkStd, 'TranShkStd': TranShkStd} def findProfile(GroFacs, Y0): From 6d521813913ee0265bb42b5ddf709a41a55fef20 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 4 Jan 2021 11:08:12 -0500 Subject: [PATCH 16/47] Fix T_age and T_cycle --- examples/Calibration/CalibTest.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/Calibration/CalibTest.py b/examples/Calibration/CalibTest.py index f529a5daf..435352054 100644 --- a/examples/Calibration/CalibTest.py +++ b/examples/Calibration/CalibTest.py @@ -42,8 +42,8 @@ min_age = birth_age, max_age = death_age) time_params = { - 'T_age': None, # Age at which agents are automatically killed - 'T_cycle': death_age - birth_age + 1, + 'T_age': death_age - birth_age, # Age at which agents are automatically killed + 'T_cycle': death_age - birth_age, 'T_retire': 0 } From 5fe659aa400fb46240775b25a7025d25536e2173 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 4 Jan 2021 11:37:25 -0500 Subject: [PATCH 17/47] More dimensionality fixes --- HARK/Calibration/Calibration.py | 6 +++--- examples/Calibration/Calib_plots.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index b80629668..36089c5c8 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -119,7 +119,7 @@ def ParseIncomeSpec(age_min, age_max, # Working period PermGroWrk, P0 = findPermGroFacs(age_min, age_ret, None, PolyCoefs, ReplRate) - PLast = findProfile(PermGroWrk, P0)[-1] + PLast = findProfile(PermGroWrk[:-1], P0)[-1] # Retirement period PermGroRet, R0 = findPermGroFacs(age_ret+1, age_max, None, @@ -154,12 +154,12 @@ def ParseIncomeSpec(age_min, age_max, else: pass - return {'PermGroFac': PermGroFac, 'pLvlInitMean': np.log(P0), + return {'PermGroFac': PermGroFac, 'P0': P0, 'pLvlInitMean': np.log(P0), 'PermShkStd': PermShkStd, 'TranShkStd': TranShkStd} def findProfile(GroFacs, Y0): - factors = np.array([Y0] + GroFacs[:-1]) + factors = np.array([Y0] + GroFacs) Y = np.cumprod(factors) return Y diff --git a/examples/Calibration/Calib_plots.py b/examples/Calibration/Calib_plots.py index 39854a735..618f51758 100644 --- a/examples/Calibration/Calib_plots.py +++ b/examples/Calibration/Calib_plots.py @@ -61,7 +61,7 @@ min_age = 21 max_age = 100 -ages = np.arange(min_age, max_age + 1) +ages = np.arange(min_age, max_age) plt.figure() for s in ['male', 'female']: From 2fddaca383f991b4b72347c2d06d2ff3e4b8e57e Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 4 Jan 2021 16:02:14 -0500 Subject: [PATCH 18/47] Add parser for time-related hark parameters --- HARK/Calibration/Calibration.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 36089c5c8..fe47fb25f 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -235,4 +235,16 @@ def parse_ssa_life_table(filename, sep, sex, min_age, max_age): # Compute survival probability LivPrb = 1 - lt['DProb'].to_numpy() - return(list(LivPrb)) \ No newline at end of file + return(list(LivPrb)) + +# %% Tools for setting time-related parameters + +def ParseTimeParams(age_birth, age_death): + + # T_cycle is the number of non-terminal periods in the agent's problem + T_cycle = age_death - age_birth + # T_age is the age at which the agents are killed with certainty in + # simulations (at the end of the T_age-th period) + T_age = age_death - age_birth + 1 + + return({'T_cycle': T_cycle, 'T_age': T_age}) \ No newline at end of file From fe620e56d9160ab012c0b00e897d14ec9dd1c4e6 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 4 Jan 2021 16:03:54 -0500 Subject: [PATCH 19/47] Use time parser --- examples/Calibration/CalibTest.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/examples/Calibration/CalibTest.py b/examples/Calibration/CalibTest.py index 435352054..507868f09 100644 --- a/examples/Calibration/CalibTest.py +++ b/examples/Calibration/CalibTest.py @@ -21,6 +21,7 @@ from HARK.Calibration.Calibration import ( ParseIncomeSpec, + ParseTimeParams, parse_ssa_life_table, CGM_income, Cagetti_income @@ -41,16 +42,12 @@ sep = ',', sex = 'male', min_age = birth_age, max_age = death_age) -time_params = { - 'T_age': death_age - birth_age, # Age at which agents are automatically killed - 'T_cycle': death_age - birth_age, - 'T_retire': 0 -} +time_params = ParseTimeParams(age_birth = birth_age, age_death = death_age) params = copy(init_lifecycle) params.update(time_params) params.update(income_params) -params.update({'LivPrb': liv_prb}) +params.update({'LivPrb': liv_prb, 'aNrmInitMean': -100}) # %% Create and solve agent Agent = IndShockConsumerType(**params) From 7509a9e439ba2ca650cb451e8860c0dc37aca61d Mon Sep 17 00:00:00 2001 From: MateoVG Date: Thu, 7 Jan 2021 09:58:38 -0500 Subject: [PATCH 20/47] Label plots --- examples/Calibration/CalibTest.py | 1 + examples/Calibration/Calib_plots.py | 9 ++++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/examples/Calibration/CalibTest.py b/examples/Calibration/CalibTest.py index 507868f09..87901488f 100644 --- a/examples/Calibration/CalibTest.py +++ b/examples/Calibration/CalibTest.py @@ -93,5 +93,6 @@ label = 'Consumption') plt.legend() plt.xlabel('Age') +plt.ylabel('Thousands of 1992 USD') plt.title('Variable Means Conditional on Survival') plt.grid() diff --git a/examples/Calibration/Calib_plots.py b/examples/Calibration/Calib_plots.py index 618f51758..4e33eb7d8 100644 --- a/examples/Calibration/Calib_plots.py +++ b/examples/Calibration/Calib_plots.py @@ -33,6 +33,8 @@ plt.plot(ages, MeanY, label = label) plt.title('CGM') +plt.xlabel('Age') +plt.ylabel('Mean permanent income (1992 USD)') plt.legend() plt.show() @@ -54,6 +56,9 @@ plt.plot(ages, MeanY, label = label) plt.title('Cagetti') +plt.xlabel('Age') + +plt.ylabel('Mean permanent income (1992 USD)') plt.legend() plt.show() @@ -71,5 +76,7 @@ min_age = min_age, max_age = max_age) plt.plot(ages, LivPrb, label = s) - + +plt.title('SSA Survival Probabilities') +plt.xlabel('Age') plt.legend() \ No newline at end of file From b04e22ea5705eaa0ebe0ebde08b40895ddddabd3 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Thu, 14 Jan 2021 14:48:34 -0500 Subject: [PATCH 21/47] Rename "calibplots" --- ...{Calib_plots.py => Income_calibrations.py} | 19 +------------------ 1 file changed, 1 insertion(+), 18 deletions(-) rename examples/Calibration/{Calib_plots.py => Income_calibrations.py} (71%) diff --git a/examples/Calibration/Calib_plots.py b/examples/Calibration/Income_calibrations.py similarity index 71% rename from examples/Calibration/Calib_plots.py rename to examples/Calibration/Income_calibrations.py index 618f51758..cf32f5710 100644 --- a/examples/Calibration/Calib_plots.py +++ b/examples/Calibration/Income_calibrations.py @@ -55,21 +55,4 @@ plt.title('Cagetti') plt.legend() -plt.show() - -# %% Life probabilities - -min_age = 21 -max_age = 100 -ages = np.arange(min_age, max_age) - -plt.figure() -for s in ['male', 'female']: - - LivPrb = parse_ssa_life_table(filename = 'LifeTables/SSA_LifeTable2017.csv', - sep = ',', sex = s, - min_age = min_age, max_age = max_age) - - plt.plot(ages, LivPrb, label = s) - -plt.legend() \ No newline at end of file +plt.show() \ No newline at end of file From ee23a6b95124d5c82b29c0d070903edeb647ad68 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Thu, 14 Jan 2021 14:59:41 -0500 Subject: [PATCH 22/47] Use new SSA module in example --- examples/Calibration/Income_calibrations.py | 1 - .../{CalibTest.py => Life_Cycle_example.py} | 12 +++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) rename examples/Calibration/{CalibTest.py => Life_Cycle_example.py} (87%) diff --git a/examples/Calibration/Income_calibrations.py b/examples/Calibration/Income_calibrations.py index 2bc2db37a..d9fc454f1 100644 --- a/examples/Calibration/Income_calibrations.py +++ b/examples/Calibration/Income_calibrations.py @@ -8,7 +8,6 @@ from HARK.Calibration.Calibration import ( ParseIncomeSpec, findProfile, - parse_ssa_life_table, Cagetti_income, CGM_income ) diff --git a/examples/Calibration/CalibTest.py b/examples/Calibration/Life_Cycle_example.py similarity index 87% rename from examples/Calibration/CalibTest.py rename to examples/Calibration/Life_Cycle_example.py index 87901488f..9bb1cf30d 100644 --- a/examples/Calibration/CalibTest.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -22,11 +22,12 @@ from HARK.Calibration.Calibration import ( ParseIncomeSpec, ParseTimeParams, - parse_ssa_life_table, CGM_income, Cagetti_income ) +from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table + import matplotlib.pyplot as plt import pandas as pd from copy import copy @@ -36,11 +37,12 @@ death_age = 90 income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, - **CGM_income['NoHS']) + **Cagetti_income['NoHS']) -liv_prb = parse_ssa_life_table(filename = 'LifeTables/SSA_LifeTable2017.csv', - sep = ',', sex = 'male', - min_age = birth_age, max_age = death_age) +# We need survival probabilities only up to death_age-1, because survival +# probability at death_age is 1. +liv_prb = parse_ssa_life_table(female = True, cross_sec = True, year = 2004, + min_age = birth_age, max_age = death_age - 1) time_params = ParseTimeParams(age_birth = birth_age, age_death = death_age) From 927afd43c291fd3a29b61671d71ae7c8cc7628a5 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Fri, 15 Jan 2021 10:17:43 -0500 Subject: [PATCH 23/47] Rm old ssa function --- HARK/Calibration/Calibration.py | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index fe47fb25f..0ffff40eb 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -214,29 +214,6 @@ def findProfile(GroFacs, Y0): 'BaseYear': 1992} } -# %% Tools for survival probabilities - -def parse_ssa_life_table(filename, sep, sex, min_age, max_age): - - lt = pd.read_csv(filename, sep = sep, header=[0,1,2]) - - # Death probability column depends on sex - if sex == 'female': - death_col = 4 - else: - death_col = 1 - - # Keep only age and death probability - lt = pd.DataFrame({'Age': lt.iloc[:,0], 'DProb': lt.iloc[:,death_col]}) - # And relevant years - lt = lt[lt['Age'] >= min_age] - lt = lt[lt['Age'] < max_age].sort_values(by = ['Age']) - - # Compute survival probability - LivPrb = 1 - lt['DProb'].to_numpy() - - return(list(LivPrb)) - # %% Tools for setting time-related parameters def ParseTimeParams(age_birth, age_death): From 665aa46b6dcad071a1fa89b2e50b47025ab77db2 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Fri, 15 Jan 2021 10:35:38 -0500 Subject: [PATCH 24/47] First test of sabelhaus song --- HARK/Calibration/Calibration.py | 42 +++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 0ffff40eb..c8262ab71 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -6,14 +6,13 @@ """ import numpy as np -import pandas as pd +from HARK.Calibration.Income.sabelhaus_song.SabelhausSongProfiles import sabelhaus_song_var_profile __all__ = [ "Cagetti_income", "CGM_income", "ParseIncomeSpec", - "findProfile", - "parse_ssa_life_table" + "findProfile" ] def AgeLogPolyToGrowthRates(coefs, age_min, age_max): @@ -96,6 +95,7 @@ def ParseIncomeSpec(age_min, age_max, PolyCoefs = None, ReplRate = None, PolyRetir = None, PermShkStd = None, TranShkStd = None, + SabelhausSong = False, **unused): # There is no income distribution for the last period, as the distributions @@ -139,20 +139,38 @@ def ParseIncomeSpec(age_min, age_max, # Volatilities - if isinstance(PermShkStd, float) and isinstance(TranShkStd, float): - + if SabelhausSong: + if age_ret is None: - PermShkStd = [PermShkStd] * N_periods - TranShkStd = [TranShkStd] * N_periods - + IncShkStds = sabelhaus_song_var_profile(cohort = 1950, + age_min = age_min, age_max= age_max) + PermShkStd = IncShkStds['PermShkStd'] + TranShkStd = IncShkStds['TranShkStd'] + else: - PermShkStd = [PermShkStd] * N_work_periods + [0.0] * N_ret_periods - TranShkStd = [TranShkStd] * N_work_periods + [0.0] * N_ret_periods - + IncShkStds = sabelhaus_song_var_profile(cohort = 1950, + age_min = age_min, age_max= age_ret) + PermShkStd = IncShkStds['PermShkStd'] + [0.0] * N_ret_periods + TranShkStd = IncShkStds['TranShkStd'] + [0.0] * N_ret_periods + else: - pass + + if isinstance(PermShkStd, float) and isinstance(TranShkStd, float): + + if age_ret is None: + + PermShkStd = [PermShkStd] * N_periods + TranShkStd = [TranShkStd] * N_periods + + else: + + PermShkStd = [PermShkStd] * N_work_periods + [0.0] * N_ret_periods + TranShkStd = [TranShkStd] * N_work_periods + [0.0] * N_ret_periods + + else: + pass return {'PermGroFac': PermGroFac, 'P0': P0, 'pLvlInitMean': np.log(P0), 'PermShkStd': PermShkStd, 'TranShkStd': TranShkStd} From d27a7eccc6d33afbad15a0ce72e31c601af5998a Mon Sep 17 00:00:00 2001 From: MateoVG Date: Fri, 15 Jan 2021 10:46:32 -0500 Subject: [PATCH 25/47] Add Sabelhaus Song test --- examples/Calibration/Life_Cycle_example.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index 9bb1cf30d..d21378b4b 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -37,7 +37,7 @@ death_age = 90 income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, - **Cagetti_income['NoHS']) + **Cagetti_income['College'], SabelhausSong=True) # We need survival probabilities only up to death_age-1, because survival # probability at death_age is 1. From d0f9ca009268d4447d83d4b5bff7e1cae4dd1cb1 Mon Sep 17 00:00:00 2001 From: MateoVG Date: Mon, 18 Jan 2021 15:21:31 -0500 Subject: [PATCH 26/47] incorporate initial distributions --- examples/Calibration/Life_Cycle_example.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index d21378b4b..4927e4d20 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -27,7 +27,7 @@ ) from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table - +from HARK.datasets.SCF.WealthIncomeDist.parser import income_wealth_dists_from_scf import matplotlib.pyplot as plt import pandas as pd from copy import copy @@ -35,9 +35,15 @@ # %% Alter calibration birth_age = 21 death_age = 90 +education = 'HS' +# Income specification income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, - **Cagetti_income['College'], SabelhausSong=True) + **CGM_income[education], SabelhausSong=True) + +# Initial distribution of wealth and permanent income +dist_params = income_wealth_dists_from_scf(age = birth_age, + education = education) # We need survival probabilities only up to death_age-1, because survival # probability at death_age is 1. @@ -49,7 +55,8 @@ params = copy(init_lifecycle) params.update(time_params) params.update(income_params) -params.update({'LivPrb': liv_prb, 'aNrmInitMean': -100}) +params.update(dist_params) +params.update({'LivPrb': liv_prb}) # %% Create and solve agent Agent = IndShockConsumerType(**params) @@ -60,7 +67,7 @@ # Setup # Number of agents and periods in the simulation. -Agent.AgentCount = 50 +Agent.AgentCount = 500 Agent.T_sim = 200 # Set up the variables we want to keep track of. From 6457cfc8feeb937ec4645dd6d8da5ae2104539a6 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 22 Jan 2021 17:52:44 -0500 Subject: [PATCH 27/47] Update examples with infl. adjustments --- examples/Calibration/Life_Cycle_example.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index 4927e4d20..0ce761d30 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -33,17 +33,17 @@ from copy import copy # %% Alter calibration -birth_age = 21 +birth_age = 25 death_age = 90 -education = 'HS' +education = 'College' # Income specification income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, **CGM_income[education], SabelhausSong=True) # Initial distribution of wealth and permanent income -dist_params = income_wealth_dists_from_scf(age = birth_age, - education = education) +dist_params = income_wealth_dists_from_scf(base_year=1992, age = birth_age, + education = education, wave = 1995) # We need survival probabilities only up to death_age-1, because survival # probability at death_age is 1. From 3ec00c3b7a35fa4b3fc6396052ee8ae2bbdc2918 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Mon, 25 Jan 2021 16:18:12 -0500 Subject: [PATCH 28/47] Adjust volatility timing --- HARK/Calibration/Calibration.py | 28 +++++++++++++-------- examples/Calibration/Income_calibrations.py | 2 +- 2 files changed, 18 insertions(+), 12 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index c8262ab71..02b06ce07 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -6,7 +6,7 @@ """ import numpy as np -from HARK.Calibration.Income.sabelhaus_song.SabelhausSongProfiles import sabelhaus_song_var_profile +from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile __all__ = [ "Cagetti_income", @@ -98,12 +98,13 @@ def ParseIncomeSpec(age_min, age_max, SabelhausSong = False, **unused): - # There is no income distribution for the last period, as the distributions - # are forward-looking and the last period is terminal. + # How many non-terminal periods are there. N_periods = age_max - age_min if age_ret is not None: + # How many non terminal periods are spent working N_work_periods = age_ret - age_min + 1 + # How many non terminal periods are spent in retirement N_ret_periods = age_max - age_ret - 1 # Growth factors @@ -136,24 +137,29 @@ def ParseIncomeSpec(age_min, age_max, else: pass - - + # Volatilities + # In this section, it is important to keep in mind that IncomeDstn[t] + # is the income distribution from period t to t+1, as perceived in period + # t. + # Therefore (assuming an annual model with agents entering at age 0), + # IncomeDstn[3] would contain the distribution of income shocks that occur + # at the start of age 4. if SabelhausSong: if age_ret is None: IncShkStds = sabelhaus_song_var_profile(cohort = 1950, - age_min = age_min, age_max= age_max) + age_min = age_min + 1, age_max= age_max) PermShkStd = IncShkStds['PermShkStd'] TranShkStd = IncShkStds['TranShkStd'] else: IncShkStds = sabelhaus_song_var_profile(cohort = 1950, - age_min = age_min, age_max= age_ret) - PermShkStd = IncShkStds['PermShkStd'] + [0.0] * N_ret_periods - TranShkStd = IncShkStds['TranShkStd'] + [0.0] * N_ret_periods + age_min = age_min + 1, age_max= age_ret) + PermShkStd = IncShkStds['PermShkStd'] + [0.0] * (N_ret_periods + 1) + TranShkStd = IncShkStds['TranShkStd'] + [0.0] * (N_ret_periods + 1) else: @@ -166,8 +172,8 @@ def ParseIncomeSpec(age_min, age_max, else: - PermShkStd = [PermShkStd] * N_work_periods + [0.0] * N_ret_periods - TranShkStd = [TranShkStd] * N_work_periods + [0.0] * N_ret_periods + PermShkStd = [PermShkStd] * (N_work_periods - 1) + [0.0] * (N_ret_periods + 1) + TranShkStd = [TranShkStd] * (N_work_periods - 1) + [0.0] * (N_ret_periods + 1) else: pass diff --git a/examples/Calibration/Income_calibrations.py b/examples/Calibration/Income_calibrations.py index cf32f5710..00b5cc443 100644 --- a/examples/Calibration/Income_calibrations.py +++ b/examples/Calibration/Income_calibrations.py @@ -8,10 +8,10 @@ from HARK.Calibration.Calibration import ( ParseIncomeSpec, findProfile, - parse_ssa_life_table, Cagetti_income, CGM_income ) + import numpy as np import matplotlib.pyplot as plt From 3eade4c49499821d47b78fc25673e3e4b288a783 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Mon, 25 Jan 2021 17:05:57 -0500 Subject: [PATCH 29/47] CPI adj --- examples/Calibration/Income_calibrations.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/Calibration/Income_calibrations.py b/examples/Calibration/Income_calibrations.py index 85ffaac98..fecdae719 100644 --- a/examples/Calibration/Income_calibrations.py +++ b/examples/Calibration/Income_calibrations.py @@ -16,7 +16,7 @@ import matplotlib.pyplot as plt # %% CGM calibration -target_year = 1992 +target_year = 2018 age_min = 21 age_max = 100 @@ -50,7 +50,8 @@ label = spec[0] - params = ParseIncomeSpec(age_min = age_min, age_max = age_max, **spec[1]) + params = ParseIncomeSpec(age_min = age_min, age_max = age_max, + TargetYear = target_year, **spec[1]) MeanY = findProfile(params['PermGroFac'], params['P0']) plt.plot(ages, MeanY, label = label) From 0ff6414b85361b89c95afbd71bef26d84001bd23 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Mon, 1 Feb 2021 20:56:52 -0500 Subject: [PATCH 30/47] Add exact Cagetti calibration and year trends --- HARK/Calibration/Calibration.py | 326 ++++++++++++-------- examples/Calibration/Income_calibrations.py | 51 +-- 2 files changed, 225 insertions(+), 152 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 97ee415a6..5f9b0c565 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -9,12 +9,8 @@ from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile from HARK.datasets.cpi.us.CPITools import cpi_deflator -__all__ = [ - "Cagetti_income", - "CGM_income", - "ParseIncomeSpec", - "findProfile" -] +__all__ = ["Cagetti_income", "CGM_income", "ParseIncomeSpec", "findProfile"] + def AgeLogPolyToGrowthRates(coefs, age_min, age_max): """ @@ -49,95 +45,126 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): """ # Figure out the degree of the polynomial deg = len(coefs) - 1 - + # Create age matrices - age_10 = np.arange(age_min, age_max + 1).reshape(age_max - age_min +1,1)/10 - age_mat = np.hstack(list(map(lambda n: age_10**n, range(deg+1)))) - + age_10 = np.arange(age_min, age_max + 1).reshape(age_max - age_min + 1, 1) / 10 + age_mat = np.hstack(list(map(lambda n: age_10 ** n, range(deg + 1)))) + # Fing the value of the polynomial lnYDet = np.dot(age_mat, np.array(coefs)) - + # Find the starting level Y0 = np.exp(lnYDet[0]) - + # Compute growth factors GrowthFac = np.exp(np.diff(lnYDet)) # The last growth factor is nan: we do not know lnYDet(age_max+1) GrowthFac = np.append(GrowthFac, np.nan) - + return GrowthFac.tolist(), Y0 -def findPermGroFacs(age_min, age_max, age_ret, PolyCoefs, ReplRate): - + +def findPermGroFacs(age_min, age_max, age_ret, AgePolyCoefs, ReplRate): + if age_ret is None: - GroFacs, Y0 = AgeLogPolyToGrowthRates(PolyCoefs, age_min, age_max) + GroFacs, Y0 = AgeLogPolyToGrowthRates(AgePolyCoefs, age_min, age_max) else: # First find working age growth rates and starting income - WrkGroFacs, Y0 = AgeLogPolyToGrowthRates(PolyCoefs, age_min, age_ret) - + WrkGroFacs, Y0 = AgeLogPolyToGrowthRates(AgePolyCoefs, age_min, age_ret) + # Replace the last item, which must be NaN, with the replacement rate WrkGroFacs[-1] = ReplRate - + # Now create the retirement phase n_ret_years = age_max - age_ret RetGroFacs = [1.0] * (n_ret_years - 1) + [np.nan] - + # Concatenate GroFacs = WrkGroFacs + RetGroFacs - + return GroFacs, Y0 - -def ParseIncomeSpec(BaseYear, age_min, age_max, - age_ret = None, - PolyCoefs = None, ReplRate = None, - PolyRetir = None, - PermShkStd = None, TranShkStd = None, - SabelhausSong = False, TargetYear = None, - **unused): - + +def ParseIncomeSpec( + base_monet_year, + age_min, + age_max, + age_ret=None, + AgePolyCoefs=None, + ReplRate=None, + AgePolyRetir=None, + YearTrend=None, + start_year=None, + PermShkStd=None, + TranShkStd=None, + SabelhausSong=False, + adjust_infl_to=None, + **unused +): + + income_params = {} # How many non-terminal periods are there. N_periods = age_max - age_min - + if age_ret is not None: # How many non terminal periods are spent working N_work_periods = age_ret - age_min + 1 # How many non terminal periods are spent in retirement - N_ret_periods = age_max - age_ret - 1 - + N_ret_periods = age_max - age_ret - 1 + # Growth factors - if PolyCoefs is not None: - - if PolyRetir is None: - - PermGroFac, P0 = findPermGroFacs(age_min, age_max, age_ret, - PolyCoefs, ReplRate) - + if AgePolyCoefs is not None: + + if AgePolyRetir is None: + + PermGroFac, P0 = findPermGroFacs( + age_min, age_max, age_ret, AgePolyCoefs, ReplRate + ) + else: - + # Working period - PermGroWrk, P0 = findPermGroFacs(age_min, age_ret, None, - PolyCoefs, ReplRate) + PermGroWrk, P0 = findPermGroFacs( + age_min, age_ret, None, AgePolyCoefs, ReplRate + ) PLast = findProfile(PermGroWrk[:-1], P0)[-1] - + # Retirement period - PermGroRet, R0 = findPermGroFacs(age_ret+1, age_max, None, - PolyRetir, ReplRate) - + PermGroRet, R0 = findPermGroFacs( + age_ret + 1, age_max, None, AgePolyRetir, ReplRate + ) + # Input the replacement rate into the Work grow factors - PermGroWrk[-1] = R0/PLast + PermGroWrk[-1] = R0 / PLast PermGroFac = PermGroWrk + PermGroRet - + # In any case, PermGroFac[-1] will be np.nan, signaling that there is # no expected growth in the terminal period. Discard it, as HARK expect # list of growth rates for non-terminal periods PermGroFac = PermGroFac[:-1] - + + # Apply the yearly trend if it is given + if YearTrend is not None: + + # Compute and apply the compounding yearly growth factor + YearGroFac = np.exp(YearTrend["Coef"]) + PermGroFac = [x * YearGroFac for x in PermGroFac] + + # Aggregate growth + income_params["PermGroFacAgg"] = YearGroFac + + # Adjust P0 with the year trend + P0 = P0 * np.power(YearGroFac, start_year - YearTrend["ZeroYear"]) + + income_params["PermGroFac"] = PermGroFac + else: - pass + + # Placeholder for future ways of storing income calibrations + raise NotImplementedError() # Volatilities # In this section, it is important to keep in mind that IncomeDstn[t] @@ -147,120 +174,157 @@ def ParseIncomeSpec(BaseYear, age_min, age_max, # IncomeDstn[3] would contain the distribution of income shocks that occur # at the start of age 4. if SabelhausSong: - + if age_ret is None: - - IncShkStds = sabelhaus_song_var_profile(cohort = 1950, - age_min = age_min + 1, age_max= age_max) - PermShkStd = IncShkStds['PermShkStd'] - TranShkStd = IncShkStds['TranShkStd'] - + + IncShkStds = sabelhaus_song_var_profile( + cohort=1950, age_min=age_min + 1, age_max=age_max + ) + PermShkStd = IncShkStds["PermShkStd"] + TranShkStd = IncShkStds["TranShkStd"] + else: - - IncShkStds = sabelhaus_song_var_profile(cohort = 1950, - age_min = age_min + 1, age_max= age_ret) - PermShkStd = IncShkStds['PermShkStd'] + [0.0] * (N_ret_periods + 1) - TranShkStd = IncShkStds['TranShkStd'] + [0.0] * (N_ret_periods + 1) - + + IncShkStds = sabelhaus_song_var_profile( + cohort=1950, age_min=age_min + 1, age_max=age_ret + ) + PermShkStd = IncShkStds["PermShkStd"] + [0.0] * (N_ret_periods + 1) + TranShkStd = IncShkStds["TranShkStd"] + [0.0] * (N_ret_periods + 1) + else: - + if isinstance(PermShkStd, float) and isinstance(TranShkStd, float): - + if age_ret is None: - + PermShkStd = [PermShkStd] * N_periods TranShkStd = [TranShkStd] * N_periods - + else: - - PermShkStd = [PermShkStd] * (N_work_periods - 1) + [0.0] * (N_ret_periods + 1) - TranShkStd = [TranShkStd] * (N_work_periods - 1) + [0.0] * (N_ret_periods + 1) - + + PermShkStd = [PermShkStd] * (N_work_periods - 1) + [0.0] * ( + N_ret_periods + 1 + ) + TranShkStd = [TranShkStd] * (N_work_periods - 1) + [0.0] * ( + N_ret_periods + 1 + ) + else: - pass - + + raise NotImplementedError() + + income_params["PermShkStd"] = PermShkStd + income_params["TranShkStd"] = TranShkStd + # Apply inflation adjustment if requested - if TargetYear is not None: - + if adjust_infl_to is not None: + # Deflate using the CPI september measurement, which is what the SCF # uses. - defl = cpi_deflator(from_year = BaseYear, to_year = TargetYear, - base_month='SEP')[0] - + defl = cpi_deflator( + from_year=base_monet_year, to_year=adjust_infl_to, base_month="SEP" + )[0] + else: - + defl = 1 - + P0 = P0 * defl - - return {'PermGroFac': PermGroFac, 'P0': P0, 'pLvlInitMean': np.log(P0), - 'PermShkStd': PermShkStd, 'TranShkStd': TranShkStd} - + income_params["P0"] = P0 + income_params["pLvlInitMean"] = np.log(P0) + + return income_params + + def findProfile(GroFacs, Y0): - + factors = np.array([Y0] + GroFacs) Y = np.cumprod(factors) - + return Y - + # Processes from Cocco, Gomes, Maenhout (2005). CGM_income = { - 'NoHS' : {'PolyCoefs': [-2.1361 + 2.6275, 0.1684*10, -0.0353*10, 0.0023*10], - 'age_ret': 65, - 'ReplRate': 0.8898, - 'PermShkStd': np.sqrt(0.0105), - 'TranShkStd': np.sqrt(0.1056), - 'BaseYear': 1992}, - - 'HS' : {'PolyCoefs': [-2.1700 + 2.7004, 0.1682*10, -0.0323*10, 0.0020*10], - 'age_ret': 65, - 'ReplRate': 0.6821, - 'PermShkStd': np.sqrt(0.0106), - 'TranShkStd': np.sqrt(0.0738), - 'BaseYear': 1992}, - - 'College' : {'PolyCoefs': [-4.3148 + 2.3831, 0.3194*10, -0.0577*10, 0.0033*10], - 'age_ret': 65, - 'ReplRate': 0.9389, - 'PermShkStd': np.sqrt(0.0169), - 'TranShkStd': np.sqrt(0.0584), - 'BaseYear': 1992} + "NoHS": { + "AgePolyCoefs": [-2.1361 + 2.6275, 0.1684 * 10, -0.0353 * 10, 0.0023 * 10], + "age_ret": 65, + "ReplRate": 0.8898, + "PermShkStd": np.sqrt(0.0105), + "TranShkStd": np.sqrt(0.1056), + "base_monet_year": 1992, + }, + "HS": { + "AgePolyCoefs": [-2.1700 + 2.7004, 0.1682 * 10, -0.0323 * 10, 0.0020 * 10], + "age_ret": 65, + "ReplRate": 0.6821, + "PermShkStd": np.sqrt(0.0106), + "TranShkStd": np.sqrt(0.0738), + "base_monet_year": 1992, + }, + "College": { + "AgePolyCoefs": [-4.3148 + 2.3831, 0.3194 * 10, -0.0577 * 10, 0.0033 * 10], + "age_ret": 65, + "ReplRate": 0.9389, + "PermShkStd": np.sqrt(0.0169), + "TranShkStd": np.sqrt(0.0584), + "base_monet_year": 1992, + }, } # Processes from Cagetti (2003). # - He uses volatilities from Carroll-Samwick (1997) Cagetti_income = { - 'NoHS' : {'PolyCoefs': [1.2430, 0.6941, 0.0361, -0.0259, 0.0018], - 'PolyRetir': [3.6872, -0.1034], - 'age_ret': 65, - 'PermShkStd': np.sqrt(0.0214), # Take 9-12 from CS - 'TranShkStd': np.sqrt(0.0658), # Take 9-12 from CS - 'BaseYear': 1992}, - - 'HS' : {'PolyCoefs': [3.0551, -0.6925, 0.4339, -0.0703, 0.0035], - 'PolyRetir': [4.1631, -0.1378], - 'age_ret': 65, - 'PermShkStd': np.sqrt(0.0277), # Take HS diploma from CS - 'TranShkStd': np.sqrt(0.0431), # Take HS diploma from CS - 'BaseYear': 1992}, - - 'College' : {'PolyCoefs': [2.1684, 0.5230, -0.0002, -0.0057, 0.0001], - 'PolyRetir': [3.7636, -0.0369], - 'age_ret': 65, - 'PermShkStd': np.sqrt(0.0146), # Take College degree from CS - 'TranShkStd': np.sqrt(0.0385), # Take College degree from CS - 'BaseYear': 1992} + "NoHS": { + "AgePolyCoefs": [7.99641616, 1.06559456, -0.14449728, 0.00048128, 0.0004096], + "AgePolyRetir": [10.84636791, -0.24562326], + "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, + "age_ret": 65, + "PermShkStd": np.sqrt(0.0214), # Take 9-12 from CS + "TranShkStd": np.sqrt(0.0658), # Take 9-12 from CS + "base_monet_year": 1992, + }, + "HS": { + "AgePolyCoefs": [ + 10.01333075, + -0.563234304, + 0.348710528, + -0.059442176, + 0.002947072, + ], + "AgePolyRetir": [11.21721558, -0.26820465], + "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, + "age_ret": 65, + "PermShkStd": np.sqrt(0.0277), # Take HS diploma from CS + "TranShkStd": np.sqrt(0.0431), # Take HS diploma from CS + "base_monet_year": 1992, + }, + "College": { + "AgePolyCoefs": [ + 9.916855488, + -0.057984416, + 0.146196992, + -0.027623424, + 0.001282048, + ], + "AgePolyRetir": [10.81011279, -0.16610233], + "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, + "age_ret": 65, + "PermShkStd": np.sqrt(0.0146), # Take College degree from CS + "TranShkStd": np.sqrt(0.0385), # Take College degree from CS + "base_monet_year": 1992, + }, } # %% Tools for setting time-related parameters + def ParseTimeParams(age_birth, age_death): - + # T_cycle is the number of non-terminal periods in the agent's problem T_cycle = age_death - age_birth # T_age is the age at which the agents are killed with certainty in # simulations (at the end of the T_age-th period) T_age = age_death - age_birth + 1 - - return({'T_cycle': T_cycle, 'T_age': T_age}) \ No newline at end of file + + return {"T_cycle": T_cycle, "T_age": T_age} diff --git a/examples/Calibration/Income_calibrations.py b/examples/Calibration/Income_calibrations.py index fecdae719..5fa82c9ac 100644 --- a/examples/Calibration/Income_calibrations.py +++ b/examples/Calibration/Income_calibrations.py @@ -9,14 +9,14 @@ ParseIncomeSpec, findProfile, Cagetti_income, - CGM_income + CGM_income, ) import numpy as np import matplotlib.pyplot as plt # %% CGM calibration -target_year = 2018 +adjust_infl_to = 1992 age_min = 21 age_max = 100 @@ -25,16 +25,17 @@ plt.figure() for spec in CGM_income.items(): - + label = spec[0] - - params = ParseIncomeSpec(age_min = age_min, age_max = age_max, - TargetYear = target_year, **spec[1]) - MeanY = findProfile(params['PermGroFac'], params['P0']) - - plt.plot(ages, MeanY, label = label) - -plt.title('CGM') + + params = ParseIncomeSpec( + age_min=age_min, age_max=age_max, adjust_infl_to=adjust_infl_to, **spec[1] + ) + MeanY = findProfile(params["PermGroFac"], params["P0"]) + + plt.plot(ages, MeanY, label=label) + +plt.title("CGM") plt.legend() plt.show() @@ -42,20 +43,28 @@ age_min = 25 age_max = 91 +# Cagetti has a year trend in his specification, so we have to say on what +# year agents enter the model. +start_year = 1980 ages = np.arange(age_min, age_max + 1) plt.figure() for spec in Cagetti_income.items(): - + label = spec[0] - - params = ParseIncomeSpec(age_min = age_min, age_max = age_max, - TargetYear = target_year, **spec[1]) - MeanY = findProfile(params['PermGroFac'], params['P0']) - - plt.plot(ages, MeanY, label = label) - -plt.title('Cagetti') + + params = ParseIncomeSpec( + age_min=age_min, + age_max=age_max, + adjust_infl_to=adjust_infl_to, + start_year=start_year, + **spec[1] + ) + MeanY = findProfile(params["PermGroFac"], params["P0"]) + + plt.plot(ages, MeanY, label=label) + +plt.title("Cagetti") plt.legend() -plt.show() \ No newline at end of file +plt.show() From c92167fe56b274eaacacdd55ab2a06840c6bc1df Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Tue, 2 Feb 2021 09:54:22 -0500 Subject: [PATCH 31/47] Fix example --- examples/Calibration/Life_Cycle_example.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index 0ce761d30..8433fdde1 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -27,7 +27,7 @@ ) from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table -from HARK.datasets.SCF.WealthIncomeDist.parser import income_wealth_dists_from_scf +from HARK.datasets.SCF.WealthIncomeDist.SCFDistTools import income_wealth_dists_from_scf import matplotlib.pyplot as plt import pandas as pd from copy import copy From 9e85a87a538b9328056e7785fb64425adc3cded0 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Thu, 4 Feb 2021 18:11:41 -0500 Subject: [PATCH 32/47] Comment minor functions --- HARK/Calibration/Calibration.py | 85 +++++++++++++++++++++++++++------ 1 file changed, 71 insertions(+), 14 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 5f9b0c565..31fbb9245 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -15,8 +15,8 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): """ The deterministic component of permanent income is often expressed as a - log-polynomial of age. In HARK, this part of the income process is - expressed in a sequence of growth factors 'PermGroFac'. + log-polynomial of age. In multiple HARK models, this part of the income + process is expressed in a sequence of growth factors 'PermGroFac'. This function computes growth factors from the coefficients of a log-polynomial specification @@ -40,8 +40,8 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): GrowthFac : [float] of length age_max - age_min + 1 List of growth factors that replicate the polynomial. - Y0 : float - Initial level of income + P0 : float + Initial level of income implied my the polynomial """ # Figure out the degree of the polynomial deg = len(coefs) - 1 @@ -54,20 +54,54 @@ def AgeLogPolyToGrowthRates(coefs, age_min, age_max): lnYDet = np.dot(age_mat, np.array(coefs)) # Find the starting level - Y0 = np.exp(lnYDet[0]) + P0 = np.exp(lnYDet[0]) # Compute growth factors GrowthFac = np.exp(np.diff(lnYDet)) # The last growth factor is nan: we do not know lnYDet(age_max+1) GrowthFac = np.append(GrowthFac, np.nan) - return GrowthFac.tolist(), Y0 + return GrowthFac.tolist(), P0 def findPermGroFacs(age_min, age_max, age_ret, AgePolyCoefs, ReplRate): + """ + Finds initial income and sequence of growth factors from a polynomial + specification of log-income, an optional retirement age and a replacement + rate. + + Retirement income will be Income_{age_ret} * ReplRate. - if age_ret is None: + Parameters + ---------- + age_min : int + Initial age at which to compute the income specification. + age_max : int + Maximum age up to which the income process must be specified. + age_ret : int + Age of retirement. Note that retirement happens after labor income is + received. For example, age_ret = 65 then the agent will receive labor + income up to age 65 and retirement benefits starting at age 66. + If age_ret is None, there will be no retirement. + AgePolyCoefs : numpy array or list of floats + Coefficients of the income log-polynomial, in ascending degree order + (starting with the constant). Income follows the specification: + ln(P)_age = \sum_{i=1}^{len(AgePolyCoefs)} (age/10)^i * AgePolyCoefs[i] + ReplRate : float + Replacement rate for retirement income. + Returns + ------- + GroFacs : list + List of income growth factors. + Y0 : float + Level of income at age_min + """ + + if age_ret is None: + + # If there is no retirement, the age polynomial applies for the whole + # lifetime GroFacs, Y0 = AgeLogPolyToGrowthRates(AgePolyCoefs, age_min, age_max) else: @@ -157,7 +191,8 @@ def ParseIncomeSpec( income_params["PermGroFacAgg"] = YearGroFac # Adjust P0 with the year trend - P0 = P0 * np.power(YearGroFac, start_year - YearTrend["ZeroYear"]) + if start_year is not None: + P0 = P0 * np.power(YearGroFac, start_year - YearTrend["ZeroYear"]) income_params["PermGroFac"] = PermGroFac @@ -274,10 +309,14 @@ def findProfile(GroFacs, Y0): # Processes from Cagetti (2003). # - He uses volatilities from Carroll-Samwick (1997) +# - He expresses income in dollars. It is more amicable to express it in +# thousands of dollars, also making it comparable to CGM. Thus, we substract +# ln(1e3) = 3*ln(10) from intercepts, which is equivalent to dividing income +# by a thousand. Cagetti_income = { "NoHS": { - "AgePolyCoefs": [7.99641616, 1.06559456, -0.14449728, 0.00048128, 0.0004096], - "AgePolyRetir": [10.84636791, -0.24562326], + "AgePolyCoefs": [7.99641616 - 3.0*np.log(10), 1.06559456, -0.14449728, 0.00048128, 0.0004096], + "AgePolyRetir": [10.84636791- 3.0*np.log(10), -0.24562326], "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, "age_ret": 65, "PermShkStd": np.sqrt(0.0214), # Take 9-12 from CS @@ -286,13 +325,13 @@ def findProfile(GroFacs, Y0): }, "HS": { "AgePolyCoefs": [ - 10.01333075, + 10.01333075- 3.0*np.log(10), -0.563234304, 0.348710528, -0.059442176, 0.002947072, ], - "AgePolyRetir": [11.21721558, -0.26820465], + "AgePolyRetir": [11.21721558- 3.0*np.log(10), -0.26820465], "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, "age_ret": 65, "PermShkStd": np.sqrt(0.0277), # Take HS diploma from CS @@ -301,13 +340,13 @@ def findProfile(GroFacs, Y0): }, "College": { "AgePolyCoefs": [ - 9.916855488, + 9.916855488- 3.0*np.log(10), -0.057984416, 0.146196992, -0.027623424, 0.001282048, ], - "AgePolyRetir": [10.81011279, -0.16610233], + "AgePolyRetir": [10.81011279- 3.0*np.log(10), -0.16610233], "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, "age_ret": 65, "PermShkStd": np.sqrt(0.0146), # Take College degree from CS @@ -320,7 +359,25 @@ def findProfile(GroFacs, Y0): def ParseTimeParams(age_birth, age_death): + """ + Converts simple statements of the age at which an agent is born and the + age at which he dies with certaintiy into the parameters that HARK needs + for figuring out the timing of the model. + Parameters + ---------- + age_birth : int + Age at which the agent enters the model, e.g., 21. + age_death : int + Age at which the agent dies with certainty, e.g., 100. + + Returns + ------- + dict + Dictionary with parameters "T_cycle" and "T_age" which HARK expects + and which map to the birth and death ages specified by the user. + + """ # T_cycle is the number of non-terminal periods in the agent's problem T_cycle = age_death - age_birth # T_age is the age at which the agents are killed with certainty in From d5c8561e45219caae75db0a0ef9ebf137147a1e4 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Thu, 4 Feb 2021 18:18:19 -0500 Subject: [PATCH 33/47] Lint --- HARK/Calibration/Calibration.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 31fbb9245..ad605c689 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -97,9 +97,9 @@ def findPermGroFacs(age_min, age_max, age_ret, AgePolyCoefs, ReplRate): Y0 : float Level of income at age_min """ - + if age_ret is None: - + # If there is no retirement, the age polynomial applies for the whole # lifetime GroFacs, Y0 = AgeLogPolyToGrowthRates(AgePolyCoefs, age_min, age_max) @@ -308,6 +308,8 @@ def findProfile(GroFacs, Y0): } # Processes from Cagetti (2003). +# - The author generously provided estimates from which the age polynomials +# and yearly trends were recovered. # - He uses volatilities from Carroll-Samwick (1997) # - He expresses income in dollars. It is more amicable to express it in # thousands of dollars, also making it comparable to CGM. Thus, we substract @@ -315,8 +317,14 @@ def findProfile(GroFacs, Y0): # by a thousand. Cagetti_income = { "NoHS": { - "AgePolyCoefs": [7.99641616 - 3.0*np.log(10), 1.06559456, -0.14449728, 0.00048128, 0.0004096], - "AgePolyRetir": [10.84636791- 3.0*np.log(10), -0.24562326], + "AgePolyCoefs": [ + 7.99641616 - 3.0 * np.log(10), + 1.06559456, + -0.14449728, + 0.00048128, + 0.0004096, + ], + "AgePolyRetir": [10.84636791 - 3.0 * np.log(10), -0.24562326], "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, "age_ret": 65, "PermShkStd": np.sqrt(0.0214), # Take 9-12 from CS @@ -325,13 +333,13 @@ def findProfile(GroFacs, Y0): }, "HS": { "AgePolyCoefs": [ - 10.01333075- 3.0*np.log(10), + 10.01333075 - 3.0 * np.log(10), -0.563234304, 0.348710528, -0.059442176, 0.002947072, ], - "AgePolyRetir": [11.21721558- 3.0*np.log(10), -0.26820465], + "AgePolyRetir": [11.21721558 - 3.0 * np.log(10), -0.26820465], "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, "age_ret": 65, "PermShkStd": np.sqrt(0.0277), # Take HS diploma from CS @@ -340,13 +348,13 @@ def findProfile(GroFacs, Y0): }, "College": { "AgePolyCoefs": [ - 9.916855488- 3.0*np.log(10), + 9.916855488 - 3.0 * np.log(10), -0.057984416, 0.146196992, -0.027623424, 0.001282048, ], - "AgePolyRetir": [10.81011279- 3.0*np.log(10), -0.16610233], + "AgePolyRetir": [10.81011279 - 3.0 * np.log(10), -0.16610233], "YearTrend": {"Coef": 0.016, "ZeroYear": 1980}, "age_ret": 65, "PermShkStd": np.sqrt(0.0146), # Take College degree from CS From ed1c466cf92dd2afa269f7647f166346bba0b7d8 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Thu, 4 Feb 2021 18:36:58 -0500 Subject: [PATCH 34/47] Small comment --- HARK/Calibration/Calibration.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index ad605c689..9ad67c1f4 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -272,7 +272,23 @@ def ParseIncomeSpec( def findProfile(GroFacs, Y0): + """ + Generates a sequence {Y_{t}}_{t=0}^N from an initial Y_0 and a sequence + of growth factors GroFac[n] = Y_{n+1}/Y_n + + Parameters + ---------- + GroFacs : list or numpy array + Growth factors in chronological order. + Y0 : float + initial value of the series. + Returns + ------- + Y : numpy array + Array with the values of the series. + + """ factors = np.array([Y0] + GroFacs) Y = np.cumprod(factors) From 1c1c050ba2b938e06dd49760ab0c2368e5cc9493 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Thu, 4 Feb 2021 18:50:16 -0500 Subject: [PATCH 35/47] Improve presentation in income example file --- examples/Calibration/Income_calibrations.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/examples/Calibration/Income_calibrations.py b/examples/Calibration/Income_calibrations.py index 5fa82c9ac..879d46618 100644 --- a/examples/Calibration/Income_calibrations.py +++ b/examples/Calibration/Income_calibrations.py @@ -15,9 +15,11 @@ import numpy as np import matplotlib.pyplot as plt -# %% CGM calibration +# What year to use as the base monetary year? +# (pick 1992 as it is used by both of the papers we are comparing) adjust_infl_to = 1992 +# %% Cocco, Gomes, Maenhout (2005) calibration age_min = 21 age_max = 100 @@ -35,11 +37,18 @@ plt.plot(ages, MeanY, label=label) -plt.title("CGM") +plt.title( + "Mean paths of permanent income calibrations in\n" + + "Cocco, Gomes & Maenhout (2005)" +) +plt.xlabel("Age") +plt.ylabel( + "Mean Permanent Income,\n" + "Thousands of {} U.S. dollars".format(adjust_infl_to) +) plt.legend() plt.show() -# %% Cagetti calibration +# %% Cagetti (2003) calibration age_min = 25 age_max = 91 @@ -65,6 +74,10 @@ plt.plot(ages, MeanY, label=label) -plt.title("Cagetti") +plt.title("Mean paths of permanent income calibrations in\n" + "Cagetti (2003)") +plt.xlabel("Age") +plt.ylabel( + "Mean Permanent Income,\n" + "Thousands of {} U.S. dollars".format(adjust_infl_to) +) plt.legend() plt.show() From f82115f37e30be396e9299febf94eafbdf0ebe65 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Thu, 4 Feb 2021 18:57:10 -0500 Subject: [PATCH 36/47] Minor changes to LC simulation example --- examples/Calibration/Life_Cycle_example.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index 8433fdde1..e7c3ab49e 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -35,11 +35,12 @@ # %% Alter calibration birth_age = 25 death_age = 90 +income_calib = CGM_income education = 'College' # Income specification income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, - **CGM_income[education], SabelhausSong=True) + **income_calib[education], SabelhausSong=True) # Initial distribution of wealth and permanent income dist_params = income_wealth_dists_from_scf(base_year=1992, age = birth_age, @@ -54,8 +55,8 @@ params = copy(init_lifecycle) params.update(time_params) -params.update(income_params) params.update(dist_params) +params.update(income_params) params.update({'LivPrb': liv_prb}) # %% Create and solve agent From 332e3ef33b0f0b4bdd996e8ee96be053b096495f Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 14:07:36 -0500 Subject: [PATCH 37/47] Add docstring to the main function --- HARK/Calibration/Calibration.py | 87 +++++++++++++++++++++++++++++++-- 1 file changed, 84 insertions(+), 3 deletions(-) diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Calibration.py index 9ad67c1f4..33fdf492b 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Calibration.py @@ -136,8 +136,79 @@ def ParseIncomeSpec( TranShkStd=None, SabelhausSong=False, adjust_infl_to=None, - **unused ): + """ + A function that produces income growth rates and income shock volatilities + + Parameters + ---------- + base_monet_year : int + Base monetary year in which the income process is specified. Answer to + "In what year's U.S. dollars was income expressed in the process that + will be parsed?". + age_min : int + Age at which agents enter the model. + age_max : int + Age at whih agents die with certainty. E.g., if age_max = 100, the + agent dies at the end of his 100th year of life. + age_ret : int, optional + Age of retirement. The default is None. + AgePolyCoefs : numpy array or list of floats + Coefficients of the income log-polynomial, in ascending degree order + (starting with the constant). Permanent income follows the specification: + ln(P)_age = \sum_{i=1}^{len(AgePolyCoefs)} (age/10)^i * AgePolyCoefs[i]. + The default is None. + ReplRate : float, optional + Replacement rate for retirement income. Retirement income will be + Income_{age_ret} * ReplRate. The default is None. + AgePolyRetir : numpy array or list of floats + Specifies a different age polynomial for income after retirement. It + follows the same convention as AgePolyCoefs. The default is None. + YearTrend : dict, optional + Dictionary with entries "Coef" (float) and "ZeroYear" (int). Allows + a time trend to be added to log-income. If provided, mean log-income at + age a and year t will be: + ln P = polynomial(a) + Coef * (t - ZeroYear) + The default is None. + start_year : int, optional + Year at which the agent enters the model. This is important only for + specifications with a time-trend for income profiles. + The default is None. + PermShkStd : float, optional + Standard deviation of log-permanent-income shocks, if it is constant. + The default is None. + TranShkStd : float, optional + Standard deviation of log-transitory-income shocks, if it is constant. + The default is None. + SabelhausSong : bool, optional + Indicates whether to use transitory and permanent income shock + volatilities from Sabelhaus & Song (2010) "The Great Moderation in + Micro Labor Earnings". The default is False. + adjust_infl_to : int, optional + Year at which nominal quantities should be expressed. Answers the + question "In what year's U.S. dollars should income be expressed". + The default is None. In such case, base_monet_year will be used. + + Returns + ------- + income_params : dict + Dictionary with entries: + - P0: initial level of permanent income. + - pLvlInitMean: mean of the distribution of log-permanent income. + np.log(P0) = pLvlInitMean + - PermGroFac : list of deterministic growth factors for permanent + income. + - PermShkStd: list of standard deviations of shocks to + log-permanent income. + - TranShkStd: list of standard deviations of transitory shocks + to income. + - PermGroFacAgg: if a yearly trend in income is provided, this will + be the aggregate level of growth in permanent incomes. + + This dictionary has the names and formats that various models in HARK + expect, so that it can be directly updated into other parameter + dictionaries. + """ income_params = {} # How many non-terminal periods are there. @@ -295,7 +366,15 @@ def findProfile(GroFacs, Y0): return Y -# Processes from Cocco, Gomes, Maenhout (2005). +# Processes from Cocco, Gomes, Maenhout (2005): +# Cocco, J. F., Gomes, F. J., & Maenhout, P. J. (2005). Consumption and +# portfolio choice over the life cycle. The Review of Financial Studies, +# 18(2), 491-533. +# - The profiles are provided as presented in the original paper. +# - It seem to us that income peaks at very young ages. +# - We suspect this might be due to the author's treatment of trends in income +# growth. +# - This can be adressed using the YearTrend and pLvlGroFacAgg options. CGM_income = { "NoHS": { "AgePolyCoefs": [-2.1361 + 2.6275, 0.1684 * 10, -0.0353 * 10, 0.0023 * 10], @@ -323,7 +402,9 @@ def findProfile(GroFacs, Y0): }, } -# Processes from Cagetti (2003). +# Processes from Cagetti (2003) +# Cagetti, M. (2003). Wealth accumulation over the life cycle and precautionary +# savings. Journal of Business & Economic Statistics, 21(3), 339-353. # - The author generously provided estimates from which the age polynomials # and yearly trends were recovered. # - He uses volatilities from Carroll-Samwick (1997) From 01fa2055013c0142594b09d7eb406df7ad5ae7c2 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 15:24:24 -0500 Subject: [PATCH 38/47] Merge SabSong and income file into single file --- .../{Calibration.py => Income/IncomeTools.py} | 414 +++++++++++++++--- .../Income/SabelhausSongProfiles.py | 303 ------------- 2 files changed, 358 insertions(+), 359 deletions(-) rename HARK/Calibration/{Calibration.py => Income/IncomeTools.py} (65%) delete mode 100644 HARK/Calibration/Income/SabelhausSongProfiles.py diff --git a/HARK/Calibration/Calibration.py b/HARK/Calibration/Income/IncomeTools.py similarity index 65% rename from HARK/Calibration/Calibration.py rename to HARK/Calibration/Income/IncomeTools.py index 33fdf492b..9aa22d604 100644 --- a/HARK/Calibration/Calibration.py +++ b/HARK/Calibration/Income/IncomeTools.py @@ -5,11 +5,58 @@ @author: Mateo """ +# %% Preamble + import numpy as np -from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile +from warnings import warn +from HARK.interpolation import LinearInterp from HARK.datasets.cpi.us.CPITools import cpi_deflator -__all__ = ["Cagetti_income", "CGM_income", "ParseIncomeSpec", "findProfile"] +__all__ = [ + "ParseTimeParams", + "Sabelhaus_Song_cohort_trend", + "Sabelhaus_Song_all_years", + "sabelhaus_song_var_profile", + "Cagetti_income", + "CGM_income", + "ParseIncomeSpec", + "findProfile", +] + + +# %% Tools for setting time-related parameters + + +def ParseTimeParams(age_birth, age_death): + """ + Converts simple statements of the age at which an agent is born and the + age at which he dies with certaintiy into the parameters that HARK needs + for figuring out the timing of the model. + + Parameters + ---------- + age_birth : int + Age at which the agent enters the model, e.g., 21. + age_death : int + Age at which the agent dies with certainty, e.g., 100. + + Returns + ------- + dict + Dictionary with parameters "T_cycle" and "T_age" which HARK expects + and which map to the birth and death ages specified by the user. + + """ + # T_cycle is the number of non-terminal periods in the agent's problem + T_cycle = age_death - age_birth + # T_age is the age at which the agents are killed with certainty in + # simulations (at the end of the T_age-th period) + T_age = age_death - age_birth + 1 + + return {"T_cycle": T_cycle, "T_age": T_age} + + +# %% Tools for finding the mean profiles of permanent income. def AgeLogPolyToGrowthRates(coefs, age_min, age_max): @@ -122,6 +169,314 @@ def findPermGroFacs(age_min, age_max, age_ret, AgePolyCoefs, ReplRate): return GroFacs, Y0 +def findProfile(GroFacs, Y0): + """ + Generates a sequence {Y_{t}}_{t=0}^N from an initial Y_0 and a sequence + of growth factors GroFac[n] = Y_{n+1}/Y_n + + Parameters + ---------- + GroFacs : list or numpy array + Growth factors in chronological order. + Y0 : float + initial value of the series. + + Returns + ------- + Y : numpy array + Array with the values of the series. + + """ + factors = np.array([Y0] + GroFacs) + Y = np.cumprod(factors) + + return Y + + +# %% Tools for life-cycle profiles of income volatility + +# The raw results shared by John Sabelhaus contain the following two +# sets of estimates (with and without cohor trends), which we will +# use for constructing the age profiles. + +# The first specification contains a cohort trend. The variance of +# (transitory or permanent) shocks to income of a person born in year +# "cohort" and who is now age "age" is +# age_dummy(age) + beta * (cohort - 1926) +# Where we have dummies for ages 27 to 54 +Sabelhaus_Song_cohort_trend = { + "Ages": np.arange(27, 55), + "AgeDummiesPrm": np.array( + [ + 0.0837941, + 0.0706855, + 0.0638561, + 0.0603879, + 0.0554693, + 0.0532388, + 0.0515262, + 0.0486079, + 0.0455297, + 0.0456573, + 0.0417433, + 0.0420146, + 0.0391508, + 0.0395776, + 0.0369826, + 0.0387158, + 0.0365356, + 0.036701, + 0.0364236, + 0.0358601, + 0.0348528, + 0.0362901, + 0.0373366, + 0.0372724, + 0.0401297, + 0.0415868, + 0.0434772, + 0.046668, + ] + ), + "AgeDummiesTrn": np.array( + [ + 0.1412842, + 0.1477754, + 0.1510265, + 0.1512203, + 0.1516837, + 0.151412, + 0.1489388, + 0.148521, + 0.1470632, + 0.143514, + 0.1411806, + 0.1378733, + 0.135245, + 0.1318365, + 0.1299689, + 0.1255799, + 0.1220823, + 0.1178995, + 0.1148793, + 0.1107577, + 0.1073337, + 0.102347, + 0.0962066, + 0.0918819, + 0.0887777, + 0.0835057, + 0.0766663, + 0.0698848, + ] + ), + "CohortCoefPrm": -0.0005966, + "CohortCoefTrn": -0.0017764 / 2, +} + +# The second specification contains no cohort trend. The variance of +# (transitory or permanent) shocks to income of a person born in year +# "cohort" and who is now age "age" is: age_dummy(age) +# Where we have dummies for ages 27 to 54. We use this "aggregate" +# specification if no cohort is provided. +Sabelhaus_Song_all_years = { + "Ages": np.arange(27, 55), + "AgeDummiesPrm": np.array( + [ + 0.0599296, + 0.0474176, + 0.0411848, + 0.0383132, + 0.0339912, + 0.0323573, + 0.0312414, + 0.0289196, + 0.0264381, + 0.0271623, + 0.0238449, + 0.0247128, + 0.0224456, + 0.0234691, + 0.0214706, + 0.0238005, + 0.0222169, + 0.0229789, + 0.0232982, + 0.0233312, + 0.0229205, + 0.0249545, + 0.0265975, + 0.02713, + 0.0305839, + 0.0326376, + 0.0351246, + 0.038912, + ] + ), + "AgeDummiesTrn": np.array( + [ + 0.1061999, + 0.1135794, + 0.1177187, + 0.1188007, + 0.1201523, + 0.1207688, + 0.1191838, + 0.1196542, + 0.1190846, + 0.1164236, + 0.1149784, + 0.1125594, + 0.1108192, + 0.1082989, + 0.1073195, + 0.1038188, + 0.1012094, + 0.0979148, + 0.0957829, + 0.0925495, + 0.0900136, + 0.0859151, + 0.0806629, + 0.0772264, + 0.0750105, + 0.0706267, + 0.0646755, + 0.0587821, + ] + ), + "CohortCoefPrm": 0, + "CohortCoefTrn": 0, +} + + +def sabelhaus_song_var_profile(age_min=27, age_max=54, cohort=None, smooth=True): + """ + This is a function to find the life-cycle profiles of the volatilities + of transitory and permanent shocks to income using the estimates in + [1] Sabelhaus and Song (2010). + + Parameters + ---------- + age_min : int, optional + Minimum age at which to construct volatilities. The default is 27. + age_max : int, optional + Maximum age at which to construct volatilities. The default is 54. + cohort : int, optional + Birth year of the hypothetical person for which the volatilities will + be constructed. The default is None, and in this case the we will + use the specification that does not have cohort trends. + smooth: bool, optional + Boolean indicating whether to smooth the variance profile estimates + using third degree polynomials for the age dummies estimated by + Sabelhaus and Song. If False, the original dummies are used. + + Returns + ------- + profiles : dict + Dictionary with entries: + - Ages: array of ages for which we found income volatilities in + ascending order + - TranShkStd: array of standard deviations of transitory income + shocks. Position n corresponds to Ages[n]. + - PermShkStd: array of standard deviations of permanent income + shocks. Position n corresponds to Ages[n]. + + Note that TransShkStd[n] and PermShkStd[n] are the volatilities of + shocks _experienced_ at age Age[n], (not those expected at Age[n+1] + from the perspective of Age[n]). + + Note that Sabelhaus and Song work in discrete time and with periods + that represent one year. Therefore, the outputs must be interpreted + at the yearly frequency. + """ + + assert age_max >= age_min, ( + "The maximum age can not be lower than the " + "minimum age." + ) + + # Determine which set of estimates to use based on wether a cohort is + # provided or not. + if cohort is None: + + spec = Sabelhaus_Song_all_years + cohort = 0 + warn("No cohort was provided. Using aggregate specification.") + + else: + + spec = Sabelhaus_Song_cohort_trend + + # Extract coefficients + beta_eps = spec["CohortCoefTrn"] + beta_eta = spec["CohortCoefPrm"] + tran_age_dummies = spec["AgeDummiesTrn"] + perm_age_dummies = spec["AgeDummiesPrm"] + + # Smooth out dummies using a 3rd degree polynomial if requested + if smooth: + + # Fit polynomials + tran_poly = np.poly1d(np.polyfit(spec["Ages"], tran_age_dummies, deg=3)) + perm_poly = np.poly1d(np.polyfit(spec["Ages"], perm_age_dummies, deg=3)) + + # Replace dummies + tran_age_dummies = tran_poly(spec["Ages"]) + perm_age_dummies = perm_poly(spec["Ages"]) + + # Make interpolators for transitory and permanent dummies. Alter to use + # flat extrapolation. + + # We use Sabelhaus and Song (2010) dummies for ages 27-54 and extrapolate + # outside of that just using the endpoints. + tran_dummy_interp = LinearInterp( + np.arange(min(spec["Ages"]) - 1, max(spec["Ages"]) + 2), + np.concatenate( + [[tran_age_dummies[0]], tran_age_dummies, [tran_age_dummies[-1]]] + ), + lower_extrap=True, + ) + + perm_dummy_interp = LinearInterp( + np.arange(min(spec["Ages"]) - 1, max(spec["Ages"]) + 2), + np.concatenate( + [[perm_age_dummies[0]], perm_age_dummies, [perm_age_dummies[-1]]] + ), + lower_extrap=True, + ) + + if age_min < 27 or age_max > 54: + warn( + "Sabelhaus and Song (2010) provide variance profiles for ages " + + "27 to 54. Extrapolating variances using the extreme points." + ) + + if cohort < 1926 or cohort > 1980: + warn( + "Sabelhaus and Song (2010) use data from birth cohorts " + + "[1926,1980]. Extrapolating variances." + ) + + cohort = max(min(cohort, 1980), 1926) + + # Construct variances + # They use 1926 as the base year for cohort effects. + ages = np.arange(age_min, age_max + 1) + tran_std = tran_dummy_interp(ages) + (cohort - 1926) * beta_eps + perm_std = perm_dummy_interp(ages) + (cohort - 1926) * beta_eta + + profiles = { + "Age": list(ages), + "TranShkStd": list(tran_std), + "PermShkStd": list(perm_std), + } + + return profiles + + +# %% Encompassing tool to parse full income specifications + + def ParseIncomeSpec( base_monet_year, age_min, @@ -342,29 +697,7 @@ def ParseIncomeSpec( return income_params -def findProfile(GroFacs, Y0): - """ - Generates a sequence {Y_{t}}_{t=0}^N from an initial Y_0 and a sequence - of growth factors GroFac[n] = Y_{n+1}/Y_n - - Parameters - ---------- - GroFacs : list or numpy array - Growth factors in chronological order. - Y0 : float - initial value of the series. - - Returns - ------- - Y : numpy array - Array with the values of the series. - - """ - factors = np.array([Y0] + GroFacs) - Y = np.cumprod(factors) - - return Y - +# %% Income specifications from various papers # Processes from Cocco, Gomes, Maenhout (2005): # Cocco, J. F., Gomes, F. J., & Maenhout, P. J. (2005). Consumption and @@ -459,34 +792,3 @@ def findProfile(GroFacs, Y0): "base_monet_year": 1992, }, } - -# %% Tools for setting time-related parameters - - -def ParseTimeParams(age_birth, age_death): - """ - Converts simple statements of the age at which an agent is born and the - age at which he dies with certaintiy into the parameters that HARK needs - for figuring out the timing of the model. - - Parameters - ---------- - age_birth : int - Age at which the agent enters the model, e.g., 21. - age_death : int - Age at which the agent dies with certainty, e.g., 100. - - Returns - ------- - dict - Dictionary with parameters "T_cycle" and "T_age" which HARK expects - and which map to the birth and death ages specified by the user. - - """ - # T_cycle is the number of non-terminal periods in the agent's problem - T_cycle = age_death - age_birth - # T_age is the age at which the agents are killed with certainty in - # simulations (at the end of the T_age-th period) - T_age = age_death - age_birth + 1 - - return {"T_cycle": T_cycle, "T_age": T_age} diff --git a/HARK/Calibration/Income/SabelhausSongProfiles.py b/HARK/Calibration/Income/SabelhausSongProfiles.py deleted file mode 100644 index a9b145dce..000000000 --- a/HARK/Calibration/Income/SabelhausSongProfiles.py +++ /dev/null @@ -1,303 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Jan 14 15:49:43 2021 - -@author: Mateo - -This scripts contains representations of the results in: - -[1] Sabelhaus, J., & Song, J. (2010). The great moderation in micro labor - earnings. Journal of Monetary Economics, 57(4), 391-403. - -It provides functions to produce life-cycle profiles of the variances of -income shocks from the raw results from the paper. - -The raw estimates were generously shared by John Sabelhaus. - -""" - -import numpy as np -from warnings import warn -from HARK.interpolation import LinearInterp - - -# The raw results shared by John Sabelhaus contain the following two -# sets of estimates (with and without cohor trends), which we will -# use for constructing the age profiles. - -# The first specification contains a cohort trend. The variance of -# (transitory or permanent) shocks to income of a person born in year -# "cohort" and who is now age "age" is -# age_dummy(age) + beta * (cohort - 1926) -# Where we have dummies for ages 27 to 54 -Sabelhaus_Song_cohort_trend = { - "Ages": np.arange(27, 55), - "AgeDummiesPrm": np.array( - [ - 0.0837941, - 0.0706855, - 0.0638561, - 0.0603879, - 0.0554693, - 0.0532388, - 0.0515262, - 0.0486079, - 0.0455297, - 0.0456573, - 0.0417433, - 0.0420146, - 0.0391508, - 0.0395776, - 0.0369826, - 0.0387158, - 0.0365356, - 0.036701, - 0.0364236, - 0.0358601, - 0.0348528, - 0.0362901, - 0.0373366, - 0.0372724, - 0.0401297, - 0.0415868, - 0.0434772, - 0.046668, - ] - ), - "AgeDummiesTrn": np.array( - [ - 0.1412842, - 0.1477754, - 0.1510265, - 0.1512203, - 0.1516837, - 0.151412, - 0.1489388, - 0.148521, - 0.1470632, - 0.143514, - 0.1411806, - 0.1378733, - 0.135245, - 0.1318365, - 0.1299689, - 0.1255799, - 0.1220823, - 0.1178995, - 0.1148793, - 0.1107577, - 0.1073337, - 0.102347, - 0.0962066, - 0.0918819, - 0.0887777, - 0.0835057, - 0.0766663, - 0.0698848, - ] - ), - "CohortCoefPrm": -0.0005966, - "CohortCoefTrn": -0.0017764 / 2, -} - -# The second specification contains no cohort trend. The variance of -# (transitory or permanent) shocks to income of a person born in year -# "cohort" and who is now age "age" is: age_dummy(age) -# Where we have dummies for ages 27 to 54. We use this "aggregate" -# specification if no cohort is provided. -Sabelhaus_Song_all_years = { - "Ages": np.arange(27, 55), - "AgeDummiesPrm": np.array( - [ - 0.0599296, - 0.0474176, - 0.0411848, - 0.0383132, - 0.0339912, - 0.0323573, - 0.0312414, - 0.0289196, - 0.0264381, - 0.0271623, - 0.0238449, - 0.0247128, - 0.0224456, - 0.0234691, - 0.0214706, - 0.0238005, - 0.0222169, - 0.0229789, - 0.0232982, - 0.0233312, - 0.0229205, - 0.0249545, - 0.0265975, - 0.02713, - 0.0305839, - 0.0326376, - 0.0351246, - 0.038912, - ] - ), - "AgeDummiesTrn": np.array( - [ - 0.1061999, - 0.1135794, - 0.1177187, - 0.1188007, - 0.1201523, - 0.1207688, - 0.1191838, - 0.1196542, - 0.1190846, - 0.1164236, - 0.1149784, - 0.1125594, - 0.1108192, - 0.1082989, - 0.1073195, - 0.1038188, - 0.1012094, - 0.0979148, - 0.0957829, - 0.0925495, - 0.0900136, - 0.0859151, - 0.0806629, - 0.0772264, - 0.0750105, - 0.0706267, - 0.0646755, - 0.0587821, - ] - ), - "CohortCoefPrm": 0, - "CohortCoefTrn": 0, -} - - -def sabelhaus_song_var_profile(age_min=27, age_max=54, cohort=None, - smooth = True): - """ - This is a function to find the life-cycle profiles of the volatilities - of transitory and permanent shocks to income using the estimates in - [1] Sabelhaus and Song (2010). - - Parameters - ---------- - age_min : int, optional - Minimum age at which to construct volatilities. The default is 27. - age_max : int, optional - Maximum age at which to construct volatilities. The default is 54. - cohort : int, optional - Birth year of the hypothetical person for which the volatilities will - be constructed. The default is None, and in this case the we will - use the specification that does not have cohort trends. - smooth: bool, optional - Boolean indicating whether to smooth the variance profile estimates - using third degree polynomials for the age dummies estimated by - Sabelhaus and Song. If False, the original dummies are used. - - Returns - ------- - profiles : dict - Dictionary with entries: - - Ages: array of ages for which we found income volatilities in - ascending order - - TranShkStd: array of standard deviations of transitory income - shocks. Position n corresponds to Ages[n]. - - PermShkStd: array of standard deviations of permanent income - shocks. Position n corresponds to Ages[n]. - - Note that TransShkStd[n] and PermShkStd[n] are the volatilities of - shocks _experienced_ at age Age[n], (not those expected at Age[n+1] - from the perspective of Age[n]). - - Note that Sabelhaus and Song work in discrete time and with periods - that represent one year. Therefore, the outputs must be interpreted - at the yearly frequency. - """ - - assert age_max >= age_min, ( - "The maximum age can not be lower than the " + "minimum age." - ) - - # Determine which set of estimates to use based on wether a cohort is - # provided or not. - if cohort is None: - - spec = Sabelhaus_Song_all_years - cohort = 0 - warn("No cohort was provided. Using aggregate specification.") - - else: - - spec = Sabelhaus_Song_cohort_trend - - # Extract coefficients - beta_eps = spec["CohortCoefTrn"] - beta_eta = spec["CohortCoefPrm"] - tran_age_dummies = spec["AgeDummiesTrn"] - perm_age_dummies = spec["AgeDummiesPrm"] - - # Smooth out dummies using a 3rd degree polynomial if requested - if smooth: - - # Fit polynomials - tran_poly = np.poly1d(np.polyfit(spec['Ages'], - tran_age_dummies, deg = 3)) - perm_poly = np.poly1d(np.polyfit(spec['Ages'], - perm_age_dummies, deg = 3)) - - # Replace dummies - tran_age_dummies = tran_poly(spec['Ages']) - perm_age_dummies = perm_poly(spec['Ages']) - - # Make interpolators for transitory and permanent dummies. Alter to use - # flat extrapolation. - - # We use Sabelhaus and Song (2010) dummies for ages 27-54 and extrapolate - # outside of that just using the endpoints. - tran_dummy_interp = LinearInterp( - np.arange(min(spec['Ages'])-1, max(spec['Ages'])+2), - np.concatenate( - [[tran_age_dummies[0]], tran_age_dummies, [tran_age_dummies[-1]]] - ), - lower_extrap=True, - ) - - perm_dummy_interp = LinearInterp( - np.arange(min(spec['Ages'])-1, max(spec['Ages'])+2), - np.concatenate( - [[perm_age_dummies[0]], perm_age_dummies, [perm_age_dummies[-1]]] - ), - lower_extrap=True, - ) - - if age_min < 27 or age_max > 54: - warn( - "Sabelhaus and Song (2010) provide variance profiles for ages " - + "27 to 54. Extrapolating variances using the extreme points." - ) - - if cohort < 1926 or cohort > 1980: - warn( - "Sabelhaus and Song (2010) use data from birth cohorts " - + "[1926,1980]. Extrapolating variances." - ) - - cohort = max(min(cohort, 1980), 1926) - - # Construct variances - # They use 1926 as the base year for cohort effects. - ages = np.arange(age_min, age_max + 1) - tran_std = tran_dummy_interp(ages) + (cohort - 1926) * beta_eps - perm_std = perm_dummy_interp(ages) + (cohort - 1926) * beta_eta - - profiles = { - "Age": list(ages), - "TranShkStd": list(tran_std), - "PermShkStd": list(perm_std), - } - - return profiles From e721e3af62d71471230e0d35785925ba55d1412d Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 15:47:47 -0500 Subject: [PATCH 39/47] Adjust examples --- examples/Calibration/Income_calibrations.py | 2 +- examples/Calibration/Sabelhaus_Song_var_profiles.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/Calibration/Income_calibrations.py b/examples/Calibration/Income_calibrations.py index 879d46618..33234764c 100644 --- a/examples/Calibration/Income_calibrations.py +++ b/examples/Calibration/Income_calibrations.py @@ -5,7 +5,7 @@ @author: Mateo """ -from HARK.Calibration.Calibration import ( +from HARK.Calibration.Income.IncomeTools import ( ParseIncomeSpec, findProfile, Cagetti_income, diff --git a/examples/Calibration/Sabelhaus_Song_var_profiles.py b/examples/Calibration/Sabelhaus_Song_var_profiles.py index 636fc9b2c..14bb8b766 100644 --- a/examples/Calibration/Sabelhaus_Song_var_profiles.py +++ b/examples/Calibration/Sabelhaus_Song_var_profiles.py @@ -14,7 +14,7 @@ """ import matplotlib.pyplot as plt -from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile +from HARK.Calibration.Income.IncomeTools import sabelhaus_song_var_profile # Set up ages and cohorts at which we will get the variances age_min = 27 From 66a8620814ba553aaaf4be9ba453d355b2918ff7 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 15:48:22 -0500 Subject: [PATCH 40/47] Add tests for mean income profiles --- ...ausSongProfiles.py => test_IncomeTools.py} | 223 +++++++++++++++++- 1 file changed, 210 insertions(+), 13 deletions(-) rename HARK/Calibration/Income/tests/{test_SabelhausSongProfiles.py => test_IncomeTools.py} (58%) diff --git a/HARK/Calibration/Income/tests/test_SabelhausSongProfiles.py b/HARK/Calibration/Income/tests/test_IncomeTools.py similarity index 58% rename from HARK/Calibration/Income/tests/test_SabelhausSongProfiles.py rename to HARK/Calibration/Income/tests/test_IncomeTools.py index 9e781d3ca..70ef8835f 100644 --- a/HARK/Calibration/Income/tests/test_SabelhausSongProfiles.py +++ b/HARK/Calibration/Income/tests/test_IncomeTools.py @@ -8,9 +8,210 @@ # Bring in modules we need import unittest import numpy as np -from HARK.Calibration.Income.SabelhausSongProfiles import sabelhaus_song_var_profile +from HARK.Calibration.Income.IncomeTools import ( + sabelhaus_song_var_profile, + ParseIncomeSpec, + findProfile, + CGM_income, + Cagetti_income, +) +# %% Mean income profile tests +class test_income_paths(unittest.TestCase): + def setUp(self): + + # Assign a result from Cocco, Gomes, Maenhout + self.cgm_hs_mean_p = np.array( + [ + 16.8338, + 17.8221, + 18.7965, + 19.7509, + 20.6796, + 21.5773, + 22.4388, + 23.2596, + 24.0359, + 24.7642, + 25.4417, + 26.0662, + 26.6362, + 27.1506, + 27.6092, + 28.0122, + 28.3603, + 28.6548, + 28.8974, + 29.0902, + 29.2357, + 29.3367, + 29.3963, + 29.4178, + 29.4046, + 29.3602, + 29.2884, + 29.1927, + 29.0771, + 28.9451, + 28.8004, + 28.6468, + 28.4876, + 28.3266, + 28.167, + 28.0122, + 27.8655, + 27.7301, + 27.6092, + 27.5059, + 27.4232, + 27.3643, + 27.3323, + 27.3304, + 27.3619, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + 18.6635, + ] + ) + + # And a result from Cagetti + self.cagetti_college_mean_p = np.array( + [ + 29.8545, + 31.0534, + 32.2894, + 33.5602, + 34.8632, + 36.1954, + 37.5537, + 38.9343, + 40.3332, + 41.746, + 43.1681, + 44.5943, + 46.0193, + 47.4376, + 48.8433, + 50.2304, + 51.5926, + 52.9237, + 54.2173, + 55.4672, + 56.6672, + 57.8109, + 58.8927, + 59.9067, + 60.8477, + 61.7106, + 62.4908, + 63.1844, + 63.7878, + 64.2979, + 64.7123, + 65.0294, + 65.2479, + 65.3675, + 65.3882, + 65.311, + 65.1372, + 64.8691, + 64.5092, + 64.0609, + 63.5278, + 31.8833, + 31.8638, + 31.8444, + 31.825, + 31.8056, + 31.7862, + 31.7668, + 31.7474, + 31.728, + 31.7087, + 31.6893, + 31.67, + 31.6507, + 31.6314, + 31.6121, + 31.5928, + 31.5735, + 31.5542, + 31.535, + 31.5158, + 31.4965, + 31.4773, + 31.4581, + 31.4389, + 31.4197, + 31.4006, + ] + ) + + def test_CGM(self): + + adjust_infl_to = 1992 + age_min = 21 + age_max = 100 + spec = CGM_income["HS"] + params = ParseIncomeSpec( + age_min=age_min, age_max=age_max, adjust_infl_to=adjust_infl_to, **spec + ) + MeanP = findProfile(params["PermGroFac"], params["P0"]) + self.assertTrue(np.allclose(self.cgm_hs_mean_p, MeanP, atol=1e-03)) + + def test_Cagetti(self): + + adjust_infl_to = 1992 + age_min = 25 + age_max = 91 + start_year = 1980 + spec = Cagetti_income["College"] + params = ParseIncomeSpec( + age_min=age_min, + age_max=age_max, + adjust_infl_to=adjust_infl_to, + start_year=start_year, + **spec + ) + MeanP = findProfile(params["PermGroFac"], params["P0"]) + + self.assertTrue(np.allclose(self.cagetti_college_mean_p, MeanP, atol=1e-03)) + + +# %% Volatility profile tests class test_SabelhausSongProfiles(unittest.TestCase): def setUp(self): @@ -283,38 +484,34 @@ def test_smoothing(self): # 1940 self.assertTrue( np.allclose( - np.array(smooth1940["TranShkStd"]), self.Fig6Coh1940Tran, rtol= rtol + np.array(smooth1940["TranShkStd"]), self.Fig6Coh1940Tran, rtol=rtol ) ) self.assertTrue( np.allclose( - np.array(smooth1940["PermShkStd"]), self.Fig6Coh1940Perm, atol= rtol + np.array(smooth1940["PermShkStd"]), self.Fig6Coh1940Perm, atol=rtol ) ) - + # 1965 self.assertTrue( np.allclose( - np.array(smooth1965["TranShkStd"]), self.Fig6Coh1965Tran, rtol= rtol + np.array(smooth1965["TranShkStd"]), self.Fig6Coh1965Tran, rtol=rtol ) ) self.assertTrue( np.allclose( - np.array(smooth1965["PermShkStd"]), self.Fig6Coh1965Perm, atol= rtol + np.array(smooth1965["PermShkStd"]), self.Fig6Coh1965Perm, atol=rtol ) ) - + # Aggregate self.assertTrue( - np.allclose( - np.array(smoothAgg["TranShkStd"]), self.AggTran, rtol= rtol - ) + np.allclose(np.array(smoothAgg["TranShkStd"]), self.AggTran, rtol=rtol) ) self.assertTrue( - np.allclose( - np.array(smoothAgg["PermShkStd"]), self.AggPerm, atol= rtol - ) + np.allclose(np.array(smoothAgg["PermShkStd"]), self.AggPerm, atol=rtol) ) From 63dc98173ad48f76cccecef23b695518751f3611 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateo=20Vel=C3=A1squez-Giraldo?= Date: Fri, 5 Feb 2021 16:00:19 -0500 Subject: [PATCH 41/47] Update README.md --- HARK/Calibration/Income/README.md | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/HARK/Calibration/Income/README.md b/HARK/Calibration/Income/README.md index e83abe6cd..d9cff42ef 100644 --- a/HARK/Calibration/Income/README.md +++ b/HARK/Calibration/Income/README.md @@ -2,8 +2,14 @@ This folder contains various tools that are used for replicating calibrations/estimates of income processes from the literature. -## `SabelhausSongProfiles.py` +## `IncomeTools.py` -This file contains tools for producing life-cycle profiles of the volatilities of transitory and permanent shocks to income using the results from [Sabelhaus & Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358). +This is the only file for now. Its main components are: -The estimates used to reproduce the results from the paper were generously provided by [John Sabelhaus](https://www.sites.google.com/view/johnsabelhaus). +- `sabelhaus_song_var_profile`: a function for producing life-cycle profiles of the volatilities of transitory and permanent shocks to income using the results from [Sabelhaus & Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358). The estimates used to reproduce the results from the paper were generously provided by [John Sabelhaus](https://www.sites.google.com/view/johnsabelhaus). + +- `Cagetti_income`: a representation of the income specification used by [Cagetti (2003)](https://www.jstor.org/stable/1392584?seq=1). The estimates used to reproduce these results were generously shared by [Marco Cagetti](https://www.federalreserve.gov/econres/marco-cagetti.htm). + +- `CGM_income`: a representation of the income specification used by [Cocco, Gomes, and Maenhout (2005)](https://academic.oup.com/rfs/article-abstract/18/2/491/1599892). The estimates used to reproduce these results come directly from the published version of the paper. + +- `ParseIncomeSpec`: a function that takes in calibrations in various formats that are common in the literature and produces the parameters that represent them in the format that HARK expects. From a6fa62ebc0e2a833b1c0a237382976a0a0603fcc Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 16:08:17 -0500 Subject: [PATCH 42/47] Improve and format LC example --- examples/Calibration/Life_Cycle_example.py | 72 ++++++++++++---------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index e7c3ab49e..b8c3c2662 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -15,15 +15,14 @@ # %% from HARK.ConsumptionSaving.ConsIndShockModel import ( -IndShockConsumerType, -init_lifecycle + IndShockConsumerType, + init_lifecycle, ) -from HARK.Calibration.Calibration import ( +from HARK.Calibration.Income.IncomeTools import ( ParseIncomeSpec, ParseTimeParams, CGM_income, - Cagetti_income ) from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table @@ -35,29 +34,39 @@ # %% Alter calibration birth_age = 25 death_age = 90 +adjust_infl_to = 1992 income_calib = CGM_income -education = 'College' +education = "College" # Income specification -income_params = ParseIncomeSpec(age_min = birth_age, age_max = death_age, - **income_calib[education], SabelhausSong=True) +income_params = ParseIncomeSpec( + age_min=birth_age, + age_max=death_age, + adjust_infl_to=adjust_infl_to, + **income_calib[education], + SabelhausSong=True +) # Initial distribution of wealth and permanent income -dist_params = income_wealth_dists_from_scf(base_year=1992, age = birth_age, - education = education, wave = 1995) +dist_params = income_wealth_dists_from_scf( + base_year=adjust_infl_to, age=birth_age, education=education, wave=1995 +) # We need survival probabilities only up to death_age-1, because survival # probability at death_age is 1. -liv_prb = parse_ssa_life_table(female = True, cross_sec = True, year = 2004, - min_age = birth_age, max_age = death_age - 1) +liv_prb = parse_ssa_life_table( + female=True, cross_sec=True, year=2004, min_age=birth_age, max_age=death_age - 1 +) -time_params = ParseTimeParams(age_birth = birth_age, age_death = death_age) +# Parameters related to the number of periods implied by the calibration +time_params = ParseTimeParams(age_birth=birth_age, age_death=death_age) +# Update all the new parameters params = copy(init_lifecycle) params.update(time_params) params.update(dist_params) params.update(income_params) -params.update({'LivPrb': liv_prb}) +params.update({"LivPrb": liv_prb}) # %% Create and solve agent Agent = IndShockConsumerType(**params) @@ -65,14 +74,12 @@ # %% Simulation -# Setup - # Number of agents and periods in the simulation. Agent.AgentCount = 500 -Agent.T_sim = 200 +Agent.T_sim = 200 # Set up the variables we want to keep track of. -Agent.track_vars = ['aNrmNow','cNrmNow', 'pLvlNow', 't_age','mNrmNow'] +Agent.track_vars = ["aNrmNow", "cNrmNow", "pLvlNow", "t_age", "mNrmNow"] # Run the simulations Agent.initializeSim() @@ -80,29 +87,28 @@ # %% Extract and format simulation results -raw_data = {'Age': Agent.history['t_age'].flatten()+birth_age - 1, - 'pIncome': Agent.history['pLvlNow'].flatten(), - 'nrmM': Agent.history['mNrmNow'].flatten(), - 'nrmC': Agent.history['cNrmNow'].flatten()} +raw_data = { + "Age": Agent.history["t_age"].flatten() + birth_age - 1, + "pIncome": Agent.history["pLvlNow"].flatten(), + "nrmM": Agent.history["mNrmNow"].flatten(), + "nrmC": Agent.history["cNrmNow"].flatten(), +} Data = pd.DataFrame(raw_data) -Data['Cons'] = Data.nrmC * Data.pIncome -Data['M'] = Data.nrmM * Data.pIncome +Data["Cons"] = Data.nrmC * Data.pIncome +Data["M"] = Data.nrmM * Data.pIncome # %% Plots # Find the mean of each variable at every age -AgeMeans = Data.groupby(['Age']).mean().reset_index() +AgeMeans = Data.groupby(["Age"]).median().reset_index() plt.figure() -plt.plot(AgeMeans.Age, AgeMeans.pIncome, - label = 'Permanent Income') -plt.plot(AgeMeans.Age, AgeMeans.M, - label = 'Market resources') -plt.plot(AgeMeans.Age, AgeMeans.Cons, - label = 'Consumption') +plt.plot(AgeMeans.Age, AgeMeans.pIncome, label="Permanent Income") +plt.plot(AgeMeans.Age, AgeMeans.M, label="Market resources") +plt.plot(AgeMeans.Age, AgeMeans.Cons, label="Consumption") plt.legend() -plt.xlabel('Age') -plt.ylabel('Thousands of 1992 USD') -plt.title('Variable Means Conditional on Survival') +plt.xlabel("Age") +plt.ylabel("Thousands of {} USD".format(adjust_infl_to)) +plt.title("Variable Medians Conditional on Survival") plt.grid() From 47c65592741d6b36ec0f88313c8a471643fe9679 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 16:10:02 -0500 Subject: [PATCH 43/47] Remove initial calibrations from other files --- HARK/Calibration/CGMRemark.py | 216 ----------------- HARK/Calibration/Pandemic.py | 334 -------------------------- HARK/Calibration/SolvingMicroDSOPs.py | 173 ------------- 3 files changed, 723 deletions(-) delete mode 100644 HARK/Calibration/CGMRemark.py delete mode 100644 HARK/Calibration/Pandemic.py delete mode 100644 HARK/Calibration/SolvingMicroDSOPs.py diff --git a/HARK/Calibration/CGMRemark.py b/HARK/Calibration/CGMRemark.py deleted file mode 100644 index 224454907..000000000 --- a/HARK/Calibration/CGMRemark.py +++ /dev/null @@ -1,216 +0,0 @@ -import numpy as np -# %% Preferences - -# Relative risk aversion -CRRA = 10 -# Discount factor -DiscFac = 0.96 - -# Survival probabilities from the author's Fortran code -n = 80 -survprob = np.zeros(n+1) -survprob[1] = 0.99845 -survprob[2] = 0.99839 -survprob[3] = 0.99833 -survprob[4] = 0.9983 -survprob[5] = 0.99827 -survprob[6] = 0.99826 -survprob[7] = 0.99824 -survprob[8] = 0.9982 -survprob[9] = 0.99813 -survprob[10] = 0.99804 -survprob[11] = 0.99795 -survprob[12] = 0.99785 -survprob[13] = 0.99776 -survprob[14] = 0.99766 -survprob[15] = 0.99755 -survprob[16] = 0.99743 -survprob[17] = 0.9973 -survprob[18] = 0.99718 -survprob[19] = 0.99707 -survprob[20] = 0.99696 -survprob[21] = 0.99685 -survprob[22] = 0.99672 -survprob[23] = 0.99656 -survprob[24] = 0.99635 -survprob[25] = 0.9961 -survprob[26] = 0.99579 -survprob[27] = 0.99543 -survprob[28] = 0.99504 -survprob[29] = 0.99463 -survprob[30] = 0.9942 -survprob[31] = 0.9937 -survprob[32] = 0.99311 -survprob[33] = 0.99245 -survprob[34] = 0.99172 -survprob[35] = 0.99091 -survprob[36] = 0.99005 -survprob[37] = 0.98911 -survprob[38] = 0.98803 -survprob[39] = 0.9868 -survprob[40] = 0.98545 -survprob[41] = 0.98409 -survprob[42] = 0.9827 -survprob[43] = 0.98123 -survprob[44] = 0.97961 -survprob[45] = 0.97786 -survprob[46] = 0.97603 -survprob[47] = 0.97414 -survprob[48] = 0.97207 -survprob[49] = 0.9697 -survprob[50] = 0.96699 -survprob[51] = 0.96393 -survprob[52] = 0.96055 -survprob[53] = 0.9569 -survprob[54] = 0.9531 -survprob[55] = 0.94921 -survprob[56] = 0.94508 -survprob[57] = 0.94057 -survprob[58] = 0.9357 -survprob[59] = 0.93031 -survprob[60] = 0.92424 -survprob[61] = 0.91717 -survprob[62] = 0.90922 -survprob[63] = 0.90089 -survprob[64] = 0.89282 -survprob[65] = 0.88503 -survprob[66] = 0.87622 -survprob[67] = 0.86576 -survprob[68] = 0.8544 -survprob[69] = 0.8423 -survprob[70] = 0.82942 -survprob[71] = 0.8154 -survprob[72] = 0.80002 -survprob[73] = 0.78404 -survprob[74] = 0.76842 -survprob[75] = 0.75382 -survprob[76] = 0.73996 -survprob[77] = 0.72464 -survprob[78] = 0.71057 -survprob[79] = 0.6961 -survprob[80] = 0.6809 - -# Fix indexing problem (fortran starts at 1, python at 0) -survprob = np.delete(survprob, [0]) -# Now we have 80 probabilities of death, -# for ages 20 to 99. - -# Labor income - -# They assume its a polinomial of age. Here are the coefficients -a=-2.170042+2.700381 -b1=0.16818 -b2=-0.0323371/10 -b3=0.0019704/100 - -time_params = {'Age_born': 20, 'Age_retire': 65, 'Age_death': 100} -t_start = time_params['Age_born'] -t_ret = time_params['Age_retire'] # We are currently interpreting this as the last period of work -t_end = time_params['Age_death'] - -# They assume retirement income is a fraction of labor income in the -# last working period -repl_fac = 0.68212 - -# Compute average income at each point in (working) life -f = np.arange(t_start, t_ret+1,1) -f = a + b1*f + b2*(f**2) + b3*(f**3) -det_work_inc = np.exp(f) - -# Retirement income -det_ret_inc = repl_fac*det_work_inc[-1]*np.ones(t_end - t_ret) - -# Get a full vector of the deterministic part of income -det_income = np.concatenate((det_work_inc, det_ret_inc)) - -# ln Gamma_t+1 = ln f_t+1 - ln f_t -gr_fac = np.exp(np.diff(np.log(det_income))) - -# Now we have growth factors for T_end-1 periods. - -# Finally define the normalization factor used by CGM, for plots. -# ### IMPORTANT ### -# We adjust this normalization factor for what we believe is a typo in the -# original article. See the REMARK jupyter notebook for details. -norm_factor = det_income * np.exp(0) - -# %% Shocks - -# Transitory and permanent shock variance from the paper -std_tran_shock = np.sqrt(0.0738) -std_perm_shock = np.sqrt(0.0106) - -# Vectorize. (HARK turns off these shocks after T_retirement) -std_tran_vec = np.array([std_tran_shock]*(t_end-t_start)) -std_perm_vec = np.array([std_perm_shock]*(t_end-t_start)) - -# %% Financial instruments - -# Risk-free factor -Rfree = 1.02 - -# Creation of risky asset return distributions - -Mu = 0.06 # Equity premium -Std = 0.157 # standard deviation of rate-of-return shocks - -RiskyAvg = Mu + Rfree -RiskyStd = Std - -# Make a dictionary to specify the rest of params -dict_portfolio = { - # Usual params - 'CRRA': CRRA, - 'Rfree': Rfree, - 'DiscFac': DiscFac, - - # Life cycle - 'T_age' : t_end-t_start+1, # Time of death - 'T_cycle' : t_end-t_start, # Number of non-terminal periods - 'T_retire':t_ret-t_start+1, - 'LivPrb': survprob.tolist(), - 'PermGroFac': gr_fac.tolist(), - 'cycles': 1, - - # Income shocks - 'PermShkStd': std_perm_vec, - 'PermShkCount': 3, - 'TranShkStd': std_tran_vec, - 'TranShkCount': 3, - 'UnempPrb': 0, - 'UnempPrbRet': 0, - 'IncUnemp': 0, - 'IncUnempRet': 0, - 'BoroCnstArt': 0, - 'tax_rate':0.0, - - # Portfolio related params - 'RiskyAvg': RiskyAvg, - 'RiskyStd': RiskyStd, - 'RiskyCount': 3, - 'RiskyShareCount': 30, - - # Grid - 'aXtraMin': 0.001, - 'aXtraMax': 400, - 'aXtraCount': 400, - 'aXtraExtra': [None], - 'aXtraNestFac': 3, - - # General - 'vFuncBool': False, - 'CubicBool': False, - - # Simulation params - 'AgentCount': 10, - 'pLvlInitMean' : np.log(det_income[0]), # Mean of log initial permanent income (only matters for simulation) - 'pLvlInitStd' : std_perm_shock, # Standard deviation of log initial permanent income (only matters for simulation) - 'T_sim': (t_end - t_start+1)*50, - - # Unused params required for simulation - 'PermGroFacAgg': 1, - 'aNrmInitMean': -50.0, # Agents start with 0 assets (this is log-mean) - 'aNrmInitStd' : 0.0 -} - -age_plot_params = [20, 30, 55, 75] diff --git a/HARK/Calibration/Pandemic.py b/HARK/Calibration/Pandemic.py deleted file mode 100644 index 84c5f3696..000000000 --- a/HARK/Calibration/Pandemic.py +++ /dev/null @@ -1,334 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -import os -import csv -from HARK.distribution import Uniform -from importlib import reload -figs_dir = '../../Figures/' - -# Import configurable parameters, and keep them updated through reload -import parameter_config -reload(parameter_config) -from parameter_config import * - -############################################################################### - -# Size of simulations -AgentCountTotal = 1000000 # Total simulated population -T_sim = 13 # Number of quarters to simulate in counterfactuals - -# Basic lifecycle length parameters (don't touch these) -init_age = 24 -working_T = 41*4 # Number of working periods -retired_T = 55*4 # Number of retired periods -T_cycle = working_T + retired_T - -# Define the distribution of the discount factor for each eduation level -DiscFacCount = 7 -DiscFacDstnD = Uniform(DiscFacMeanD-DiscFacSpread, DiscFacMeanD+DiscFacSpread).approx(DiscFacCount) -DiscFacDstnH = Uniform(DiscFacMeanH-DiscFacSpread, DiscFacMeanH+DiscFacSpread).approx(DiscFacCount) -DiscFacDstnC = Uniform(DiscFacMeanC-DiscFacSpread, DiscFacMeanC+DiscFacSpread).approx(DiscFacCount) -DiscFacDstns = [DiscFacDstnD, DiscFacDstnH, DiscFacDstnC] - -# Define permanent income growth rates for each education level (from Cagetti 2003) -PermGroRte_d_ann = [5.2522391e-002, 5.0039782e-002, 4.7586132e-002, 4.5162424e-002, 4.2769638e-002, 4.0408757e-002, 3.8080763e-002, 3.5786635e-002, 3.3527358e-002, 3.1303911e-002, 2.9117277e-002, 2.6968437e-002, 2.4858374e-002, 2.2788068e-002, 2.0758501e-002, 1.8770655e-002, 1.6825511e-002, 1.4924052e-002, 1.3067258e-002, 1.1256112e-002, 9.4915947e-003, 7.7746883e-003, 6.1063742e-003, 4.4876340e-003, 2.9194495e-003, 1.4028022e-003, -6.1326258e-005, -1.4719542e-003, -2.8280999e-003, -4.1287819e-003, -5.3730185e-003, -6.5598280e-003, -7.6882288e-003, -8.7572392e-003, -9.7658777e-003, -1.0713163e-002, -1.1598112e-002, -1.2419745e-002, -1.3177079e-002, -1.3869133e-002, -4.3985368e-001, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003, -8.5623256e-003] -PermGroRte_h_ann = [4.1102173e-002, 4.1194381e-002, 4.1117402e-002, 4.0878307e-002, 4.0484168e-002, 3.9942056e-002, 3.9259042e-002, 3.8442198e-002, 3.7498596e-002, 3.6435308e-002, 3.5259403e-002, 3.3977955e-002, 3.2598035e-002, 3.1126713e-002, 2.9571062e-002, 2.7938153e-002, 2.6235058e-002, 2.4468848e-002, 2.2646594e-002, 2.0775369e-002, 1.8862243e-002, 1.6914288e-002, 1.4938576e-002, 1.2942178e-002, 1.0932165e-002, 8.9156095e-003, 6.8995825e-003, 4.8911556e-003, 2.8974003e-003, 9.2538802e-004, -1.0178097e-003, -2.9251214e-003, -4.7894755e-003, -6.6038005e-003, -8.3610250e-003, -1.0054077e-002, -1.1675886e-002, -1.3219380e-002, -1.4677487e-002, -1.6043137e-002, -5.5864350e-001, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002, -1.0820465e-002] -PermGroRte_c_ann = [3.9375106e-002, 3.9030288e-002, 3.8601230e-002, 3.8091011e-002, 3.7502710e-002, 3.6839406e-002, 3.6104179e-002, 3.5300107e-002, 3.4430270e-002, 3.3497746e-002, 3.2505614e-002, 3.1456953e-002, 3.0354843e-002, 2.9202363e-002, 2.8002591e-002, 2.6758606e-002, 2.5473489e-002, 2.4150316e-002, 2.2792168e-002, 2.1402124e-002, 1.9983263e-002, 1.8538663e-002, 1.7071404e-002, 1.5584565e-002, 1.4081224e-002, 1.2564462e-002, 1.1037356e-002, 9.5029859e-003, 7.9644308e-003, 6.4247695e-003, 4.8870812e-003, 3.3544449e-003, 1.8299396e-003, 3.1664424e-004, -1.1823620e-003, -2.6640003e-003, -4.1251914e-003, -5.5628564e-003, -6.9739162e-003, -8.3552918e-003, -6.8938447e-001, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004, -6.1023256e-004] -PermGroRte_d_ann += 31*[PermGroRte_d_ann[-1]] # Add 31 years of the same permanent income growth rate to the end of the sequence -PermGroRte_h_ann += 31*[PermGroRte_h_ann[-1]] -PermGroRte_c_ann += 31*[PermGroRte_c_ann[-1]] -PermGroRte_d_retire = PermGroRte_d_ann[40] # Store the big shock to permanent income at retirement -PermGroRte_h_retire = PermGroRte_h_ann[40] -PermGroRte_c_retire = PermGroRte_c_ann[40] -PermGroRte_d_ann[40] = PermGroRte_d_ann[39] # Overwrite the "retirement drop" with the adjacent growth rate -PermGroRte_h_ann[40] = PermGroRte_h_ann[39] -PermGroRte_c_ann[40] = PermGroRte_c_ann[39] -PermGroFac_d = [] -PermGroFac_h = [] -PermGroFac_c = [] -for j in range(len(PermGroRte_d_ann)): # Make sequences of quarterly permanent income growth factors from annual permanent income growth rates - PermGroFac_d += 4*[(1 + PermGroRte_d_ann[j])**0.25] - PermGroFac_h += 4*[(1 + PermGroRte_h_ann[j])**0.25] - PermGroFac_c += 4*[(1 + PermGroRte_c_ann[j])**0.25] -PermGroFac_d[working_T-1] = 1 + PermGroRte_d_retire # Put the big shock at retirement back into the sequence -PermGroFac_h[working_T-1] = 1 + PermGroRte_h_retire -PermGroFac_c[working_T-1] = 1 + PermGroRte_c_retire -SkillRot = 0.00125 -for t in range(T_cycle): - PermGroFac_d[t] = PermGroFac_d[t]*np.ones(6) - PermGroFac_d[t][1:3] = PermGroFac_d[t][0] - SkillRot - PermGroFac_d[t][3:6] = PermGroFac_d[t][0] - SkillRot - PermGroFac_h[t] = PermGroFac_h[t]*np.ones(6) - PermGroFac_h[t][1:3] = PermGroFac_h[t][0] - SkillRot - PermGroFac_h[t][3:6] = PermGroFac_h[t][0] - SkillRot - PermGroFac_c[t] = PermGroFac_c[t]*np.ones(6) - PermGroFac_c[t][1:3] = PermGroFac_c[t][0] - SkillRot - PermGroFac_c[t][3:6] = PermGroFac_c[t][0] - SkillRot -PermGroFac_d_small = [PermGroFac_d[t][:2] for t in range(T_cycle)] -PermGroFac_h_small = [PermGroFac_h[t][:2] for t in range(T_cycle)] -PermGroFac_c_small = [PermGroFac_c[t][:2] for t in range(T_cycle)] - -# Define the paths of permanent and transitory shocks (from Sabelhaus and Song) -TranShkStd = (np.concatenate((np.linspace(0.1,0.12,17), 0.12*np.ones(17), np.linspace(0.12,0.075,61), np.linspace(0.074,0.007,68), np.zeros(retired_T+1)))*4)**0.5 -TranShkStd = np.ndarray.tolist(TranShkStd) -PermShkStd = np.concatenate((((0.00011342*(np.linspace(24,64.75,working_T-1)-47)**2 + 0.01)/(11.0/4.0))**0.5,np.zeros(retired_T+1))) -PermShkStd[123:162] = PermShkStd[122] # Don't extrapolate -PermShkStd = np.ndarray.tolist(PermShkStd) - -# Import survival probabilities from SSA data -data_location = os.path.dirname(os.path.abspath(__file__)) -f = open(data_location + '/' + 'USactuarial.txt','r') -actuarial_reader = csv.reader(f,delimiter='\t') -raw_actuarial = list(actuarial_reader) -base_death_probs = [] -for j in range(len(raw_actuarial)): - base_death_probs += [float(raw_actuarial[j][4])] # This effectively assumes that everyone is a white woman -f.close() - -# Import adjustments for education and apply them to the base mortality rates -f = open(data_location + '/' + 'EducMortAdj.txt','r') -adjustment_reader = csv.reader(f,delimiter=' ') -raw_adjustments = list(adjustment_reader) -d_death_probs = [] -h_death_probs = [] -c_death_probs = [] -for j in range(76): - d_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[j][1])] - h_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[j][2])] - c_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[j][3])] -for j in range(76,96): - d_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[75][1])] - h_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[75][2])] - c_death_probs += [base_death_probs[j + init_age]*float(raw_adjustments[75][3])] -LivPrb_d = [] -LivPrb_h = [] -LivPrb_c = [] -for j in range(len(d_death_probs)): # Convert annual mortality rates to quarterly survival rates - LivPrb_d += 4*[(1 - d_death_probs[j])**0.25] - LivPrb_h += 4*[(1 - h_death_probs[j])**0.25] - LivPrb_c += 4*[(1 - c_death_probs[j])**0.25] -LivPrb_d_small = [LivPrb_d[t]*np.ones(2) for t in range(T_cycle)] -LivPrb_h_small = [LivPrb_h[t]*np.ones(2) for t in range(T_cycle)] -LivPrb_c_small = [LivPrb_c[t]*np.ones(2) for t in range(T_cycle)] -LivPrb_d = [LivPrb_d[t]*np.ones(6) for t in range(T_cycle)] -LivPrb_h = [LivPrb_h[t]*np.ones(6) for t in range(T_cycle)] -LivPrb_c = [LivPrb_c[t]*np.ones(6) for t in range(T_cycle)] - - -def makeMrkvArray(Urate, Uspell, Dspell, Lspell, Dexit=0): - ''' - Make an age-varying list of Markov transition matrices given the unemployment - rate (in normal times), the average length of an unemployment spell, the average - length of a deep unemployment spell, and the average length of a lockdown. - Parameters - ---------- - Urate: float - Erogodic unemployment rate - Uspell: float - Expected length of unemployment spell - Dspell: float - Expected length of deep unemployment spell - Lspell: float - Expected length of pandemic - Dexit: float - How likely to leave deep unemployment when pandemic ends - 0 => No correlation between exiting deep unemployment and pandemic ending - 1 => Certain to exit deep unemployment when pandemic ends - ''' - U_persist = 1.-1./Uspell - E_persist = 1.-Urate*(1.-U_persist)/(1.-Urate) - D_persist = 1.-1./Dspell - u = U_persist - e = E_persist - d = D_persist - r = 1-1./Lspell - dr = r*(1-Dexit) - - MrkvArray_working = np.array([[e, 1-e, 0.0, 0.0, 0.0, 0.0], - [1-u, u, 0.0, 0.0, 0.0, 0.0], - [0.0, 1-d, d, 0.0, 0.0, 0.0], - [e*(1-r), (1-e)*(1-r), 0.0, e*r, (1-e)*r, 0.0], - [(1-u)*(1-r), u*(1-r), 0.0, (1-u)*r, u*r, 0.0], - [0.0, (1-d)*(1-dr), d*(1-dr), 0.0, (1-d)*dr, d*dr] - ]) - MrkvArray_retired = np.array([[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [1.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [1.0, 0.0, 0.0, 0.0, 0.0, 0.0], - [(1-r), 0.0, 0.0, r, 0.0, 0.0], - [(1-r), 0.0, 0.0, r, 0.0, 0.0], - [(1-r), 0.0, 0.0, r, 0.0, 0.0] - ]) - MrkvArray = (working_T-1)*[MrkvArray_working] + (retired_T+1)*[MrkvArray_retired] - return MrkvArray - - -# Make Markov transition arrays among discrete states in each period of the lifecycle (ACTUAL / SIMULATION) -MrkvArray_real = makeMrkvArray(Urate_normal, Uspell, Dspell_real, Lspell_real) - -# Make Markov transition arrays among discrete states in each period of the lifecycle (PERCEIVED / SIMULATION) -MrkvArray_pcvd = makeMrkvArray(Urate_normal, Uspell, Dspell_pcvd, Lspell_pcvd) - -# Make a two state Markov array ("small") that is only used when generating the initial distribution of states -U_logit_base = np.log(1./(1.-Urate_normal)-1.) -U_persist = 1.-1./Uspell -E_persist = 1.-Urate_normal*(1.-U_persist)/(1.-Urate_normal) -u = U_persist -e = E_persist -MrkvArray_working_small = np.array([[e, 1-e], - [1-u, u] - ]) -MrkvArray_retired_small = np.array([[1., 0.], - [1., 0.] - ]) -MrkvArray_small = (working_T-1)*[MrkvArray_working_small] + (retired_T+1)*[MrkvArray_retired_small] - -# Define a parameter dictionaries for three education levels -init_dropout = {"cycles" : 1, - "T_cycle": T_cycle, - "T_retire": working_T-1, - 'T_sim': T_cycle, - 'T_age': T_cycle+1, - 'AgentCount': 10000, - "PermGroFacAgg": PermGroFacAgg, - "PopGroFac": PopGroFac, - "CRRA": CRRA, - "DiscFac": 0.98, # This will be overwritten at type construction - "Rfree_big" : np.array(6*[1.01]), - "PermGroFac_big": PermGroFac_d, - "LivPrb_big": LivPrb_d, - "uPfac_big" : np.array(3*[1.0] + 3*[uPfac_L]), - "MrkvArray_big" : MrkvArray_pcvd, - "Rfree" : np.array(2*[1.01]), - "PermGroFac": PermGroFac_d_small, - "LivPrb": LivPrb_d_small, - "uPfac" : np.array(2*[1.0]), - "MrkvArray" : MrkvArray_small, # Yes, this is intentional - "MrkvArray_pcvd" : MrkvArray_small, # Yes, this is intentional - "MrkvArray_sim" : MrkvArray_real, - "BoroCnstArt": 0.0, - "PermShkStd": PermShkStd, - "PermShkCount": PermShkCount, - "TranShkStd": TranShkStd, - "TranShkCount": TranShkCount, - "UnempPrb": 0.0, # Unemployment is modeled as a Markov state - "UnempPrbRet": UnempPrbRet, - "IncUnemp": IncUnemp, - "IncUnempRet": IncUnempRet, - "aXtraMin": aXtraMin, - "aXtraMax": aXtraMax, - "aXtraCount": aXtraCount, - "aXtraExtra": aXtraExtra, - "aXtraNestFac": aXtraNestFac, - "CubicBool": False, - "vFuncBool": False, - 'aNrmInitMean': np.log(0.00001), # Initial assets are zero - 'aNrmInitStd': 0.0, - 'pLvlInitMean': pLvlInitMeanD, - 'pLvlInitStd': pLvlInitStd, - "MrkvPrbsInit" : np.array([1-Urate_normal, Urate_normal] + 4*[0.0]), - 'Urate' : Urate_normal, - 'Uspell' : Uspell, - 'L_shared' : L_shared, - 'Lspell_pcvd' : Lspell_pcvd, - 'Lspell_real' : Lspell_real, - 'Dspell_pcvd' : Dspell_pcvd, - 'Dspell_real' : Dspell_real, - 'EducType': 0, - 'UpdatePrb': UpdatePrb, - 'track_vars' : [] - } -if L_shared: - init_dropout['T_lockdown'] = int(Lspell_real) - -adj_highschool = { - "PermGroFac" : PermGroFac_h_small, - "LivPrb" : LivPrb_h_small, - "PermGroFac_big" : PermGroFac_h, - "LivPrb_big" : LivPrb_h, - 'pLvlInitMean' : pLvlInitMeanH, - 'EducType' : 1} -init_highschool = init_dropout.copy() -init_highschool.update(adj_highschool) - -adj_college = { - "PermGroFac" : PermGroFac_c_small, - "LivPrb" : LivPrb_c_small, - "PermGroFac_big" : PermGroFac_c, - "LivPrb_big" : LivPrb_c, - 'pLvlInitMean' : pLvlInitMeanC, - 'EducType' : 2} -init_college = init_dropout.copy() -init_college.update(adj_college) - -# Define a dictionary to represent the baseline scenario -base_dict = { - 'PanShock' : False, - 'StimMax' : 0., - 'StimCut0' : None, - 'StimCut1' : None, - 'BonusUnemp' : 0., - 'BonusDeep' : 0., - 'T_ahead' : 0, - 'UnempD' : U_logit_base, - 'UnempH' : U_logit_base, - 'UnempC' : U_logit_base, - 'UnempP' : 0., - 'UnempA1' : 0., - 'UnempA2' : 0., - 'DeepD' : -np.inf, - 'DeepH' : -np.inf, - 'DeepC' : -np.inf, - 'DeepP' : 0., - 'DeepA1' : 0., - 'DeepA2' : 0., - 'Dspell_pcvd' : Dspell_pcvd, # These five parameters don't do anything in baseline scenario - 'Dspell_real' : Dspell_real, - 'Lspell_pcvd' : Lspell_pcvd, - 'Lspell_real' : Lspell_real, - 'L_shared' : L_shared - } - -# Define a dictionary to mutate baseline for the pandemic -pandemic_changes = { - 'PanShock' : True, - 'UnempD' : UnempD + Unemp0, - 'UnempH' : UnempH + Unemp0, - 'UnempC' : UnempC + Unemp0, - 'UnempP' : UnempP, - 'UnempA1' : UnempA1, - 'UnempA2' : UnempA2, - 'DeepD' : DeepD + Deep0, - 'DeepH' : DeepH + Deep0, - 'DeepC' : DeepC + Deep0, - 'DeepP' : DeepP, - 'DeepA1' : DeepA1, - 'DeepA2' : DeepA2, - } - -# Define a dictionary to mutate baseline for the fiscal stimulus -stimulus_changes = { - 'StimMax' : StimMax, - 'StimCut0' : StimCut0, - 'StimCut1' : StimCut1, - 'BonusUnemp' : BonusUnemp, - 'BonusDeep' : BonusDeep, - 'T_ahead' : T_ahead, - } - -# Define a dictionary to mutate baseline for a deep unemployment pandemic -deep_pandemic_changes = { - 'PanShock' : True, - 'UnempD' : UnempD + DeepPanAdj1, - 'UnempH' : UnempH + DeepPanAdj1, - 'UnempC' : UnempC + DeepPanAdj1, - 'UnempP' : UnempP + DeepPanAdj3, - 'UnempA1' : UnempA1, - 'UnempA2' : UnempA2, - 'DeepD' : DeepD + DeepPanAdj2, - 'DeepH' : DeepH + DeepPanAdj2, - 'DeepC' : DeepC + DeepPanAdj2, - 'DeepP' : DeepP, - 'DeepA1' : DeepA1, - 'DeepA2' : DeepA2, - } - diff --git a/HARK/Calibration/SolvingMicroDSOPs.py b/HARK/Calibration/SolvingMicroDSOPs.py deleted file mode 100644 index c684145fb..000000000 --- a/HARK/Calibration/SolvingMicroDSOPs.py +++ /dev/null @@ -1,173 +0,0 @@ -''' -Specifies the full set of calibrated values required to estimate the SolvingMicroDSOPs -model. The empirical data is stored in a separate csv file and is loaded in SetupSCFdata. -''' -from __future__ import print_function -import numpy as np -# --------------------------------------------------------------------------------- -# Debugging flags -# --------------------------------------------------------------------------------- - -show_PermGroFacAgg_error = False -# Error Notes: -# This sets a "quick fix" to the error, AttributeError: 'TempConsumerType' object has no attribute 'PermGroFacAgg' -# If you set this flag to "True" you will see the error. A more thorough fix is to -# fix the place where this error was introduced (Set to "True" and it will appear; -# this was almost certainly introduced when the code was extended to be used in the -# GE setting). An even more thorough solution, which moves beyond the scope of -# fixing this error, is adding unit tests to ID when changes to some code will -# break things elsewhere. -# Note: alternatively, decide that the "init_consumer_objects['PermGroFacAgg'] = 1.0" -# line below fixes it properly ('feature not a bug') and remove all this text. - -# --------------------------------------------------------------------------------- -# - 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 = np.array([0.17, 0.5, 0.83]) # Three point discrete distribution of initial w -initial_wealth_income_ratio_probs = np.array([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 - -# ----------------------------------------------------------------------------- -# -- Set up the dictionary "container" for making a basic lifecycle type ------ -# ----------------------------------------------------------------------------- - -# 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_rate':0.0, - 'vFuncBool':vFuncBool, - 'CubicBool':CubicBool - } - -if show_PermGroFacAgg_error: - pass # do nothing -else: - print("***NOTE: using a 'quick fix' for an attribute error. See 'Error Notes' in EstimationParameter.py for further discussion.***") - init_consumer_objects['PermGroFacAgg'] = 1.0 - - -if __name__ == '__main__': - print("Sorry, EstimationParameters doesn't actually do anything on its own.") - print("This module is imported by StructEstimation, providing calibrated ") - print("parameters for the example estimation. Please see that module if you ") - print("want more interesting output.") From 286f8d6da16e5d362ad39db34b6ab763d0ad3296 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 16:11:27 -0500 Subject: [PATCH 44/47] Rm test notebook --- examples/Calibration/CalibTest.ipynb | 55 ---------------------------- 1 file changed, 55 deletions(-) delete mode 100644 examples/Calibration/CalibTest.ipynb diff --git a/examples/Calibration/CalibTest.ipynb b/examples/Calibration/CalibTest.ipynb deleted file mode 100644 index 25af34153..000000000 --- a/examples/Calibration/CalibTest.ipynb +++ /dev/null @@ -1,55 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType" - ] - } - ], - "metadata": { - "jupytext": { - "formats": "ipynb,py:percent" - }, - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.9" - }, - "latex_envs": { - "LaTeX_envs_menu_present": true, - "autoclose": false, - "autocomplete": true, - "bibliofile": "biblio.bib", - "cite_by": "apalike", - "current_citInitial": 1, - "eqLabelWithNumbers": true, - "eqNumInitial": 1, - "hotkeys": { - "equation": "Ctrl-E", - "itemize": "Ctrl-I" - }, - "labels_anchors": false, - "latex_user_defs": false, - "report_style_numbering": false, - "user_envs_cfg": false - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From d42f1fbfd6a7bb9fdcd5cf8f7919394ea5ff5d40 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Fri, 5 Feb 2021 16:13:52 -0500 Subject: [PATCH 45/47] Comment --- HARK/Calibration/Income/IncomeTools.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/HARK/Calibration/Income/IncomeTools.py b/HARK/Calibration/Income/IncomeTools.py index 9aa22d604..d89455818 100644 --- a/HARK/Calibration/Income/IncomeTools.py +++ b/HARK/Calibration/Income/IncomeTools.py @@ -671,7 +671,8 @@ def ParseIncomeSpec( ) else: - + + # Placeholder for future ways of specifying volatilities raise NotImplementedError() income_params["PermShkStd"] = PermShkStd From 7edb8451d1009d12324a810db12238d79cc5e677 Mon Sep 17 00:00:00 2001 From: Mateo Velasquez-Giraldo Date: Mon, 8 Feb 2021 10:24:42 -0500 Subject: [PATCH 46/47] Update example --- examples/Calibration/Life_Cycle_example.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/Calibration/Life_Cycle_example.py b/examples/Calibration/Life_Cycle_example.py index b8c3c2662..adab6fcf2 100644 --- a/examples/Calibration/Life_Cycle_example.py +++ b/examples/Calibration/Life_Cycle_example.py @@ -79,7 +79,7 @@ Agent.T_sim = 200 # Set up the variables we want to keep track of. -Agent.track_vars = ["aNrmNow", "cNrmNow", "pLvlNow", "t_age", "mNrmNow"] +Agent.track_vars = ["aNrm", "cNrm", "pLvl", "t_age", "mNrm"] # Run the simulations Agent.initializeSim() @@ -89,9 +89,9 @@ raw_data = { "Age": Agent.history["t_age"].flatten() + birth_age - 1, - "pIncome": Agent.history["pLvlNow"].flatten(), - "nrmM": Agent.history["mNrmNow"].flatten(), - "nrmC": Agent.history["cNrmNow"].flatten(), + "pIncome": Agent.history["pLvl"].flatten(), + "nrmM": Agent.history["mNrm"].flatten(), + "nrmC": Agent.history["cNrm"].flatten(), } Data = pd.DataFrame(raw_data) From c429270ea5c40f93d04ebfa679d34ddaf9315631 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateo=20Vel=C3=A1squez-Giraldo?= Date: Thu, 11 Feb 2021 10:04:54 -0500 Subject: [PATCH 47/47] Update CHANGELOG.md --- Documentation/CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index 6d625370e..17999c3f4 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -27,7 +27,7 @@ interpolation for problems with CRRA utility. See [#888](https://github.com/econ * Adds a module for extracting initial distributions of permanent income (`pLvl`) and normalized assets (`aNrm`) from the SCF [#932](https://github.com/econ-ark/HARK/pull/932). * Fix the return fields of `dcegm/calcCrossPoints`[#909](https://github.com/econ-ark/HARK/pull/909). * Corrects location of constructor documentation to class string for Sphinx rendering [#908](https://github.com/econ-ark/HARK/pull/908) -* Adds a module for producing life-cycle profiles of income shock variances from [Sabelhaus and Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358). See [#921](https://github.com/econ-ark/HARK/pull/921). +* Adds a module with tools for parsing and using various income calibrations from the literature. It includes the option of using life-cycle profiles of income shock variances from [Sabelhaus and Song (2010)](https://www.sciencedirect.com/science/article/abs/pii/S0304393210000358). See [#921](https://github.com/econ-ark/HARK/pull/921), [#941](https://github.com/econ-ark/HARK/pull/941). * remove "Now" from model variable names [#936](https://github.com/econ-ark/HARK/pull/936) * Moves state MrkvNow to shocks['Mrkv'] in AggShockMarkov and KrusellSmith models [#935](https://github.com/econ-ark/HARK/pull/935)