target shares

.gitignore

@@ -148,3 +148,4 @@ cython_debug/
 
 .ipynb_checkpoints
 _build
+code/data/S&P Target Date glidepath.xlsx

code/estimark/estimation.py

@@ -52,6 +52,7 @@
     scf_mapping,
     scf_weights,
 )
+from estimark.snp import snp
 
 # Pathnames to the other files:
 # Relative directory for primitive parameter files
@@ -115,7 +116,10 @@ def make_estimation_agent(
     # 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"]
+    track_vars = ["bNrm"]
+    if "Portfolio" in agent_name:
+        track_vars += ["Share"]
+    estimation_agent.track_vars = track_vars
     # Draw initial assets for each consumer
     estimation_agent.aNrmInit = DiscreteDistribution(
         options["prob_w_to_y"],
@@ -157,6 +161,14 @@ def get_targeted_moments(
     return target_moments
 
 
+share_moments = get_targeted_moments(
+    data=snp["S&P Target Date Equity allocation"].to_numpy(),
+    weights=np.ones(len(snp)),
+    groups=snp["age_groups"].to_numpy(),
+    mapping=scf_mapping,
+)
+
+
 def get_initial_guess(agent_name):
     # start from previous estimation results if available
 
@@ -172,35 +184,33 @@ def get_initial_guess(agent_name):
 
 
 # Define the objective function for the simulated method of moments estimation
-# TODO: params, bounds, agent
-def simulate_moments(params, agent):
+def simulate_moments(params, agent, agent_name):
     """A quick check to make sure that the parameter values are within bounds.
     Far flung falues of DiscFacAdj or CRRA might cause an error during solution or
     simulation, so the objective function doesn't even bother with them.
     """
     DiscFacAdj, CRRA = params
 
-    # # TODO: bounds should be handled by the optimizer
-    # bounds_DiscFacAdj = options["bounds_DiscFacAdj"]
-    # bounds_CRRA = options["bounds_CRRA"]
-
-    # if (
-    #     DiscFacAdj < bounds_DiscFacAdj[0]
-    #     or DiscFacAdj > bounds_DiscFacAdj[1]
-    #     or CRRA < bounds_CRRA[0]
-    #     or CRRA > bounds_CRRA[1]
-    # ):
-    #     return 1e30 * np.ones(len(scf_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["timevary_DiscFac"]]
     agent.CRRA = CRRA
-    if hasattr(agent, "BeqCRRA"):
-        agent.BeqCRRA = CRRA
+    # if hasattr(agent, "BeqCRRA"):
+    #     agent.BeqCRRA = [CRRA] * len(options["timevary_DiscFac"])
     # 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")
+
+    # simulate with true parameters (override subjective beliefs)
+    if "(Stock)" in agent_name and "Portfolio" in agent_name:
+        agent.RiskyAvg = agent.RiskyAvgTrue
+        agent.RiskyStd = agent.RiskyStdTrue
+        agent.update_RiskyDstn()
+    if "(Labor)" in agent_name:
+        agent.TranShkStd = init_consumer_objects["TranShkStd"]
+        agent.PermShkStd = init_consumer_objects["PermShkStd"]
+        agent.update_income_process()
+
     max_sim_age = max([max(ages) for ages in scf_mapping]) + 1
     # Initialize the simulation by clearing histories, resetting initial values
     agent.initialize_sim()
@@ -219,14 +229,19 @@ def simulate_moments(params, agent):
 
     sim_moments = np.array(sim_moments)
 
-    # # TODO: too many of these, check if solving/simulating has bug
-    # if np.isnan(sim_moments).any():
-    #     return 1e30 * np.ones(len(scf_mapping))
+    if "Portfolio" in agent_name:
+        sim_share_history = agent.history["Share"]
+        share_moments = []
+        for g in range(group_count):
+            cohort_indices = scf_mapping[g]
+            share_moments += [np.median(sim_share_history[cohort_indices])]
+
+        sim_moments = np.append(sim_moments, share_moments)
 
     return sim_moments
 
 
-def smm_obj_func(params, agent, moments):
+def smm_obj_func(params, agent, moments, agent_name):
     """The objective function for the SMM estimation.  Given values of discount factor
     adjuster DiscFacAdj, coeffecient of relative risk aversion CRRA, a base consumer
     agent type, empirical data, and calibrated parameters, this function calculates
@@ -273,15 +288,25 @@ def smm_obj_func(params, agent, moments):
         median wealth-to-permanent-income ratio in the simulation.
 
     """
-    sim_moments = simulate_moments(params, agent)
+    if "Portfolio" in agent_name:
+        moments = np.append(moments, share_moments)
+
+    sim_moments = simulate_moments(params, agent, agent_name)
     errors = moments - sim_moments
     loss = np.dot(errors, errors)
 
     return {"value": loss, "root_contributions": errors}
 
 
 # Define the bootstrap procedure
-def calculate_se_bootstrap(initial_estimate, N, agent, seed=0, verbose=False):
+def calculate_se_bootstrap(
+    initial_estimate,
+    N,
+    agent,
+    agent_name,
+    seed=0,
+    verbose=False,
+):
     """Calculates standard errors by repeatedly re-estimating the model with datasets
     resampled from the actual data.
 
@@ -333,6 +358,7 @@ def smm_obj_func_bootstrap(params):
                 params,
                 agent=agent,
                 moments=bootstrap_moments,
+                agent_name=agent_name,
             )
 
         # Estimate the model with the bootstrap data and add to list of estimates
@@ -376,6 +402,7 @@ def smm_obj_func_redux(params):
             params,
             agent=estimation_agent,
             moments=target_moments,
+            agent_name=agent_name,
         )
 
     t_start_estimate = time()
@@ -389,7 +416,7 @@ def smm_obj_func_redux(params):
         lower_bounds=np.array(
             [options["bounds_DiscFacAdj"][0], options["bounds_CRRA"][0]],
         ),
-        # multistart=True,
+        multistart=True,
         error_handling="continue",
     )
     t_end_estimate = time()
@@ -549,6 +576,7 @@ def do_make_contour_plot(agent_name, estimation_agent, model_estimate, target_mo
                 np.array([DiscFacAdj, CRRA]),
                 agent=estimation_agent,
                 moments=target_moments,
+                agent_name=agent_name,
             )
     smm_contour = plt.contourf(CRRA_mesh, DiscFacAdj_mesh, smm_obj_levels, level_count)
     t_end_contour = time()

code/estimark/parameters.py

@@ -16,7 +16,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
-aXtraCount = 100  # Number of points in the grid of "assets above minimum"
+aXtraCount = 20  # Number of points in the grid of "assets above minimum"
 
 # Artificial borrowing constraint; imposed minimum level of end-of period assets
 BoroCnstArt = 0.0
@@ -34,20 +34,22 @@
 IncUnemp = 0.3  # Unemployment benefits replacement rate
 IncUnempRet = 0.0  # "Unemployment" benefits when retired
 
-final_age = 90  # Age at which the problem ends (die with certainty)
+final_age = 120  # 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
 terminal_t = final_age - initial_age  # Total number of periods in the model
 retirement_t = retirement_age - initial_age - 1
 
+final_age_data = 95  # Age at which the data ends
+
 # Initial guess of the coefficient of relative risk aversion during estimation (rho)
 init_CRRA = 5.0
 # Initial guess of the adjustment to the discount factor during estimation (beth)
 init_DiscFacAdj = 0.99
 # Bounds for beth; if violated, objective function returns "penalty value"
 bounds_DiscFacAdj = [0.5, 1.5]
 # Bounds for rho; if violated, objective function returns "penalty value"
-bounds_CRRA = [1.1, 10.0]
+bounds_CRRA = [1.1, 20.0]
 
 # Income
 ss_variances = True
@@ -64,8 +66,9 @@
 # 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"
-timevary_DiscFac = np.genfromtxt(csv_file_path) * 0.0 + 1.0  # TODO
-constant_DiscFac = np.ones_like(timevary_DiscFac)
+# timevary_DiscFac = np.genfromtxt(csv_file_path) * 0.0 + 1.0  # TODO
+# constant_DiscFac = np.ones_like(timevary_DiscFac)
+timevary_DiscFac = np.ones_like(inc_calib["PermShkStd"])
 
 # Survival probabilities over the lifecycle
 liv_prb = parse_ssa_life_table(
@@ -82,7 +85,7 @@
 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
+seed = 1132023  # Just an integer to seed the estimation
 
 options = {
     "init_w_to_y": init_w_to_y,

code/estimark/scf.py

@@ -6,7 +6,7 @@
 import numpy as np  # Numerical Python
 import pandas as pd
 
-from estimark.parameters import initial_age
+from estimark.parameters import final_age_data, 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"
@@ -21,7 +21,10 @@
 
 # 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)]
+age_groups = [
+    list(range(start, start + 5))
+    for start in range(initial_age + 1, final_age_data + 1, 5)
+]
 
 # Initialize empty lists for the data
 scf_data = []  # Ratio of wealth to permanent income

code/estimark/snp.py

@@ -0,0 +1,25 @@
+from pathlib import Path
+
+import numpy as np  # Numerical Python
+import pandas as pd
+
+from estimark.parameters import final_age_data, initial_age
+
+file_path = (
+    Path(__file__).resolve().parent / ".." / "data" / "S&P Target Date glidepath.xlsx"
+)
+
+# Load data
+snp = pd.read_excel(file_path)
+
+# Filter data using loc
+snp = snp.loc[
+    (snp["Current Age"] >= initial_age) & (snp["Current Age"] <= final_age_data)
+]
+
+# Create age groups and code NaNs as 0
+bins = range(initial_age, final_age_data + 1, 5)
+labels = np.arange(1, len(bins))
+snp["age_groups"] = pd.cut(snp["Current Age"], bins=bins, labels=labels)
+
+# Get targeted moments