Merge branch 'master' into master

Documentation/CHANGELOG.md

@@ -17,6 +17,7 @@ Release Date: TBD
 - Adds a discretize method to DBlocks and RBlocks (#1460)[https://github.com/econ-ark/HARK/pull/1460]
 - Allows structural equations in model files to be provided in string form [#1427](https://github.com/econ-ark/HARK/pull/1427)
 - Introduces `HARK.parser' module for parsing configuration files into models [#1427](https://github.com/econ-ark/HARK/pull/1427)
+- Allows construction of shocks with arguments based on mathematical expressions [#1464](https://github.com/econ-ark/HARK/pull/1464)
 
 #### Minor Changes
 

HARK/distribution.py

@@ -311,6 +311,8 @@ def __new__(
         mu: Union[float, np.ndarray] = 0.0,
         sigma: Union[float, np.ndarray] = 1.0,
         seed: Optional[int] = 0,
+        mean=None,
+        std=None,
     ):
         """
         Create a new Lognormal distribution. If sigma is zero, return a
@@ -324,6 +326,10 @@ def __new__(
             Standard deviation of underlying normal distribution, by default 1.0
         seed : Optional[int], optional
             Seed for random number generator, by default None
+        mean: optional
+            For alternative mean/std parameterization, the mean of the lognormal distribution
+        std: optional
+            For alternative mean/std parameterization, the standard deviation of the lognormal distribution
 
         Returns
         -------
@@ -342,7 +348,15 @@ def __init__(
         mu: Union[float, np.ndarray] = 0.0,
         sigma: Union[float, np.ndarray] = 1.0,
         seed: Optional[int] = 0,
+        mean=None,
+        std=None,
     ):
+        if mean is not None and sigma is not None:
+            mean_squared = mean**2
+            variance = std**2
+            mu = np.log(mean / (np.sqrt(1.0 + variance / mean_squared)))
+            sigma = np.sqrt(np.log(1.0 + variance / mean_squared))
+
         self.mu = np.asarray(mu)
         self.sigma = np.asarray(sigma)
 

HARK/model.py

@@ -78,6 +78,57 @@ def discretized_shock_dstn(shocks, disc_params):
     return all_shock_dstn
 
 
+def construct_shocks(shock_data, scope):
+    """
+    Returns a dictionary from shock labels to Distributions.
+
+    When the corresponding value in shock_data contains
+    a distribution constructor and input information,
+    any symbolic expressions used in the inputs are
+    evaluated in the provided scope.
+
+    Parameters
+    ------------
+
+     shock_data: Mapping(str, Distribution or tuple)
+        A mapping from variable names to Distribution objects,
+        representing exogenous shocks.
+
+        Optionally, the mapping can be to tuples of Distribution
+        constructors and dictionary of input arguments.
+        In this case, the dictionary can map argument names to
+        numbers, or to strings. The strings are parsed as
+        mathematical expressions and evaluated in the scope
+        of a calibration dictionary.
+
+    scope: dict(str, values)
+        Variables assigned to numerical values.
+        The scope in which expressions will be evaluated
+    """
+    sd = deepcopy(shock_data)
+
+    for v in sd:
+        if isinstance(sd[v], tuple):
+            dist_class = sd[v][0]
+
+            dist_args = sd[v][1]  # should be a dictionary
+
+            for a in dist_args:
+                if isinstance(dist_args[a], str):
+                    arg_lambda = math_text_to_lambda(dist_args[a])
+                    arg_value = arg_lambda(
+                        *[scope[var] for var in signature(arg_lambda).parameters]
+                    )
+
+                    dist_args[a] = arg_value
+
+            dist = dist_class(**dist_args)
+
+            sd[v] = dist
+
+    return sd
+
+
 def simulate_dynamics(
     dynamics: Mapping[str, Union[Callable, Control]],
     pre: Mapping[str, Any],
@@ -149,10 +200,17 @@ class DBlock(Block):
 
     Parameters
     ----------
-    shocks: Mapping(str, Distribution)
+    shocks: Mapping(str, Distribution or tuple)
         A mapping from variable names to Distribution objects,
         representing exogenous shocks.
 
+        Optionally, the mapping can be to tuples of Distribution
+        constructors and dictionary of input arguments.
+        In this case, the dictionary can map argument names to
+        numbers, or to strings. The strings are parsed as
+        mathematical expressions and evaluated in the scope
+        of a calibration dictionary.
+
     dynamics: Mapping(str, str or callable)
         A dictionary mapping variable names to mathematical expressions.
         These expressions can be simple functions, in which case the
@@ -167,6 +225,13 @@ class DBlock(Block):
     dynamics: dict = field(default_factory=dict)
     reward: dict = field(default_factory=dict)
 
+    def construct_shocks(self, calibration):
+        """
+        Constructs all shocks given calibration.
+        This method mutates the DBlock.
+        """
+        self.shocks = construct_shocks(self.shocks, calibration)
+
     def discretize(self, disc_params):
         """
         Returns a new DBlock which is a copy of this one, but with shock discretized.
@@ -299,6 +364,13 @@ class RBlock(Block):
     description: str = ""
     blocks: List[Block] = field(default_factory=list)
 
+    def construct_shocks(self, calibration):
+        """
+        Construct all shocks given a calibration dictionary.
+        """
+        for b in self.blocks:
+            b.construct_shocks(calibration)
+
     def discretize(self, disc_params):
         """
         Recursively discretizes all the blocks.

HARK/models/consumer.py

@@ -7,28 +7,25 @@
 Microeconomic Dynamic Stochastic Optimization Problems"
 """
 
-# TODO: Include these in calibration, then construct shocks
-LivPrb = 0.98
-TranShkStd = 0.1
-RiskyStd = 0.1
-
 calibration = {
     "DiscFac": 0.96,
     "CRRA": 2.0,
     "R": 1.03,  # note: this can be overriden by the portfolio dynamics
     "Rfree": 1.03,
     "EqP": 0.02,
-    "LivPrb": LivPrb,
+    "LivPrb": 0.98,
     "PermGroFac": 1.01,
     "BoroCnstArt": None,
+    "TranShkStd": 0.1,
+    "RiskyStd": 0.1,
 }
 
 consumption_block = DBlock(
     **{
         "name": "consumption",
         "shocks": {
-            "live": Bernoulli(p=LivPrb),  # Move to tick or mortality block?
-            "theta": MeanOneLogNormal(sigma=TranShkStd),
+            "live": (Bernoulli, {"p": "LivPrb"}),  # Move to tick or mortality block?
+            "theta": (MeanOneLogNormal, {"sigma": "TranShkStd"}),
         },
         "dynamics": {
             "b": lambda k, R: k * R,
@@ -46,8 +43,8 @@
     **{
         "name": "consumption normalized",
         "shocks": {
-            "live": Bernoulli(p=LivPrb),  # Move to tick or mortality block?
-            "theta": MeanOneLogNormal(sigma=TranShkStd),
+            "live": (Bernoulli, {"p": "LivPrb"}),  # Move to tick or mortality block?
+            "theta": (MeanOneLogNormal, {"sigma": "TranShkStd"}),
         },
         "dynamics": {
             "b": lambda k, R, PermGroFac: k * R / PermGroFac,
@@ -63,9 +60,7 @@
     **{
         "name": "portfolio",
         "shocks": {
-            "risky_return": Lognormal.from_mean_std(
-                calibration["Rfree"] + calibration["EqP"], RiskyStd
-            )
+            "risky_return": (Lognormal, {"mean": "Rfree + EqP", "std": "RiskyStd"})
         },
         "dynamics": {
             "stigma": Control(["a"]),

HARK/simulation/monte_carlo.py

@@ -14,7 +14,7 @@
 )
 from HARK.model import Aggregate
 from HARK.model import DBlock
-from HARK.model import simulate_dynamics
+from HARK.model import construct_shocks, simulate_dynamics
 
 
 def draw_shocks(shocks: Mapping[str, Distribution], conditions: Sequence[int]):
@@ -139,7 +139,11 @@ def __init__(
 
         self.calibration = calibration
         self.block = block
-        self.shocks = block.get_shocks()
+
+        # shocks are exogenous (but for age) but can depend on calibration
+        raw_shocks = block.get_shocks()
+        self.shocks = construct_shocks(raw_shocks, calibration)
+
         self.dynamics = block.get_dynamics()
         self.dr = dr
         self.initial = initial

HARK/tests/test_model.py

@@ -34,6 +34,7 @@ class test_DBlock(unittest.TestCase):
     def setUp(self):
         self.test_block_A = model.DBlock(**test_block_A_data)
         self.cblock = cons.consumption_block_normalized
+        self.cblock.construct_shocks(cons.calibration)
 
         # prior states relative to the decision, so with realized shocks.
         self.dpre = {"k": 2, "R": 1.05, "PermGroFac": 1.1, "theta": 1, "CRRA": 2}
@@ -98,6 +99,7 @@ def test_init(self):
         self.assertEqual(len(r_block_tree.get_shocks()), 3)
 
     def test_discretize(self):
+        self.cpp.construct_shocks(cons.calibration)
         cppd = self.cpp.discretize({"theta": {"N": 5}, "risky_return": {"N": 6}})
 
         self.assertEqual(len(cppd.get_shocks()["theta"].pmv), 5)