From d529a44c8eddc1bbb60a183dc02f5cd6ca3634e0 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 3 Jun 2024 09:05:05 -0400 Subject: [PATCH 01/16] sketch of a general value backwards induction algorithm built on DBlocks --- Documentation/CHANGELOG.md | 1 + Documentation/reference/tools/algos.rst | 7 + Documentation/reference/tools/algos/vbi.rst | 7 + Documentation/reference/tools/index.rst.orig | 27 ++++ HARK/algos/vbi.py | 153 +++++++++++++++++++ HARK/model.py | 8 + 6 files changed, 203 insertions(+) create mode 100644 Documentation/reference/tools/algos.rst create mode 100644 Documentation/reference/tools/algos/vbi.rst create mode 100644 Documentation/reference/tools/index.rst.orig create mode 100644 HARK/algos/vbi.py diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index a1b7f9604..87f33a401 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -24,6 +24,7 @@ Note: Due to major changes on this release, you may need to adjust how AgentType - All methods that construct inputs for solvers are now functions that are specified in the dictionary attribute `constructors`. [#1410](https://github.com/econ-ark/HARK/pull/1410) - Such constructed inputs can use alternate parameterizations / formats by changing the `constructor` function and providing its arguments in `parameters`. - Move `HARK.datasets` to `HARK.Calibration` for better organization of data and calibration tools. [#1430](https://github.com/econ-ark/HARK/pull/1430) +- Adds `HARK.algos.vbi` as a general algorithm for solving the optimization step of an agent's problem. ### Minor Changes diff --git a/Documentation/reference/tools/algos.rst b/Documentation/reference/tools/algos.rst new file mode 100644 index 000000000..70f2b822f --- /dev/null +++ b/Documentation/reference/tools/algos.rst @@ -0,0 +1,7 @@ +algos +-------- + +.. toctree:: + :maxdepth: 3 + + algos/vbi diff --git a/Documentation/reference/tools/algos/vbi.rst b/Documentation/reference/tools/algos/vbi.rst new file mode 100644 index 000000000..918b8020d --- /dev/null +++ b/Documentation/reference/tools/algos/vbi.rst @@ -0,0 +1,7 @@ +algos.vbi +---------- + +.. automodule:: HARK.algos.vbi + :members: + :undoc-members: + :show-inheritance: diff --git a/Documentation/reference/tools/index.rst.orig b/Documentation/reference/tools/index.rst.orig new file mode 100644 index 000000000..b1bb19aef --- /dev/null +++ b/Documentation/reference/tools/index.rst.orig @@ -0,0 +1,27 @@ +:orphan: + +Tools +===== + +<<<<<<< HEAD +.. toctree:: + :maxdepth: 3 + + algos + core + dcegm + distribution + econforgeinterp + estimation + frame + helpers + interpolation + numba + parallel + rewards + simulation + utilities + validators +======= +See :doc:`../index`. +>>>>>>> master diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py new file mode 100644 index 000000000..54efa2753 --- /dev/null +++ b/HARK/algos/vbi.py @@ -0,0 +1,153 @@ +""" +Use backwards induction to derive the arrival value function +from a continuation value function and stage dynamics. +""" + +from HARK.model import DBlock +import itertools +from scipy.optimize import minimize, brentq +from typing import Mapping + + +def get_action_rule(action): + """ + Produce a function from any inputs to a given value. + This is useful for constructing decision rules with fixed actions. + """ + + def ar(*args): + return action + + return ar + + +def ar_from_data(da): + """ + Produce a function from any inputs to a given value. + This is useful for constructing decision rules with fixed actions. + """ + + def ar(**args): + return da.interp(**args) + + return ar + + +Grid = Mapping[str, Sequence] + + +def grid_to_data_array( + grid: Grid = {}, ## TODO: Better data structure here. +): + """ + Construct a zero-valued DataArray with the coordinates + based on the Grid passed in. + + Parameters + ---------- + grid: Grid + A mapping from variable labels to a sequence of numerical values. + + Returns + -------- + da xarray.DataArray + An xarray.DataArray with coordinates given by both grids. + """ + + coords = {**grid} + + da = xr.DataArray( + np.empty([len(v) for v in coords.values()]), dims=coords.keys(), coords=coords + ) + + return da + + +def vbi_solve( + block: DBlock, continuation, state_grid: Grid, disc_params, calibration={} +): + """ + Solve a DBlock using backwards induction on the value function. + + Parameters + ----------- + block + continuation + + state_grid: Grid + This is a grid over all variables that the optimization will range over. + This should be just the information set of the decision variables. + + disc_params + calibration + """ + + # state-rule value function + srv_function = block.get_state_rule_value_function_from_continuation(continuation) + + controls = block.get_controls() + + # pseudo + policy_data = grid_to_data_array(state_grid) + value_data = grid_to_data_array(state_grid) + + # loop through every point in the state grid + for state_point in itertools.product(*state_grid.values()): + # build a dictionary from these states, as scope for the optimization + state_vals = {k: v for k, v in zip(state_grid.keys(), state_point)} + + # copy calibration + # update with state_vals + # this becomes the before-action states + pre_states = state_vals + parameters + + # prepare function to optimize + def negated_value(a): # old! (should be negative) + dr = {c: get_action_rule(a[i]) for i, c in enumerate(controls)} + + # negative, for minimization later + return -srv_function(pre_states, dr) + + ## get lower bound. + ## not yet implemented + lower_bound = np.array([-1e-12] * len(controls)) ## a really low number! + + ## get upper bound + ## not yet implemented + upper_bound = np.array([1e12] * len(controls)) + + # pseudo + # optimize_action(pre_states, srv_function) + + a_best, root_res = minimize( # choice of + negated_value, full_output=True + ) + + dr_best = {c: get_action_rule(a_best[i]) for i, c in enumerate(controls)} + + if root_res.converged: + policy_data.sel(**state_vals).variable.data.put(0, a_best) + value_data.sel(**state_vals).variable.data.put( + 0, srv_function(pre_states, dr_best) + ) + else: + print(f"Optimization failure at {state_vals}.") + print(root_res) + + dr_best = {c: get_action_rule(root_res[i]) for i, c in enumerate(controls)} + + policy_data.sel(**state_vals).variable.data.put(0, res.root) # ? + value_data.sel(**state_vals).variable.data.put( + 0, srv_function(pre_states, dr_best) + ) + + # use the xarray interpolator to create a decision rule. + dr_from_data = { + c: ar_from_data(da) # maybe needs to be more sensitive to the information set + for i, c in enumerate(controls) + } + + dec_vf = block.get_decision_value_function(dr_from_data, continuation) + arr_vf = block.get_arrival_value_function(disc_params, dr_from_data, continuation) + + return dr_from_data, dec_vf, arr_vf \ No newline at end of file diff --git a/HARK/model.py b/HARK/model.py index c97a2ad76..a36bc57fa 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -161,6 +161,14 @@ def get_dynamics(self): def get_vars(self): return list(self.shocks.keys()) + list(self.dynamics.keys()) + def get_controls(self): + """ + TODO: Repeated in RBlock. Move to higher order class. + """ + dyn = self.get_dynamics() + + return [varn for varn in dyn if isinstance(dyn[varn], Control)] + def transition(self, pre, dr): """ Returns variable values given previous values and decision rule for all controls. From 431fc80b39b955c359ec45d28be72a104ed186e5 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 3 Jun 2024 09:14:00 -0400 Subject: [PATCH 02/16] ruff --- HARK/algos/vbi.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index 54efa2753..0120b3ea6 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -5,7 +5,7 @@ from HARK.model import DBlock import itertools -from scipy.optimize import minimize, brentq +from scipy.optimize import minimize from typing import Mapping @@ -150,4 +150,4 @@ def negated_value(a): # old! (should be negative) dec_vf = block.get_decision_value_function(dr_from_data, continuation) arr_vf = block.get_arrival_value_function(disc_params, dr_from_data, continuation) - return dr_from_data, dec_vf, arr_vf \ No newline at end of file + return dr_from_data, dec_vf, arr_vf From fd603b5e184b218ed5c4bb286ce65cbce19a35af Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 3 Jun 2024 09:46:34 -0400 Subject: [PATCH 03/16] fixes for doc errors --- HARK/algos/__init__.py | 0 HARK/algos/vbi.py | 3 ++- 2 files changed, 2 insertions(+), 1 deletion(-) create mode 100644 HARK/algos/__init__.py diff --git a/HARK/algos/__init__.py b/HARK/algos/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index 0120b3ea6..4b6049700 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -6,7 +6,8 @@ from HARK.model import DBlock import itertools from scipy.optimize import minimize -from typing import Mapping +from typing import Mapping, Sequence +import xarray as xr def get_action_rule(action): From b1daaff215aa9dc6d46dec903265912519e82069 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 3 Jun 2024 09:54:10 -0400 Subject: [PATCH 04/16] adding algos to docs table of contents --- Documentation/reference/index.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/Documentation/reference/index.rst b/Documentation/reference/index.rst index de65cad62..9d0f223b3 100644 --- a/Documentation/reference/index.rst +++ b/Documentation/reference/index.rst @@ -5,6 +5,7 @@ API Reference :caption: Tools :maxdepth: 1 + tools/algos tools/core tools/dcegm tools/distribution From 510acaf0c5bb9e677344493c0f8f2d9d6cfc5ba8 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 4 Jun 2024 15:26:11 -0400 Subject: [PATCH 05/16] automated test for vbi solution algorithm --- HARK/algos/tests/test_vbi.py | 34 ++++++++++++++++++++++++++++++++++ HARK/algos/vbi.py | 34 +++++++++++++++++++++------------- HARK/model.py | 13 ++++++++++++- 3 files changed, 67 insertions(+), 14 deletions(-) create mode 100644 HARK/algos/tests/test_vbi.py diff --git a/HARK/algos/tests/test_vbi.py b/HARK/algos/tests/test_vbi.py new file mode 100644 index 000000000..a698bacec --- /dev/null +++ b/HARK/algos/tests/test_vbi.py @@ -0,0 +1,34 @@ +import HARK.algos.vbi as vbi +from HARK.distribution import Bernoulli +from HARK.model import Control, DBlock +import numpy as np + +import unittest + + +block_1 = DBlock( + **{ + "name": "vbi_test_1", + "shocks": { + "coin": Bernoulli(p=0.5), + }, + "dynamics": { + "m": lambda y, coin: y + coin, + "c": Control(["m"]), + "a": lambda m, c: m - c, + }, + "reward": {"u": lambda c: 1 - (c - 1) ** 2}, + } +) + + +class test_vbi(unittest.TestCase): + # def setUp(self): + # pass + + def test_solve_block_1(self): + state_grid = {"m": np.linspace(0, 2, 10)} + + dr, dec_vf, arr_vf = vbi.vbi_solve(block_1, lambda a: a, state_grid) + + self.assertAlmostEqual(dr["c"](**{"m": 1}), 0.5) diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index 4b6049700..e56e04618 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -5,6 +5,7 @@ from HARK.model import DBlock import itertools +import numpy as np from scipy.optimize import minimize from typing import Mapping, Sequence import xarray as xr @@ -16,7 +17,7 @@ def get_action_rule(action): This is useful for constructing decision rules with fixed actions. """ - def ar(*args): + def ar(): return action return ar @@ -29,7 +30,7 @@ def ar_from_data(da): """ def ar(**args): - return da.interp(**args) + return da.interp(**args).values.tolist() return ar @@ -65,7 +66,7 @@ def grid_to_data_array( def vbi_solve( - block: DBlock, continuation, state_grid: Grid, disc_params, calibration={} + block: DBlock, continuation, state_grid: Grid, disc_params={}, calibration={} ): """ Solve a DBlock using backwards induction on the value function. @@ -97,10 +98,11 @@ def vbi_solve( # build a dictionary from these states, as scope for the optimization state_vals = {k: v for k, v in zip(state_grid.keys(), state_point)} - # copy calibration - # update with state_vals - # this becomes the before-action states - pre_states = state_vals + parameters + # The value of the action is computed given + # the problem calibration and the states for the current point on the + # state-grid. + pre_states = calibration.copy() + pre_states.update(state_vals) # prepare function to optimize def negated_value(a): # old! (should be negative) @@ -120,14 +122,18 @@ def negated_value(a): # old! (should be negative) # pseudo # optimize_action(pre_states, srv_function) - a_best, root_res = minimize( # choice of - negated_value, full_output=True + res = minimize( # choice of + negated_value, + 0, # x0 is starting guess, here arbitrary. ) + print(res) - dr_best = {c: get_action_rule(a_best[i]) for i, c in enumerate(controls)} + dr_best = {c: get_action_rule(res.x[i]) for i, c in enumerate(controls)} - if root_res.converged: - policy_data.sel(**state_vals).variable.data.put(0, a_best) + if res.success: + policy_data.sel(**state_vals).variable.data.put( + 0, res.x[0] + ) # will only work for scalar actions value_data.sel(**state_vals).variable.data.put( 0, srv_function(pre_states, dr_best) ) @@ -144,7 +150,9 @@ def negated_value(a): # old! (should be negative) # use the xarray interpolator to create a decision rule. dr_from_data = { - c: ar_from_data(da) # maybe needs to be more sensitive to the information set + c: ar_from_data( + policy_data + ) # maybe needs to be more sensitive to the information set for i, c in enumerate(controls) } diff --git a/HARK/model.py b/HARK/model.py index a36bc57fa..7fea2c991 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -173,7 +173,18 @@ def transition(self, pre, dr): """ Returns variable values given previous values and decision rule for all controls. """ - return simulate_dynamics(self.dynamics, pre, dr) + dyn = self.dynamics.copy() + + # don't simulate values that have already been given. + # this will break if there's a directly recursive label, + # i.e. if dynamics at time t for variable 'a' + # depend on state of 'a' at time t-1 + # This is a forbidden case in CDC's design. + for varn in pre: + if varn in dyn: + del dyn[varn] + + return simulate_dynamics(dyn, pre, dr) def calc_reward(self, vals): """ From 7b0bf4f43e1bb4448dbc9e52e429900087d23429 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 4 Jun 2024 15:53:24 -0400 Subject: [PATCH 06/16] changelog for vbi solver in 0.16.0 release now --- Documentation/CHANGELOG.md | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index e636e8b0b..a02479600 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -8,6 +8,17 @@ For more information on HARK, see [our Github organization](https://github.com/e ## Changes +### 0.16.0 + +Under development + +#### Major Changes + +- Adds `HARK.algos.vbi` as a general algorithm for solving the optimization step of an agent's problem. [#1438](https://github.com/econ-ark/HARK/pull/1438) + +#### Minor Changes + + ### 0.15.0 Release Date: June 4, 2024 @@ -16,7 +27,7 @@ Note: Due to major changes on this release, you may need to adjust how AgentType This release drops support for Python 3.8 and 3.9, consistent with SPEC 0, and adds support for Python 3.11 and 3.12. We expect that all HARK features still work with the older versions, but they are no longer part of our testing regimen. -### Major Changes +#### Major Changes - Drop official support for Python 3.8 and 3.9, add support for 3.11 and 3.12. [#1415](https://github.com/econ-ark/HARK/pull/1415) - Replace object-oriented solvers with single function versions. [#1394](https://github.com/econ-ark/HARK/pull/1394) @@ -27,9 +38,8 @@ This release drops support for Python 3.8 and 3.9, consistent with SPEC 0, and a - All methods that construct inputs for solvers are now functions that are specified in the dictionary attribute `constructors`. [#1410](https://github.com/econ-ark/HARK/pull/1410) - Such constructed inputs can use alternate parameterizations / formats by changing the `constructor` function and providing its arguments in `parameters`. - Move `HARK.datasets` to `HARK.Calibration` for better organization of data and calibration tools. [#1430](https://github.com/econ-ark/HARK/pull/1430) -- Adds `HARK.algos.vbi` as a general algorithm for solving the optimization step of an agent's problem. -### Minor Changes +#### Minor Changes - Add option to pass pre-built grid to `LinearFast`. [#1388](https://github.com/econ-ark/HARK/pull/1388) - Moves calculation of stable points out of ConsIndShock solver, into method called by post_solve [#1349](https://github.com/econ-ark/HARK/pull/1349) From 0e9300c0047fe337cf1cad5b8320459032087c57 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 4 Jun 2024 17:14:06 -0400 Subject: [PATCH 07/16] first try at adding upper bound information to consumer model --- HARK/model.py | 7 ++++--- HARK/models/consumer.py | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/HARK/model.py b/HARK/model.py index 7fea2c991..30f31bc59 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -31,12 +31,13 @@ class Control: Parameters ---------- - args : list of str + iset : list of str The labels of the variables that are in the information set of this control. """ - def __init__(self, args): - pass + def __init__(self, iset, upper_bound=None): + self.iset = iset + self.upper_bound = upper_bound def discretized_shock_dstn(shocks, disc_params): diff --git a/HARK/models/consumer.py b/HARK/models/consumer.py index bee7b102a..245405403 100644 --- a/HARK/models/consumer.py +++ b/HARK/models/consumer.py @@ -52,7 +52,7 @@ "dynamics": { "b": lambda k, R, PermGroFac: k * R / PermGroFac, "m": lambda b, theta: b + theta, - "c": Control(["m"]), + "c": Control(["m"], upper_bound=lambda m: m), "a": lambda m, c: m - c, }, "reward": {"u": lambda c, CRRA: c ** (1 - CRRA) / (1 - CRRA)}, From 1770ab0f2e13f292e46eaf3e178df9ab5dc6b9a8 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 5 Jun 2024 16:40:53 -0400 Subject: [PATCH 08/16] using lower and upper bounds given to Control in vbi algorithm; passing test on normalized consumption block --- HARK/algos/tests/test_vbi.py | 18 ++++++++++--- HARK/algos/vbi.py | 30 ++++++++++++++------- HARK/model.py | 24 ++++++++++++----- HARK/models/consumer.py | 2 +- HARK/models/perfect_foresight.py | 2 +- HARK/models/perfect_foresight_normalized.py | 2 +- 6 files changed, 57 insertions(+), 21 deletions(-) diff --git a/HARK/algos/tests/test_vbi.py b/HARK/algos/tests/test_vbi.py index a698bacec..becdac17a 100644 --- a/HARK/algos/tests/test_vbi.py +++ b/HARK/algos/tests/test_vbi.py @@ -1,8 +1,8 @@ import HARK.algos.vbi as vbi from HARK.distribution import Bernoulli from HARK.model import Control, DBlock +import HARK.models.consumer as cons import numpy as np - import unittest @@ -14,7 +14,7 @@ }, "dynamics": { "m": lambda y, coin: y + coin, - "c": Control(["m"]), + "c": Control(["m"], lower_bound=lambda m: 0, upper_bound=lambda m: m), "a": lambda m, c: m - c, }, "reward": {"u": lambda c: 1 - (c - 1) ** 2}, @@ -29,6 +29,18 @@ class test_vbi(unittest.TestCase): def test_solve_block_1(self): state_grid = {"m": np.linspace(0, 2, 10)} - dr, dec_vf, arr_vf = vbi.vbi_solve(block_1, lambda a: a, state_grid) + dr, dec_vf, arr_vf = vbi.solve(block_1, lambda a: a, state_grid) self.assertAlmostEqual(dr["c"](**{"m": 1}), 0.5) + + def test_solve_consumption_problem(self): + state_grid = {"m": np.linspace(0, 5, 10)} + + dr, dec_vf, arr_vf = vbi.solve( + cons.consumption_block_normalized, + lambda a: 0, + state_grid, + calibration=cons.calibration, + ) + + self.assertAlmostEqual(dr["c"](**{"m": 1.5}), 1.5) diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index e56e04618..ce5e18ceb 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -4,6 +4,7 @@ """ from HARK.model import DBlock +from inspect import signature import itertools import numpy as np from scipy.optimize import minimize @@ -65,7 +66,7 @@ def grid_to_data_array( return da -def vbi_solve( +def solve( block: DBlock, continuation, state_grid: Grid, disc_params={}, calibration={} ): """ @@ -112,19 +113,30 @@ def negated_value(a): # old! (should be negative) return -srv_function(pre_states, dr) ## get lower bound. - ## not yet implemented - lower_bound = np.array([-1e-12] * len(controls)) ## a really low number! + ## assumes only one control currently + lower_bound = -1e-6 ## a really low number! + feq = block.dynamics[controls[0]].lower_bound + if feq is not None: + lower_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) ## get upper bound - ## not yet implemented - upper_bound = np.array([1e12] * len(controls)) + ## assumes only one control currently + upper_bound = 1e-12 # a very high number + feq = block.dynamics[controls[0]].upper_bound + if feq is not None: + upper_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) # pseudo # optimize_action(pre_states, srv_function) + bounds = ((lower_bound, upper_bound),) + + print(bounds) + res = minimize( # choice of negated_value, - 0, # x0 is starting guess, here arbitrary. + 1, # x0 is starting guess, here arbitrary. + bounds=bounds, ) print(res) @@ -139,11 +151,11 @@ def negated_value(a): # old! (should be negative) ) else: print(f"Optimization failure at {state_vals}.") - print(root_res) + print(res) - dr_best = {c: get_action_rule(root_res[i]) for i, c in enumerate(controls)} + dr_best = {c: get_action_rule(res.x[i]) for i, c in enumerate(controls)} - policy_data.sel(**state_vals).variable.data.put(0, res.root) # ? + policy_data.sel(**state_vals).variable.data.put(0, res.x[0]) # ? value_data.sel(**state_vals).variable.data.put( 0, srv_function(pre_states, dr_best) ) diff --git a/HARK/model.py b/HARK/model.py index 30f31bc59..fd5a4a4b3 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -33,10 +33,14 @@ class Control: ---------- iset : list of str The labels of the variables that are in the information set of this control. + + upper_bound : function + An 'equation function' which evaluates to the upper bound of the control variable. """ - def __init__(self, iset, upper_bound=None): + def __init__(self, iset, lower_bound=None, upper_bound=None): self.iset = iset + self.lower_bound = lower_bound self.upper_bound = upper_bound @@ -172,18 +176,26 @@ def get_controls(self): def transition(self, pre, dr): """ - Returns variable values given previous values and decision rule for all controls. + Computes the state variables following pre-given states, + given a decision rule for all controls. """ dyn = self.dynamics.copy() - # don't simulate values that have already been given. + # don't simulate any states that are logically prior + # to those that have already been given. + met_pre = False # this is a hack; really should use dependency graph + for varn in list(dyn.keys()): + if not met_pre: + if varn in pre: + met_pre = True + del dyn[varn] + elif varn not in pre and varn not in dr: + del dyn[varn] + # this will break if there's a directly recursive label, # i.e. if dynamics at time t for variable 'a' # depend on state of 'a' at time t-1 # This is a forbidden case in CDC's design. - for varn in pre: - if varn in dyn: - del dyn[varn] return simulate_dynamics(dyn, pre, dr) diff --git a/HARK/models/consumer.py b/HARK/models/consumer.py index 245405403..747bd6957 100644 --- a/HARK/models/consumer.py +++ b/HARK/models/consumer.py @@ -34,7 +34,7 @@ "b": lambda k, R: k * R, "y": lambda p, theta: p * theta, "m": lambda b, y: b + y, - "c": Control(["m"]), + "c": Control(["m"], upper_bound=lambda m: m), "p": lambda PermGroFac, p: PermGroFac * p, "a": lambda m, c: m - c, }, diff --git a/HARK/models/perfect_foresight.py b/HARK/models/perfect_foresight.py index 22b02ad83..149b80d27 100644 --- a/HARK/models/perfect_foresight.py +++ b/HARK/models/perfect_foresight.py @@ -24,7 +24,7 @@ "dynamics": { "y": lambda p: p, "m": lambda Rfree, a, y: Rfree * a + y, - "c": Control(["m"]), + "c": Control(["m"], upper_bound=lambda m: m), "p": lambda PermGroFac, p: PermGroFac * p, "a": lambda m, c: m - c, }, diff --git a/HARK/models/perfect_foresight_normalized.py b/HARK/models/perfect_foresight_normalized.py index 544a14fbc..841ad6fbe 100644 --- a/HARK/models/perfect_foresight_normalized.py +++ b/HARK/models/perfect_foresight_normalized.py @@ -26,7 +26,7 @@ "r_eff": lambda Rfree, PermGroFac: Rfree / PermGroFac, "b_nrm": lambda r_eff, a_nrm: r_eff * a_nrm, "m_nrm": lambda b_nrm: b_nrm + 1, - "c_nrm": Control(["m_nrm"]), + "c_nrm": Control(["m_nrm"], upper_bound=lambda m_nrm: m_nrm), "a_nrm": lambda m_nrm, c_nrm: m_nrm - c_nrm, }, "reward": {"u": lambda c: c ** (1 - CRRA) / (1 - CRRA)}, From aba2c556beffd6e28c23fe84f6331b9cb53ae392 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 5 Jun 2024 16:47:30 -0400 Subject: [PATCH 09/16] passing disc_params to consumption block vbi solver test --- HARK/algos/tests/test_vbi.py | 1 + 1 file changed, 1 insertion(+) diff --git a/HARK/algos/tests/test_vbi.py b/HARK/algos/tests/test_vbi.py index becdac17a..87e55d436 100644 --- a/HARK/algos/tests/test_vbi.py +++ b/HARK/algos/tests/test_vbi.py @@ -40,6 +40,7 @@ def test_solve_consumption_problem(self): cons.consumption_block_normalized, lambda a: 0, state_grid, + disc_params = {"theta": {"N" : 7}}, calibration=cons.calibration, ) From e67662795a5063e08e191ac9d7eaff8fca0bba46 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 5 Jun 2024 16:51:02 -0400 Subject: [PATCH 10/16] ruff fixes --- HARK/algos/tests/test_vbi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HARK/algos/tests/test_vbi.py b/HARK/algos/tests/test_vbi.py index 87e55d436..7bf77e9bb 100644 --- a/HARK/algos/tests/test_vbi.py +++ b/HARK/algos/tests/test_vbi.py @@ -40,7 +40,7 @@ def test_solve_consumption_problem(self): cons.consumption_block_normalized, lambda a: 0, state_grid, - disc_params = {"theta": {"N" : 7}}, + disc_params={"theta": {"N": 7}}, calibration=cons.calibration, ) From 3e90f64de33279dba9e227fc64f4892efdf92781 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 5 Jun 2024 17:03:18 -0400 Subject: [PATCH 11/16] fixing the general model tests -- 'screening' transition is optional --- HARK/algos/vbi.py | 4 +++- HARK/model.py | 39 +++++++++++++++++++++------------------ 2 files changed, 24 insertions(+), 19 deletions(-) diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index ce5e18ceb..d453d5a00 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -86,7 +86,9 @@ def solve( """ # state-rule value function - srv_function = block.get_state_rule_value_function_from_continuation(continuation) + srv_function = block.get_state_rule_value_function_from_continuation( + continuation, screen=True + ) controls = block.get_controls() diff --git a/HARK/model.py b/HARK/model.py index fd5a4a4b3..6fb862684 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -174,28 +174,29 @@ def get_controls(self): return [varn for varn in dyn if isinstance(dyn[varn], Control)] - def transition(self, pre, dr): + def transition(self, pre, dr, screen=False): """ Computes the state variables following pre-given states, given a decision rule for all controls. """ dyn = self.dynamics.copy() - # don't simulate any states that are logically prior - # to those that have already been given. - met_pre = False # this is a hack; really should use dependency graph - for varn in list(dyn.keys()): - if not met_pre: - if varn in pre: - met_pre = True - del dyn[varn] - elif varn not in pre and varn not in dr: - del dyn[varn] - - # this will break if there's a directly recursive label, - # i.e. if dynamics at time t for variable 'a' - # depend on state of 'a' at time t-1 - # This is a forbidden case in CDC's design. + if screen: + # don't simulate any states that are logically prior + # to those that have already been given. + met_pre = False # this is a hack; really should use dependency graph + for varn in list(dyn.keys()): + if not met_pre: + if varn in pre: + met_pre = True + del dyn[varn] + elif varn not in pre and varn not in dr: + del dyn[varn] + + # this will break if there's a directly recursive label, + # i.e. if dynamics at time t for variable 'a' + # depend on state of 'a' at time t-1 + # This is a forbidden case in CDC's design. return simulate_dynamics(dyn, pre, dr) @@ -211,7 +212,9 @@ def calc_reward(self, vals): return rvals - def get_state_rule_value_function_from_continuation(self, continuation): + def get_state_rule_value_function_from_continuation( + self, continuation, screen=False + ): """ Given a continuation value function, returns a state-rule value function: the value for each state and decision rule. @@ -220,7 +223,7 @@ def get_state_rule_value_function_from_continuation(self, continuation): """ def state_rule_value_function(pre, dr): - vals = self.transition(pre, dr) + vals = self.transition(pre, dr, screen=screen) r = list(self.calc_reward(vals).values())[0] # a hack; to be improved cv = continuation( *[vals[var] for var in signature(continuation).parameters] From 32eb4c11e9004dbe2ca8fb4d95e276ab03d47044 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 5 Jun 2024 17:04:55 -0400 Subject: [PATCH 12/16] documentating optional screening of transition function --- HARK/model.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/HARK/model.py b/HARK/model.py index 6fb862684..00fe2b2f1 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -178,6 +178,15 @@ def transition(self, pre, dr, screen=False): """ Computes the state variables following pre-given states, given a decision rule for all controls. + + Parameters + ----------- + pre + dr + + screen: Boolean + If True, the remove any dynamics that are prior to the first given state. + Defaults to False. """ dyn = self.dynamics.copy() From 6ec5d5eb7fa220375cfccfe8637910cf39969e25 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 18 Jun 2024 15:59:53 -0400 Subject: [PATCH 13/16] touch ups --- HARK/algos/vbi.py | 3 --- HARK/models/consumer.py | 6 ++++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index d453d5a00..a0a215b06 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -133,14 +133,11 @@ def negated_value(a): # old! (should be negative) bounds = ((lower_bound, upper_bound),) - print(bounds) - res = minimize( # choice of negated_value, 1, # x0 is starting guess, here arbitrary. bounds=bounds, ) - print(res) dr_best = {c: get_action_rule(res.x[i]) for i, c in enumerate(controls)} diff --git a/HARK/models/consumer.py b/HARK/models/consumer.py index 747bd6957..3900d08a1 100644 --- a/HARK/models/consumer.py +++ b/HARK/models/consumer.py @@ -84,5 +84,7 @@ } ) -cons_problem = RBlock(blocks=[consumption_block, tick_block]) -cons_portfolio_problem = RBlock(blocks=[consumption_block, portfolio_block, tick_block]) +cons_problem = RBlock(blocks=[consumption_block_normalized, tick_block]) +cons_portfolio_problem = RBlock( + blocks=[consumption_block_normalized, portfolio_block, tick_block] +) From 7040c44d172eeee616f58eba9defac010bde2fc0 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 18 Jun 2024 18:21:32 -0400 Subject: [PATCH 14/16] don't need to initialize income in normalized models --- HARK/models/test_models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/HARK/models/test_models.py b/HARK/models/test_models.py index 9cb9152d1..1a51d2cde 100644 --- a/HARK/models/test_models.py +++ b/HARK/models/test_models.py @@ -103,7 +103,7 @@ def setUp(self): { # initial states "k": Lognormal(-6, 0), #'live' : 1, - "p": 1.0, + # "p": 1.0, - not needed, using normalized problem }, agent_count=2, T_sim=5, @@ -121,7 +121,7 @@ def setUp(self): { # initial states "k": Lognormal(-6, 0), #'live' : 1, - "p": 1.0, + # "p": 1.0,- not needed, using normalized problem "R": 1.03, }, agent_count=2, From 61fd2f73df611c3b8d704f5ad387840a8f849c3c Mon Sep 17 00:00:00 2001 From: sb Date: Thu, 14 Nov 2024 10:42:13 -0500 Subject: [PATCH 15/16] vbi algo handling cases where stage has 0 controls --- HARK/algos/tests/test_vbi.py | 24 +++++++++++ HARK/algos/vbi.py | 83 +++++++++++++++++++----------------- HARK/model.py | 1 + 3 files changed, 68 insertions(+), 40 deletions(-) diff --git a/HARK/algos/tests/test_vbi.py b/HARK/algos/tests/test_vbi.py index 7bf77e9bb..ff15a76be 100644 --- a/HARK/algos/tests/test_vbi.py +++ b/HARK/algos/tests/test_vbi.py @@ -21,6 +21,21 @@ } ) +block_2 = DBlock( # has no control variable + **{ + "name": "vbi_test_1", + "shocks": { + "coin": Bernoulli(p=0.5), + }, + "dynamics": { + "m": lambda y, coin: y + coin, + "a": lambda m: m - 1, + }, + "reward": {"u": lambda m: 0}, + } +) + + class test_vbi(unittest.TestCase): # def setUp(self): @@ -33,6 +48,15 @@ def test_solve_block_1(self): self.assertAlmostEqual(dr["c"](**{"m": 1}), 0.5) + def test_solve_block_2(self): + # no control variable case. + state_grid = {"m": np.linspace(0, 2, 10)} + + dr, dec_vf, arr_vf = vbi.solve(block_2, lambda a: a, state_grid) + + # arrival value function gives the correct expect value of continuation + self.assertAlmostEqual(arr_vf({"y": 10}), 9.5) + def test_solve_consumption_problem(self): state_grid = {"m": np.linspace(0, 5, 10)} diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index a0a215b06..9d06712d3 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -114,50 +114,53 @@ def negated_value(a): # old! (should be negative) # negative, for minimization later return -srv_function(pre_states, dr) - ## get lower bound. - ## assumes only one control currently - lower_bound = -1e-6 ## a really low number! - feq = block.dynamics[controls[0]].lower_bound - if feq is not None: - lower_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) - - ## get upper bound - ## assumes only one control currently - upper_bound = 1e-12 # a very high number - feq = block.dynamics[controls[0]].upper_bound - if feq is not None: - upper_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) - - # pseudo - # optimize_action(pre_states, srv_function) - - bounds = ((lower_bound, upper_bound),) - - res = minimize( # choice of - negated_value, - 1, # x0 is starting guess, here arbitrary. - bounds=bounds, - ) - - dr_best = {c: get_action_rule(res.x[i]) for i, c in enumerate(controls)} - - if res.success: - policy_data.sel(**state_vals).variable.data.put( - 0, res.x[0] - ) # will only work for scalar actions - value_data.sel(**state_vals).variable.data.put( - 0, srv_function(pre_states, dr_best) + if len(controls) == 0: + # if no controls, no optimization is necessary + pass + elif len(controls) == 1: + ## get lower bound. + ## assumes only one control currently + lower_bound = -1e-6 ## a really low number! + feq = block.dynamics[controls[0]].lower_bound + if feq is not None: + lower_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) + + ## get upper bound + ## assumes only one control currently + upper_bound = 1e-12 # a very high number + feq = block.dynamics[controls[0]].upper_bound + if feq is not None: + upper_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) + + bounds = ((lower_bound, upper_bound),) + + res = minimize( # choice of + negated_value, + 1, # x0 is starting guess, here arbitrary. + bounds=bounds, ) - else: - print(f"Optimization failure at {state_vals}.") - print(res) dr_best = {c: get_action_rule(res.x[i]) for i, c in enumerate(controls)} - policy_data.sel(**state_vals).variable.data.put(0, res.x[0]) # ? - value_data.sel(**state_vals).variable.data.put( - 0, srv_function(pre_states, dr_best) - ) + if res.success: + policy_data.sel(**state_vals).variable.data.put( + 0, res.x[0] + ) # will only work for scalar actions + value_data.sel(**state_vals).variable.data.put( + 0, srv_function(pre_states, dr_best) + ) + else: + print(f"Optimization failure at {state_vals}.") + print(res) + + dr_best = {c: get_action_rule(res.x[i]) for i, c in enumerate(controls)} + + policy_data.sel(**state_vals).variable.data.put(0, res.x[0]) # ? + value_data.sel(**state_vals).variable.data.put( + 0, srv_function(pre_states, dr_best) + ) + elif len(controls) > 1: + raise Exception(f"Value backup iteration is not yet implemented for stages with {len(controls)} > 1 control variables.") # use the xarray interpolator to create a decision rule. dr_from_data = { diff --git a/HARK/model.py b/HARK/model.py index 191546801..c281f42f0 100644 --- a/HARK/model.py +++ b/HARK/model.py @@ -343,6 +343,7 @@ def get_state_rule_value_function_from_continuation( def state_rule_value_function(pre, dr): vals = self.transition(pre, dr, screen=screen) r = list(self.calc_reward(vals).values())[0] # a hack; to be improved + # this assumes a single reward variable; instead, a named could be passed in. cv = continuation( *[vals[var] for var in signature(continuation).parameters] ) From 8f45b797b6a40eec45c59af16dc822d62ce8d906 Mon Sep 17 00:00:00 2001 From: sb Date: Thu, 14 Nov 2024 10:43:22 -0500 Subject: [PATCH 16/16] rufff fixes --- HARK/algos/tests/test_vbi.py | 3 +-- HARK/algos/vbi.py | 14 ++++++++++---- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/HARK/algos/tests/test_vbi.py b/HARK/algos/tests/test_vbi.py index ff15a76be..fa9b5399d 100644 --- a/HARK/algos/tests/test_vbi.py +++ b/HARK/algos/tests/test_vbi.py @@ -21,7 +21,7 @@ } ) -block_2 = DBlock( # has no control variable +block_2 = DBlock( # has no control variable **{ "name": "vbi_test_1", "shocks": { @@ -36,7 +36,6 @@ ) - class test_vbi(unittest.TestCase): # def setUp(self): # pass diff --git a/HARK/algos/vbi.py b/HARK/algos/vbi.py index 9d06712d3..3deba378b 100644 --- a/HARK/algos/vbi.py +++ b/HARK/algos/vbi.py @@ -114,7 +114,7 @@ def negated_value(a): # old! (should be negative) # negative, for minimization later return -srv_function(pre_states, dr) - if len(controls) == 0: + if len(controls) == 0: # if no controls, no optimization is necessary pass elif len(controls) == 1: @@ -123,14 +123,18 @@ def negated_value(a): # old! (should be negative) lower_bound = -1e-6 ## a really low number! feq = block.dynamics[controls[0]].lower_bound if feq is not None: - lower_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) + lower_bound = feq( + *[pre_states[var] for var in signature(feq).parameters] + ) ## get upper bound ## assumes only one control currently upper_bound = 1e-12 # a very high number feq = block.dynamics[controls[0]].upper_bound if feq is not None: - upper_bound = feq(*[pre_states[var] for var in signature(feq).parameters]) + upper_bound = feq( + *[pre_states[var] for var in signature(feq).parameters] + ) bounds = ((lower_bound, upper_bound),) @@ -160,7 +164,9 @@ def negated_value(a): # old! (should be negative) 0, srv_function(pre_states, dr_best) ) elif len(controls) > 1: - raise Exception(f"Value backup iteration is not yet implemented for stages with {len(controls)} > 1 control variables.") + raise Exception( + f"Value backup iteration is not yet implemented for stages with {len(controls)} > 1 control variables." + ) # use the xarray interpolator to create a decision rule. dr_from_data = {