use estimagic

code/estimark/estimation.py

@@ -9,6 +9,8 @@
 income as defined in ConsIndShockModel.
 """
 
+import pandas as pd
+
 import csv
 from time import time  # Timing utility
 
@@ -51,6 +53,8 @@
     scf_weights,
 )
 
+import estimagic as em
+
 # Pathnames to the other files:
 # Relative directory for primitive parameter files
 tables_dir = "content/tables/"
@@ -153,12 +157,12 @@ def get_targeted_moments(
 def get_initial_guess(agent_name):
     # start from previous estimation results if available
 
-    csv_file_path = tables_dir + agent_name + "_estimate_results.csv"
+    csv_file_path = f"{tables_dir}{agent_name}_estimate_results.csv"
 
     try:
-        with open(csv_file_path, "r") as file:
-            initial_guess = np.loadtxt(file, skiprows=1, delimiter=",")
-    except FileNotFoundError:
+        res = pd.read_csv(csv_file_path, header=None)
+        initial_guess = res.iloc[:2, 1].astype(float).tolist()
+    except (FileNotFoundError, IndexError):
         initial_guess = [options["init_DiscFacAdj"], options["init_CRRA"]]
 
     return initial_guess
@@ -175,6 +179,7 @@ def simulate_moments(params, agent):
 
     DiscFacAdj, CRRA = params
 
+    # todo: bounds should be handled by the optimizer
     bounds_DiscFacAdj = options["bounds_DiscFacAdj"]
     bounds_CRRA = options["bounds_CRRA"]
 
@@ -211,6 +216,10 @@ 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))
+
     return sim_moments
 
 
@@ -367,13 +376,23 @@ def smm_obj_func_redux(params):
         )
 
     t_start_estimate = time()
-    model_estimate = minimize_nelder_mead(
-        smm_obj_func_redux, initial_guess, verbose=True
+    res = em.minimize(
+        smm_obj_func_redux,
+        initial_guess,
+        algorithm="scipy_neldermead",
+        upper_bounds=np.array(
+            [options["bounds_DiscFacAdj"][1], options["bounds_CRRA"][1]]
+        ),
+        lower_bounds=np.array(
+            [options["bounds_DiscFacAdj"][0], options["bounds_CRRA"][0]]
+        ),
+        multistart=True,
     )
     t_end_estimate = time()
-
     time_to_estimate = t_end_estimate - t_start_estimate
 
+    model_estimate = res.params
+
     # Calculate minutes and remaining seconds
     minutes, seconds = divmod(time_to_estimate, 60)
     print(f"Time to estimate: {int(minutes)} min, {int(seconds)} sec.")
@@ -385,8 +404,19 @@ def smm_obj_func_redux(params):
 
     with open(estimate_results_file, "wt") as f:
         writer = csv.writer(f)
-        writer.writerow(["DiscFacAdj", "CRRA"])
-        writer.writerow([model_estimate[0], model_estimate[1]])
+
+        writer.writerow(["DiscFacAdj", model_estimate[0]])
+        writer.writerow(["CRRA", model_estimate[1]])
+        writer.writerow(["time_to_estimate", time_to_estimate])
+
+        for key in vars(res):
+            if key not in [
+                "history",
+                "convergence_report",
+                "multistart_info",
+                "algorithm_output",
+            ]:
+                writer.writerow([key, getattr(res, key)])
 
     return model_estimate, time_to_estimate
 

code/estimark/parameters.py

@@ -18,7 +18,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 = 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
@@ -47,9 +47,9 @@
 # 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.0001, 15.0]
+bounds_DiscFacAdj = [0.5, 1.5]
 # Bounds for rho; if violated, objective function returns "penalty value"
-bounds_CRRA = [0.0001, 15.0]
+bounds_CRRA = [1.1, 10.0]
 
 # Income
 ss_variances = True

content/tables/IndShockSub(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9826952831981901,3.768367697583021
+DiscFacAdj,0.9826417687605296
+CRRA,3.8206605179089514
+time_to_estimate,273.235800743103
+params,[0.986417687605296, 3.8206605179089514]
+criterion,0.9501149506157349
+start_criterion,0.9501149506157349
+start_params,[0.9826417687605296, 3.8206605179089514]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,348
+n_derivative_evaluations,
+n_iterations,177

content/tables/IndShockSub(Stock)(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.990059541559503,4.128239153345464
+DiscFacAdj,0.9900891154527984
+CRRA,4.100899179990721
+time_to_estimate,229.75330018997192
+params,[0.990891154527984, 4.100899179990721]
+criterion,1.133972554008415
+start_criterion,1.133972554014542
+start_params,[0.990089115470004, 4.100899171954525]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,298
+n_derivative_evaluations,
+n_iterations,157

content/tables/IndShockSub(Stock)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9628905096680382,1.442713138370085
+DiscFacAdj,0.9626659339838747
+CRRA,1.450655649628727
+time_to_estimate,278.38573718070984
+params,[0.966659339838747, 1.450655649628727]
+criterion,5.728180713297819
+start_criterion,5.728180713297819
+start_params,[0.9626659339838747, 1.450655649628727]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,336
+n_derivative_evaluations,
+n_iterations,177

content/tables/IndShock_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9564067036912898,1.3327634296841504
+DiscFacAdj,0.9562461537688056
+CRRA,1.3387132918993245
+time_to_estimate,280.9034581184387
+params,[0.952461537688056, 1.3387132918993245]
+criterion,5.3374265466181985
+start_criterion,5.337426546939854
+start_params,[0.9562461533319309, 1.3387133067909955]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,326
+n_derivative_evaluations,
+n_iterations,168

content/tables/PortfolioSub(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9909618480820204,5.986965941730887
+DiscFacAdj,0.9809297742019563
+CRRA,5.618112957475339
+time_to_estimate,603.8376204967499
+params,[0.989297742019563, 5.618112957475339]
+criterion,2.665832181614686
+start_criterion,2.666329553441798
+start_params,[0.9809051537052127, 5.51663750438404]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,335
+n_derivative_evaluations,
+n_iterations,171

content/tables/PortfolioSub(Stock)(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9911792224594733,6.034080046465693
+DiscFacAdj,0.9911785704835548
+CRRA,6.034152052728421
+time_to_estimate,641.9356205463409
+params,[0.991785704835548, 6.034152052728421]
+criterion,4.963398992709653
+start_criterion,4.963398992709666
+start_params,[0.9911785704835547, 6.034152052728421]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,375
+n_derivative_evaluations,
+n_iterations,185

content/tables/PortfolioSub(Stock)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9913555579379695,6.010819574001086
+DiscFacAdj,0.601523327991334
+CRRA,8.317976231224396
+time_to_estimate,581.5355100631714
+params,[0.60523327991334, 8.317976231224396]
+criterion,13.408766511615521
+start_criterion,13.409113025891063
+start_params,[0.6018344625104538, 8.31331685373428]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,322
+n_derivative_evaluations,
+n_iterations,164

content/tables/Portfolio_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.927363265767304,1.4738986792731867
+DiscFacAdj,0.92750018606101
+CRRA,1.4600877205227354
+time_to_estimate,643.3507697582245
+params,[0.9250018606101, 1.4600877205227354]
+criterion,8.947670723898293
+start_criterion,8.947670724237371
+start_params,[0.9275001864339323, 1.4600876974824124]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,377
+n_derivative_evaluations,
+n_iterations,193

content/tables/WarmGlowPortfolioSub(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9756576356903954,8.845926577241347
+DiscFacAdj,0.9780667056620715
+CRRA,5.686287061969352
+time_to_estimate,560.7655048370361
+params,[0.970667056620715, 5.686287061969352]
+criterion,2.51020591715548
+start_criterion,2.51026533163372
+start_params,[0.9780780923691006, 5.681434389740495]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,296
+n_derivative_evaluations,
+n_iterations,157

content/tables/WarmGlowPortfolioSub(Stock)(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9756558075814038,8.846028566854438
+DiscFacAdj,0.9760267917068915
+CRRA,8.720502188660854
+time_to_estimate,797.6410620212555
+params,[0.970267917068915, 8.720502188660854]
+criterion,3.1554655208054507
+start_criterion,3.155465521671051
+start_params,[0.9760267914875025, 8.720502190077596]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,503
+n_derivative_evaluations,
+n_iterations,266

content/tables/WarmGlowPortfolioSub(Stock)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9756576356903954,8.845926577241347
+DiscFacAdj,0.5989021241004722
+CRRA,8.372281839569975
+time_to_estimate,571.794606924057
+params,[0.599021241004722, 8.372281839569975]
+criterion,13.411282146299964
+start_criterion,13.411306700788627
+start_params,[0.6008994808072222, 8.34917075904358]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,303
+n_derivative_evaluations,
+n_iterations,156

content/tables/WarmGlowPortfolio_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9756570720304689,8.846018147853837
+DiscFacAdj,0.9276703440995966
+CRRA,1.454260290778273
+time_to_estimate,548.6548085212708
+params,[0.926703440995966, 1.454260290778273]
+criterion,8.922913379322182
+start_criterion,8.922913379322186
+start_params,[0.9276703440995965, 1.454260290778273]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,301
+n_derivative_evaluations,
+n_iterations,154

content/tables/WarmGlowSub(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9791782282367655,4.231669394199844
+DiscFacAdj,0.9789789171789187
+CRRA,4.301704752722339
+time_to_estimate,278.50627040863037
+params,[0.979789171789187, 4.301704752722339]
+criterion,0.8849773434503059
+start_criterion,0.8849773434503059
+start_params,[0.9789789171789187, 4.301704752722339]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,351
+n_derivative_evaluations,
+n_iterations,180

content/tables/WarmGlowSub(Stock)(Labor)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9861730357696492,4.425749877454312
+DiscFacAdj,0.9860432110010826
+CRRA,4.476854619240293
+time_to_estimate,261.37478494644165
+params,[0.980432110010826, 4.476854619240293]
+criterion,1.0401767350001887
+start_criterion,1.040176735000193
+start_params,[0.9860432110010825, 4.476854619240292]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,322
+n_derivative_evaluations,
+n_iterations,168

content/tables/WarmGlowSub(Stock)Market_estimate_results.csv

@@ -1,2 +1,15 @@
-DiscFacAdj,CRRA
-0.9861734987743034,4.425678459697778
+DiscFacAdj,0.9623327944560536
+CRRA,1.3852355775601488
+time_to_estimate,283.0823395252228
+params,[0.963327944560536, 1.3852355775601488]
+criterion,5.410433480030225
+start_criterion,5.410433480030225
+start_params,[0.9623327944560536, 1.3852355775601488]
+algorithm,multistart_scipy_neldermead
+direction,minimize
+n_free,2
+message,Optimization terminated successfully.
+success,True
+n_criterion_evaluations,326
+n_derivative_evaluations,
+n_iterations,166