From 41d31684f1bb8a898531f4d6dc8db4b8721b1c80 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 12 Jun 2023 13:52:35 -0400 Subject: [PATCH 01/14] SolvingMicroDSOPs vector data and stub tests for FOC and EGM algorithms --- HARK/algos/tests/smdsops_aVec.npy | Bin 0 -> 168 bytes HARK/algos/tests/smdsops_cVec2.npy | Bin 0 -> 168 bytes HARK/algos/tests/smdsops_cVec_egm.npy | Bin 0 -> 168 bytes HARK/algos/tests/smdsops_mVec.npy | Bin 0 -> 168 bytes HARK/algos/tests/smdsops_mVec_egm.npy | Bin 0 -> 168 bytes HARK/algos/tests/test_algos.py | 38 ++++++++++++++++++++++++++ 6 files changed, 38 insertions(+) create mode 100644 HARK/algos/tests/smdsops_aVec.npy create mode 100644 HARK/algos/tests/smdsops_cVec2.npy create mode 100644 HARK/algos/tests/smdsops_cVec_egm.npy create mode 100644 HARK/algos/tests/smdsops_mVec.npy create mode 100644 HARK/algos/tests/smdsops_mVec_egm.npy create mode 100644 HARK/algos/tests/test_algos.py diff --git a/HARK/algos/tests/smdsops_aVec.npy b/HARK/algos/tests/smdsops_aVec.npy new file mode 100644 index 0000000000000000000000000000000000000000..4db49809a1da1fffa95fccca99216184063901e6 GIT binary patch literal 168 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= jXCxM+0{I%II+{8PwF(pfE(R3v!5+$WfY2OJTEGDS!Sx+y literal 0 HcmV?d00001 diff --git a/HARK/algos/tests/smdsops_cVec2.npy b/HARK/algos/tests/smdsops_cVec2.npy new file mode 100644 index 0000000000000000000000000000000000000000..723cb842cdb932bd26d7b74af12c43bcb7bc18ef GIT binary patch literal 168 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I%II+{8PwF(pfuDKGvk$r8a?R)YMKVUic%3ef?Uu{CrH~Z3O3iX2D|Jlb0 LFSU8OkJ$kL5XmiX literal 0 HcmV?d00001 diff --git a/HARK/algos/tests/smdsops_cVec_egm.npy b/HARK/algos/tests/smdsops_cVec_egm.npy new file mode 100644 index 0000000000000000000000000000000000000000..ece9ee4eb9475fa45f587ce939f784876c792a4e GIT binary patch literal 168 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I%II+{8PwF(pfF5^63<9?0j_9B%Lx*uHr+1Ez$Ju%wH?jW@AsnlOF0f*9k LAJezEiZ}oO<$cIAm9K1Pg*Sh literal 0 HcmV?d00001 diff --git a/HARK/algos/tests/smdsops_mVec_egm.npy b/HARK/algos/tests/smdsops_mVec_egm.npy new file mode 100644 index 0000000000000000000000000000000000000000..03ac5bf443d1245f98d54862ecfa1537fe9cd437 GIT binary patch literal 168 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I%II+{8PwF(pfF5^63<9?0j_FXn>#44279ro$47fIe1b`V 0)) + +class egm_test(unittest.TestCase): + """ + cVec_egm = egm(aVec) (machine precision) + pi(mVec_egm) == cVec_egm + """ + def setUp(self): + self.aVec = np.load("smdsops_aVec.npy") + self.cVec_egm = np.load("smdsops_cVec_egm.npy") + self.mVec_egm = np.load("smdsops_mVec_egm.npy") + + def test_egm(self): + + self.assertTrue(np.all(self.aVec + self.cVec_egm == self.mVec_egm)) \ No newline at end of file From a1390fe739e4730c93e0ba47e612123040b028d5 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 13 Jun 2023 14:11:24 -0400 Subject: [PATCH 02/14] porting in code from stages --- HARK/algos/egm.py | 125 ++++++++++++++++++ HARK/algos/foc.py | 234 +++++++++++++++++++++++++++++++++ HARK/algos/tests/test_algos.py | 1 + 3 files changed, 360 insertions(+) create mode 100644 HARK/algos/egm.py create mode 100644 HARK/algos/foc.py diff --git a/HARK/algos/egm.py b/HARK/algos/egm.py new file mode 100644 index 000000000..93ce77c1c --- /dev/null +++ b/HARK/algos/egm.py @@ -0,0 +1,125 @@ + + +def analytic_pi_star_y(self, + y, + v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0 + ): + """ + The optimal action that results in output values y. + + Available only with reward_der_inv and transition_inv + are well-defined. + """ + if self.reward_der_inv is None: + raise Exception("No inverse marginal reward function found.") + + if self.transition_inv is None: + raise Exception("No inverse transition function found. ") + + if self.transition_der_a is None or not( + isinstance(self.transition_der_a, float) or isinstance(self.transition_der_a, int) + ): + raise Exception(f"No constant transition derivative found. transition_der_a is {self.transition_der_a}") + + if not isinstance(self.discount, float) or isinstance(self.discount, int): + raise Exception("Analytic pi_star_y requires constant discount factor (rendering B' = 0).") + + ### TEST: available T_der as constant. + + v_y_der_at_y = v_y_der(y) + + if isinstance(v_y_der_at_y, xr.DataArray): + v_y_der_at_y = v_y_der_at_y.values # np.atleast1darray() ? + + + if 0 > v_y_der_at_y: + raise Exception(f"Negative marginal value {v_y_der_at_y} computes at y value of {y}. Reward is {- self.discount * self.transition_der_a * v_y_der_at_y}") + + return self.reward_der_inv(- self.discount * self.transition_der_a * v_y_der_at_y) + +def optimal_policy_egm(self, + y_grid : Mapping[str, Sequence] = {}, + v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, + ): + """ + Given a grid over output + and marginal output value function, + compute the optimal action. + + ## NO SHOCKS ALLOWED ## + + This depends on the stage having an + *inverse marginal reward function* + and *inverse transition function*. + + Does not use rootfinding! + The 'grid' over the input + + ### ASSUMES: No discounting this phase, + ### and... T' = -1 ??? + """ + + if self.reward_der_inv is None: + raise Exception("No inverse marginal reward function found. EGM requires reward_der_inv defined for this stage.") + + if self.transition_inv is None: + raise Exception("No inverse transition function found. EGM requires transition_inv defined for this stage.") + + pi_y_data = xr.DataArray( + np.zeros([len(v) for v in y_grid.values()]), + dims = y_grid.keys(), + coords = y_grid + ) + + # Collecting data for the real optimal policy with respect to inputs + x_val_data = [] + a_val_data = [] + + for y_point in itertools.product(*y_grid.values()): + y_vals = {k : v for k, v in zip(y_grid.keys() , y_point)} + + acts = self.analytic_pi_star_y(y_vals, v_y_der) + + pi_y_data.sel(**y_vals).variable.data.put(0, acts) + + x_vals = self.transition_inv(y_vals, self.action_zip(acts)) + + x_val_data.append(x_vals) + a_val_data.append(self.action_zip(acts)) + + ## TODO is this dealing with repeated values? + x_coords = { + x : np.array([xv[x] for xv in x_val_data]) + for x + in self.inputs + } + + for x_label in self.solution_points.coords: + for sol_x_val in self.solution_points.coords[x_label]: + if sol_x_val.values not in x_coords[x_label]: + ii = np.searchsorted(x_coords[x_label], sol_x_val) + x_coords[x_label] = np.insert(x_coords[x_label], ii, sol_x_val) + + pi_data = xr.DataArray( + np.zeros([len(v) for v in x_coords.values()]), + dims = x_coords.keys(), + coords = x_coords + ) + + if 'pi*' in self.solution_points: + for x_point in itertools.product(*x_coords.values()): + x_vals = {k : v for k, v in zip(x_coords.keys() , x_point)} + + if label_index_in_dataset(x_vals, self.solution_points['pi*']): + acts = np.atleast_1d(self.solution_points['pi*'].sel(x_vals)) + + pi_data.sel(**x_vals, **{}).variable.data.put(0, acts) + + for i, x_vals in enumerate(x_val_data): + x_vals = x_val_data[i] + a_vals = a_val_data[i] + acts = [a_vals[a] for a in a_vals] + + pi_data.sel(**x_vals).variable.data.put(0, acts) + + return pi_data, pi_y_data diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py new file mode 100644 index 000000000..9a882602d --- /dev/null +++ b/HARK/algos/foc.py @@ -0,0 +1,234 @@ +import numpy as np +from scipy.optimize import minimize, brentq +import xarray as xr + +""" +Sargent + Stachurski: +x - states +z - shocks +a - actions + +New: +y - post-states + +French notation: +g - transtion function x,z, a -> y + + + +Question: + - Do we want to have a canonical 'grid' object that includes variable names and allows for non-cartesian grids? + +""" + +def xndindex(ds, dims=None): + """ + There is currently no integrated way to iterate over an xarray.DataArray with its coordinate labels. + + This method is a workaround from: + https://github.com/pydata/xarray/issues/2805#issuecomment-1255029201 + """ + if dims is None: + dims = ds.dims + elif type(dims) is str: + dims=[dims] + else: + pass + + for d in dims: + if d not in ds.dims: + raise ValueError("Invalid dimension '{}'. Available dimensions {}".format(d, ds.dims)) + + iter_dict = {k:v for k,v in ds.sizes.items() if k in dims} + for d,k in zip(repeat(tuple(iter_dict.keys())),zip(np.ndindex(tuple(iter_dict.values())))): + yield {k:l for k,l in zip(d,k[0])} + + +def xz_grids_to_data_array( + x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. + k_grid : Mapping[str, Sequence] = {} + ): + + coords = {**x_grid, **k_grid} + + da = xr.DataArray( + np.zeros([len(v) for v in coords.values()]), + dims = coords.keys(), + coords = coords + + return + + +def optimal_policy_foc(self, + x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. + z_grid : Mapping[str, Sequence] = {}, + v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, + optimizer_args = None # TODO: For brettq. + ): + """ + Given a grid over input and shock state values, + and marginal output value function, + compute the optimal action. + + Uses root finding and the first order condition. + + This is written with the expectation that: + - Actions are scalar + - the Q function is concave over the action space. + + Functionality is not guaranteed otherwise. + + + NOTE: This does not put pre-defined solution points into the solution data. + That needs to be added in a different step. + + + TODO: K -> Z + TODO: Transition function as an argument + + Parameters + ----------- + + x_grid : Mapping[str, Sequence] + A mapping of state variable names to a sequence of grid values + + z_grid: Mapping[str, Sequence] + A mapping of shock variable names to grid values + + + """ + # Set up data arrays with coordinates based on the grid. + pi_data = xz_grids_to_data_array(x_grid, z_grid) + q_der_data = xz_grids_to_data_array(x_grid, z_grid) + ## May need to expand this for multivariate y + y_data = xz_grids_to_data_array(x_grid, z_grid) + + xz_grid_size = np.prod([len(xv) for xv in x_grid.values()]) * np.prod([len(zv) for zv in z_grid.values()]) + + "dq_da(x,z,a) = dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da(x,z,a)" + + def foc(a): + a_vals = {an : av for an,av in zip(self.actions, (a,))} + return self.dq_da(x_vals, k_vals, a_vals, v_y_der) + + + + # TODO: replace these with iterators.... + #xz_iterator = xndindex(pi_data) + #import pdb; pdb.set_trace() + + + + for x_point in itertools.product(*x_grid.values()): + x_vals = {k : v for k, v in zip(x_grid.keys() , x_point)} + + for z_point in itertools.product(*z_grid.values()): + z_vals = {k : v for k, v in zip(k_grid.keys() , k_point)} + + # repeated code with the other optimizer -- can be functionalized out? + if len(self.actions) == 0: + q_der_xz = self.d_q_d_a( + x_vals, + z_vals, + {}, # no actions + v_y_der + ) + + # this should just be the default value for pi + pi_data.sel(**x_vals, **z_vals).variable.data.put(0, np.nan) + + # then this can be looped in differently. + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) + + ## TODO: define g + y = g(x_vals, z_vals, self.action_zip([np.nan])) + y_n = np.array([y[k] for k in y]) + y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) + + + else: + + # these lower bounds as arugments to the + lower_bound = self.action_lower_bound(x_vals, k_vals) + upper_bound = self.action_upper_bound(x_vals, k_vals) + + ## what if no lower bound? + q_der_lower = None + if lower_bound[0] is not None: + q_der_lower = self.d_q_d_a( + x_vals, + k_vals, + self.action_zip(lower_bound), + v_y_der + ) + else: + lower_bound[0] = 1e-12 ## a really high number! + + q_der_upper = None + if upper_bound[0] is not None: + q_der_upper = self.d_q_d_a( + x_vals, + k_vals, + self.action_zip(upper_bound), + v_y_der + ) + else: + upper_bound[0] = 1e12 ## a really high number! + + ## TODO: Better handling of case when there is a missing bound? + if q_der_lower is not None and q_der_upper is not None and q_der_lower < 0 and q_der_upper < 0: + a0 = lower_bound + + pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) + q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_lower) + elif q_der_lower is not None and q_der_upper is not None and q_der_lower > 0 and q_der_upper > 0: + a0 = upper_bound + + pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) + q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_upper) + else: + ## Better exception handling here + ## asserting that Q is concave + if q_der_lower is not None and q_der_upper is not None and not(q_der_lower > 0 and q_der_upper < 0): + raise Exception("Cannot solve for optimal policy with FOC if Q is not concave!") + + a0, root_res = brentq( + foc, + lower_bound[0], # only works with scalar actions + upper_bound[0], # only works with scalar actions + full_output = True + ) + + if root_res.converged: + pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) + + q_der_xz = self.dq_da( + x_vals, + z_vals, + self.action_zip((a0,)), # actions are scalar + v_y_der + ) + + q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_xk) + else: + print(f"Rootfinding failure at {x_vals}, {z_vals}.") + print(root_res) + + #raise Exception("Failed to optimize.") + pi_data.sel(**x_vals, **z_vals).variable.data.put(0, root_res.root) + + q_der_xz = dq_da( + x_vals, + z_vals, + self.action_zip((root_res.root,)), # actions are scalar + v_y_der + ) + + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) + + acts = np.atleast_1d(pi_data.sel(**x_vals, **z_vals).values) + y = g(x_vals, z_vals, self.action_zip(acts)) + y_n = np.array([y[k] for k in y]) + y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) + + return pi_data, q_der_data, y_data \ No newline at end of file diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 510bb0c43..67dd82ba0 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -4,6 +4,7 @@ # Bring in modules we need import unittest +import HARK.algos.foc import numpy as np class foc_test(unittest.TestCase): From 72d9e755faff4b1c43623fb17fb6bd807e878592 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 13 Jun 2023 16:38:06 -0400 Subject: [PATCH 03/14] adding gothic class for test and basic problem definition code (from 'stage' earlier) --- HARK/algos/foc.py | 4 +- HARK/algos/tests/gothic_class.py | 289 +++++++++++++++++++++++++++++++ HARK/algos/tests/test_algos.py | 82 ++++++++- 3 files changed, 372 insertions(+), 3 deletions(-) create mode 100644 HARK/algos/tests/gothic_class.py diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index 9a882602d..f8bdbfccd 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -1,5 +1,6 @@ import numpy as np from scipy.optimize import minimize, brentq +from typing import Callable, Mapping, Sequence import xarray as xr """ @@ -55,8 +56,9 @@ def xz_grids_to_data_array( np.zeros([len(v) for v in coords.values()]), dims = coords.keys(), coords = coords + ) - return + return da def optimal_policy_foc(self, diff --git a/HARK/algos/tests/gothic_class.py b/HARK/algos/tests/gothic_class.py new file mode 100644 index 000000000..0340d7a65 --- /dev/null +++ b/HARK/algos/tests/gothic_class.py @@ -0,0 +1,289 @@ +# --- +# ## A Description of the Gothic class +# This is the Gothic class from SolvingMicroDSOPs. +# It contains gold standard functions that can be used for testing. +# +# The Gothic class to define the end-of-period "gothic" functions: $\mathfrak{v}$, $\mathfrak{v}'$, and $\mathfrak{c}$, as well as the interpolations of each of these functions. +# +# Defining these in one class allows us to bundle the parameters for the problem in one place, and then hide them from the user. We have likewise bundled the parameters for the utility function and the discrete distribution approximation in their own classes. The class structure additionally allows us to bundle useful fucnitonality with the utility function and discrete distribution, such as the marginal utility in the utility class, and the expectation operator associated with the discrete distribution. The layers of abstraction provided by the object-oriented framework allow us to use the bare minimum additional parameters for each level of the code. See the notebook regarding these classes for further explanation. +# +# We define a Gothic object with a utility function $u, \beta,$ the risk parameter $\rho, \gamma,$ R, and $\theta$. +# +# Once initilized, we will have access to these methods in the Gothic class: +# +# +# V_Tminus1: GothicV at {T-1}, in levels +# +# VP_Tminus1: GothicV\' at {T-1}, in marginal values +# +# V_Tminus1_interp, and +# VP_Tminus1_interp: Both above, interpolated on an a-grid +# +# Usage: +# +# Gothic.V_Tminus1(a): Return the gothicV value for a at T-1. +# Gothic.VP_Tminus1(a): Return the gothicV\' value for a at T-1. +# +# Gothic.V_Tminus1_interp(a_grid): Return gothicV(a) as an interpolated +# function, interpolated on the a_grid +# provided. +# Gothic.VP_Tminus1_interp(a_grid): As above, return gothicV\'(a) as +# an interpolation function on the +# a_grid. +# +# Gothic.C_Tminus1(a): Return the gothicC value for a at T-1. +# +# Gothic.C_Tminus1_interp(a_grid): Return gothicC(a) as an interpolated +# function, interpolated on the a_grid +# provided. +# +# ## The Gothic class: + +from __future__ import division +from scipy.interpolate import InterpolatedUnivariateSpline +import numpy as np + +class Gothic: + + def __init__(self, u, beta, rho, Gamma, R, Income, variable_variance=False): + """ + Initialize a Gothic object. + + Args: + u (object): Utility function. Should accept a real number & have + a "prime" method which is the first derivative. + beta (float): Time discount factor. + rho (float): Risk aversion. + gamma (array): Array of gamma values, time series indexed by t. + R (float): The real return factor. Fixed in time. + Income (object): Approximated distribution for a two-shock method. + Must have method: "Income.E()." + NOTE: The convention is that permanent shock to + incom (psi) comes first, and the temporary shock + (eta) comes second in the ordered pair of the + shocks to income. Any function of which we need + to find an expectation, with respect to income, + should be defined as such. + variable_variance (boolean): If true, the Income is a list of + income objects. + Returns: + Nothing. + Raises: + [] + """ + self.u = u + self.beta = beta + self.rho = rho + self.Gamma = Gamma + self.Gamma_to_1minusRho = Gamma ** (1.0 - rho) # Define here once. + self.Gamma_to_minusRho = Gamma ** (-rho) # Define here once. + self.R = R + self.Income = Income + self.variable_variance = variable_variance + + + + def V(self, a, t=-1, v_prime=None): + """ + Given an end-of-period a value, return the GothicV_{T-1} value. + For t = None, implements equation (22) from MicroDSOP: value function at T-1 + For t != None, v_prime != None, implements equation (17) from MicroDSOP. + """ + + # Define function describing tomorrow: + if t == -1: + tp1 = -1 # Selects final value in a vector. + t = -2 + V_func = lambda tinc_shk: self.u(self.R/(self.Gamma[tp1]) * a + tinc_shk) + elif v_prime is not None: + tp1 = t + 1 + V_func = lambda tinc_shk: v_prime(self.R/(self.Gamma[tp1]) * a + tinc_shk) + else: + raise Exception("Please either specify that t=-1 (indicating solution for period T-1) or specify *both* t and v_prime.") + + if self.variable_variance: + gothicV = self.beta * self.Gamma_to_1minusRho[tp1] * self.Income[tp1].E(V_func) + # TODO: confirm that + else: + gothicV = self.beta * self.Gamma_to_1minusRho[tp1] * self.Income.E(V_func) + + return(gothicV) + + + def V_prime(self, a, t=-1, c_prime=None): + """ + Given an end-of-period a-value, return the GothicV_prime value. + If t=-1, return T-1 value; else return the t-value. + + This implements equation (19) and (30) from MicroDSOP for T-1, and + equation (18) for all previous time periods. + """ + + if t == -1: + tp1 = -1 # Selects final value in a vector. + t = -2 + Vp_func = lambda tinc_shk: psi**(-self.rho) * self.u.prime(self.R/(self.Gamma[tp1]) * a + tinc_shk) + elif c_prime is not None: + tp1 = t+1 + + #mtp1 = self.R/(self.Gamma[tp1]*psi) * a + eta + #print "mtp1", mtp1 + #g = lambda psi, eta: psi**(-self.rho) * self.u.prime(c_prime(mtp1)) + # one possible solution: + Vp_func = lambda tinc_shk, R=self.R, gamma=self.Gamma[tp1],aa=a, rho=self.rho, uP=self.u.prime, ctp1=c_prime: uP(ctp1(R/(gamma) * aa + tinc_shk)) + + else: + raise Exception("Please either specify that t=-1 (indicating solution for period T-1) or specify *both* t and c_prime.") + + if self.variable_variance: + gothicV_prime = self.beta * self.R * self.Gamma_to_minusRho[tp1] * self.Income[tp1].E(Vp_func) + else: + gothicV_prime = self.beta * self.R * self.Gamma_to_minusRho[tp1] * self.Income.E(Vp_func) + + return(gothicV_prime) + + + def C(self, a, t=-1, c_prime=None): + """ + Return the gothicC value for a. If t=-1, return the value for T-1. + + Implements equation (34) in MicroDSOP for T-1; implements equation (20) + for all other periods. + """ + + if t == -1: + scriptC = self.V_prime(a,t=-1)**(-1.0/self.rho) + elif c_prime is not None: + scriptC = self.V_prime(a, t=t, c_prime=c_prime)**(-1.0/self.rho) + else: + raise Exception("Please either specify that t=-1 (indicating solution for period T-1) or specify *both* t and c_prime.") + + return(scriptC) + + # copied from ./Code/Python/active_development/archive/Gothic Class 1shock.ipynb + def C_Tminus1(self, a): + """ + Return the gothicC value for a at T-1. Equation (34) in MicroDSOP. + """ + return self.VP_Tminus1(a)**(-1.0/self.rho) + + # copied from ./Code/Python/active_development/archive/Gothic Class 1shock.ipynb + # changed Theta -> Income + def VP_Tminus1(self, a): + """ + Given an end-of-period a-value, return the GothicV_prime_Tminus1 value. + Vectorize to work on a grid. + + This implements function (30) from MicroDSOP. + """ + # Convenience definitions. Note we take the last value of Gamma: + fancyR_T = self.R/self.Gamma[-1] + + Vp_func = lambda tinc_shk: self.u.prime(fancyR_T * a + tinc_shk) + # The value: + GVTm1P = self.beta * self.R * self.Gamma_to_minusRho[-1] * self.Income.E(Vp_func) + + return GVTm1P + + # copied from ./Code/Python/active_development/archive/Gothic Class 1shock.ipynb + # changed Theta -> Income + def C_t(self, a, c_prime, t=None): + """ + Return the gothicC value for a at t. + + This employs Equation (20) in MicroDSOP. + """ + # Quick comparison test against hand-coded equation (76): + + if t is None: + t = -1 + + E_sum = 0.0 + for theta in self.Income.X: + fancyR_tp1 = self.R/self.Gamma[t+1] + c_tp1 = c_prime(fancyR_tp1*a + theta) + + E_sum += c_tp1**(-self.rho) + + alt_scriptC = (self.beta * self.R * (self.Gamma[t+1] ** (-self.rho)) * (1.0/self.Income.N) * E_sum) ** (-1.0/self.rho) + + scriptC = self.VP_t(a, c_prime, t)**(-1.0/self.rho) + + #print "alt_scriptC", alt_scriptC + #print "scriptC", scriptC + + tempdiff = alt_scriptC - scriptC + assert np.abs(tempdiff) < 1e-10, "in Gothic.C_t, manually calculated scriptC(a) != computed scriptC, by this much: " + str(tempdiff) + " values: alt_scriptC: " + str(alt_scriptC) + " scriptC: " + str(scriptC) + + return scriptC + + # copied from ./Code/Python/active_development/archive/Gothic Class 1shock.ipynb + # changed Theta -> Income + def VP_t(self, a, c_prime, t=None): + """ + Given a next-period consumption function, find the Vprime function for this period. + + This implements function (__) from MicroDSOP. + """ + + if t is None: + Gamma_to_mRho = self.Gamma_to_minusRho[0] + scriptR_tp1 = self.R/self.Gamma[0] + else: + Gamma_to_mRho = self.Gamma_to_minusRho[t+1] + scriptR_tp1 = self.R/self.Gamma[t+1] + + Vp_func = lambda tinc_shk: self.u.prime(c_prime(scriptR_tp1 * a + tinc_shk)) + # The value: + GVPt = self.beta * self.R * Gamma_to_mRho * self.Income.E(Vp_func) + + return GVPt +# - + +import numpy as np +import scipy.stats as stats +import pylab as plt +from scipy.optimize import brentq +## TODO Need these resources here also? +from resources import DiscreteApproximation, Utility, DiscreteApproximationTwoMeanOneIndependentLognormalDistribs, DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event + +# General parameters: +rho = 2.0 +beta = 0.96 +gamma = np.array([1.0,1.0,1.0]) # A three-element "time series;" this + # structure needed for gothic class below +R = 1.02 + +# Define discrete approximation: +sigma = 1.0 +#mu = -0.5*(sigma**2) +#z = stats.lognorm(sigma, 0, np.exp(mu)) # Create "frozen" distribution instance +theta_grid_N = 7 + +sigma2 = 1.0 +N2 = 7 +p0 = 0.01 + +# Create a discrete approximation instance: +#theta = DiscreteApproximation(N=theta_grid_N, cdf=z.cdf, pdf=z.pdf, invcdf=z.ppf) +income = DiscreteApproximationTwoMeanOneIndependentLognormalDistribs( + N1=theta_grid_N, sigma1=sigma, N2=N2, sigma2=sigma2) +#DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event( +# N1=theta_grid_N, sigma1=sigma, N2=N2, sigma2=sigma2, pZevent=p0, z=0.0) + +# M grid: +m_min, m_max, m_size = 0.0, 4.0, 5 # Assign multiple values on a line +m_grid = np.linspace(m_min, m_max, m_size) + +# Set up a-grid: +a_min, a_max, a_size = 0.0, 4.0, 5 +a_grid = np.linspace(a_min, a_max, a_size) + +self_a_min = -min(income.X2) * R/gamma[0] # Self-imposed minimum a +self_c_min = min(m_grid) - self_a_min # Self-imposed min c + +# Define utility: +u = Utility(gamma=rho) + +# Create a Gothic object with these specific parameters: +gothic = Gothic(u, beta, rho, gamma, R, income) \ No newline at end of file diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 67dd82ba0..b2cf1fa5a 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -4,9 +4,42 @@ # Bring in modules we need import unittest -import HARK.algos.foc +from HARK.algos.foc import optimal_policy_foc +from gothic_class import gothic +from rewards import CRRAutilityP_inv import numpy as np + +""" + + +g = lambda x, k, a : {'a' : x['m'] - a['c']}, +dg_dx = 1, ## Used in FOC method, step 5 +dg_da = -1, ## Used in FOC method, step 5 +g_inv = lambda y, a : {'m' : y['a'] + a['c']}, ## Used in EGM method, step 8 +r = lambda x, k, a : u(a['c']), +dr_da = lambda x, k, a: u.prime(a['c']), +dr_inv = lambda uP : (CRRAutilityP_inv(uP, rho),), + +x = ['m'], +a = ['c'], +y = ['a'], + +TODO: Where are the constants from? +action_upper_bound = lambda x, k: (x['m'] + gamma[0] * theta.X[0] / R,), + +#discount = beta, <-- Removed because beta is in gothic V! + +##### Inputs to optimizers, interpolators, solvers... +optimizer_args = { + 'method' : 'Nelder-Mead', + 'options' : { + 'maxiter': 1e3, + #'disp' : True + } +}, +""" + class foc_test(unittest.TestCase): """ FOC test: @@ -22,7 +55,38 @@ def setUp(self): def test_x(self): - self.assertTrue(np.all(self.cVec2 > 0)) + g = lambda x, k, a : {'a' : x['m'] - a['c']}, + dg_dx = 1, ## Used in FOC method, step 5 + dg_da = -1, ## Used in FOC method, step 5 + g_inv = lambda y, a : {'m' : y['a'] + a['c']}, ## Used in EGM method, step 8 + r = lambda x, k, a : u(a['c']), + dr_da = lambda x, k, a: u.prime(a['c']), + dr_inv = lambda uP : (CRRAutilityP_inv(uP, rho),), + + action_upper_bound = lambda x, k: (x['m'] + gamma[0] * theta.X[0] / R,), + + #discount = beta, <-- Removed because beta is in gothic V! + + ##### Inputs to optimizers, interpolators, solvers... + ## TODO: Is this used? + optimizer_args = { + 'method' : 'Nelder-Mead', + 'options' : { + 'maxiter': 1e3, + #'disp' : True + } + } + + def consumption_v_y_der(y : Mapping[str,Any]): + return gothic.VP_Tminus1(y['a']) + + + pi_star, q_der, y_data = optimal_policy_foc( + {'m' : self.mVec}, + v_y_der = consumption_v_y_der + ) + + self.assertTrue(np.all(self.cVec2 == pi_star.values)) class egm_test(unittest.TestCase): """ @@ -36,4 +100,18 @@ def setUp(self): def test_egm(self): + """ + + ### EGM test from SolvingMicroDSOPs with stages + def consumption_v_y_der(y : Mapping[str,Any]): + return gothic.VP_Tminus1(y['a']) + + pi, pi_y = stage.optimal_policy_egm( + y_grid = {'a' : aVec}, + v_y_der = consumption_v_y_der + ) + + pi.interp({'m' : mVec_egm}) - cFunc_egm(mVec_egm) == 0 + """ + self.assertTrue(np.all(self.aVec + self.cVec_egm == self.mVec_egm)) \ No newline at end of file From 5b1ca7fe341bf984e096f86d02ed68f5b6201d45 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 13 Jun 2023 17:09:47 -0400 Subject: [PATCH 04/14] setting up use of gothic_class in algos tests --- HARK/algos/__init__.py | 0 HARK/algos/foc.py | 150 +-- HARK/algos/tests/__init__.py | 0 HARK/algos/tests/test_algos.py | 18 +- HARK/gothic/__init__.py | 0 HARK/{algos/tests => gothic}/gothic_class.py | 2 +- HARK/gothic/resources.py | 1021 ++++++++++++++++++ 7 files changed, 1108 insertions(+), 83 deletions(-) create mode 100644 HARK/algos/__init__.py create mode 100644 HARK/algos/tests/__init__.py create mode 100644 HARK/gothic/__init__.py rename HARK/{algos/tests => gothic}/gothic_class.py (98%) create mode 100644 HARK/gothic/resources.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/foc.py b/HARK/algos/foc.py index f8bdbfccd..a1d219288 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -1,3 +1,4 @@ +import itertools import numpy as np from scipy.optimize import minimize, brentq from typing import Callable, Mapping, Sequence @@ -47,10 +48,10 @@ def xndindex(ds, dims=None): def xz_grids_to_data_array( x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. - k_grid : Mapping[str, Sequence] = {} + z_grid : Mapping[str, Sequence] = {} ): - coords = {**x_grid, **k_grid} + coords = {**x_grid, **z_grid} da = xr.DataArray( np.zeros([len(v) for v in coords.values()]), @@ -61,7 +62,7 @@ def xz_grids_to_data_array( return da -def optimal_policy_foc(self, +def optimal_policy_foc( x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. z_grid : Mapping[str, Sequence] = {}, v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, @@ -100,7 +101,7 @@ def optimal_policy_foc(self, """ # Set up data arrays with coordinates based on the grid. - pi_data = xz_grids_to_data_array(x_grid, z_grid) + pi_data = xz_grids_to_data_array(x_grid, z_grid) ## TODO: Default value for pi should be nan. q_der_data = xz_grids_to_data_array(x_grid, z_grid) ## May need to expand this for multivariate y y_data = xz_grids_to_data_array(x_grid, z_grid) @@ -125,9 +126,10 @@ def foc(a): x_vals = {k : v for k, v in zip(x_grid.keys() , x_point)} for z_point in itertools.product(*z_grid.values()): - z_vals = {k : v for k, v in zip(k_grid.keys() , k_point)} + z_vals = {k : v for k, v in zip(z_grid.keys() , z_point)} # repeated code with the other optimizer -- can be functionalized out? + """ if len(self.actions) == 0: q_der_xz = self.d_q_d_a( x_vals, @@ -146,87 +148,85 @@ def foc(a): y = g(x_vals, z_vals, self.action_zip([np.nan])) y_n = np.array([y[k] for k in y]) y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) - - + """ + + # these lower bounds as arugments to the + lower_bound = self.action_lower_bound(x_vals, k_vals) + upper_bound = self.action_upper_bound(x_vals, k_vals) + + ## what if no lower bound? + q_der_lower = None + if lower_bound[0] is not None: + q_der_lower = self.d_q_d_a( + x_vals, + k_vals, + self.action_zip(lower_bound), + v_y_der + ) else: + lower_bound[0] = 1e-12 ## a really high number! + + q_der_upper = None + if upper_bound[0] is not None: + q_der_upper = self.d_q_d_a( + x_vals, + k_vals, + self.action_zip(upper_bound), + v_y_der + ) + else: + upper_bound[0] = 1e12 ## a really high number! - # these lower bounds as arugments to the - lower_bound = self.action_lower_bound(x_vals, k_vals) - upper_bound = self.action_upper_bound(x_vals, k_vals) + ## TODO: Better handling of case when there is a missing bound? + if q_der_lower is not None and q_der_upper is not None and q_der_lower < 0 and q_der_upper < 0: + a0 = lower_bound - ## what if no lower bound? - q_der_lower = None - if lower_bound[0] is not None: - q_der_lower = self.d_q_d_a( - x_vals, - k_vals, - self.action_zip(lower_bound), - v_y_der - ) - else: - lower_bound[0] = 1e-12 ## a really high number! - - q_der_upper = None - if upper_bound[0] is not None: - q_der_upper = self.d_q_d_a( - x_vals, - k_vals, - self.action_zip(upper_bound), - v_y_der - ) - else: - upper_bound[0] = 1e12 ## a really high number! - - ## TODO: Better handling of case when there is a missing bound? - if q_der_lower is not None and q_der_upper is not None and q_der_lower < 0 and q_der_upper < 0: - a0 = lower_bound + pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) + q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_lower) + elif q_der_lower is not None and q_der_upper is not None and q_der_lower > 0 and q_der_upper > 0: + a0 = upper_bound + pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) + q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_upper) + else: + ## Better exception handling here + ## asserting that Q is concave + if q_der_lower is not None and q_der_upper is not None and not(q_der_lower > 0 and q_der_upper < 0): + raise Exception("Cannot solve for optimal policy with FOC if Q is not concave!") + + a0, root_res = brentq( + foc, + lower_bound[0], # only works with scalar actions + upper_bound[0], # only works with scalar actions + full_output = True + ) + + if root_res.converged: pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) - q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_lower) - elif q_der_lower is not None and q_der_upper is not None and q_der_lower > 0 and q_der_upper > 0: - a0 = upper_bound - pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) - q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_upper) - else: - ## Better exception handling here - ## asserting that Q is concave - if q_der_lower is not None and q_der_upper is not None and not(q_der_lower > 0 and q_der_upper < 0): - raise Exception("Cannot solve for optimal policy with FOC if Q is not concave!") - - a0, root_res = brentq( - foc, - lower_bound[0], # only works with scalar actions - upper_bound[0], # only works with scalar actions - full_output = True + q_der_xz = self.dq_da( + x_vals, + z_vals, + self.action_zip((a0,)), # actions are scalar + v_y_der ) - if root_res.converged: - pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) - - q_der_xz = self.dq_da( - x_vals, - z_vals, - self.action_zip((a0,)), # actions are scalar - v_y_der - ) - - q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_xk) - else: - print(f"Rootfinding failure at {x_vals}, {z_vals}.") - print(root_res) + q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_xk) + else: + print(f"Rootfinding failure at {x_vals}, {z_vals}.") + print(root_res) - #raise Exception("Failed to optimize.") - pi_data.sel(**x_vals, **z_vals).variable.data.put(0, root_res.root) + #raise Exception("Failed to optimize.") + pi_data.sel(**x_vals, **z_vals).variable.data.put(0, root_res.root) - q_der_xz = dq_da( - x_vals, - z_vals, - self.action_zip((root_res.root,)), # actions are scalar - v_y_der - ) + q_der_xz = dq_da( + x_vals, + z_vals, + self.action_zip((root_res.root,)), # actions are scalar + v_y_der + ) - q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) acts = np.atleast_1d(pi_data.sel(**x_vals, **z_vals).values) y = g(x_vals, z_vals, self.action_zip(acts)) diff --git a/HARK/algos/tests/__init__.py b/HARK/algos/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index b2cf1fa5a..877745dba 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -5,10 +5,14 @@ import unittest from HARK.algos.foc import optimal_policy_foc -from gothic_class import gothic -from rewards import CRRAutilityP_inv +from HARK.gothic.gothic_class import gothic +from HARK.rewards import CRRAutilityP_inv import numpy as np +import os +from typing import Any, Mapping +__location__ = os.path.realpath( + os.path.join(os.getcwd(), os.path.dirname(__file__))) """ @@ -49,8 +53,8 @@ class foc_test(unittest.TestCase): def setUp(self): - self.mVec = np.load("smdsops_mVec.npy") - self.cVec2 = np.load("smdsops_cVec2.npy") + self.mVec = np.load(os.path.join(__location__, "smdsops_mVec.npy")) + self.cVec2 = np.load(os.path.join(__location__, "smdsops_cVec2.npy")) def test_x(self): @@ -94,9 +98,9 @@ class egm_test(unittest.TestCase): pi(mVec_egm) == cVec_egm """ def setUp(self): - self.aVec = np.load("smdsops_aVec.npy") - self.cVec_egm = np.load("smdsops_cVec_egm.npy") - self.mVec_egm = np.load("smdsops_mVec_egm.npy") + self.aVec = np.load(os.path.join(__location__, "smdsops_aVec.npy")) + self.cVec_egm = np.load(os.path.join(__location__, "smdsops_cVec_egm.npy")) + self.mVec_egm = np.load(os.path.join(__location__, "smdsops_mVec_egm.npy")) def test_egm(self): diff --git a/HARK/gothic/__init__.py b/HARK/gothic/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/HARK/algos/tests/gothic_class.py b/HARK/gothic/gothic_class.py similarity index 98% rename from HARK/algos/tests/gothic_class.py rename to HARK/gothic/gothic_class.py index 0340d7a65..d9d889e0c 100644 --- a/HARK/algos/tests/gothic_class.py +++ b/HARK/gothic/gothic_class.py @@ -245,7 +245,7 @@ def VP_t(self, a, c_prime, t=None): import pylab as plt from scipy.optimize import brentq ## TODO Need these resources here also? -from resources import DiscreteApproximation, Utility, DiscreteApproximationTwoMeanOneIndependentLognormalDistribs, DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event +from HARK.gothic.resources import DiscreteApproximation, Utility, DiscreteApproximationTwoMeanOneIndependentLognormalDistribs, DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event # General parameters: rho = 2.0 diff --git a/HARK/gothic/resources.py b/HARK/gothic/resources.py new file mode 100644 index 000000000..2983266b4 --- /dev/null +++ b/HARK/gothic/resources.py @@ -0,0 +1,1021 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:light +# text_representation: +# extension: .py +# format_name: light +# format_version: '1.5' +# jupytext_version: 1.14.4 +# kernelspec: +# display_name: Python 3 (ipykernel) +# language: python +# name: python3 +# --- + +# # Resources + +# This is a library of classes and functions which I often employ in simulations and numerical work. The classes (and occasional function) are organized alphabetically by titled section (within each section, the classes themselves are organized by importance/inheritance, as needed). +# +# The classes included in Resources are used by many different projects and are gathered one place for convenience. Specific significant classes (such as the Gothic class) get their own library files. +# +# ***NOTE:*** If at all possible, please add to the "resources.py" file by updating *this* Resources.ipynb file and then manually exporting to the "resources.py" file in the same directory (and then, of course, opening a Python command prompt and executing "import resources" to overwrite and update the "resource.pyc" file, which is ultimately used by import commands in other Python programs). +# +# Note, finally, that this is a cleaned-up and pared-down version of my original "resources.py" file, which may still be found in "resources_backup.py." + +# + +# Import required Python libraries +from __future__ import division # In Python 3.x this will not be needed. +import math +import numpy as np +from warnings import warn +import scipy.stats as stats +from scipy.integrate import quad +from numpy.random import RandomState +from scipy.interpolate import interp1d +import pylab as plt + + +from numpy import log, exp # For ease of reading in get_improved_grid + + +# - + +# ## Discrete Approximation to a Continuous Distribution + +# This implements a Python version of the discrete approximation code used in Carroll's Microeconomic DSOP notes. + +# + code_folding=[2, 49, 65] +class DiscreteApproximation(object): + + def __init__(self, N, cdf, pdf, invcdf, precise_summation=False): + """ + Given a cdf, pdf, and invererse cdf, and number of evenly-spaced bins + (a partition) of a state space, find the conditional expected values + of each evenly-sized partition, E[X|partition] + + An equidistant grid is imposed upon the [0,1] range, and the inverse cdf + is employed to determine the appropriate "bin cutoffs." + + The cdf, pdf, and inverse cdf need to be "frozen" distributions -- they + need to have their parameters already set and only accept one value to + produce output. + + Note that the MATLAB version has a manually-entered minpoint and maxpoint. + These numbers, used "raw" as an upper bound in the Fortran QUADPACK + automatic integration function (which is used by both Numpy and Octave, + as well as many others) cause numerical instability in the solution. + + The solution is to use np.inf as the upper bound, which happens + automatically when the state cutoffs are determined. The np.inf values + forces correct evaluation of the integral. + For more details please see discussion in this stackoverflow: + http://stackoverflow.com/questions/24003694/discontinuity-in-results-when-using-scipy-integrate-quad + """ + self.cdf = cdf + self.pdf = pdf + self.invcdf = invcdf + self.N = N # Number of bins + probs_cutoffs = np.arange(N+1)/N # Includes 0 and 1 + state_cutoffs = invcdf(probs_cutoffs) # State cutoff values, each bin + #bin_probs = np.diff(probs_cutoffs) # Total prob mass, each bin + bin_probs = np.zeros(N) + + # Find the E[X|bin] values: + F = lambda x: x*self.pdf(x) + Ebins = [] + for i, (x0, x1) in enumerate(zip(state_cutoffs[:-1], state_cutoffs[1:])): + bin_probs[i] = cdf(x1) - cdf(x0) ## pdf between x0 and x1 + cond_mean, err = quad(F, x0, x1) + Ebins.append(cond_mean/bin_probs[i]) + + self.X = np.array(Ebins) + self.pmf = bin_probs + + self.precise_summation = precise_summation + + + def E(self, f=None): + """ + Find expectation of f(X) over the discrete space. + """ + if f is None: + if self.precise_summation: + return math.fsum(self.pmf * self.X) # High precision + else: + return self.pmf.dot(self.X) # ~100x faster + else: + if self.precise_summation: + return math.fsum(np.multiply(self.pmf, f(self.X))) # High precision + else: + return self.pmf.dot(f(self.X)) # ~100x faster + + + def plot(self, x0, x1, gridsize=100): + """ + Plot the discrete approximation against the actual distribution. + """ + grid = np.linspace(x0, x1, gridsize) + + plt.plot(self.X, self.cdf(self.X), 'ok') + plt.plot(grid, self.cdf(grid), 'k-') + plt.axis([x0, x1, 0, 1]) + plt.hlines(self.pmf.cumsum(),x0, x1,color ="black",linestyles='--') ## added by Tao to show equiprobable ranges + plt.title('Discrete Approximation to Lognormal Distribution') + plt.xlabel('theta') + plt.ylabel('CDF') + plt.grid() + plt.show() + + +# - + +# ## A Specific Discrete Approximation to the Mean-one lognormal: +# +# Now we use simple inheritance to create a mean-one lognormal instance quickly and easily. + +# + code_folding=[] +class DiscreteApproximationToMeanOneLogNormal(DiscreteApproximation): + """ + Extension of the DiscreteApproximation class, which creates a mean-one + lognormal approximation, given standard deviation 'sigma' and number + of grid-points 'N'. + """ + + def __init__(self, N, sigma, precise_summation=False): + """ + Given N and sigma, create a lognormal "frozen" distribution object, + and initilize an appropriate DiscreteApproximation. + + N: integer, number of grid points to interpolate over + sigma: double, standard deviation of the lognormal distribution + """ + + # Save the one bit of new information (N saved in super()): + self.sigma = sigma + self.mu = -0.5*(self.sigma**2) + + # Set up the "frozen" lognormal: + distrib = stats.lognorm(self.sigma, 0, np.exp(self.mu)) + + # Create the class: + super(DiscreteApproximationToMeanOneLogNormal, self).__init__( + N, distrib.cdf, distrib.pdf, distrib.ppf, + precise_summation=precise_summation) + + # Everything is done now. + + + +# + code_folding=[] +# Testing: +if __name__ == "__main__": + import numpy as np + import pylab as plt + from warnings import warn + from scipy.integrate import quad + from numpy.random import RandomState + from scipy.interpolate import interp1d + import scipy.stats as stats + from copy import deepcopy + + # Create a 2-D discrete approximation instance: + sigma = 0.1 + N = 7 + + LNDiscApprox = DiscreteApproximationToMeanOneLogNormal(N, sigma) + + print("LNDiscApprox.E():"+ str(LNDiscApprox.E())) + # Check that "raw expectation" is correct: + assert np.max(np.abs(LNDiscApprox.E() - 1.0)) < 1e-16, "Unconditional expectation is not correct!" + print("Success: Unconditional expectation is correct!") + + # Check that the expectation applied to a misc function is correct: + # Manually, the 7-point approximation for each gridpoint should be: + manual_X = np.array([ 0.85043016002691718125, 0.91862318529875575113, + 0.95908470592906813756, 0.99506598629571241243, + 1.03241349447674446438, 1.0779763032188010019 , + 1.16640616475400205054]) + manual_p = np.repeat(1.0/7, 7) + + # Check for a random function: + g = lambda x1: x1**2 + 1.5 + + # Manually construct the expectation of g: + manual_EgX = 0.0 + for x, p in zip(manual_X, manual_p): + manual_EgX += g(x) * p + + # Now compare manual value against vectorized: + EgX = LNDiscApprox.E(g) + + #print "manual_EgXY:", manual_EgXY + #print "EgXY:", EgXY + + assert np.abs(manual_EgX - EgX) < 1e-12,"Eg(X) not equal between the values that have expectations." + print("Success: Eg(X) = manually calculated values.") + print("All tests passed successfully.") + + LNDiscApprox.plot(0.7, 1.4) + + # Comparing to the 0630 Matlab definition: + matsig = np.array([0.509520331925153, 0.667826497278589, 0.776380649071803, 0.879396570886877, 0.989752424342583, 1.121403448305962, 1.305157824866014, 1.750562243647017]) + matN = 8 + + +# - + +# ## Discrete Approximation to Two Independent Continuous Distributions +# +# A discrete approximation which neatly bundles two independent continuous distributions into a single discrete approximation object. This "buys" the user a simple expectations operator. +# +# Further below we will include a simple extension to automatically set up two mean-one lognormal approximations, and two log-normal approximations with discrete probability of a 0-valued event. +# + +# + code_folding=[164] +class DiscreteApproximationTwoIndependentDistribs(object): + + def __init__(self, N1, cdf1, pdf1, invcdf1, N2, cdf2, pdf2, invcdf2, precise_summation=False): + """ + Given a cdf, pdf, and invererse cdf, and number of evenly-spaced bins + (a partition) of a state space, find the conditional expected values + of each evenly-sized partition, E[X|partition] + + An equidistant grid is imposed upon the [0,1] range, and the inverse cdf + is employed to determine the appropriate "bin cutoffs." + + The cdf, pdf, and inverse cdf need to be "frozen" distributions -- they + need to have their parameters already set and only accept one value to + produce output. + + Note that the "first" distribution, X1, will be represented as running + along the first axis of the joined space, while the second distribution, + X2, will run along the horizontal axis. + + We will take: + X1 = xrow = [a, b] + X2 = ycol = [d, e, f] + + we want: + x_mesh = [[a, a, a], + [b, b, b]] + + y_mesh = [[d, e, f], + [d, e, f]] + + to represent the joint distribution. + + + Note 2: the MATLAB version has a manually-entered minpoint and maxpoint. + These numbers, used "raw" as an upper bound in the Fortran QUADPACK + automatic integration function (which is used by both Numpy and Octave, + as well as many others) cause numerical instability in the solution. + + The solution is to use np.inf as the upper bound, which happens + automatically when the state cutoffs are determined. The np.inf values + forces correct evaluation of the integral. + + + """ + + self.precise_summation = precise_summation + # Used to implement very high-precision calculation of the expected value. Good + # for summations of very big, very small, or very close-to-zero numbers. + + # ------ Set up first discrete approx ------ + self.cdf1 = cdf1 + self.pdf1 = pdf1 + self.invcdf1 = invcdf1 + self.N1 = N1 # Number of bins + + probs_cutoffs1 = np.arange(N1+1.0)/N1 # Includes 0 and 1 + state_cutoffs1 = invcdf1(probs_cutoffs1) # State cutoff values, each bin + bin_probs1 = np.zeros(N1) + + # Find the E[X|bin] values: + F1 = lambda x: x*self.pdf1(x) + Ebins1 = [] + for i, (x0, x1) in enumerate(zip(state_cutoffs1[:-1], state_cutoffs1[1:])): + bin_probs1[i] = cdf1(x1) - cdf1(x0) + cond_mean1, err1 = quad(F1, x0, x1) + Ebins1.append(cond_mean1/bin_probs1[i]) + + self.X1 = np.array(Ebins1) + self.pmf1 = bin_probs1 + + + # ------ Set up second discrete approx ------ + self.cdf2 = cdf2 + self.pdf2 = pdf2 + self.invcdf2 = invcdf2 + self.N2 = N2 # Number of bins + + probs_cutoffs2 = np.arange(N2+1.0)/N2 # Includes 0 and 1 + state_cutoffs2 = invcdf2(probs_cutoffs2) # State cutoff values, each bin + bin_probs2 = np.zeros(N2) + + # Find the E[X|bin] values: + F2 = lambda x: x*self.pdf2(x) + Ebins2 = [] + for i, (x0, x1) in enumerate(zip(state_cutoffs2[:-1], state_cutoffs2[1:])): + bin_probs2[i] = cdf2(x1) - cdf2(x0) + cond_mean2, err2 = quad(F2, x0, x1) + Ebins2.append(cond_mean2/bin_probs2[i]) + + self.X2 = np.array(Ebins2) + self.pmf2 = bin_probs2 + + ''' + For + xrow = [a, b] + ycol = [d, e, f] + + we want: + x_mesh = [[a, a, a], + [b, b, b]] + + y_mesh = [[d, e, f], + [d, e, f]] + + ''' + + # Now construct the X1/X2 mesh values: + nrow = len(self.X1) + ncol = len(self.X2) + + # We don't use the np.meshgrid commands, because they do not + # easily support non-symmetric grids. + self.X1mesh,self.X2mesh = np.meshgrid(self.X1,self.X2) + self.X1mesh = self.X1mesh.T + self.X2mesh = self.X2mesh.T + #self.X1mesh = np.tile(np.transpose([self.X1]), (1, ncol)) + #self.X2mesh = np.tile(self.X2, (nrow, 1)) + + # Construct the appropriate probability for each point: + self.pmf = np.zeros_like(self.X1mesh) + # Loop to fill in probs: + for i, p1 in enumerate(self.pmf1): + for j, p2 in enumerate(self.pmf2): + self.pmf[i,j] = p1*p2 + + # Create flat versions of all these: + self.flatX1 = self.X1mesh.ravel() + self.flatX2 = self.X1mesh.ravel() + self.flatpmf = self.pmf.ravel() + + # Check that sums to 1: + #print "np.sum(self.pmf):", np.sum(self.pmf) + assert np.abs(np.sum(self.pmf) - 1.0) < 1e-10, "Total 2D pmf doesn't sum to 1." + + + def condE1(self, x2, f=None): + pass + + + def E(self, f=None): + """ + Find the expectation of f over X1, X2. + + Note f must work for being applied to a grid of values. + """ + if f is None: + # Return simple conditional expectation. + + if self.precise_summation: + return( np.array( [math.fsum(self.pmf1 * self.X1), math.fsum(self.pmf2 * self.X2)] ) ) + else: + return( np.array( [self.pmf1.dot(self.X1), self.pmf2.dot(self.X2)] ) ) + + else: + fval = f(self.X1mesh, self.X2mesh) + #fval = f(self.flatX1, self.flatX2) + a = np.multiply(self.pmf, fval) + if self.precise_summation: + return(math.fsum(a)) # High-floating-point precision summation. + else: + return( np.sum(a) ) # 100x faster C summation + # np.multiply works the same for both numpy arrays and matrices; + # potentially important if for some reason, f takes an array and + # returns a matrix. Likely uncommon, but just want to be safe. + + + def plot(self, x0, x1, x20, x21, gridsize=100): + """ + Plot the discrete approximation against the actual distribution. + """ + grid = np.linspace(x0, x1, gridsize) + + plt.plot(self.X1, self.cdf1(self.X1), 'ok') + plt.plot(grid, self.cdf1(grid), 'k-') + plt.hlines(self.pmf1.cumsum(),x0, x1,color ="black",linestyles='--') ## added by Tao to show equiprobable ranges + plt.axis([x0, x1, 0, 1]) + plt.title('Discrete Approximation to Lognormal Distribution') + plt.xlabel('theta') + plt.ylabel('CDF') + plt.grid() + plt.show() + + grid = np.linspace(x20, x21, gridsize) + + plt.plot(self.X2, self.cdf2(self.X2), 'ok') + plt.plot(grid, self.cdf2(grid), 'k-') + plt.hlines(self.pmf2.cumsum(),x0, x1,color ="black",linestyles='--') ## added by Tao to show equiprobable ranges + plt.axis([x20, x21, 0, 1]) + plt.title('Second Discrete Approximation to Lognormal Distribution') + plt.xlabel('theta') + plt.ylabel('CDF') + plt.grid() + plt.show() + + +# + code_folding=[] +# Immediately run some tests. +# Note that this will not run unless this is executed as the "main" file. +if __name__ == "__main__": + import numpy as np + import pylab as plt + from warnings import warn + from scipy.integrate import quad + from numpy.random import RandomState + from scipy.interpolate import interp1d + import scipy.stats as stats + from copy import deepcopy + + # Create a 2-D discrete approximation instance: + sigma1 = 0.1 + mu1 = -0.5*(sigma1**2) + N1 = 7 + z1 = stats.lognorm(sigma1, 0, np.exp(mu1)) # Create "frozen" distribution instance + + sigma2 = 0.1 + mu2 = -0.5*(sigma2**2) + N2 = 8 + z2 = stats.lognorm(sigma2, 0, np.exp(mu2)) # Create "frozen" distribution instance + + TwoDimDiscApprox = DiscreteApproximationTwoIndependentDistribs(N1, z1.cdf, z1.pdf, z1.ppf, N2, z2.cdf, z2.pdf, z2.ppf) + + + # Check that "raw expectation" is correct: + assert np.max(np.abs(TwoDimDiscApprox.E() - 1.0)) < 1e-16, "Unconditional expectation is not correct!" + + # Check that the expectation applied to a misc function is correct: + # Manually, the 7-point approximation for each gridpoint should be: + manual_X1 = np.array([ 0.85043016002691718125, 0.91862318529875575113, + 0.95908470592906813756, 0.99506598629571241243, + 1.03241349447674446438, 1.0779763032188010019 , + 1.16640616475400205054]) + manual_p1 = np.repeat(1.0/7, 7) + + manual_X2 = deepcopy(manual_X1) + manual_p2 = np.repeat(1.0/7, 7) + + # Check for a random function: + g = lambda x1, x2: x1**2 + 1.5*x2 + + # Manually construct the expectation of g: + manual_EgXY = 0.0 + for x1, p1 in zip(manual_X1, manual_p1): + for x2, p2 in zip(manual_X2, manual_p2): + manual_EgXY += g(x1, x2) * p1 * p2 + + # Now compare manual value against vectorized: + EgXY = TwoDimDiscApprox.E(g) + + #print "manual_EgXY:", manual_EgXY + #print "EgXY:", EgXY + + assert np.abs(manual_EgXY - EgXY) < 1e-12, "Eg(X,Y) not equal between the values that have expectations." + + TwoDimDiscApprox.plot(0.7, 1.4, 0.7, 1.4) + + +# + code_folding=[0, 8] +class DiscreteApproximationTwoMeanOneIndependentLognormalDistribs(DiscreteApproximationTwoIndependentDistribs): + """ + Extend the "DiscreteApproximationTwoIndependentDistribs" to automatically define a + lognormal, mean 1 distribution. + + Nathan Palmer + """ + + def __init__(self, N1, sigma1, N2, sigma2, precise_summation=False): + """ + Given N1, sigma1, N2, sigma2, create two lognormal "frozen" distributions, + and initilize an appropriate DiscreteApproximationTwoIndependentDistribs. + + Very simple inheritance exercise. + + N1: integer, number of grid points to interpolate over for distrib 1 + sigma1: double, standard deviation of the first lognormal distribution + N2: integer, number of grid points to interpolate over for distrib 2 + sigma2: double, standard deviation of the second lognormal distribution + """ + + # Save the one bit of new information (N saved in super()): + self.sigma1 = sigma1 + self.sigma2 = sigma2 + + # Set up the "frozen" lognormals: + self.mu1 = -0.5*(self.sigma1**2) + distrib1 = stats.lognorm(self.sigma1, 0, np.exp(self.mu1)) + + self.mu2 = -0.5*(self.sigma2**2) + distrib2 = stats.lognorm(self.sigma2, 0, np.exp(self.mu2)) + + # Create the class: + super(DiscreteApproximationTwoMeanOneIndependentLognormalDistribs, self).__init__( + N1=N1, cdf1=distrib1.cdf, pdf1=distrib1.pdf, invcdf1=distrib1.ppf, + N2=N2, cdf2=distrib2.cdf, pdf2=distrib2.pdf, invcdf2=distrib2.ppf, + precise_summation=precise_summation) + + +#DiscreteApproximationTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_0_Event + +# - + +# Immediately run some tests. +# Note that this will not run unless this is executed as the "main" file. +if __name__ == "__main__": + import numpy as np + import pylab as plt + from warnings import warn + from scipy.integrate import quad + from numpy.random import RandomState + from scipy.interpolate import interp1d + import scipy.stats as stats + from copy import deepcopy + + # Create a 2-D discrete approximation instance: + sigma1 = 0.1 + mu1 = -0.5*(sigma1**2) + N1 = 7 + z1 = stats.lognorm(sigma1, 0, np.exp(mu1)) # Create "frozen" distribution instance + + sigma2 = 0.1 + mu2 = -0.5*(sigma2**2) + N2 = 7 + z2 = stats.lognorm(sigma2, 0, np.exp(mu2)) # Create "frozen" distribution instance + + TwoDimDiscApprox = DiscreteApproximationTwoMeanOneIndependentLognormalDistribs( + N1=N1, sigma1=sigma1, N2=N2, sigma2=sigma2) + + + # Check that mu calculated correctly: + assert np.max(np.abs(TwoDimDiscApprox.mu1 - mu1)) < 1e-16, "Mu1 not calculated correctly!" + assert np.max(np.abs(TwoDimDiscApprox.mu2 - mu2)) < 1e-16, "Mu2 not calculated correctly!" + print("M1 and Mu2 were both calculated correctly!") + + # Check that "raw expectation" is correct: + assert np.max(np.abs(TwoDimDiscApprox.E() - 1.0)) < 1e-16, "Unconditional expectation is not correct!" + print("Unconditional expectation *is* correct: E(X1, X2) ="+str(TwoDimDiscApprox.E())) + + # Check that the expectation applied to a misc function is correct: + # Manually, the 7-point approximation for each gridpoint should be: + manual_X1 = np.array([ 0.85043016002691718125, 0.91862318529875575113, + 0.95908470592906813756, 0.99506598629571241243, + 1.03241349447674446438, 1.0779763032188010019 , + 1.16640616475400205054]) + manual_p1 = np.repeat(1.0/7, 7) + + manual_X2 = deepcopy(manual_X1) + manual_p2 = np.repeat(1.0/7, 7) + + # Check for a random function: + g = lambda x1, x2: x1**2 + 1.5*x2 + + # Manually construct the expectation of g: + manual_EgXY = 0.0 + for x1, p1 in zip(manual_X1, manual_p1): + for x2, p2 in zip(manual_X2, manual_p2): + manual_EgXY += g(x1, x2) * p1 * p2 + + # Now compare manual value against vectorized: + EgXY = TwoDimDiscApprox.E(g) + + #print "manual_EgXY:", manual_EgXY + #print "EgXY:", EgXY + + assert np.abs(manual_EgXY - EgXY) < 1e-12, "Eg(X,Y) not equal manually calculated joint value." + print("Eg(X,Y) *does* equal manually calculated joint value:") + print("\tE[g(X1, X2)] from class:"+ str(EgXY)+ "\n\tand from manual calc:"+ str(manual_EgXY) + + "\n\tnp.abs(manual_EgXY - EgXY):"+str(np.abs(manual_EgXY - EgXY))) + + TwoDimDiscApprox.plot(0.7, 1.4, 0.7, 1.4) + + +# + code_folding=[] +class DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event(DiscreteApproximationTwoMeanOneIndependentLognormalDistribs): + """ + ---------------------------------------------------------------------------- + Extend the "DiscreteApproximationTwoMeanOneIndependentLognormalDistribs" + class by introducing a discrete probability of a discrete-valued event occuring + for the second random variable, with value "z." + + The second variable experiences a z-event with probability pZevent. If the + state vector for the second random variable, X2, already includes a z value, + no changes are made (and a warning is raised). + + Otherwise the state vector X2 is prepended with z, and all *other* + values of X2 are set to X2/(1-pZevent). + + The probability vector pmf2 is prepended with pZevent, and all other + pmf2 are set to pmf2*(1-pZevent). + + Finally, the state-space is re-determined and the total pmf matrix, + "pmf", is re-calcualted with the new values. + + All other methods should function the same as before. + This is still a relatively simple example of multiple inheritance. + """ + + def __init__(self, N1, sigma1, N2, sigma2, pZevent, z=0.0, precise_summation=False): + + # Execute the initilizer for parent own class: + super(DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event, self).__init__( + N1=N1, sigma1=sigma1, N2=N2, sigma2=sigma2, precise_summation=precise_summation) + + self.pZevent = pZevent + + if not 0 <= pZevent <= 1: + raise Exception("The probability that discrete event z = "+ z +" has a probability not in the range [0,1]: pZevent = "+str(pZevent)) + + # Update X2 vector: + if z in self.X2: + warn("Discrete shock value "+ z +" already exists in the RV discrete space self.X2. Please confirm that inputs are correct") + else: + self.X2 = np.append(z, self.X2/(1.0 - pZevent)) + + # Update pmf2: + self.pmf2 = np.append(pZevent, self.pmf2 * (1.0 - pZevent)) + + # Update total state-space: + nrow = len(self.X1) + ncol = len(self.X2) + + # We don't use the np.meshgrid commands, because they do not + # easily support non-symmetric grids. + self.X1mesh = np.tile(np.transpose([self.X1]), (1, ncol)) + self.X2mesh = np.tile(self.X2, (nrow, 1)) + + # Construct the appropriate probability for each point: + self.pmf = np.zeros_like(self.X1mesh) + + # Loop to fill in probs: + for i, p1 in enumerate(self.pmf1): + for j, p2 in enumerate(self.pmf2): + self.pmf[i,j] = p1*p2 + + # Flat versions: + self.flatX1 = self.X1mesh.ravel() + self.flatX2 = self.X2mesh.ravel() + self.flatpmf = self.pmf.ravel() + + # Check that sums to 1: + #print "np.sum(self.pmf):", np.sum(self.pmf) + assert np.abs(np.sum(self.pmf) - 1.0) < 1e-10, "Total 2D pmf doesn't sum to 1." + + + + +# + code_folding=[] +# Immediately run some tests. +# Note that this will not run unless this is executed as the "main" file. +if __name__ == "__main__": + import numpy as np + import pylab as plt + from warnings import warn + from scipy.integrate import quad + from numpy.random import RandomState + from scipy.interpolate import interp1d + import scipy.stats as stats + from copy import deepcopy + + # Create a 2-D discrete approximation instance: + sigma1 = 0.1 + mu1 = -0.5*(sigma1**2) + N1 = 7 + z1 = stats.lognorm(sigma1, 0, np.exp(mu1)) # Create "frozen" distribution instance + + sigma2 = 0.1 + mu2 = -0.5*(sigma2**2) + N2 = 7 + z2 = stats.lognorm(sigma2, 0, np.exp(mu2)) # Create "frozen" distribution instance + + # Define a "z" event: + prob_Z_event = 0.1 # If we push this to 0.5 or 0.999, we'll see the "error" between manaual and numpy calc grow ~1e-15 to 1e-13 or 1e-14 + z_value = 500.0 # Same here; as we increase order of magnitude, we'll see "error" between manual and numpy calc grow. + # See user warning below. + + # UPDATE: NOTE that this problem was entirely solved by using math.fsum() to conduct high-precision summation. + + TwoDimDiscApproxZ = DiscreteApproximationToTwoMeanOneIndependentLognormalDistribsWithDiscreteProb_Z_Event( + N1=N1, sigma1=sigma1, N2=N2, sigma2=sigma2, + pZevent=prob_Z_event, z=z_value, precise_summation=True) + # Try precise_summation=False to see errors emerge. + + # Check that mu calculated correctly: + assert np.max(np.abs(TwoDimDiscApproxZ.mu1 - mu1)) < 1e-16, "Mu1 not calculated correctly!" + assert np.max(np.abs(TwoDimDiscApproxZ.mu2 - mu2)) < 1e-16, "Mu2 not calculated correctly!" + print("M1 and Mu2 were both calculated correctly!") + + # Check that the expectation applied to a misc function is correct: + # Manually, the 7-point approximation for each gridpoint should be: + manual_X1 = np.array([ 0.85043016002691718125, 0.91862318529875575113, + 0.95908470592906813756, 0.99506598629571241243, + 1.03241349447674446438, 1.0779763032188010019 , + 1.16640616475400205054]) + manual_p1 = np.repeat(1.0/7, 7) + + manual_X2 = deepcopy(manual_X1) + manual_p2 = np.repeat(1.0/7, 7) + + # Manually adjust X2 for the 0-valued event: + manual_X2 = np.append(z_value, manual_X2/(1.0-prob_Z_event)) + manual_p2 = np.append(prob_Z_event, manual_p2*(1.0-prob_Z_event)) + + # Manually calculate the unconditional expectation: + #manual_EX1 = np.dot(manual_X1, manual_p1) + #manual_EX2 = np.dot(manual_X2, manual_p2) + manual_EX1 = math.fsum(manual_X1*manual_p1) + manual_EX2 = math.fsum(manual_X2*manual_p2) + manual_EX = np.array([manual_EX1, manual_EX2]) + + # Check that "raw expectation" is correct: + print("TwoDimDiscApprox.E()"+ str(TwoDimDiscApproxZ.E())) + print("Manual E[(X1, X2)]"+str(manual_EX)) + print("max(abs(diff)):"+str(np.max(np.abs(TwoDimDiscApproxZ.E() - manual_EX)))) + # This is the value we "know" it should be, for a shock of 0.0: + #assert np.max(np.abs(TwoDimDiscApprox.E() - 1.0)) < 1e-16, "Unconditional expectation is not correct for a shock value of 0!" + #print "Unconditional expectation *is* correct: E(X1, X2) =", TwoDimDiscApprox.E() + + # With the manually calculated value: + assert np.max(np.abs(TwoDimDiscApproxZ.E() - manual_EX)) < 1e-12, "Unconditional expectation is not correct!" + print("Unconditional expectation *is* correct: E(X1, X2) =\n\tFrom class:"+ str(TwoDimDiscApproxZ.E())+ "\n\tManual calc:"+ str(manual_EX)) + + + # Check for a random function: + g = lambda x1, x2: x1**2 + 1.5*x2 + + # Manually construct the expectation of g: + #manual_EgXY = 0.0 + temp_manual_EgXY = [] + for x1, p1 in zip(manual_X1, manual_p1): + for x2, p2 in zip(manual_X2, manual_p2): + #manual_EgXY += g(x1, x2) * p1 * p2 + temp_manual_EgXY.append(g(x1, x2) * p1 * p2) + + manual_EgXY = math.fsum(temp_manual_EgXY) + + # Now compare manual value against vectorized: + EgXY = TwoDimDiscApproxZ.E(g) + + #print "manual_EgXY:", manual_EgXY + #print "EgXY:", EgXY + + print("TwoDimDiscApprox.E(g)"+str(TwoDimDiscApproxZ.E(g))) + print("Manual E[g(X1, X2)]"+str(manual_EgXY)) + print("max(abs(diff)):"+str(np.max(np.abs(TwoDimDiscApproxZ.E(g) - manual_EgXY)))) + + + assert np.abs(manual_EgXY - EgXY) < 1e-16, "Eg(X,Y) not equal manually calculated joint value." + print("Eg(X,Y) *does* equal manually calculated joint value:") + print("\tE[g(X1, X2)] from class:"+ str(EgXY)+"s\n\tand from manual calc:"+str(manual_EgXY)) + + + warn("\n\nNOTE: There is a very small difference between the manual and Numpy-summed calculation of the e[g(X1,X2)]. 1e-15, but still there, and not 1e-16. Very mild concern -- but only if this grows with increasing numbers.\n\n") + # NOTE: This may very well be due to how Numpy can "correct" for precision errors in addition of very small or very large numbers. + warn("\n\nUPDATE NOTE: this problem was entirely solved by using math.fsum() to conduct high-precision summation. The trouble: math.fsum is 100x slower than np.sum(). Use precise_summation=True to see this version work.") + + TwoDimDiscApproxZ.plot(0.7, 1.4, 0.7, 1.4) + + +# - + +# ## Discrete Random Variable +# +# A very simple wrapper to provide a discerete random variable with only two values -- "employed" and "unemployed" income. + +# + code_folding=[9] +class SimpleDiscreteRandomVariable(stats.rv_discrete): + + def __init__(self, values, probs, certain_value, name="discrete random variable", precise_summation=False): + discrete_RV_value_prob = [values, probs] + self.precise_summation = precise_summation + super(SimpleDiscreteRandomVariable, self).__init__(min(0,min(values)), np.inf, name, None, 1e-10, discrete_RV_value_prob) + # stats.rv_discrete(self, a=0, b=inf, name=None, badvalue=None, moment_tol=1e-08, values=None, ...) + self.certain_value=certain_value + + def E(self, f=None): + """ + Find expectation of f(X) over the discrete space. + """ + if f is None: + if self.precise_summation: + return math.fsum(np.multiply(self.pk, self.xk)) # High precision + else: + return self.mean() + else: + fX = f(self.certain_value, self.xk) + if self.precise_summation: + return math.fsum(np.multiply(self.pk, fX)) # High precision + else: + return np.dot(self.pk, fX) # ~100x faster + + +# - + +# ### Testing +# +# Simple testing. + +if __name__ == "__main__": + + p=0.0005 + vals = np.array([0.0, 1.0/(1.0-p)]) + probs = np.array([p, 1.0-p]) + + test = SimpleDiscreteRandomVariable(values=vals, probs=probs, certain_value=1.0, name="discrete random variable", precise_summation=True) + + manual_EX = math.fsum(vals*probs) + + # Check that "raw expectation" is correct: + print("test.E()"+str(test.E())) + print("Manual E[X]"+ str(manual_EX)) + print("max(abs(diff)):"+str(np.max(np.abs(test.E() - manual_EX)))) + + # With the manually calculated value: + assert np.max(np.abs(test.E() - manual_EX)) < 1e-12, "Unconditional expectation is not correct!" + print("Unconditional expectation *is* correct: E(X) =\n\tFrom class:"+str(test.E())+"\n\tManual calc:"+str(manual_EX)) + + # Check for a random function: + g = lambda x1, x2: x1**2 + pi*x2 #g = lambda x1: x1**2 + 1.5 + + # Manually construct the expectation of g: + temp_manual_EgX = [] + for x1, p1 in zip(vals, probs): + temp_manual_EgX.append(g(1.0, x1) * p1) + + manual_EgX = math.fsum(temp_manual_EgX) + + # Now compare manual value against vectorized: + EgX = test.E(g) + + #print "manual_EgXY:", manual_EgXY + #print "EgXY:", EgXY + + print("test.E(g)"+str(test.E(g))) + print("Manual E[g(X)]"+ str(manual_EgX)) + print("max(abs(diff)):"+str(np.max(np.abs(test.E(g) - manual_EgX)))) + + + assert np.abs(manual_EgX - EgX) < 1e-16, "Eg(X) not equal manually calculated joint value." + print("Eg(X) *does* equal manually calculated joint value:") + print("\tE[g(X)] from class:"+str(EgX)+"\n\tand from manual calc:"+str(manual_EgX)) + + +# ## Improved Grid + +# Unlike most other definitions in this library, this is a simple function, not a class. This is implements the multi-exponential grid discussed in Carroll's Microeconomioc DSOP notes. + +# Note: this is taken almost verbatim from Chris Carroll's original Matlab code +# Define function: +def get_improved_grid(minval, maxval, size): + """ + This function finds values on a 1D grid such that the multi-exponential + growth rate* from each point to the next is constant (instead of, eg., + imposing a constant absolute gap between points). + + This is a nearly verbatim translation of Chris Carroll's DSOP code from + MATLAB to Python. + + *That is, exp(exp(exp(...))) for some number of exponentiations n. + + Args: + minval (float): Minimum grid value. + maxval (float): Maximum grid value. + size (int): Number of points in the grid. + Returns: + new_a_grid (list): eee-spaced grid. + Raises: + [] + """ + gMinMin = 0.01*minval; + gMaxMax = 10*maxval; + gridMin = log(1+log(1+log(1+gMaxMax))) + gridMax = (log(1+log(1+log(1+gMaxMax)))-log(1+log(1+log(1+gMinMin))))/size + index = log(1+log(1+log(1+gMinMin))) + (log(1+log(1+log(1+gMaxMax)))-log(1+log(1+log(1+gMinMin))))/size + i = 1 + point = 0 + points = [] + while point < gridMin: + point = point + indexThis + points.append(point) + i += 1 + + new_a_grid = exp(exp(exp(points)-1)-1)-1 + return new_a_grid + + + +# ## Utility functions + +# Python's object-oriented class structure makes it convenient to define a utility function and bundle with it properties such as the first and second derivative. "Resources" includes definitions for CRRA utility ("Utility") and exponential utility ("UtilityExponential"). +# +# (Note that I may need to eventually refactor "Utility" to "UtilityCRRA.") + +# + code_folding=[0, 5, 33, 51] +class Utility(object): + """ + Define a CRRA utility function in levels and + in first and second derivative. + """ + def __init__(self, gamma): + self.gamma = gamma + + def __call__(self, c): + if self.gamma == 1: + # Log case: + return( np.log(c) ) + else: + # Regular case: + return( c**(1.0 - self.gamma) / (1.0 - self.gamma) ) + + def prime(self, c): + if self.gamma == 1: + # Log case: + return(1.0/c) + else: + # Regular case: + return(c**-self.gamma) + + def prime2(self, c): + if self.gamma == 1: + # Log case: + return(-c**(-2.0)) + else: + # Regular case: + return(-self.gamma*c**(-self.gamma-1.0)) + + +class UtilityExponential: + """ + Define an exponential utility function in levels and + in first and second derivative. + """ + def __init__(self, gamma): + self.gamma = gamma + + def __call__(self, c): + return( c**self.gamma ) + + def prime(self, c): + return(self.gamma*c**(self.gamma-1.0)) + + def prime2(self, c): + return((self.gamma-1.0)*self.gamma*c**(self.gamma-2.0)) + + +class UtilityWithCMin: + """ + Define a CRRA utility function in levels and in + first derivative, with an imposed lower bound on + utility (via a lower bound on consumption.) + + Note that this is not the preferred way to impose + bounded utility; this is however included here + for backwards compatability with some simulations. + + """ + def __init__(self, gamma, cmin): + self.gamma = gamma + self.cmin = cmin + if gamma == 1: + # Log case: + self.__call__ = self.call_log + self.prime = self.prime_log + else: + # Regular case: + self.__call__ = self.call_regular + self.prime = self.prime_regular + + def prime_regular(self, c): + """ + The derivative when gamma != 1. + """ + return(c**-self.gamma) + + def prime_log(self, c): + """ + The derivative when gamma == 1. + """ + return(1.0/c) + + def call_regular(self, c): + return( np.maximum(c,self.cmin)**(1.0 - self.gamma) / (1.0 - self.gamma) ) + + def call_log(self, c): + return( np.log(np.maximum(c, self.cmin)) ) + +# - + + From 11c9a9fac127de9cc54d5523601c7714d6574ece Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 14 Jun 2023 08:08:25 -0400 Subject: [PATCH 05/14] working on generalizing FOC code --- HARK/algos/foc.py | 74 +++++++++++++++++++--------------- HARK/algos/tests/test_algos.py | 10 ++++- 2 files changed, 51 insertions(+), 33 deletions(-) diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index a1d219288..2c921c211 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -23,6 +23,8 @@ """ +## TODO: Action handling + def xndindex(ds, dims=None): """ There is currently no integrated way to iterate over an xarray.DataArray with its coordinate labels. @@ -63,11 +65,20 @@ def xz_grids_to_data_array( def optimal_policy_foc( - x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. - z_grid : Mapping[str, Sequence] = {}, - v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, - optimizer_args = None # TODO: For brettq. - ): + g : Callable[[Mapping, Mapping, Mapping], float], #= lambda x, z, a : {'a' : x['m'] - a['c']}, + r, # = lambda x, z, a : u(a['c']), + dr_da, # = lambda x, z, a: u.prime(a['c']), + dr_inv, # = lambda uP : (CRRAutilityP_inv(uP, rho),), + dg_dx = 1, ## Used in FOC method, step 5 + dg_da = -1, ## Used in FOC method, step 5 + x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. + z_grid : Mapping[str, Sequence] = {}, + v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, + discount = 1, + action_upper_bound = None, # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), + action_lower_bound = None, ## TODO: What are the default bounds? + optimizer_args = None # TODO: For brettq. + ): """ Given a grid over input and shock state values, and marginal output value function, @@ -108,13 +119,8 @@ def optimal_policy_foc( xz_grid_size = np.prod([len(xv) for xv in x_grid.values()]) * np.prod([len(zv) for zv in z_grid.values()]) - "dq_da(x,z,a) = dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da(x,z,a)" - - def foc(a): - a_vals = {an : av for an,av in zip(self.actions, (a,))} - return self.dq_da(x_vals, k_vals, a_vals, v_y_der) - - + def dq_da(x,z,a): + return dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da(x,z,a) # TODO: replace these with iterators.... #xz_iterator = xndindex(pi_data) @@ -150,30 +156,34 @@ def foc(a): y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) """ + def foc(a): + a_vals = {an : av for an,av in zip(self.actions, (a,))} + return dq_da(x_vals, z_vals, a_vals, v_y_der) + # these lower bounds as arugments to the - lower_bound = self.action_lower_bound(x_vals, k_vals) - upper_bound = self.action_upper_bound(x_vals, k_vals) + lower_bound = action_lower_bound(x_vals, z_vals) + upper_bound = action_upper_bound(x_vals, z_vals) ## what if no lower bound? q_der_lower = None if lower_bound[0] is not None: - q_der_lower = self.d_q_d_a( - x_vals, - k_vals, - self.action_zip(lower_bound), - v_y_der - ) + q_der_lower = dq_da( + x_vals, + k_vals, + self.action_zip(lower_bound), + v_y_der + ) else: lower_bound[0] = 1e-12 ## a really high number! q_der_upper = None if upper_bound[0] is not None: - q_der_upper = self.d_q_d_a( - x_vals, - k_vals, - self.action_zip(upper_bound), - v_y_der - ) + q_der_upper = dq_da( + x_vals, + z_vals, + self.action_zip(upper_bound), + v_y_der + ) else: upper_bound[0] = 1e12 ## a really high number! @@ -181,13 +191,13 @@ def foc(a): if q_der_lower is not None and q_der_upper is not None and q_der_lower < 0 and q_der_upper < 0: a0 = lower_bound - pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) - q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_lower) + pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_lower) elif q_der_lower is not None and q_der_upper is not None and q_der_lower > 0 and q_der_upper > 0: a0 = upper_bound - pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) - q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_upper) + pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_upper) else: ## Better exception handling here ## asserting that Q is concave @@ -202,7 +212,7 @@ def foc(a): ) if root_res.converged: - pi_data.sel(**x_vals, **k_vals).variable.data.put(0, a0) + pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) q_der_xz = self.dq_da( x_vals, @@ -211,7 +221,7 @@ def foc(a): v_y_der ) - q_der_data.sel(**x_vals, **k_vals).variable.data.put(0, q_der_xk) + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) else: print(f"Rootfinding failure at {x_vals}, {z_vals}.") print(root_res) diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 877745dba..dff55e675 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -86,8 +86,16 @@ def consumption_v_y_der(y : Mapping[str,Any]): pi_star, q_der, y_data = optimal_policy_foc( + g, + r, + dr_da, + dr_inv, + dg_dx, + dg_da, {'m' : self.mVec}, - v_y_der = consumption_v_y_der + v_y_der = consumption_v_y_der, + action_upper_bound = None, # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), + action_lower_bound = None, ) self.assertTrue(np.all(self.cVec2 == pi_star.values)) From 8fbca335732a29184edf71827b61cf8865007434 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 14 Jun 2023 08:50:11 -0400 Subject: [PATCH 06/14] fixing references to self.actions and self.action_zip --- HARK/algos/foc.py | 29 ++++++++++++++++++----------- HARK/algos/tests/test_algos.py | 3 ++- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index 2c921c211..a8fdd9f8c 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -1,7 +1,8 @@ +from dataclasses import field import itertools import numpy as np from scipy.optimize import minimize, brentq -from typing import Callable, Mapping, Sequence +from typing import Callable, Mapping, Sequence, Tuple import xarray as xr """ @@ -66,6 +67,7 @@ def xz_grids_to_data_array( def optimal_policy_foc( g : Callable[[Mapping, Mapping, Mapping], float], #= lambda x, z, a : {'a' : x['m'] - a['c']}, + actions, r, # = lambda x, z, a : u(a['c']), dr_da, # = lambda x, z, a: u.prime(a['c']), dr_inv, # = lambda uP : (CRRAutilityP_inv(uP, rho),), @@ -75,8 +77,8 @@ def optimal_policy_foc( z_grid : Mapping[str, Sequence] = {}, v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, discount = 1, - action_upper_bound = None, # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), - action_lower_bound = None, ## TODO: What are the default bounds? + action_upper_bound : Callable[[Mapping, Mapping], Sequence[float]] = field(default = None), # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), + action_lower_bound : Callable[[Mapping, Mapping], Sequence[float]] = field(default = None), optimizer_args = None # TODO: For brettq. ): """ @@ -126,7 +128,12 @@ def dq_da(x,z,a): #xz_iterator = xndindex(pi_data) #import pdb; pdb.set_trace() - + def action_zip(self, a : Tuple): + """ + Wraps a tuple of values for an action in a dictionary with labels. + Useful for converting between forms of model equations. + """ + return {an : av for an,av in zip(self.actions, a)} for x_point in itertools.product(*x_grid.values()): x_vals = {k : v for k, v in zip(x_grid.keys() , x_point)} @@ -157,7 +164,7 @@ def dq_da(x,z,a): """ def foc(a): - a_vals = {an : av for an,av in zip(self.actions, (a,))} + a_vals = action_zip((a,)) return dq_da(x_vals, z_vals, a_vals, v_y_der) # these lower bounds as arugments to the @@ -169,8 +176,8 @@ def foc(a): if lower_bound[0] is not None: q_der_lower = dq_da( x_vals, - k_vals, - self.action_zip(lower_bound), + z_vals, + action_zip(lower_bound), v_y_der ) else: @@ -181,7 +188,7 @@ def foc(a): q_der_upper = dq_da( x_vals, z_vals, - self.action_zip(upper_bound), + action_zip(upper_bound), v_y_der ) else: @@ -214,10 +221,10 @@ def foc(a): if root_res.converged: pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) - q_der_xz = self.dq_da( + q_der_xz = dq_da( x_vals, z_vals, - self.action_zip((a0,)), # actions are scalar + action_zip((a0,)), # actions are scalar v_y_der ) @@ -232,7 +239,7 @@ def foc(a): q_der_xz = dq_da( x_vals, z_vals, - self.action_zip((root_res.root,)), # actions are scalar + action_zip((root_res.root,)), # actions are scalar v_y_der ) diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index dff55e675..58cf7374e 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -87,6 +87,7 @@ def consumption_v_y_der(y : Mapping[str,Any]): pi_star, q_der, y_data = optimal_policy_foc( g, + ['a'], r, dr_da, dr_inv, @@ -95,7 +96,7 @@ def consumption_v_y_der(y : Mapping[str,Any]): {'m' : self.mVec}, v_y_der = consumption_v_y_der, action_upper_bound = None, # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), - action_lower_bound = None, + action_lower_bound = lambda x, z: 0, ) self.assertTrue(np.all(self.cVec2 == pi_star.values)) From 0860539f2482d6c139fd165faa8ce4c8f7054344 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 14 Jun 2023 09:22:40 -0400 Subject: [PATCH 07/14] more fixes --- HARK/algos/foc.py | 11 ++++--- HARK/algos/tests/test_algos.py | 55 ++++++++++++++++++++++++++-------- 2 files changed, 50 insertions(+), 16 deletions(-) diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index a8fdd9f8c..2490c4641 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -121,19 +121,22 @@ def optimal_policy_foc( xz_grid_size = np.prod([len(xv) for xv in x_grid.values()]) * np.prod([len(zv) for zv in z_grid.values()]) - def dq_da(x,z,a): - return dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da(x,z,a) + def dq_da(x,z,a, v_y_der): + print(dr_da, x, z, a, discount, v_y_der) + return dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da # could optionally be dg_da(x,z,a) # TODO: replace these with iterators.... #xz_iterator = xndindex(pi_data) #import pdb; pdb.set_trace() - def action_zip(self, a : Tuple): + def action_zip(a : Tuple): """ Wraps a tuple of values for an action in a dictionary with labels. Useful for converting between forms of model equations. + + References 'actions' argument of optimal_policy_foc() """ - return {an : av for an,av in zip(self.actions, a)} + return {an : av for an,av in zip(actions, a)} for x_point in itertools.product(*x_grid.values()): x_vals = {k : v for k, v in zip(x_grid.keys() , x_point)} diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 58cf7374e..b3e260c95 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -6,14 +6,45 @@ from HARK.algos.foc import optimal_policy_foc from HARK.gothic.gothic_class import gothic -from HARK.rewards import CRRAutilityP_inv +from HARK.gothic.resources import ( + Utility, + DiscreteApproximation, + DiscreteApproximationTwoIndependentDistribs, +) + +from HARK.rewards import CRRAutilityP, CRRAutilityP_inv import numpy as np import os +from scipy import stats from typing import Any, Mapping __location__ = os.path.realpath( os.path.join(os.getcwd(), os.path.dirname(__file__))) + +# From SolvingMicroDSOPS +# Set up general parameters: + +rho = 2.0 ### CRRA coefficient +beta = 0.96 ### discount factor +gamma = Gamma = np.array([1.0]) # permanent income growth factor +# A one-element "time series" array +# (the array structure needed for gothic class below) +R = 1.02 ## Risk free interest factor + +# Define utility: +u = Utility(rho) + +theta_sigma = 0.5 +theta_mu = -0.5 * (theta_sigma**2) +theta_z = stats.lognorm( + theta_sigma, 0, np.exp(theta_mu) +) # Create "frozen" distribution instance +theta_grid_N = 7 ### how many grid points to approximate this continuous distribution + +theta = DiscreteApproximation( + N=theta_grid_N, cdf=theta_z.cdf, pdf=theta_z.pdf, invcdf=theta_z.ppf +) """ @@ -57,15 +88,15 @@ def setUp(self): self.cVec2 = np.load(os.path.join(__location__, "smdsops_cVec2.npy")) - def test_x(self): + def test_optimal_policy_foc(self): - g = lambda x, k, a : {'a' : x['m'] - a['c']}, - dg_dx = 1, ## Used in FOC method, step 5 - dg_da = -1, ## Used in FOC method, step 5 - g_inv = lambda y, a : {'m' : y['a'] + a['c']}, ## Used in EGM method, step 8 - r = lambda x, k, a : u(a['c']), - dr_da = lambda x, k, a: u.prime(a['c']), - dr_inv = lambda uP : (CRRAutilityP_inv(uP, rho),), + g = lambda x, k, a : {'a' : x['m'] - a['c']} + dg_dx = 1 ## Used in FOC method, step 5 + dg_da = -1 ## Used in FOC method, step 5 + g_inv = lambda y, a : {'m' : y['a'] + a['c']} ## Used in EGM method, step 8 + r = lambda x, k, a : u(a['c']) + dr_da = lambda x, k, a: (CRRAutilityP(a['c'], rho),) # u.prime(a['c']) + dr_inv = lambda uP : (CRRAutilityP_inv(uP, rho),) action_upper_bound = lambda x, k: (x['m'] + gamma[0] * theta.X[0] / R,), @@ -87,7 +118,7 @@ def consumption_v_y_der(y : Mapping[str,Any]): pi_star, q_der, y_data = optimal_policy_foc( g, - ['a'], + ['c'], r, dr_da, dr_inv, @@ -95,8 +126,8 @@ def consumption_v_y_der(y : Mapping[str,Any]): dg_da, {'m' : self.mVec}, v_y_der = consumption_v_y_der, - action_upper_bound = None, # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), - action_lower_bound = lambda x, z: 0, + action_upper_bound = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), + action_lower_bound = lambda x, z: (0,), ) self.assertTrue(np.all(self.cVec2 == pi_star.values)) From 01347c7cf3c29a80035085e31c4788570bd34c61 Mon Sep 17 00:00:00 2001 From: sb Date: Wed, 14 Jun 2023 09:30:51 -0400 Subject: [PATCH 08/14] passing FOC test! --- HARK/algos/foc.py | 3 +-- HARK/algos/tests/test_algos.py | 7 ++++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index 2490c4641..7eba581d5 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -122,7 +122,6 @@ def optimal_policy_foc( xz_grid_size = np.prod([len(xv) for xv in x_grid.values()]) * np.prod([len(zv) for zv in z_grid.values()]) def dq_da(x,z,a, v_y_der): - print(dr_da, x, z, a, discount, v_y_der) return dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da # could optionally be dg_da(x,z,a) # TODO: replace these with iterators.... @@ -249,7 +248,7 @@ def foc(a): q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) acts = np.atleast_1d(pi_data.sel(**x_vals, **z_vals).values) - y = g(x_vals, z_vals, self.action_zip(acts)) + y = g(x_vals, z_vals, action_zip(acts)) y_n = np.array([y[k] for k in y]) y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index b3e260c95..677196cd0 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -5,7 +5,7 @@ import unittest from HARK.algos.foc import optimal_policy_foc -from HARK.gothic.gothic_class import gothic +from HARK.gothic.gothic_class import Gothic from HARK.gothic.resources import ( Utility, DiscreteApproximation, @@ -45,9 +45,10 @@ theta = DiscreteApproximation( N=theta_grid_N, cdf=theta_z.cdf, pdf=theta_z.pdf, invcdf=theta_z.ppf ) -""" +gothic = Gothic(u, beta, rho, gamma, R, theta) +""" g = lambda x, k, a : {'a' : x['m'] - a['c']}, dg_dx = 1, ## Used in FOC method, step 5 dg_da = -1, ## Used in FOC method, step 5 @@ -130,7 +131,7 @@ def consumption_v_y_der(y : Mapping[str,Any]): action_lower_bound = lambda x, z: (0,), ) - self.assertTrue(np.all(self.cVec2 == pi_star.values)) + self.assertTrue(np.all(abs(self.cVec2 - pi_star.values) < 1e-12)) class egm_test(unittest.TestCase): """ From 00a5a7fcff73ea1885bdf4bd1991b54bee5f2ad9 Mon Sep 17 00:00:00 2001 From: sb Date: Fri, 16 Jun 2023 10:52:44 -0400 Subject: [PATCH 09/14] working egm function and test --- HARK/algos/egm.py | 117 +++++++++++++++++---------------- HARK/algos/tests/test_algos.py | 30 +++++++-- 2 files changed, 87 insertions(+), 60 deletions(-) diff --git a/HARK/algos/egm.py b/HARK/algos/egm.py index 93ce77c1c..184fcd4e4 100644 --- a/HARK/algos/egm.py +++ b/HARK/algos/egm.py @@ -1,53 +1,60 @@ - - -def analytic_pi_star_y(self, - y, - v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0 - ): +from dataclasses import field +import itertools +import numpy as np +from scipy.optimize import minimize, brentq +from typing import Callable, Mapping, Sequence, Tuple +import xarray as xr + +def analytic_pi_y_star( + y, + v_y_der : Callable[[Mapping, Mapping, Mapping], float], + dr_da_inv : Callable[[float], float], + dg_da = -1.0, + discount = 1.0 + ): """ The optimal action that results in output values y. - Available only with reward_der_inv and transition_inv - are well-defined. - """ - if self.reward_der_inv is None: - raise Exception("No inverse marginal reward function found.") - - if self.transition_inv is None: - raise Exception("No inverse transition function found. ") + Assumes: + - dg_da is a constant + - discount factor is constant - if self.transition_der_a is None or not( - isinstance(self.transition_der_a, float) or isinstance(self.transition_der_a, int) + Params + ------ + """ + if dg_da is None or not( + isinstance(dg_da, float) or isinstance(dg_da, int) ): - raise Exception(f"No constant transition derivative found. transition_der_a is {self.transition_der_a}") + raise Exception(f"No constant transition derivative found. transition_der_a is {dg_da}") - if not isinstance(self.discount, float) or isinstance(self.discount, int): + if not (isinstance(discount, float) or isinstance(discount, int)): raise Exception("Analytic pi_star_y requires constant discount factor (rendering B' = 0).") - ### TEST: available T_der as constant. - v_y_der_at_y = v_y_der(y) if isinstance(v_y_der_at_y, xr.DataArray): v_y_der_at_y = v_y_der_at_y.values # np.atleast1darray() ? - if 0 > v_y_der_at_y: - raise Exception(f"Negative marginal value {v_y_der_at_y} computes at y value of {y}. Reward is {- self.discount * self.transition_der_a * v_y_der_at_y}") - - return self.reward_der_inv(- self.discount * self.transition_der_a * v_y_der_at_y) - -def optimal_policy_egm(self, - y_grid : Mapping[str, Sequence] = {}, - v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, - ): + raise Exception(f"Negative marginal value {v_y_der_at_y} computes at y value of {y}. Reward is {- discount * dg_da * v_y_der_at_y}") + + return dr_da_inv(- discount * dg_da * v_y_der_at_y) + +def egm( + inputs, + actions, + g_inv : Callable[[Mapping, Mapping, Mapping], float], + dr_da_inv, # = lambda uP : (CRRAutilityP_inv(uP, rho),), + dg_da = -1, + y_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. + v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, + discount = 1, + ): """ Given a grid over output and marginal output value function, compute the optimal action. - ## NO SHOCKS ALLOWED ## - This depends on the stage having an *inverse marginal reward function* and *inverse transition function*. @@ -59,12 +66,8 @@ def optimal_policy_egm(self, ### and... T' = -1 ??? """ - if self.reward_der_inv is None: - raise Exception("No inverse marginal reward function found. EGM requires reward_der_inv defined for this stage.") - - if self.transition_inv is None: - raise Exception("No inverse transition function found. EGM requires transition_inv defined for this stage.") - + ## can be functionalized out once we settle + ## on grid-to-DataArray conversions pi_y_data = xr.DataArray( np.zeros([len(v) for v in y_grid.values()]), dims = y_grid.keys(), @@ -75,46 +78,48 @@ def optimal_policy_egm(self, x_val_data = [] a_val_data = [] + ## duplicated from foc.py; move to a shared helper library? + def action_zip(a : Tuple): + """ + Wraps a tuple of values for an action in a dictionary with labels. + Useful for converting between forms of model equations. + + References 'actions' argument of optimal_policy_foc() + """ + return {an : av for an,av in zip(actions, a)} + + # can I pass in the grid, rather than iterate? for y_point in itertools.product(*y_grid.values()): y_vals = {k : v for k, v in zip(y_grid.keys() , y_point)} - acts = self.analytic_pi_star_y(y_vals, v_y_der) + acts = analytic_pi_y_star( + y_vals, + v_y_der, + dr_da_inv, + dg_da, + discount + ) pi_y_data.sel(**y_vals).variable.data.put(0, acts) - x_vals = self.transition_inv(y_vals, self.action_zip(acts)) + x_vals = g_inv(y_vals, action_zip(acts)) x_val_data.append(x_vals) - a_val_data.append(self.action_zip(acts)) + a_val_data.append(action_zip(acts)) ## TODO is this dealing with repeated values? x_coords = { x : np.array([xv[x] for xv in x_val_data]) for x - in self.inputs + in inputs } - for x_label in self.solution_points.coords: - for sol_x_val in self.solution_points.coords[x_label]: - if sol_x_val.values not in x_coords[x_label]: - ii = np.searchsorted(x_coords[x_label], sol_x_val) - x_coords[x_label] = np.insert(x_coords[x_label], ii, sol_x_val) - pi_data = xr.DataArray( np.zeros([len(v) for v in x_coords.values()]), dims = x_coords.keys(), coords = x_coords ) - if 'pi*' in self.solution_points: - for x_point in itertools.product(*x_coords.values()): - x_vals = {k : v for k, v in zip(x_coords.keys() , x_point)} - - if label_index_in_dataset(x_vals, self.solution_points['pi*']): - acts = np.atleast_1d(self.solution_points['pi*'].sel(x_vals)) - - pi_data.sel(**x_vals, **{}).variable.data.put(0, acts) - for i, x_vals in enumerate(x_val_data): x_vals = x_val_data[i] a_vals = a_val_data[i] diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 677196cd0..00dd4a278 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -4,6 +4,7 @@ # Bring in modules we need import unittest +from HARK.algos.egm import egm from HARK.algos.foc import optimal_policy_foc from HARK.gothic.gothic_class import Gothic from HARK.gothic.resources import ( @@ -97,7 +98,7 @@ def test_optimal_policy_foc(self): g_inv = lambda y, a : {'m' : y['a'] + a['c']} ## Used in EGM method, step 8 r = lambda x, k, a : u(a['c']) dr_da = lambda x, k, a: (CRRAutilityP(a['c'], rho),) # u.prime(a['c']) - dr_inv = lambda uP : (CRRAutilityP_inv(uP, rho),) + dr_da_inv = lambda uP : (CRRAutilityP_inv(uP, rho),) action_upper_bound = lambda x, k: (x['m'] + gamma[0] * theta.X[0] / R,), @@ -115,14 +116,13 @@ def test_optimal_policy_foc(self): def consumption_v_y_der(y : Mapping[str,Any]): return gothic.VP_Tminus1(y['a']) - pi_star, q_der, y_data = optimal_policy_foc( g, ['c'], r, dr_da, - dr_inv, + dr_da_inv, dg_dx, dg_da, {'m' : self.mVec}, @@ -159,4 +159,26 @@ def consumption_v_y_der(y : Mapping[str,Any]): pi.interp({'m' : mVec_egm}) - cFunc_egm(mVec_egm) == 0 """ - self.assertTrue(np.all(self.aVec + self.cVec_egm == self.mVec_egm)) \ No newline at end of file + g = lambda x, k, a : {'a' : x['m'] - a['c']} + dg_dx = 1 ## Used in FOC method, step 5 + dg_da = -1 ## Used in FOC method, step 5 + g_inv = lambda y, a : {'m' : y['a'] + a['c']} ## Used in EGM method, step 8 + r = lambda x, k, a : u(a['c']) + dr_da = lambda x, k, a: (CRRAutilityP(a['c'], rho),) # u.prime(a['c']) + dr_da_inv = lambda uP : (CRRAutilityP_inv(uP, rho),) + + def consumption_v_y_der(y : Mapping[str,Any]): + return gothic.VP_Tminus1(y['a']) + + pi_data, pi_y_data = egm( + ['m'], + ['c'], + g_inv, + dr_da_inv, + dg_da, + y_grid = {'a' : self.aVec}, + v_y_der = consumption_v_y_der, + ) + + self.assertTrue(np.all(self.cVec_egm == pi_y_data.values)) + self.assertTrue(np.all(self.mVec_egm == pi_data.coords['m'].values)) From 1a9846e6545b6dde465997b3abb5b7fe09bfa6e3 Mon Sep 17 00:00:00 2001 From: sb Date: Fri, 16 Jun 2023 11:21:26 -0400 Subject: [PATCH 10/14] prints to test the test for failing GitHUb executed tests --- HARK/algos/tests/test_algos.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 00dd4a278..18dfdc9e9 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -180,5 +180,9 @@ def consumption_v_y_der(y : Mapping[str,Any]): v_y_der = consumption_v_y_der, ) + if not np.all(self.cVec_egm == pi_y_data.values): + print(self.cVec_egm) + print(pi_y_data.values) + self.assertTrue(np.all(self.cVec_egm == pi_y_data.values)) self.assertTrue(np.all(self.mVec_egm == pi_data.coords['m'].values)) From 5e8f5cb415bc1ef7d4033a90b446f91dbda28a86 Mon Sep 17 00:00:00 2001 From: sb Date: Fri, 16 Jun 2023 11:41:17 -0400 Subject: [PATCH 11/14] assert messages > print statements; up to 1e-12 precision not equality? --- HARK/algos/tests/test_algos.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 18dfdc9e9..3007559e9 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -180,9 +180,11 @@ def consumption_v_y_der(y : Mapping[str,Any]): v_y_der = consumption_v_y_der, ) - if not np.all(self.cVec_egm == pi_y_data.values): - print(self.cVec_egm) - print(pi_y_data.values) - - self.assertTrue(np.all(self.cVec_egm == pi_y_data.values)) - self.assertTrue(np.all(self.mVec_egm == pi_data.coords['m'].values)) + self.assertTrue( + np.all(abs(self.cVec_egm - pi_y_data.values < 1e-12)), + f"Target values: {[x for x in self.cVec_egm]}; Test values: {[x for x in pi_y_data.values]}" + ) + self.assertTrue( + np.all(abs(self.mVec_egm - pi_data.coords['m'].values < 1e-12 )), + f"Target values: {self.mVec_egm}; Test values: {pi_data.coords['m'].values}" + ) From 4935438ee94f67b884e4217032c32bffdffd9f30 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 19 Jun 2023 17:02:54 -0400 Subject: [PATCH 12/14] better documentation for algos --- Documentation/reference/tools/algos.rst | 8 + Documentation/reference/tools/algos/egm.rst | 7 + Documentation/reference/tools/algos/foc.rst | 7 + Documentation/reference/tools/index.rst | 1 + HARK/algos/egm.py | 79 +++++-- HARK/algos/foc.py | 215 ++++++++++---------- 6 files changed, 192 insertions(+), 125 deletions(-) create mode 100644 Documentation/reference/tools/algos.rst create mode 100644 Documentation/reference/tools/algos/egm.rst create mode 100644 Documentation/reference/tools/algos/foc.rst diff --git a/Documentation/reference/tools/algos.rst b/Documentation/reference/tools/algos.rst new file mode 100644 index 000000000..b6451c4b3 --- /dev/null +++ b/Documentation/reference/tools/algos.rst @@ -0,0 +1,8 @@ +algos +-------- + +.. toctree:: + :maxdepth: 3 + + algos/egm + algos/foc \ No newline at end of file diff --git a/Documentation/reference/tools/algos/egm.rst b/Documentation/reference/tools/algos/egm.rst new file mode 100644 index 000000000..f9daee533 --- /dev/null +++ b/Documentation/reference/tools/algos/egm.rst @@ -0,0 +1,7 @@ +algos.egm +-------- + +.. automodule:: HARK.algos.egm + :members: + :undoc-members: + :show-inheritance: diff --git a/Documentation/reference/tools/algos/foc.rst b/Documentation/reference/tools/algos/foc.rst new file mode 100644 index 000000000..cd80d3aee --- /dev/null +++ b/Documentation/reference/tools/algos/foc.rst @@ -0,0 +1,7 @@ +algos.foc +-------- + +.. automodule:: HARK.algos.foc + :members: + :undoc-members: + :show-inheritance: diff --git a/Documentation/reference/tools/index.rst b/Documentation/reference/tools/index.rst index 660c606cd..851a0a3ec 100644 --- a/Documentation/reference/tools/index.rst +++ b/Documentation/reference/tools/index.rst @@ -4,6 +4,7 @@ Tools .. toctree:: :maxdepth: 3 + algos core dcegm distribution diff --git a/HARK/algos/egm.py b/HARK/algos/egm.py index 184fcd4e4..46d5deb81 100644 --- a/HARK/algos/egm.py +++ b/HARK/algos/egm.py @@ -13,14 +13,32 @@ def analytic_pi_y_star( discount = 1.0 ): """ - The optimal action that results in output values y. + The action which, when chosen, results in the post-state value Y. - Assumes: - - dg_da is a constant - - discount factor is constant + Parameters + ----------- + y: + The post-state value Y. + TODO: type signature of input y - Params - ------ + v_y_der: + Derivative post-value function over post-states Y. + + dr_da_inv: + Derivative of the reward function r with respect to actions A, inverted. + + dg_da: + Derivateive of the transition function g with respect to actions A. Must be a constant. + + discount: + Discount factor. Must be a constant. + + Returns + ------- + + The actions chosen that result in post-state Y. + + TODO: Type signature of output """ if dg_da is None or not( isinstance(dg_da, float) or isinstance(dg_da, int) @@ -41,8 +59,8 @@ def analytic_pi_y_star( return dr_da_inv(- discount * dg_da * v_y_der_at_y) def egm( - inputs, - actions, + inputs : Sequence[str], + actions : Sequence[str], g_inv : Callable[[Mapping, Mapping, Mapping], float], dr_da_inv, # = lambda uP : (CRRAutilityP_inv(uP, rho),), dg_da = -1, @@ -51,19 +69,44 @@ def egm( discount = 1, ): """ - Given a grid over output - and marginal output value function, - compute the optimal action. + Given post-states Y and information about an agent's problem, + compute the actions which, if chosen, result in those post-states, + and the corresponding states X. + + This is a method of computing data about the optimal policy of an + agent that does not use rootfinding. + + Parameters + ------------ + inputs: + Ordered labels of the state variables. + + actions: + Ordered labels of the action variables. + + g_inv: + Inverse of the transition function g. + + dr_da_inv: + Derivative of the reward function r with respect to actions A, inverted. + + dg_da: + Derivative of the transition function g with respect to action A. + + y_grid: + + v_y_der: + Derivative of post-value function v_y over post states Y + + discount: + A discount fact. Scalar. - This depends on the stage having an - *inverse marginal reward function* - and *inverse transition function*. + Returns + -------- - Does not use rootfinding! - The 'grid' over the input + pi_data: - ### ASSUMES: No discounting this phase, - ### and... T' = -1 ??? + pi_y_data: """ ## can be functionalized out once we settle diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index 7eba581d5..7b906daea 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -45,13 +45,15 @@ def xndindex(ds, dims=None): raise ValueError("Invalid dimension '{}'. Available dimensions {}".format(d, ds.dims)) iter_dict = {k:v for k,v in ds.sizes.items() if k in dims} - for d,k in zip(repeat(tuple(iter_dict.keys())),zip(np.ndindex(tuple(iter_dict.values())))): + for d,k in zip(itertools.repeat(tuple(iter_dict.keys())),zip(np.ndindex(tuple(iter_dict.values())))): yield {k:l for k,l in zip(d,k[0])} +Grid = Mapping[str, Sequence] + def xz_grids_to_data_array( - x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. - z_grid : Mapping[str, Sequence] = {} + x_grid : Grid = {}, ## TODO: Better data structure here. + z_grid : Grid = {} ): coords = {**x_grid, **z_grid} @@ -66,20 +68,19 @@ def xz_grids_to_data_array( def optimal_policy_foc( - g : Callable[[Mapping, Mapping, Mapping], float], #= lambda x, z, a : {'a' : x['m'] - a['c']}, - actions, - r, # = lambda x, z, a : u(a['c']), - dr_da, # = lambda x, z, a: u.prime(a['c']), - dr_inv, # = lambda uP : (CRRAutilityP_inv(uP, rho),), - dg_dx = 1, ## Used in FOC method, step 5 - dg_da = -1, ## Used in FOC method, step 5 - x_grid : Mapping[str, Sequence] = {}, ## TODO: Better data structure here. - z_grid : Mapping[str, Sequence] = {}, - v_y_der : Callable[[Mapping, Mapping, Mapping], float] = lambda x : 0, + g : Callable[[Grid, Grid, Grid], float], + actions: Sequence[str], + r : Callable[[Grid, Grid, Grid], float], + dr_da, + dr_inv, + dg_dx = 1, + dg_da = -1, + x_grid : Grid = {}, + z_grid : Grid = {}, + v_y_der : Callable[[Grid, Grid, Grid], float] = lambda x : 0, discount = 1, - action_upper_bound : Callable[[Mapping, Mapping], Sequence[float]] = field(default = None), # = lambda x, z: (x['m'] + gamma[0] * theta.X[0] / R,), - action_lower_bound : Callable[[Mapping, Mapping], Sequence[float]] = field(default = None), - optimizer_args = None # TODO: For brettq. + action_upper_bound : Callable[[Grid, Grid], Sequence[float]] = field(default = None), + action_lower_bound : Callable[[Grid, Grid], Sequence[float]] = field(default = None), ): """ Given a grid over input and shock state values, @@ -105,23 +106,62 @@ def optimal_policy_foc( Parameters ----------- - x_grid : Mapping[str, Sequence] - A mapping of state variable names to a sequence of grid values - - z_grid: Mapping[str, Sequence] - A mapping of shock variable names to grid values + g : + Transition functiong g. + + actions: + Ordered labels of action variables. + + r : + Reward function r. + + dr_da : + Derivative of reward function r with respect to actions A. + + dr_inv : + Inverse of derivative of reward function + + dg_dx: + Derivative of transition function g with respect to states X. + + dg_da: + Derivative of transition function g with respect to actions A. + + x_grid: + + z_grid: + + v_y_der: + Derivative of the post-value function v_y with respect to post-states Y. + + discount: + + action_upper_bound: + action_lower_bound: + + Returns + -------- + + pi_star_data: + Data for the optimal policy pi_star. + + q_der_data: + Data for the derivative of the action value function function q. + + y_data: """ # Set up data arrays with coordinates based on the grid. - pi_data = xz_grids_to_data_array(x_grid, z_grid) ## TODO: Default value for pi should be nan. + pi_star_data = xz_grids_to_data_array(x_grid, z_grid) ## TODO: Default value for pi should be nan. q_der_data = xz_grids_to_data_array(x_grid, z_grid) ## May need to expand this for multivariate y y_data = xz_grids_to_data_array(x_grid, z_grid) - xz_grid_size = np.prod([len(xv) for xv in x_grid.values()]) * np.prod([len(zv) for zv in z_grid.values()]) - def dq_da(x,z,a, v_y_der): + """ + Derivative of the action-value function q with respect to actions a + """ return dr_da(x,z,a) + discount * v_y_der(g(x,z,a)) * dg_da # could optionally be dg_da(x,z,a) # TODO: replace these with iterators.... @@ -142,40 +182,15 @@ def action_zip(a : Tuple): for z_point in itertools.product(*z_grid.values()): z_vals = {k : v for k, v in zip(z_grid.keys() , z_point)} - - # repeated code with the other optimizer -- can be functionalized out? - """ - if len(self.actions) == 0: - q_der_xz = self.d_q_d_a( - x_vals, - z_vals, - {}, # no actions - v_y_der - ) - - # this should just be the default value for pi - pi_data.sel(**x_vals, **z_vals).variable.data.put(0, np.nan) - - # then this can be looped in differently. - q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) - - ## TODO: define g - y = g(x_vals, z_vals, self.action_zip([np.nan])) - y_n = np.array([y[k] for k in y]) - y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) - """ - + def foc(a): a_vals = action_zip((a,)) return dq_da(x_vals, z_vals, a_vals, v_y_der) - # these lower bounds as arugments to the - lower_bound = action_lower_bound(x_vals, z_vals) - upper_bound = action_upper_bound(x_vals, z_vals) - ## what if no lower bound? q_der_lower = None - if lower_bound[0] is not None: + if action_lower_bound is not None: + lower_bound = action_lower_bound(x_vals, z_vals) q_der_lower = dq_da( x_vals, z_vals, @@ -183,10 +198,11 @@ def foc(a): v_y_der ) else: - lower_bound[0] = 1e-12 ## a really high number! + lower_bound = np.array([-1e-12]) ## a really low number! q_der_upper = None - if upper_bound[0] is not None: + if action_upper_bound is not None: + upper_bound = action_upper_bound(x_vals, z_vals) q_der_upper = dq_da( x_vals, z_vals, @@ -194,62 +210,47 @@ def foc(a): v_y_der ) else: - upper_bound[0] = 1e12 ## a really high number! - - ## TODO: Better handling of case when there is a missing bound? - if q_der_lower is not None and q_der_upper is not None and q_der_lower < 0 and q_der_upper < 0: - a0 = lower_bound - - pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) - q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_lower) - elif q_der_lower is not None and q_der_upper is not None and q_der_lower > 0 and q_der_upper > 0: - a0 = upper_bound - - pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) - q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_upper) - else: - ## Better exception handling here - ## asserting that Q is concave - if q_der_lower is not None and q_der_upper is not None and not(q_der_lower > 0 and q_der_upper < 0): - raise Exception("Cannot solve for optimal policy with FOC if Q is not concave!") - - a0, root_res = brentq( - foc, - lower_bound[0], # only works with scalar actions - upper_bound[0], # only works with scalar actions - full_output = True + upper_bound = np.array([1e12]) ## a really high number! + + if q_der_lower is not None and q_der_upper is not None and not(q_der_lower > 0 and q_der_upper < 0): + raise Exception("Cannot solve for optimal policy with FOC if Q is not concave!") + + a0, root_res = brentq( + foc, + lower_bound[0], # only works with scalar actions + upper_bound[0], # only works with scalar actions + full_output = True + ) + + if root_res.converged: + pi_star_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) + + q_der_xz = dq_da( + x_vals, + z_vals, + action_zip((a0,)), # actions are scalar + v_y_der ) + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) + else: + print(f"Rootfinding failure at {x_vals}, {z_vals}.") + print(root_res) - if root_res.converged: - pi_data.sel(**x_vals, **z_vals).variable.data.put(0, a0) - - q_der_xz = dq_da( - x_vals, - z_vals, - action_zip((a0,)), # actions are scalar - v_y_der - ) - - q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) - else: - print(f"Rootfinding failure at {x_vals}, {z_vals}.") - print(root_res) - - #raise Exception("Failed to optimize.") - pi_data.sel(**x_vals, **z_vals).variable.data.put(0, root_res.root) + #raise Exception("Failed to optimize.") + pi_star_data.sel(**x_vals, **z_vals).variable.data.put(0, root_res.root) - q_der_xz = dq_da( - x_vals, - z_vals, - action_zip((root_res.root,)), # actions are scalar - v_y_der - ) + q_der_xz = dq_da( + x_vals, + z_vals, + action_zip((root_res.root,)), # actions are scalar + v_y_der + ) - q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) + q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) - acts = np.atleast_1d(pi_data.sel(**x_vals, **z_vals).values) - y = g(x_vals, z_vals, action_zip(acts)) - y_n = np.array([y[k] for k in y]) - y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) + acts = np.atleast_1d(pi_star_data.sel(**x_vals, **z_vals).values) + y = g(x_vals, z_vals, action_zip(acts)) + y_n = np.array([y[k] for k in y]) + y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) - return pi_data, q_der_data, y_data \ No newline at end of file + return pi_star_data, q_der_data, y_data \ No newline at end of file From 8d095e9f5fdb348e97fff5d942d96bb4a4d25a62 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 26 Jun 2023 11:56:12 -0400 Subject: [PATCH 13/14] mathematical documentation for HARK.algos --- Documentation/conf.py | 6 + Documentation/reference/tools/algos/egm.rst | 2 +- Documentation/reference/tools/algos/foc.rst | 2 +- HARK/algos/egm.py | 55 ++++++- HARK/algos/foc.py | 151 +++++++++++--------- HARK/algos/tests/test_algos.py | 2 +- 6 files changed, 147 insertions(+), 71 deletions(-) diff --git a/Documentation/conf.py b/Documentation/conf.py index 27af3addd..3ee0caeb0 100644 --- a/Documentation/conf.py +++ b/Documentation/conf.py @@ -59,3 +59,9 @@ # nbsphinx configuration nbsphinx_execute = "never" # This is currently not working + +latex_elements = { + 'preamble': r''' +\usepackage{amsmath} +''', +} \ No newline at end of file diff --git a/Documentation/reference/tools/algos/egm.rst b/Documentation/reference/tools/algos/egm.rst index f9daee533..1aedcfe39 100644 --- a/Documentation/reference/tools/algos/egm.rst +++ b/Documentation/reference/tools/algos/egm.rst @@ -1,5 +1,5 @@ algos.egm --------- +---------- .. automodule:: HARK.algos.egm :members: diff --git a/Documentation/reference/tools/algos/foc.rst b/Documentation/reference/tools/algos/foc.rst index cd80d3aee..612e4466f 100644 --- a/Documentation/reference/tools/algos/foc.rst +++ b/Documentation/reference/tools/algos/foc.rst @@ -1,5 +1,5 @@ algos.foc --------- +---------- .. automodule:: HARK.algos.foc :members: diff --git a/HARK/algos/egm.py b/HARK/algos/egm.py index 46d5deb81..3d224ff6f 100644 --- a/HARK/algos/egm.py +++ b/HARK/algos/egm.py @@ -1,7 +1,56 @@ -from dataclasses import field +""" +This algorithm solves for an agent's policy given: + +- a state space :math:`X` +- a shock space :math:`Z` +- an action space :math:`A` +- a post-state space :math:`Y` +- a deterministic transition function :math:`g` + +The agent chooses their action after shocks have been realized: + +.. math:: g: X \\times Z \\times A \\rightarrow Y + +- a reward function :math:`r` +- a post-value function :math:`v_y` +- constraints on the action :math:`\Gamma` +- a scalar discount factor :math:`\\beta` + +The problem it solves is of the form: + +.. math:: \pi^*(x, z) = \\textrm{argmax}_{a \in \Gamma(x, z)} r(x, z, a) + \\beta v_y(g(x, z, a))] + +However, this algorithm uses the Endogenous Gridpoints Method (EGM) [1]_, +which solves this problem only indirectly. + +It starts with a grid over post-states :math:`\\bar{Y}`. + +For each value of :math:`\hat{y}` in the grid, it analytically determines +the action which, when chosen as an optimal solution to the problem, +results in that post-state. + +.. math:: + + \pi_y(y) = \\frac{\partial r}{\partial a}^{-1}(\\beta \\frac{\partial g}{\partial a} \\frac{\partial v_y}{\partial y}(\hat{y}) + +It then computes the state that corresponds to that action using +the inverse transition function: + +.. math:: + + \hat{x} = g^{-1}(\hat{y}, \pi_y(\hat{y})) + +The pairs :math:`(\hat{x}, \pi_y(\hat{y}))` then are data points for +the grid for the optimal policy :math:`\pi^*`. + +.. [1] Carroll, C. D. (2006). The method of endogenous gridpoints for solving dynamic stochastic + optimization problems. Economics letters, 91(3), 312-320. + +""" + + import itertools import numpy as np -from scipy.optimize import minimize, brentq from typing import Callable, Mapping, Sequence, Tuple import xarray as xr @@ -37,8 +86,6 @@ def analytic_pi_y_star( ------- The actions chosen that result in post-state Y. - - TODO: Type signature of output """ if dg_da is None or not( isinstance(dg_da, float) or isinstance(dg_da, int) diff --git a/HARK/algos/foc.py b/HARK/algos/foc.py index 7b906daea..542615cb5 100644 --- a/HARK/algos/foc.py +++ b/HARK/algos/foc.py @@ -1,37 +1,70 @@ -from dataclasses import field -import itertools -import numpy as np -from scipy.optimize import minimize, brentq -from typing import Callable, Mapping, Sequence, Tuple -import xarray as xr - """ -Sargent + Stachurski: -x - states -z - shocks -a - actions +This algorithm solves for an agent's policy given: + +- a state space :math:`X` +- a shock space :math:`Z` +- an action space :math:`A` +- a post-state space :math:`Y` +- a deterministic transition function :math:`g` + +The agent chooses their action after shocks have been realized: + +.. math:: g: X \\times Z \\times A \\rightarrow Y + +- a reward function :math:`r` +- a post-value function :math:`v_y` +- constraints on the action :math:`\Gamma` +- a scalar discount factor :math:`\\beta` + +The problem it solves is of the form: + +.. math:: \pi^*(x, z) = \\textrm{argmax}_{a \in \Gamma(x, z)} r(x, z, a) + \\beta v_y(g(x, z, a))] + +It solves this on each point on grids over :math:`\\bar{X}` and :math:`\\bar{Z}`, + +It does this using the first order condition (FOC). -New: -y - post-states +Given a state-action value function: -French notation: -g - transtion function x,z, a -> y +.. math:: + q(x, z, a) = r(x, z, a) + \\beta v_y(g(x, z, a)) +And the derivative of this with respect to actions: -Question: - - Do we want to have a canonical 'grid' object that includes variable names and allows for non-cartesian grids? +.. math:: + \\frac{\partial q}{\partial a}(x, z, a) = \\frac{\partial r}{\partial a}(x, z, a) + \\beta \\frac{\partial v_y}{\partial y}(g(x, z, a) \\frac{\partial g}{\partial a}(x, z, a)) + +The first order condition is that + +.. math:: + + 0 = \\frac{\partial q}{\partial a}(x, z, a) """ +from dataclasses import field +import itertools +import numpy as np +from scipy.optimize import minimize, brentq +from typing import Callable, Mapping, Sequence, Tuple +import xarray as xr + ## TODO: Action handling def xndindex(ds, dims=None): """ + This function returns an iterator over the values in a DataArray or DataSet's coordinates. + There is currently no integrated way to iterate over an xarray.DataArray with its coordinate labels. This method is a workaround from: https://github.com/pydata/xarray/issues/2805#issuecomment-1255029201 + + Parameters + ------- + + ds: xarray.DataArray or xarray.DataSet """ if dims is None: dims = ds.dims @@ -55,11 +88,28 @@ def xz_grids_to_data_array( x_grid : Grid = {}, ## TODO: Better data structure here. z_grid : Grid = {} ): + """ + Construct a zero-valued DataArray with the coordinates + based on the two Grids passed in. + + Parameters + ---------- + x_grid: Grid + A mapping from state variable labels to a sequence of numerical values. + z_grid: Grid + A mapping from shock variable labels to a sequence of numerical values. + + Returns + -------- + + da xarray.DataArray + An xarray.DataArray with coordinates given by both grids. + """ coords = {**x_grid, **z_grid} da = xr.DataArray( - np.zeros([len(v) for v in coords.values()]), + np.empty([len(v) for v in coords.values()]), dims = coords.keys(), coords = coords ) @@ -72,12 +122,12 @@ def optimal_policy_foc( actions: Sequence[str], r : Callable[[Grid, Grid, Grid], float], dr_da, - dr_inv, + dr_inv, dg_dx = 1, dg_da = -1, x_grid : Grid = {}, z_grid : Grid = {}, - v_y_der : Callable[[Grid, Grid, Grid], float] = lambda x : 0, + v_y_der : Callable[[Grid], float] = lambda y : 0, discount = 1, action_upper_bound : Callable[[Grid, Grid], Sequence[float]] = field(default = None), action_lower_bound : Callable[[Grid, Grid], Sequence[float]] = field(default = None), @@ -89,74 +139,53 @@ def optimal_policy_foc( Uses root finding and the first order condition. - This is written with the expectation that: - - Actions are scalar - - the Q function is concave over the action space. + This is written with the expectation that actions are scalar + and the q function is concave over the action space. - Functionality is not guaranteed otherwise. - - - NOTE: This does not put pre-defined solution points into the solution data. - That needs to be added in a different step. - - - TODO: K -> Z - TODO: Transition function as an argument + The functionality is not guaranteed otherwise. Parameters ----------- - - g : + g: Transition functiong g. - actions: Ordered labels of action variables. - - r : + r: Reward function r. - - dr_da : + dr_da: Derivative of reward function r with respect to actions A. - - dr_inv : + dr_inv: Inverse of derivative of reward function - dg_dx: Derivative of transition function g with respect to states X. - dg_da: Derivative of transition function g with respect to actions A. - x_grid: - + Grid of points x in state space X. A mapping between variable labels + and sequences of values z_grid: - + Grid of points z in state space Z. A mapping between variable labels + and sequences of values v_y_der: Derivative of the post-value function v_y with respect to post-states Y. - discount: - + A scalar discount factor. action_upper_bound: - + Function of x and z, giving upper bound on action a. action_lower_bound: + Function of x and z, giving lower bound on action a. Returns -------- pi_star_data: - Data for the optimal policy pi_star. - + Data for the optimal policy :math:`\pi^*`. q_der_data: - Data for the derivative of the action value function function q. - - y_data: - + Data for the derivative of the action value function function :math:`\partial q/\partial a`. """ # Set up data arrays with coordinates based on the grid. pi_star_data = xz_grids_to_data_array(x_grid, z_grid) ## TODO: Default value for pi should be nan. q_der_data = xz_grids_to_data_array(x_grid, z_grid) - ## May need to expand this for multivariate y - y_data = xz_grids_to_data_array(x_grid, z_grid) def dq_da(x,z,a, v_y_der): """ @@ -236,7 +265,6 @@ def foc(a): print(f"Rootfinding failure at {x_vals}, {z_vals}.") print(root_res) - #raise Exception("Failed to optimize.") pi_star_data.sel(**x_vals, **z_vals).variable.data.put(0, root_res.root) q_der_xz = dq_da( @@ -248,9 +276,4 @@ def foc(a): q_der_data.sel(**x_vals, **z_vals).variable.data.put(0, q_der_xz) - acts = np.atleast_1d(pi_star_data.sel(**x_vals, **z_vals).values) - y = g(x_vals, z_vals, action_zip(acts)) - y_n = np.array([y[k] for k in y]) - y_data.sel(**x_vals, **z_vals).variable.data.put(0, y_n) - - return pi_star_data, q_der_data, y_data \ No newline at end of file + return pi_star_data, q_der_data \ No newline at end of file diff --git a/HARK/algos/tests/test_algos.py b/HARK/algos/tests/test_algos.py index 3007559e9..b636ab02a 100644 --- a/HARK/algos/tests/test_algos.py +++ b/HARK/algos/tests/test_algos.py @@ -117,7 +117,7 @@ def test_optimal_policy_foc(self): def consumption_v_y_der(y : Mapping[str,Any]): return gothic.VP_Tminus1(y['a']) - pi_star, q_der, y_data = optimal_policy_foc( + pi_star, q_der = optimal_policy_foc( g, ['c'], r, From 5774394f3ff088fd7fccd19dc141648601161287 Mon Sep 17 00:00:00 2001 From: sb Date: Mon, 26 Jun 2023 11:56:32 -0400 Subject: [PATCH 14/14] CHANGELOG update --- Documentation/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/Documentation/CHANGELOG.md b/Documentation/CHANGELOG.md index e4e57967c..2bb7b2155 100644 --- a/Documentation/CHANGELOG.md +++ b/Documentation/CHANGELOG.md @@ -14,6 +14,7 @@ Release Date: TBD ### Major Changes +- Adds `HARK.algos.foc` and `HARK.algos.egm` modules for general algorithms for solving the optimization step of an agent's problem. Tests validate behavior against gold-standard implemented, contained in `HARK.gothic`. [#1283](https://github.com/econ-ark/HARK/pull/1283) - Adds `HARK.core.AgentPopulation` class to represent a population of agents with ex-ante heterogeneous parametrizations as distributions. [#1237](https://github.com/econ-ark/HARK/pull/1237) ### Minor Changes