diff --git a/code/do_all.py b/code/do_all.py index 2b20c52..981a32d 100644 --- a/code/do_all.py +++ b/code/do_all.py @@ -46,7 +46,7 @@ still run. """ -from estimark.calibration.options import ( +from estimark.options import ( all_replications, high_resource, low_resource, @@ -137,11 +137,11 @@ def run_replication(): return if subjective_markets == "2" or subjective_markets == "4": - replication_specs["subjective_stock_market"] = True + replication_specs["subjective_stock"] = True print("Adding subjective stock market beliefs...") if subjective_markets == "3" or subjective_markets == "4": - replication_specs["subjective_labor_market"] = True + replication_specs["subjective_labor"] = True print("Adding subjective labor market beliefs...") estimate(**replication_specs) diff --git a/code/do_figs.py b/code/do_figs.py deleted file mode 100644 index d2f124a..0000000 --- a/code/do_figs.py +++ /dev/null @@ -1,4 +0,0 @@ -from estimark.estimation import estimate_all - -if __name__ == "__main__": - estimate_all() diff --git a/code/estimark/calibration/estimation_parameters.py b/code/estimark/calibration/parameters.py similarity index 63% rename from code/estimark/calibration/estimation_parameters.py rename to code/estimark/calibration/parameters.py index 7ebaf25..00ff394 100644 --- a/code/estimark/calibration/estimation_parameters.py +++ b/code/estimark/calibration/parameters.py @@ -3,27 +3,12 @@ model. The empirical data is stored in a separate csv file and is loaded in setup_scf_data. """ +from pathlib import Path + import numpy as np from HARK.Calibration.Income.IncomeTools import CGM_income, parse_income_spec -from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table -from pathlib import Path from HARK.ConsumptionSaving.ConsIndShockModel import init_lifecycle - -# --------------------------------------------------------------------------------- -# 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. +from HARK.datasets.life_tables.us_ssa.SSATools import parse_ssa_life_table # --------------------------------------------------------------------------------- # - Define all of the model parameters for EstimatingMicroDSOPs and ConsumerExamples - @@ -32,9 +17,7 @@ exp_nest = 1 # Number of times to "exponentially nest" when constructing a_grid aXtraMin = 0.001 # Minimum end-of-period "assets above minimum" value aXtraMax = 100 # 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 = 200 # Number of points in the grid of "assets above minimum" +aXtraCount = 100 # Number of points in the grid of "assets above minimum" # Artificial borrowing constraint; imposed minimum level of end-of period assets BoroCnstArt = 0.0 @@ -55,17 +38,17 @@ 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 +terminal_t = final_age - initial_age # Total number of periods in the model retirement_t = retirement_age - initial_age - 1 # Initial guess of the coefficient of relative risk aversion during estimation (rho) -CRRA_start = 5.0 +init_CRRA = 5.0 # Initial guess of the adjustment to the discount factor during estimation (beth) -DiscFacAdj_start = 0.99 +init_DiscFacAdj = 0.99 # Bounds for beth; if violated, objective function returns "penalty value" -DiscFacAdj_bound = [0.0001, 15.0] +bounds_DiscFacAdj = [0.0001, 15.0] # Bounds for rho; if violated, objective function returns "penalty value" -CRRA_bound = [0.0001, 15.0] +bounds_CRRA = [0.0001, 15.0] # Income ss_variances = True @@ -73,47 +56,46 @@ # Replace retirement age income_spec["age_ret"] = retirement_age inc_calib = parse_income_spec( - age_min=initial_age, age_max=final_age, **income_spec, SabelhausSong=ss_variances + age_min=initial_age, + age_max=final_age, + **income_spec, + SabelhausSong=ss_variances, ) # Age-varying discount factors over the lifecycle, lifted from Cagetti (2003) - - # Get the directory containing the current file and construct the full path to the CSV file csv_file_path = Path(__file__).resolve().parent / ".." / "data" / "Cagetti2003.csv" -DiscFac_timevary = np.genfromtxt(csv_file_path) * 0.0 + 1.0 +timevary_DiscFac = np.genfromtxt(csv_file_path) * 0.0 + 1.0 # todo +constant_DiscFac = np.ones_like(timevary_DiscFac) # Survival probabilities over the lifecycle liv_prb = parse_ssa_life_table( - female=False, min_age=initial_age, max_age=final_age - 1, cohort=1960 + female=False, + min_age=initial_age, + max_age=final_age - 1, + cohort=1960, ) -# 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 = [ - [age for age in range(start, start + 5)] for start in range(26, 61, 5) -] - # Three point discrete distribution of initial w -initial_wealth_income_ratio_vals = np.array([0.17, 0.5, 0.83]) +init_w_to_y = np.array([0.17, 0.5, 0.83]) # Equiprobable discrete distribution of initial w -initial_wealth_income_ratio_probs = np.array([0.33333, 0.33333, 0.33334]) +prob_w_to_y = np.array([0.33333, 0.33333, 0.33334]) 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 options = { - "initial_wealth_income_ratio_vals": initial_wealth_income_ratio_vals, - "initial_wealth_income_ratio_probs": initial_wealth_income_ratio_probs, + "init_w_to_y": init_w_to_y, + "prob_w_to_y": prob_w_to_y, "num_agents": num_agents, "bootstrap_size": bootstrap_size, "seed": seed, - "DiscFacAdj_start": DiscFacAdj_start, - "CRRA_start": CRRA_start, - "DiscFacAdj_bound": DiscFacAdj_bound, - "CRRA_bound": CRRA_bound, - "DiscFac_timevary": DiscFac_timevary, + "init_DiscFacAdj": init_DiscFacAdj, + "init_CRRA": init_CRRA, + "bounds_DiscFacAdj": bounds_DiscFacAdj, + "bounds_CRRA": bounds_CRRA, + "timevary_DiscFac": timevary_DiscFac, } # ----------------------------------------------------------------------------- @@ -124,28 +106,28 @@ init_consumer_objects = { **init_lifecycle, **{ - "CRRA": CRRA_start, + "CRRA": init_CRRA, "Rfree": Rfree, "PermGroFac": inc_calib["PermGroFac"], + "PermGroFacAgg": 1.0, "BoroCnstArt": BoroCnstArt, "PermShkStd": inc_calib["PermShkStd"], "PermShkCount": PermShkCount, "TranShkStd": inc_calib["TranShkStd"], "TranShkCount": TranShkCount, - "T_cycle": TT, + "T_cycle": terminal_t, "UnempPrb": UnempPrb, "UnempPrbRet": UnempPrbRet, "T_retire": retirement_t, - "T_age": TT, + "T_age": terminal_t, "IncUnemp": IncUnemp, "IncUnempRet": IncUnempRet, "aXtraMin": aXtraMin, "aXtraMax": aXtraMax, "aXtraCount": aXtraCount, - "aXtraExtra": [aXtraExtra, aXtraHuge], "aXtraNestFac": exp_nest, "LivPrb": liv_prb, - "DiscFac": DiscFac_timevary, + "DiscFac": timevary_DiscFac, "AgentCount": num_agents, "seed": seed, "tax_rate": 0.0, @@ -158,24 +140,22 @@ ElnR = 0.020 VlnR = 0.424**2 -init_subjective_stock_market = { +TrueElnR = 0.085 +TrueVlnR = 0.170**2 + +init_subjective_stock = { "Rfree": 1.019, # from Mateo's JMP "RiskyAvg": np.exp(ElnR + 0.5 * VlnR), "RiskyStd": np.sqrt(np.exp(2 * ElnR + VlnR) * (np.exp(VlnR) - 1)), + "RiskyAvgTrue": np.exp(TrueElnR + 0.5 * TrueVlnR), + "RiskyStdTrue": np.sqrt(np.exp(2 * TrueElnR + TrueVlnR) * (np.exp(TrueVlnR) - 1)), } -init_subjective_labor_market = { # from Tao's JMP +init_subjective_labor = { # from Tao's JMP "TranShkStd": [0.03] * len(inc_calib["TranShkStd"]), "PermShkStd": [0.03] * len(inc_calib["PermShkStd"]), } -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, estimation_parameters doesn't actually do anything on its own.") diff --git a/code/estimark/calibration/setup_scf_data.py b/code/estimark/calibration/scf.py similarity index 50% rename from code/estimark/calibration/setup_scf_data.py rename to code/estimark/calibration/scf.py index 9cab37c..1ef2ba0 100644 --- a/code/estimark/calibration/setup_scf_data.py +++ b/code/estimark/calibration/scf.py @@ -2,54 +2,53 @@ Sets up the SCF data for use in the EstimatingMicroDSOPs estimation. """ -from estimark.calibration.estimation_parameters import ( - empirical_cohort_age_groups, - initial_age, -) +from pathlib import Path import numpy as np # Numerical Python import pandas as pd - - -from pathlib import Path +from estimark.calibration.parameters import initial_age # Get the directory containing the current file and construct the full path to the CSV file csv_file_path = Path(__file__).resolve().parent / ".." / "data" / "SCFdata.csv" # Read the CSV file -scf_data = pd.read_csv(csv_file_path) +scf_full_data = pd.read_csv(csv_file_path) # Keep only observations with normal incomes > 0, # otherwise wealth/income is not well defined -scf_data = scf_data[scf_data.norminc > 0.0] +scf_full_data = scf_full_data[scf_full_data.norminc > 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 +age_groups = [list(range(start, start + 5)) for start in range(26, 61, 5)] # Initialize empty lists for the data -w_to_y_data = [] # Ratio of wealth to permanent income -empirical_weights = [] # Weighting for this observation -empirical_groups = [] # Which age group this observation belongs to (1-7) +scf_data = [] # Ratio of wealth to permanent income +scf_weights = [] # Weighting for this observation +scf_groups = [] # Which age group this observation belongs to (1-7) # Only extract the data required from the SCF dataset -search_ages = [ages[-1] for ages in empirical_cohort_age_groups] +search_ages = [ages[-1] for ages in age_groups] for idx, age in enumerate(search_ages): - age_data = scf_data[scf_data.age_group.str.contains(f"{age}]")] - w_to_y_data.append(age_data.wealth_income_ratio.values) - empirical_weights.append(age_data.weight.values) + age_data = scf_full_data[scf_full_data.age_group.str.contains(f"{age}]")] + scf_data.append(age_data.wealth_income_ratio.values) + scf_weights.append(age_data.weight.values) # create a group id 1-7 for every observation - empirical_groups.append([idx + 1] * len(age_data)) + scf_groups.append([idx + 1] * len(age_data)) # Convert SCF data to numpy's array format for easier math -w_to_y_data = np.concatenate(w_to_y_data) -empirical_weights = np.concatenate(empirical_weights) -empirical_groups = np.concatenate(empirical_groups) +scf_data = np.concatenate(scf_data) +scf_weights = np.concatenate(scf_weights) +scf_groups = np.concatenate(scf_groups) # Generate a single array of SCF data, useful for resampling for bootstrap -scf_data_array = np.array([w_to_y_data, empirical_groups, empirical_weights]).T +scf_array = np.array([scf_data, scf_groups, scf_weights]).T # Generate a mapping between the real ages in the groups and the indices of simulated data -simulation_map_cohorts_to_age_indices = [] -for ages in empirical_cohort_age_groups: - simulation_map_cohorts_to_age_indices.append(np.array(ages) - initial_age) +scf_mapping = [] +for ages in age_groups: + scf_mapping.append(np.array(ages) - initial_age) if __name__ == "__main__": print("Sorry, setup_scf_data doesn't actually do anything on its own.") diff --git a/code/estimark/estimation.py b/code/estimark/estimation.py index 47696a5..4f6189a 100644 --- a/code/estimark/estimation.py +++ b/code/estimark/estimation.py @@ -9,25 +9,9 @@ income as defined in ConsIndShockModel. """ -# Parameters for the consumer type and the estimation -from estimark.calibration.estimation_parameters import ( - init_consumer_objects, - init_subjective_labor_market, - init_subjective_stock_market, - options, -) -import estimark.calibration.setup_scf_data as scf_data # SCF 2004 data on household wealth import csv -from estimark.agents import ( - BequestWarmGlowLifeCycleConsumerType, - BequestWarmGlowLifeCyclePortfolioType, - IndShkLifeCycleConsumerType, - PortfolioLifeCycleConsumerType, - WealthPortfolioLifeCycleConsumerType, -) -from time import time # Timing utility - from pathlib import Path +from time import time # Timing utility import matplotlib.pyplot as plt import numpy as np # Numeric Python @@ -41,17 +25,37 @@ # Import modules from core HARK libraries: # The consumption-saving micro model +from estimark.agents import ( + BequestWarmGlowLifeCycleConsumerType, + BequestWarmGlowLifeCyclePortfolioType, + IndShkLifeCycleConsumerType, + PortfolioLifeCycleConsumerType, + WealthPortfolioLifeCycleConsumerType, +) + +# Parameters for the consumer type and the estimation +from estimark.calibration.parameters import ( + init_consumer_objects, + init_subjective_labor, + init_subjective_stock, + options, +) +# SCF 2004 data on household wealth +from estimark.calibration.scf import ( + age_groups, + scf_array, + scf_data, + scf_groups, + scf_mapping, + scf_weights, +) # Pathnames to the other files: # Relative directory for primitive parameter files -calibration_dir = "code/calibration/" -# Relative directory for primitive parameter files tables_dir = "code/tables/" # Relative directory for primitive parameter files figures_dir = "content/figures/" -# Relative directory for primitive parameter files -code_dir = "code/" # Set booleans to determine which tasks should be done @@ -65,8 +69,8 @@ # Whether to make a contour map of the objective function local_make_contour_plot = True # Whether to use subjective beliefs -local_subjective_stock_market = False -local_subjective_labor_market = False +local_subjective_stock = False +local_subjective_labor = False # ===================================================== @@ -76,8 +80,8 @@ def make_estimation_agent( agent_name=local_agent_name, - subjective_stock_market=local_subjective_stock_market, - subjective_labor_market=local_subjective_labor_market, + subjective_stock=local_subjective_stock, + subjective_labor=local_subjective_labor, ): if agent_name == "IndShock": agent_type = IndShkLifeCycleConsumerType @@ -90,14 +94,14 @@ def make_estimation_agent( elif agent_name == "WealthPortfolio": agent_type = WealthPortfolioLifeCycleConsumerType - if subjective_stock_market or subjective_labor_market: + if subjective_stock or subjective_labor: agent_name += "Sub" - if subjective_stock_market: + if subjective_stock: agent_name += "(Stock)" - init_consumer_objects.update(init_subjective_stock_market) - if subjective_labor_market: + init_consumer_objects.update(init_subjective_stock) + if subjective_labor: agent_name += "(Labor)" - init_consumer_objects.update(init_subjective_labor_market) + init_consumer_objects.update(init_subjective_labor) agent_name += "Market" # Make a lifecycle consumer to be used for estimation, including simulated @@ -110,8 +114,8 @@ def make_estimation_agent( estimation_agent.track_vars = ["bNrm"] # Draw initial assets for each consumer estimation_agent.aNrmInit = DiscreteDistribution( - options["initial_wealth_income_ratio_probs"], - options["initial_wealth_income_ratio_vals"], + options["prob_w_to_y"], + options["init_w_to_y"], seed=options["seed"], ).draw(N=options["num_agents"]) estimation_agent.make_shock_history() @@ -132,23 +136,18 @@ def weighted_median(values, weights): return median -def get_targeted_moments( - empirical_data=scf_data.w_to_y_data, - empirical_weights=scf_data.empirical_weights, - empirical_groups=scf_data.empirical_groups, - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, +def get_target_moments( + data=scf_data, weights=scf_weights, groups=scf_groups, mapping=scf_mapping ): # Initialize - group_count = len(map_simulated_to_empirical_cohorts) - tgt_moments = np.zeros(group_count) + group_count = len(mapping) + target_moments = np.zeros(group_count) for g in range(group_count): - group_indices = empirical_groups == (g + 1) # groups are numbered from 1 - tgt_moments[g] = weighted_median( - empirical_data[group_indices], empirical_weights[group_indices] - ) + group_indices = groups == (g + 1) # groups are numbered from 1 + target_moments[g] = weighted_median(data[group_indices], weights[group_indices]) - return tgt_moments + return target_moments def get_initial_guess(agent_name): @@ -165,7 +164,7 @@ def get_initial_guess(agent_name): with open(csv_file_path, "r") as file: initial_guess = np.loadtxt(file, skiprows=1, delimiter=",") except FileNotFoundError: - initial_guess = [options["DiscFacAdj_start"], options["CRRA_start"]] + initial_guess = [options["init_DiscFacAdj"], options["init_CRRA"]] return initial_guess @@ -175,9 +174,9 @@ def simulate_moments( DiscFacAdj, CRRA, agent, - DiscFacAdj_bound=options["DiscFacAdj_bound"], - CRRA_bound=options["CRRA_bound"], - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, + bounds_DiscFacAdj=options["bounds_DiscFacAdj"], + bounds_CRRA=options["bounds_CRRA"], + mapping=scf_mapping, ): """ A quick check to make sure that the parameter values are within bounds. @@ -185,21 +184,21 @@ def simulate_moments( simulation, so the objective function doesn't even bother with them. """ if ( - DiscFacAdj < DiscFacAdj_bound[0] - or DiscFacAdj > DiscFacAdj_bound[1] - or CRRA < CRRA_bound[0] - or CRRA > CRRA_bound[1] + DiscFacAdj < bounds_DiscFacAdj[0] + or DiscFacAdj > bounds_DiscFacAdj[1] + or CRRA < bounds_CRRA[0] + or CRRA > bounds_CRRA[1] ): - return 1e30 * np.ones(len(map_simulated_to_empirical_cohorts)) + return 1e30 * np.ones(len(mapping)) # Update the agent with a new path of DiscFac based on this DiscFacAdj (and a new CRRA) - agent.DiscFac = [b * DiscFacAdj for b in options["DiscFac_timevary"]] + agent.DiscFac = [b * DiscFacAdj for b in options["timevary_DiscFac"]] agent.CRRA = CRRA # Solve the model for these parameters, then simulate wealth data agent.solve() # Solve the microeconomic model # "Unpack" the consumption function for convenient access # agent.unpack("cFunc") - max_sim_age = max([max(ages) for ages in map_simulated_to_empirical_cohorts]) + 1 + max_sim_age = max([max(ages) for ages in mapping]) + 1 # Initialize the simulation by clearing histories, resetting initial values agent.initialize_sim() agent.simulate(max_sim_age) # Simulate histories of consumption and wealth @@ -207,19 +206,13 @@ def simulate_moments( sim_w_history = agent.history["bNrm"] # Find the distance between empirical data and simulated medians for each age group - group_count = len(map_simulated_to_empirical_cohorts) + group_count = len(mapping) sim_moments = [] for g in range(group_count): # The simulated time indices corresponding to this age group - cohort_indices = map_simulated_to_empirical_cohorts[g] + cohort_indices = mapping[g] # The median of simulated wealth-to-income for this age group - sim_moments += [ - np.median( - sim_w_history[ - cohort_indices, - ] - ) - ] + sim_moments += [np.median(sim_w_history[cohort_indices])] sim_moments = np.array(sim_moments) @@ -230,10 +223,10 @@ def smm_obj_func( DiscFacAdj, CRRA, agent, - tgt_moments, - DiscFacAdj_bound=options["DiscFacAdj_bound"], - CRRA_bound=options["CRRA_bound"], - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, + target_moments, + bounds_DiscFacAdj=options["bounds_DiscFacAdj"], + bounds_CRRA=options["bounds_CRRA"], + mapping=scf_mapping, ): """ The objective function for the SMM estimation. Given values of discount factor @@ -252,16 +245,16 @@ def smm_obj_func( ---------- DiscFacAdj : float An adjustment factor to a given age-varying sequence of discount factors. - I.e. DiscFac[t] = DiscFacAdj*DiscFac_timevary[t]. + I.e. DiscFac[t] = DiscFacAdj*timevary_DiscFac[t]. CRRA : float Coefficient of relative risk aversion. agent : ConsumerType The consumer type to be used in the estimation, with all necessary para- meters defined except the discount factor and CRRA. - DiscFacAdj_bound : (float,float) + bounds_DiscFacAdj : (float,float) Lower and upper bounds on DiscFacAdj; if outside these bounds, the function simply returns a "penalty value". - DiscFacAdj_bound : (float,float) + bounds_DiscFacAdj : (float,float) Lower and upper bounds on CRRA; if outside these bounds, the function simply returns a "penalty value". empirical_data : np.array @@ -270,7 +263,7 @@ def smm_obj_func( Weights for each observation in empirical_data. empirical_groups : np.array Array of integers listing the age group for each observation in empirical_data. - map_simulated_to_empirical_cohorts : [np.array] + mapping : [np.array] List of arrays of "simulation ages" for each age grouping. E.g. if the 0th element is [1,2,3,4,5], then these time indices from the simulation correspond to the 0th empirical age group. @@ -286,11 +279,11 @@ def smm_obj_func( DiscFacAdj, CRRA, agent, - DiscFacAdj_bound, - CRRA_bound, - map_simulated_to_empirical_cohorts, + bounds_DiscFacAdj, + bounds_CRRA, + mapping, ) - errors = tgt_moments - sim_moments + errors = target_moments - sim_moments loss = np.dot(errors, errors) return loss @@ -331,25 +324,17 @@ def calculate_std_err_bootstrap(initial_estimate, N, agent, seed=0, verbose=Fals t_start = time() # Bootstrap a new dataset by resampling from the original data - bootstrap_data = ( - bootstrap_sample_from_data(scf_data.scf_data_array, seed=seed_list[n]) - ).T - w_to_y_data_bootstrap = bootstrap_data[ - 0, - ] - empirical_groups_bootstrap = bootstrap_data[ - 1, - ] - empirical_weights_bootstrap = bootstrap_data[ - 2, - ] + bootstrap_data = (bootstrap_sample_from_data(scf_array, seed=seed_list[n])).T + w_to_y_data_bootstrap = bootstrap_data[0] + empirical_groups_bootstrap = bootstrap_data[1] + empirical_weights_bootstrap = bootstrap_data[2] # Find moments with bootstrapped sample - bstrap_tgt_moments = get_targeted_moments( - empirical_data=w_to_y_data_bootstrap, - empirical_weights=empirical_weights_bootstrap, - empirical_groups=empirical_groups_bootstrap, - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, + bstrap_target_moments = get_target_moments( + data=w_to_y_data_bootstrap, + weights=empirical_weights_bootstrap, + groups=empirical_groups_bootstrap, + mapping=scf_mapping, ) # Make a temporary function for use in this estimation run @@ -358,8 +343,8 @@ def smm_obj_func_bootstrap(parameters_to_estimate): DiscFacAdj=parameters_to_estimate[0], CRRA=parameters_to_estimate[1], agent=agent, - tgt_moments=bstrap_tgt_moments, - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, + target_moments=bstrap_target_moments, + mapping=scf_mapping, ) # Estimate the model with the bootstrap data and add to list of estimates @@ -386,9 +371,7 @@ def smm_obj_func_bootstrap(parameters_to_estimate): # ================================================================= -def estimate_model_nelder_mead( - agent_name, estimation_agent, targeted_moments, initial_guess -): +def estimate_model_opt(agent_name, estimation_agent, target_moments, initial_guess): print("----------------------------------------------------------------------") print( f"Now estimating the model using Nelder-Mead from an initial guess of {initial_guess}..." @@ -406,7 +389,7 @@ def smm_obj_func_redux(parameters_to_estimate): DiscFacAdj=parameters_to_estimate[0], CRRA=parameters_to_estimate[1], agent=estimation_agent, - tgt_moments=targeted_moments, + target_moments=target_moments, ) t_start_estimate = time() @@ -496,14 +479,14 @@ def simulate_moments_redux(x): x[0], x[1], agent=estimation_agent, - DiscFacAdj_bound=options["DiscFacAdj_bound"], - CRRA_bound=options["CRRA_bound"], - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, + bounds_DiscFacAdj=options["bounds_DiscFacAdj"], + bounds_CRRA=options["bounds_CRRA"], + mapping=scf_mapping, ) return moments - n_moments = len(scf_data.simulation_map_cohorts_to_age_indices) + n_moments = len(scf_mapping) jac = np.array( [ approx_fprime( @@ -519,10 +502,7 @@ def simulate_moments_redux(x): sensitivity = np.dot(np.linalg.inv(np.dot(jac.T, jac)), jac.T) # Create lables for moments in the plots - moment_labels = [ - "[" + str(min(x)) + "," + str(max(x)) + "]" - for x in scf_data.empirical_cohort_age_groups - ] + moment_labels = ["[" + str(min(x)) + "," + str(max(x)) + "]" for x in age_groups] # Plot fig, axs = plt.subplots(len(initial_guess)) @@ -546,7 +526,7 @@ def simulate_moments_redux(x): def make_contour_plot_obj_func( - agent_name, estimation_agent, model_estimate, targeted_moments + agent_name, estimation_agent, model_estimate, target_moments ): print("``````````````````````````````````````````````````````````````````````") print("Creating the contour plot.") @@ -569,7 +549,7 @@ def make_contour_plot_obj_func( DiscFacAdj, CRRA, agent=estimation_agent, - tgt_moments=targeted_moments, + target_moments=target_moments, ) smm_contour = plt.contourf(CRRA_mesh, DiscFacAdj_mesh, smm_obj_levels, level_count) t_end_contour = time() @@ -595,8 +575,8 @@ def estimate( compute_standard_errors=local_compute_standard_errors, compute_sensitivity=local_compute_sensitivity, make_contour_plot=local_make_contour_plot, - subjective_stock_market=local_subjective_stock_market, - subjective_labor_market=local_subjective_labor_market, + subjective_stock=local_subjective_stock, + subjective_labor=local_subjective_labor, ): """ Run the main estimation procedure for SolvingMicroDSOP. @@ -619,18 +599,18 @@ def estimate( estimation_agent, agent_name = make_estimation_agent( agent_name=agent_name, - subjective_stock_market=subjective_stock_market, - subjective_labor_market=subjective_labor_market, + subjective_stock=subjective_stock, + subjective_labor=subjective_labor, ) - targeted_moments = get_targeted_moments() + target_moments = get_target_moments() initial_guess = get_initial_guess(agent_name) # Estimate the model using Nelder-Mead if estimate_model: - model_estimate, time_to_estimate = estimate_model_nelder_mead( - agent_name, estimation_agent, targeted_moments, initial_guess + model_estimate, time_to_estimate = estimate_model_opt( + agent_name, estimation_agent, target_moments, initial_guess ) # Compute standard errors by bootstrap @@ -648,212 +628,9 @@ def estimate( # Make a contour plot of the objective function if make_contour_plot: make_contour_plot_obj_func( - agent_name, estimation_agent, model_estimate, targeted_moments - ) - - -def estimate_all(): - agent_types = [ - PortfolioLifeCycleConsumerType, - BequestWarmGlowLifeCycleConsumerType, - BequestWarmGlowLifeCyclePortfolioType, - WealthPortfolioLifeCycleConsumerType, - ] - - agent_names = [ - "(a) Portfolio Choice", - "(b) Separable Wealth in Utility", - "(c) Portfolio Separable Wealth in Utility", - "(d) Non-Separable Portfolio Wealth in Utility", - ] - - fig_sensitivity = plt.figure(layout="constrained", figsize=(14, 12)) - fig_sensitivity.suptitle("Sensitivity of Moments to Parameters") - - subfigs_sensitivity = fig_sensitivity.subfigures(2, 2) - - fig_contour = plt.figure(layout="constrained", figsize=(14, 12)) - fig_contour.suptitle("Contour Plots") - - subfigs_contour = fig_contour.subplots(2, 2) - - for count, agent_type in enumerate(agent_types): - # Make a lifecycle consumer to be used for estimation, including simulated - # shocks (plus an initial distribution of wealth) - # Make a TempConsumerType for estimation - estimation_agent = agent_type(**init_consumer_objects) - # Set the number of periods to simulate - estimation_agent.T_sim = estimation_agent.T_cycle + 1 - # Choose to track bank balances as wealth - estimation_agent.track_vars = ["bNrm"] - # Draw initial assets for each consumer - estimation_agent.aNrmInit = DiscreteDistribution( - options["initial_wealth_income_ratio_probs"], - options["initial_wealth_income_ratio_vals"], - seed=options["seed"], - ).draw(N=options["num_agents"]) - estimation_agent.make_shock_history() - - targeted_moments = get_targeted_moments() - - idx = np.unravel_index(count, (2, 2)) - - initial_guess = [options["DiscFacAdj_start"], options["CRRA_start"]] - print("----------------------------------------------------------------------") - print( - f"Now estimating the model using Nelder-Mead from an initial guess of {initial_guess}..." - ) - print("----------------------------------------------------------------------") - - # Make a single-input lambda function for use in the optimizer - def smm_obj_func_redux(parameters_to_estimate): - """ - A "reduced form" of the SMM objective function, compatible with the optimizer. - Identical to smmObjectiveFunction, but takes only a single input as a length-2 - list representing [DiscFacAdj,CRRA]. - """ - return smm_obj_func( - DiscFacAdj=parameters_to_estimate[0], - CRRA=parameters_to_estimate[1], - agent=estimation_agent, - tgt_moments=targeted_moments, - ) - - t_start_estimate = time() - model_estimate = minimize_nelder_mead( - smm_obj_func_redux, initial_guess, verbose=True - ) - t_end_estimate = time() - - time_to_estimate = t_end_estimate - t_start_estimate - - # Calculate minutes and remaining seconds - minutes, seconds = divmod(time_to_estimate, 60) - print(f"Time to estimate: {int(minutes)} min, {int(seconds)} sec.") - - print( - f"Estimated values: DiscFacAdj={model_estimate[0]}, CRRA={model_estimate[1]}" - ) - - print("``````````````````````````````````````````````````````````````````````") - print("Computing sensitivity measure.") - print("``````````````````````````````````````````````````````````````````````") - - # Find the Jacobian of the function that simulates moments - def simulate_moments_redux(x): - moments = simulate_moments( - x[0], - x[1], - agent=estimation_agent, - DiscFacAdj_bound=options["DiscFacAdj_bound"], - CRRA_bound=options["CRRA_bound"], - map_simulated_to_empirical_cohorts=scf_data.simulation_map_cohorts_to_age_indices, + agent_name, estimation_agent, model_estimate, target_moments ) - return moments - - n_moments = len(scf_data.simulation_map_cohorts_to_age_indices) - jac = np.array( - [ - approx_fprime( - model_estimate, - lambda x: simulate_moments_redux(x)[j], - epsilon=0.01, - ) - for j in range(n_moments) - ] - ) - - # Compute sensitivity measure. (all moments weighted equally) - sensitivity = np.dot(np.linalg.inv(np.dot(jac.T, jac)), jac.T) - - # Create lables for moments in the plots - moment_labels = [ - "[" + str(min(x)) + "," + str(max(x)) + "]" - for x in scf_data.empirical_cohort_age_groups - ] - - # Plot - subfigs_sensitivity[idx].suptitle(agent_names[count]) - subfigsnest_sensitivity = subfigs_sensitivity[idx].subplots(2, 1) - - subfigsnest_sensitivity[0].bar( - range(n_moments), sensitivity[0, :], tick_label=moment_labels - ) - subfigsnest_sensitivity[0].set_title("DiscFacAdj") - subfigsnest_sensitivity[0].set_ylabel("Sensitivity") - subfigsnest_sensitivity[0].set_xlabel("Median W/Y Ratio") - - subfigsnest_sensitivity[1].bar( - range(n_moments), sensitivity[1, :], tick_label=moment_labels - ) - subfigsnest_sensitivity[1].set_title("CRRA") - subfigsnest_sensitivity[1].set_ylabel("Sensitivity") - subfigsnest_sensitivity[1].set_xlabel("Median W/Y Ratio") - - print("``````````````````````````````````````````````````````````````````````") - print("Creating the contour plot.") - print("``````````````````````````````````````````````````````````````````````") - t_start_contour = time() - DiscFac_star, CRRA_star = model_estimate - grid_density = 20 # Number of parameter values in each dimension - level_count = 100 # Number of contour levels to plot - DiscFacAdj_list = np.linspace( - max(DiscFac_star - 0.25, 0.5), min(DiscFac_star + 0.25, 1.05), grid_density - ) - CRRA_list = np.linspace( - max(CRRA_star - 5, 2), min(CRRA_star + 5, 8), grid_density - ) - CRRA_mesh, DiscFacAdj_mesh = np.meshgrid(CRRA_list, DiscFacAdj_list) - smm_obj_levels = np.empty([grid_density, grid_density]) - for j in range(grid_density): - DiscFacAdj = DiscFacAdj_list[j] - for k in range(grid_density): - CRRA = CRRA_list[k] - smm_obj_levels[j, k] = smm_obj_func( - DiscFacAdj, - CRRA, - agent=estimation_agent, - tgt_moments=targeted_moments, - ) - - # Create figure and axes objects - - # Plot the contour - subfigs_contour[idx].set_title(agent_names[count]) - contour = subfigs_contour[idx].contourf( - CRRA_mesh, DiscFacAdj_mesh, smm_obj_levels, level_count - ) - fig_sensitivity.colorbar(contour) - - # Plot the model estimate - subfigs_contour[idx].plot(model_estimate[1], model_estimate[0], "*r", ms=15) - - # Set axis labels and title - subfigs_contour[idx].set_xlabel( - r"coefficient of relative risk aversion $\rho$", fontsize=14 - ) - subfigs_contour[idx].set_ylabel( - r"discount factor adjustment $\beth$", fontsize=14 - ) - - # Print time to execute - t_end_contour = time() - time_to_contour = t_end_contour - t_start_contour - - # Calculate minutes and remaining seconds - minutes, seconds = divmod(time_to_contour, 60) - print(f"Time to contour: {int(minutes)} min, {int(seconds)} sec.") - - fig_sensitivity.savefig(figures_dir + "/AllSensitivity.pdf") - fig_sensitivity.savefig(figures_dir + "/AllSensitivity.png") - fig_sensitivity.savefig(figures_dir + "/AllSensitivity.svg") - - # Save and show the plot - fig_contour.savefig(figures_dir + "/AllSMMcontour.pdf") - fig_contour.savefig(figures_dir + "/AllSMMcontour.png") - fig_contour.savefig(figures_dir + "/AllSMMcontour.svg") - if __name__ == "__main__": estimate() diff --git a/code/estimark/calibration/options.py b/code/estimark/options.py similarity index 82% rename from code/estimark/calibration/options.py rename to code/estimark/options.py index c2378a4..32051bb 100644 --- a/code/estimark/calibration/options.py +++ b/code/estimark/options.py @@ -6,8 +6,8 @@ "make_contour_plot": False, "compute_standard_errors": False, "compute_sensitivity": False, - "subjective_stock_market": False, - "subjective_labor_market": False, + "subjective_stock": False, + "subjective_labor": False, } # Author note: # This takes approximately 90 seconds on a laptop with the following specs: @@ -18,8 +18,8 @@ "make_contour_plot": True, "compute_standard_errors": False, "compute_sensitivity": True, - "subjective_stock_market": False, - "subjective_labor_market": False, + "subjective_stock": False, + "subjective_labor": False, } # Author note: # This takes approximately 7 minutes on a laptop with the following specs: @@ -30,8 +30,8 @@ "make_contour_plot": False, "compute_standard_errors": True, "compute_sensitivity": True, - "subjective_stock_market": False, - "subjective_labor_market": False, + "subjective_stock": False, + "subjective_labor": False, } # Author note: # This takes approximately 30 minutes on a laptop with the following specs: @@ -42,8 +42,8 @@ "make_contour_plot": True, "compute_standard_errors": True, "compute_sensitivity": True, - "subjective_stock_market": False, - "subjective_labor_market": False, + "subjective_stock": False, + "subjective_labor": False, } # Author note: # This takes approximately 40 minutes on a laptop with the following specs: diff --git a/code/run_all.py b/code/run_all.py index 64e531b..30f2578 100644 --- a/code/run_all.py +++ b/code/run_all.py @@ -1,4 +1,4 @@ -from estimark.calibration.options import ( +from estimark.options import ( all_replications, high_resource, low_resource, @@ -63,11 +63,11 @@ def run_replication(): return if subjective_markets == "2" or subjective_markets == "4": - replication_specs["subjective_stock_market"] = True + replication_specs["subjective_stock"] = True print("Adding subjective stock market beliefs...") if subjective_markets == "3" or subjective_markets == "4": - replication_specs["subjective_labor_market"] = True + replication_specs["subjective_labor"] = True print("Adding subjective labor market beliefs...") lazy_result = dask.delayed(estimate)(**replication_specs) diff --git a/code/tables/IndShock_estimate_results.csv b/code/tables/IndShock_estimate_results.csv index b9394bb..576f746 100644 --- a/code/tables/IndShock_estimate_results.csv +++ b/code/tables/IndShock_estimate_results.csv @@ -1,2 +1,2 @@ DiscFacAdj,CRRA -0.9561845217685123,1.3419829066995432 +0.9562452937085999,1.3387414303918428 diff --git a/code/tests.py b/code/tests.py index 0c7c35d..ac288a7 100644 --- a/code/tests.py +++ b/code/tests.py @@ -1,4 +1,4 @@ -from estimark.calibration.options import low_resource, medium_resource +from estimark.options import low_resource, medium_resource from estimark.estimation import estimate