From 54db00b0a347cf8129e5bd4f9d2823f6471eb502 Mon Sep 17 00:00:00 2001 From: sb Date: Fri, 24 Apr 2020 12:14:22 -0400 Subject: [PATCH 1/3] from continuous to discretized distributions for Normal and Uniform --- HARK/distribution.py | 70 +++++++++++++++++----------------- HARK/tests/test_approxDstns.py | 2 +- 2 files changed, 35 insertions(+), 37 deletions(-) diff --git a/HARK/distribution.py b/HARK/distribution.py index 4f0a431ad..ca722bbc9 100644 --- a/HARK/distribution.py +++ b/HARK/distribution.py @@ -92,8 +92,8 @@ def __init__(self, mu = 0.0, sigma = 1.0): One or more standard deviations. Number of elements T in sigma determines number of rows of output. ''' - self.mu = 0.0 - self.sigma = 1.0 + self.mu = mu + self.sigma = sigma def draw(self, N, seed=0): ''' @@ -127,6 +127,17 @@ def draw(self, N, seed=0): draws.append(self.sigma[t]*RNG.randn(N) + self.mu[t]) return draws + def approx(self, N): + """ + Returns a discrete approximation of this distribution. + """ + x, w = np.polynomial.hermite.hermgauss(N) + # normalize w + pmf = w*np.pi**-0.5 + # correct x + X = math.sqrt(2.0)*self.sigma*x + self.mu + return DiscreteDistribution(pmf, X) + class Weibull(): ''' A Weibull distribution @@ -245,6 +256,26 @@ def draw(self, N, seed=0): draws.append(self.bot[t] + (self.top[t] - self.bot[t])*RNG.rand(N)) return draws + def approx(self, N): + ''' + Makes a discrete approximation to this uniform distribution. + + Parameters + ---------- + N : int + The number of points in the discrete approximation + + Returns + ------- + d : DiscreteDistribution + Probability associated with each point in array of discrete + points for discrete probability mass function. + ''' + pmf = np.ones(N)/float(N) + center = (self.top+self.bot)/2.0 + width = (self.top-self.bot)/2.0 + X = center + width*np.linspace(-(N-1.0)/2.0,(N-1.0)/2.0,N)/(N/2.0) + return DiscreteDistribution(pmf,X) def drawMeanOneLognormal(N, sigma=1.0, seed=0): ''' @@ -537,16 +568,9 @@ def approxMeanOneLognormal(N, sigma=1.0, **kwargs): dist = approxLognormal(N=N, mu=mu_adj, sigma=sigma, **kwargs) return dist -def approxNormal(N, mu=0.0, sigma=1.0): - x, w = np.polynomial.hermite.hermgauss(N) - # normalize w - pmf = w*np.pi**-0.5 - # correct x - X = math.sqrt(2.0)*sigma*x + mu - return DiscreteDistribution(pmf, X) def approxLognormalGaussHermite(N, mu=0.0, sigma=1.0): - d = approxNormal(N, mu, sigma) + d = Normal(mu, sigma).approx(N) return DiscreteDistribution(d.pmf, np.exp(d.X)) @@ -591,32 +615,6 @@ def approxBeta(N,a=1.0,b=1.0): pmf = np.ones(N)/float(N) return DiscreteDistribution(pmf, X) -def approxUniform(N,bot=0.0,top=1.0): - ''' - Makes a discrete approximation to a uniform distribution, given its bottom - and top limits and number of points. - - Parameters - ---------- - N : int - The number of points in the discrete approximation - bot : float - The bottom of the uniform distribution - top : float - The top of the uniform distribution - - Returns - ------- - d : DiscreteDistribution - Probability associated with each point in array of discrete - points for discrete probability mass function. - ''' - pmf = np.ones(N)/float(N) - center = (top+bot)/2.0 - width = (top-bot)/2.0 - X = center + width*np.linspace(-(N-1.0)/2.0,(N-1.0)/2.0,N)/(N/2.0) - return DiscreteDistribution(pmf,X) - def makeMarkovApproxToNormal(x_grid,mu,sigma,K=351,bound=3.5): ''' diff --git a/HARK/tests/test_approxDstns.py b/HARK/tests/test_approxDstns.py index 0a24f6326..fbe1adf98 100644 --- a/HARK/tests/test_approxDstns.py +++ b/HARK/tests/test_approxDstns.py @@ -16,7 +16,7 @@ def setUp(self): def test_mu_normal(self): for muNormal in self.muNormals: for stdNormal in self.stdNormals: - d = distribution.approxNormal(40, muNormal) + d = distribution.Normal(muNormal).approx(40) self.assertTrue(sum(d.pmf*d.X)-muNormal<1e-12) def test_mu_lognormal_from_normal(self): From d43c35eb89e3792b1222f53902ffcd6e29bd9a49 Mon Sep 17 00:00:00 2001 From: sb Date: Fri, 24 Apr 2020 13:01:30 -0400 Subject: [PATCH 2/3] approx_ stand-alone functions now approx() methods on continuous distributions, fixes #647 --- HARK/#core.py# | 1409 +++++++++++++++++ HARK/#scratch# | 0 HARK/ConsumptionSaving/ConsAggShockModel.py | 22 +- HARK/ConsumptionSaving/ConsIndShockModel.py | 28 +- .../ConsumptionSaving/ConsMarkovModel.py.orig | 1023 ++++++++++++ HARK/ConsumptionSaving/ConsPortfolioModel.py | 6 +- HARK/ConsumptionSaving/ConsPrefShockModel.py | 14 +- HARK/distribution.py | 218 +-- HARK/lifecycle.py | 29 + HARK/policy.py | 11 + HARK/solver.py | 10 + HARK/tests/test_distribution.py | 2 +- 12 files changed, 2614 insertions(+), 158 deletions(-) create mode 100644 HARK/#core.py# create mode 100644 HARK/#scratch# create mode 100644 HARK/ConsumptionSaving/ConsMarkovModel.py.orig create mode 100644 HARK/lifecycle.py create mode 100644 HARK/policy.py create mode 100644 HARK/solver.py diff --git a/HARK/#core.py# b/HARK/#core.py# new file mode 100644 index 000000000..b4b8c2a03 --- /dev/null +++ b/HARK/#core.py# @@ -0,0 +1,1409 @@ +''' +High-level functions and classes for solving a wide variety of economic models. +The "core" of HARK is a framework for "microeconomic" and "macroeconomic" +models. A micro model concerns the dynamic optimization problem for some type +of agents, where agents take the inputs to their problem as exogenous. A macro +model adds an additional layer, endogenizing some of the inputs to the micro +problem by finding a general equilibrium dynamic rule. +''' +from __future__ import print_function, division +from __future__ import absolute_import + +from builtins import str +from builtins import range +from builtins import object +import sys +import os +from distutils.dir_util import copy_tree +from .utilities import getArgNames, NullFunc +from copy import copy, deepcopy +import numpy as np +from time import time +from .parallel import multiThreadCommands, multiThreadCommandsFake + + +def distanceMetric(thing_A, thing_B): + ''' + A "universal distance" metric that can be used as a default in many settings. + + Parameters + ---------- + thing_A : object + A generic object. + thing_B : object + Another generic object. + + Returns: + ------------ + distance : float + The "distance" between thing_A and thing_B. + ''' + # Get the types of the two inputs + typeA = type(thing_A) + typeB = type(thing_B) + + if typeA is list and typeB is list: + lenA = len(thing_A) # If both inputs are lists, then the distance between + lenB = len(thing_B) # them is the maximum distance between corresponding + if lenA == lenB: # elements in the lists. If they differ in length, + distance_temp = [] # the distance is the difference in lengths. + for n in range(lenA): + distance_temp.append(distanceMetric(thing_A[n], thing_B[n])) + distance = max(distance_temp) + else: + distance = float(abs(lenA - lenB)) + # If both inputs are numbers, return their difference + elif isinstance(thing_A, (int, float)) and isinstance(thing_B, (int, float)): + distance = float(abs(thing_A - thing_B)) + # If both inputs are array-like, return the maximum absolute difference b/w + # corresponding elements (if same shape); return largest difference in dimensions + # if shapes do not align. + elif hasattr(thing_A, 'shape') and hasattr(thing_B, 'shape'): + if thing_A.shape == thing_B.shape: + distance = np.max(abs(thing_A - thing_B)) + else: + # Flatten arrays so they have the same dimensions + distance = np.max(abs(thing_A.flatten().shape[0] - thing_B.flatten().shape[0])) + # If none of the above cases, but the objects are of the same class, call + # the distance method of one on the other + elif thing_A.__class__.__name__ == thing_B.__class__.__name__: + if thing_A.__class__.__name__ == 'function': + distance = 0.0 + else: + distance = thing_A.distance(thing_B) + else: # Failsafe: the inputs are very far apart + distance = 1000.0 + return distance + + +class HARKobject(object): + ''' + A superclass for object classes in HARK. Comes with two useful methods: + a generic/universal distance method and an attribute assignment method. + ''' + def distance(self, other): + ''' + A generic distance method, which requires the existence of an attribute + called distance_criteria, giving a list of strings naming the attributes + to be considered by the distance metric. + + Parameters + ---------- + other : object + Another object to compare this instance to. + + Returns + ------- + (unnamed) : float + The distance between this object and another, using the "universal + distance" metric. + ''' + distance_list = [0.0] + for attr_name in self.distance_criteria: + try: + obj_A = getattr(self, attr_name) + obj_B = getattr(other, attr_name) + distance_list.append(distanceMetric(obj_A, obj_B)) + except AttributeError: + distance_list.append(1000.0) # if either object lacks attribute, they are not the same + return max(distance_list) + + def assignParameters(self, **kwds): + ''' + Assign an arbitrary number of attributes to this agent. + + Parameters + ---------- + **kwds : keyword arguments + Any number of keyword arguments of the form key=value. Each value + will be assigned to the attribute named in self. + + Returns + ------- + none + ''' + for key in kwds: + setattr(self, key, kwds[key]) + + def __call__(self, **kwds): + ''' + Assign an arbitrary number of attributes to this agent, as a convenience. + See assignParameters. + ''' + self.assignParameters(**kwds) + + def getAvg(self, varname, **kwds): + ''' + Calculates the average of an attribute of this instance. Returns NaN if no such attribute. + + Parameters + ---------- + varname : string + The name of the attribute whose average is to be calculated. This attribute must be an + np.array or other class compatible with np.mean. + + Returns + ------- + avg : float or np.array + The average of this attribute. Might be an array if the axis keyword is passed. + ''' + if hasattr(self, varname): + return np.mean(getattr(self, varname), **kwds) + else: + return np.nan + + +class Solution(HARKobject): + ''' + A superclass for representing the "solution" to a single period problem in a + dynamic microeconomic model. + + NOTE: This can be deprecated now that HARKobject exists, but this requires + replacing each instance of Solution with HARKobject in the other modules. + ''' + + +class AgentType(HARKobject): + ''' + A superclass for economic agents in the HARK framework. Each model should + specify its own subclass of AgentType, inheriting its methods and overwriting + as necessary. Critically, every subclass of AgentType should define class- + specific static values of the attributes time_vary and time_inv as lists of + strings. Each element of time_vary is the name of a field in AgentSubType + that varies over time in the model. Each element of time_inv is the name of + a field in AgentSubType that is constant over time in the model. + ''' + def __init__(self, solution_terminal=None, cycles=1, time_flow=True, pseudo_terminal=True, + tolerance=0.000001, seed=0, **kwds): + ''' + Initialize an instance of AgentType by setting attributes. + + Parameters + ---------- + solution_terminal : Solution + A representation of the solution to the terminal period problem of + this AgentType instance, or an initial guess of the solution if this + is an infinite horizon problem. + cycles : int + The number of times the sequence of periods is experienced by this + AgentType in their "lifetime". cycles=1 corresponds to a lifecycle + model, with a certain sequence of one period problems experienced + once before terminating. cycles=0 corresponds to an infinite horizon + model, with a sequence of one period problems repeating indefinitely. + time_flow : boolean + Whether time is currently "flowing" forward(True) or backward(False) for this + instance. Used to flip between solving (using backward iteration) + and simulating (etc). + pseudo_terminal : boolean + Indicates whether solution_terminal isn't actually part of the + solution to the problem (as a known solution to the terminal period + problem), but instead represents a "scrap value"-style termination. + When True, solution_terminal is not included in the solution; when + false, solution_terminal is the last element of the solution. + tolerance : float + Maximum acceptable "distance" between successive solutions to the + one period problem in an infinite horizon (cycles=0) model in order + for the solution to be considered as having "converged". Inoperative + when cycles>0. + seed : int + A seed for this instance's random number generator. + + Returns + ------- + None + ''' + if solution_terminal is None: + solution_terminal = NullFunc() + self.solution_terminal = solution_terminal # NOQA + self.cycles = cycles # NOQA + self.time_flow = time_flow # NOQA + self.pseudo_terminal = pseudo_terminal # NOQA + self.solveOnePeriod = NullFunc() # NOQA + self.tolerance = tolerance # NOQA + self.seed = seed # NOQA + self.track_vars = [] # NOQA + self.poststate_vars = [] # NOQA + self.read_shocks = False # NOQA + self.assignParameters(**kwds) # NOQA + self.resetRNG() # NOQA + + def timeReport(self): + ''' + Report to the user the direction that time is currently "flowing" for + this instance. Only exists as a reminder of how time_flow works. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + if self.time_flow: + print('Time varying objects are listed in ordinary chronological order.') + else: + print('Time varying objects are listed in reverse chronological order.') + + def timeFlip(self): + ''' + Reverse the flow of time for this instance. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + for name in self.time_vary: + self.__dict__[name].reverse() + self.time_flow = not self.time_flow + + def timeFwd(self): + ''' + Make time flow forward for this instance. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + if not self.time_flow: + self.timeFlip() + + def timeRev(self): + ''' + Make time flow backward for this instance. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + if self.time_flow: + self.timeFlip() + + def addToTimeVary(self, *params): + ''' + Adds any number of parameters to time_vary for this instance. + + Parameters + ---------- + params : string + Any number of strings naming attributes to be added to time_vary + + Returns + ------- + None + ''' + for param in params: + if param not in self.time_vary: + self.time_vary.append(param) + + def addToTimeInv(self, *params): + ''' + Adds any number of parameters to time_inv for this instance. + + Parameters + ---------- + params : string + Any number of strings naming attributes to be added to time_inv + + Returns + ------- + None + ''' + for param in params: + if param not in self.time_inv: + self.time_inv.append(param) + + def delFromTimeVary(self, *params): + ''' + Removes any number of parameters from time_vary for this instance. + + Parameters + ---------- + params : string + Any number of strings naming attributes to be removed from time_vary + + Returns + ------- + None + ''' + for param in params: + if param in self.time_vary: + self.time_vary.remove(param) + + def delFromTimeInv(self, *params): + ''' + Removes any number of parameters from time_inv for this instance. + + Parameters + ---------- + params : string + Any number of strings naming attributes to be removed from time_inv + + Returns + ------- + None + ''' + for param in params: + if param in self.time_inv: + self.time_inv.remove(param) + + def solve(self, verbose=False): + ''' + Solve the model for this instance of an agent type by backward induction. + Loops through the sequence of one period problems, passing the solution + from period t+1 to the problem for period t. + + Parameters + ---------- + verbose : boolean + If True, solution progress is printed to screen. + + Returns + ------- + none + ''' + + # Ignore floating point "errors". Numpy calls it "errors", but really it's excep- + # tions with well-defined answers such as 1.0/0.0 that is np.inf, -1.0/0.0 that is + # -np.inf, np.inf/np.inf is np.nan and so on. + with np.errstate(divide='ignore', over='ignore', under='ignore', invalid='ignore'): + self.preSolve() # Do pre-solution stuff + self.solution = solveAgent(self, verbose) # Solve the model by backward induction + if self.time_flow: # Put the solution in chronological order if this instance's time flow runs that way + self.solution.reverse() + self.addToTimeVary('solution') # Add solution to the list of time-varying attributes + self.postSolve() # Do post-solution stuff + + def resetRNG(self): + ''' + Reset the random number generator for this type. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + self.RNG = np.random.RandomState(self.seed) + + def checkElementsOfTimeVaryAreLists(self): + """ + A method to check that elements of time_vary are lists. + """ + for param in self.time_vary: + assert type(getattr(self, param)) == list, param + ' is not a list, but should be' + \ + ' because it is in time_vary' + + def checkRestrictions(self): + """ + A method to check that various restrictions are met for the model class. + """ + return + + def preSolve(self): + ''' + A method that is run immediately before the model is solved, to check inputs or to prepare + the terminal solution, perhaps. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + self.checkRestrictions() + self.checkElementsOfTimeVaryAreLists() + return None + + def postSolve(self): + ''' + A method that is run immediately after the model is solved, to finalize + the solution in some way. Does nothing here. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + return None + + def initializeSim(self): + ''' + Prepares this AgentType for a new simulation. Resets the internal random number generator, + makes initial states for all agents (using simBirth), clears histories of tracked variables. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + if not hasattr(self, 'T_sim'): + raise Exception('To initialize simulation variables it is necessary to first ' + + 'set the attribute T_sim to the largest number of observations ' + + 'you plan to simulate for each agent including re-births.') + elif self.T_sim <= 0: + raise Exception('T_sim represents the largest number of observations ' + + 'that can be simulated for an agent, and must be a positive number.') + + self.resetRNG() + self.t_sim = 0 + all_agents = np.ones(self.AgentCount, dtype=bool) + blank_array = np.zeros(self.AgentCount) + for var_name in self.poststate_vars: + setattr(self, var_name, copy(blank_array)) + self.t_age = np.zeros(self.AgentCount, dtype=int) # Number of periods since agent entry + self.t_cycle = np.zeros(self.AgentCount, dtype=int) # Which cycle period each agent is on + self.simBirth(all_agents) + self.clearHistory() + return None + + def simOnePeriod(self): + ''' + Simulates one period for this type. Calls the methods getMortality(), getShocks() or + readShocks, getStates(), getControls(), and getPostStates(). These should be defined for + AgentType subclasses, except getMortality (define its components simDeath and simBirth + instead) and readShocks. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + if not hasattr(self, 'solution'): + raise Exception('Model instance does not have a solution stored. To simulate, it is necessary' + ' to run the `solve()` method of the class first.') + + self.getMortality() # Replace some agents with "newborns" + if self.read_shocks: # If shock histories have been pre-specified, use those + self.readShocks() + else: # Otherwise, draw shocks as usual according to subclass-specific method + self.getShocks() + self.getStates() # Determine each agent's state at decision time + self.getControls() # Determine each agent's choice or control variables based on states + self.getPostStates() # Determine each agent's post-decision / end-of-period states using states and controls + + # Advance time for all agents + self.t_age = self.t_age + 1 # Age all consumers by one period + self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle + self.t_cycle[self.t_cycle == self.T_cycle] = 0 # Resetting to zero for those who have reached the end + + def makeShockHistory(self): + ''' + Makes a pre-specified history of shocks for the simulation. Shock variables should be named + in self.shock_vars, a list of strings that is subclass-specific. This method runs a subset + of the standard simulation loop by simulating only mortality and shocks; each variable named + in shock_vars is stored in a T_sim x AgentCount array in an attribute of self named X_hist. + Automatically sets self.read_shocks to True so that these pre-specified shocks are used for + all subsequent calls to simulate(). + + Parameters + ---------- + None + + Returns + ------- + None + ''' + # Make sure time is flowing forward and re-initialize the simulation + orig_time = self.time_flow + self.timeFwd() + self.initializeSim() + + # Make blank history arrays for each shock variable (and mortality) + for var_name in self.shock_vars: + setattr(self, var_name+'_hist', np.zeros((self.T_sim, self.AgentCount)) + np.nan) + setattr(self, 'who_dies_hist', np.zeros((self.T_sim, self.AgentCount), dtype=bool)) + + # Make and store the history of shocks for each period + for t in range(self.T_sim): + self.getMortality() + self.who_dies_hist[t,:] = self.who_dies + self.getShocks() + for var_name in self.shock_vars: + getattr(self, var_name + '_hist')[self.t_sim,:] = getattr(self, var_name) + self.t_sim += 1 + self.t_age = self.t_age + 1 # Age all consumers by one period + self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle + self.t_cycle[self.t_cycle == self.T_cycle] = 0 # Resetting to zero for those who have reached the end + + # Restore the flow of time and flag that shocks can be read rather than simulated + self.read_shocks = True + if not orig_time: + self.timeRev() + + def getMortality(self): + ''' + Simulates mortality or agent turnover according to some model-specific rules named simDeath + and simBirth (methods of an AgentType subclass). simDeath takes no arguments and returns + a Boolean array of size AgentCount, indicating which agents of this type have "died" and + must be replaced. simBirth takes such a Boolean array as an argument and generates initial + post-decision states for those agent indices. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + if self.read_shocks: + who_dies = self.who_dies_hist[self.t_sim,:] + else: + who_dies = self.simDeath() + self.simBirth(who_dies) + self.who_dies = who_dies + return None + + def simDeath(self): + ''' + Determines which agents in the current population "die" or should be replaced. Takes no + inputs, returns a Boolean array of size self.AgentCount, which has True for agents who die + and False for those that survive. Returns all False by default, must be overwritten by a + subclass to have replacement events. + + Parameters + ---------- + None + + Returns + ------- + who_dies : np.array + Boolean array of size self.AgentCount indicating which agents die and are replaced. + ''' + who_dies = np.zeros(self.AgentCount, dtype=bool) + return who_dies + + def simBirth(self, which_agents): + ''' + Makes new agents for the simulation. Takes a boolean array as an input, indicating which + agent indices are to be "born". Does nothing by default, must be overwritten by a subclass. + + Parameters + ---------- + which_agents : np.array(Bool) + Boolean array of size self.AgentCount indicating which agents should be "born". + + Returns + ------- + None + ''' + print('AgentType subclass must define method simBirth!') + return None + + def getShocks(self): + ''' + Gets values of shock variables for the current period. Does nothing by default, but can + be overwritten by subclasses of AgentType. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + return None + + def readShocks(self): + ''' + Reads values of shock variables for the current period from history arrays. For each var- + iable X named in self.shock_vars, this attribute of self is set to self.X_hist[self.t_sim,:]. + + This method is only ever called if self.read_shocks is True. This can be achieved by using + the method makeShockHistory() (or manually after storing a "handcrafted" shock history). + + Parameters + ---------- + None + + Returns + ------- + None + ''' + for var_name in self.shock_vars: + setattr(self, var_name, getattr(self, var_name + '_hist')[self.t_sim, :]) + + def getStates(self): + ''' + Gets values of state variables for the current period, probably by using post-decision states + from last period, current period shocks, and maybe market-level events. Does nothing by + default, but can be overwritten by subclasses of AgentType. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + return None + + def getControls(self): + ''' + Gets values of control variables for the current period, probably by using current states. + Does nothing by default, but can be overwritten by subclasses of AgentType. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + return None + + def getPostStates(self): + ''' + Gets values of post-decision state variables for the current period, probably by current + states and controls and maybe market-level events or shock variables. Does nothing by + default, but can be overwritten by subclasses of AgentType. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + return None + + def simulate(self, sim_periods=None): + ''' + Simulates this agent type for a given number of periods. Defaults to self.T_sim if no input. + Records histories of attributes named in self.track_vars in attributes named varname_hist. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + if not hasattr(self, 't_sim'): + raise Exception('It seems that the simulation variables were not initialize before calling ' + + 'simulate(). Call initializeSim() to initialize the variables before calling simulate() again.') + + if not hasattr(self, 'T_sim'): + raise Exception('This agent type instance must have the attribute T_sim set to a positive integer.' + + 'Set T_sim to match the largest dataset you might simulate, and run this agent\'s' + + 'initalizeSim() method before running simulate() again.') + + if sim_periods is not None and self.T_sim < sim_periods: + raise Exception('To simulate, sim_periods has to be larger than the maximum data set size ' + + 'T_sim. Either increase the attribute T_sim of this agent type instance ' + + 'and call the initializeSim() method again, or set sim_periods <= T_sim.') + + + # Ignore floating point "errors". Numpy calls it "errors", but really it's excep- + # tions with well-defined answers such as 1.0/0.0 that is np.inf, -1.0/0.0 that is + # -np.inf, np.inf/np.inf is np.nan and so on. + with np.errstate(divide='ignore', over='ignore', under='ignore', invalid='ignore'): + orig_time = self.time_flow + self.timeFwd() + if sim_periods is None: + sim_periods = self.T_sim + + for t in range(sim_periods): + self.simOnePeriod() + for var_name in self.track_vars: + getattr(self, var_name + '_hist')[self.t_sim,:] = getattr(self,var_name) + self.t_sim += 1 + + if not orig_time: + self.timeRev() + + def clearHistory(self): + ''' + Clears the histories of the attributes named in self.track_vars. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + for var_name in self.track_vars: + setattr(self, var_name + '_hist', np.zeros((self.T_sim,self.AgentCount)) + np.nan) + + +def solveAgent(agent, verbose): + ''' + Solve the dynamic model for one agent type. This function iterates on "cycles" + of an agent's model either a given number of times or until solution convergence + if an infinite horizon model is used (with agent.cycles = 0). + + Parameters + ---------- + agent : AgentType + The microeconomic AgentType whose dynamic problem is to be solved. + verbose : boolean + If True, solution progress is printed to screen (when cycles != 1). + + Returns + ------- + solution : [Solution] + A list of solutions to the one period problems that the agent will + encounter in his "lifetime". Returns in reverse chronological order. + ''' + # Record the flow of time when the Agent began the process, and make sure time is flowing backwards + original_time_flow = agent.time_flow + agent.timeRev() + + # Check to see whether this is an (in)finite horizon problem + cycles_left = agent.cycles # NOQA + infinite_horizon = cycles_left == 0 # NOQA + # Initialize the solution, which includes the terminal solution if it's not a pseudo-terminal period + solution = [] + if not agent.pseudo_terminal: + solution.append(deepcopy(agent.solution_terminal)) + + + # Initialize the process, then loop over cycles + solution_last = agent.solution_terminal # NOQA + go = True # NOQA + completed_cycles = 0 # NOQA + max_cycles = 5000 # NOQA - escape clause + if verbose: + t_last = time() + while go: + # Solve a cycle of the model, recording it if horizon is finite + solution_cycle = solveOneCycle(agent, solution_last) + if not infinite_horizon: + solution += solution_cycle + + # Check for termination: identical solutions across cycle iterations or run out of cycles + solution_now = solution_cycle[-1] + if infinite_horizon: + if completed_cycles > 0: + solution_distance = solution_now.distance(solution_last) + agent.solution_distance = solution_distance # Add these attributes so users can + agent.completed_cycles = completed_cycles # query them to see if solution is ready + go = (solution_distance > agent.tolerance and completed_cycles < max_cycles) + else: # Assume solution does not converge after only one cycle + solution_distance = 100.0 + go = True + else: + cycles_left += -1 + go = cycles_left > 0 + + # Update the "last period solution" + solution_last = solution_now + completed_cycles += 1 + + # Display progress if requested + if verbose: + t_now = time() + if infinite_horizon: + print('Finished cycle #' + str(completed_cycles) + ' in ' + str(t_now-t_last) + + ' seconds, solution distance = ' + str(solution_distance)) + else: + print('Finished cycle #' + str(completed_cycles) + ' of ' + str(agent.cycles) + + ' in ' + str(t_now-t_last) + ' seconds.') + t_last = t_now + + # Record the last cycle if horizon is infinite (solution is still empty!) + if infinite_horizon: + solution = solution_cycle # PseudoTerminal=False impossible for infinite horizon + + # Restore the direction of time to its original orientation, then return the solution + if original_time_flow: + agent.timeFwd() + return solution + + +def solveOneCycle(agent, solution_last): + ''' + Solve one "cycle" of the dynamic model for one agent type. This function + iterates over the periods within an agent's cycle, updating the time-varying + parameters and passing them to the single period solver(s). + + Parameters + ---------- + agent : AgentType + The microeconomic AgentType whose dynamic problem is to be solved. + solution_last : Solution + A representation of the solution of the period that comes after the + end of the sequence of one period problems. This might be the term- + inal period solution, a "pseudo terminal" solution, or simply the + solution to the earliest period from the succeeding cycle. + + Returns + ------- + solution_cycle : [Solution] + A list of one period solutions for one "cycle" of the AgentType's + microeconomic model. Returns in reverse chronological order. + ''' + # Calculate number of periods per cycle, defaults to 1 if all variables are time invariant + if len(agent.time_vary) > 0: + # name = agent.time_vary[0] + # T = len(eval('agent.' + name)) + T = len(agent.__dict__[agent.time_vary[0]]) + else: + T = 1 + + # Construct a dictionary to be passed to the solver + # time_inv_string = '' + # for name in agent.time_inv: + # time_inv_string += ' \'' + name + '\' : agent.' + name + ',' + # time_vary_string = '' + # for name in agent.time_vary: + # time_vary_string += ' \'' + name + '\' : None,' + # solve_dict = eval('{' + time_inv_string + time_vary_string + '}') + solve_dict = {parameter: agent.__dict__[parameter] for parameter in agent.time_inv} + solve_dict.update({parameter: None for parameter in agent.time_vary}) + + # Initialize the solution for this cycle, then iterate on periods + solution_cycle = [] + solution_next = solution_last + for t in range(T): + # Update which single period solver to use (if it depends on time) + if hasattr(agent.solveOnePeriod, "__getitem__"): + solveOnePeriod = agent.solveOnePeriod[t] + else: + solveOnePeriod = agent.solveOnePeriod + + these_args = getArgNames(solveOnePeriod) + + # Update time-varying single period inputs + for name in agent.time_vary: + if name in these_args: + # solve_dict[name] = eval('agent.' + name + '[t]') + solve_dict[name] = agent.__dict__[name][t] + solve_dict['solution_next'] = solution_next + + # Make a temporary dictionary for this period + temp_dict = {name: solve_dict[name] for name in these_args} + + # Solve one period, add it to the solution, and move to the next period + solution_t = solveOnePeriod(**temp_dict) + solution_cycle.append(solution_t) + solution_next = solution_t + + # Return the list of per-period solutions + return solution_cycle + + +# ======================================================================== +# ======================================================================== + +class Market(HARKobject): + ''' + A superclass to represent a central clearinghouse of information. Used for + dynamic general equilibrium models to solve the "macroeconomic" model as a + layer on top of the "microeconomic" models of one or more AgentTypes. + ''' + def __init__(self, agents=[], sow_vars=[], reap_vars=[], const_vars=[], track_vars=[], dyn_vars=[], + millRule=None, calcDynamics=None, act_T=1000, tolerance=0.000001,**kwds): + ''' + Make a new instance of the Market class. + + Parameters + ---------- + agents : [AgentType] + A list of all the AgentTypes in this market. + sow_vars : [string] + Names of variables generated by the "aggregate market process" that should + be "sown" to the agents in the market. Aggregate state, etc. + reap_vars : [string] + Names of variables to be collected ("reaped") from agents in the market + to be used in the "aggregate market process". + const_vars : [string] + Names of attributes of the Market instance that are used in the "aggregate + market process" but do not come from agents-- they are constant or simply + parameters inherent to the process. + track_vars : [string] + Names of variables generated by the "aggregate market process" that should + be tracked as a "history" so that a new dynamic rule can be calculated. + This is often a subset of sow_vars. + dyn_vars : [string] + Names of variables that constitute a "dynamic rule". + millRule : function + A function that takes inputs named in reap_vars and returns an object + with attributes named in sow_vars. The "aggregate market process" that + transforms individual agent actions/states/data into aggregate data to + be sent back to agents. + calcDynamics : function + A function that takes inputs named in track_vars and returns an object + with attributes named in dyn_vars. Looks at histories of aggregate + variables and generates a new "dynamic rule" for agents to believe and + act on. + act_T : int + The number of times that the "aggregate market process" should be run + in order to generate a history of aggregate variables. + tolerance: float + Minimum acceptable distance between "dynamic rules" to consider the + Market solution process converged. Distance is a user-defined metric. + + Returns + ------- + None + ''' + self.agents = agents # NOQA + self.reap_vars = reap_vars # NOQA + self.sow_vars = sow_vars # NOQA + self.const_vars = const_vars # NOQA + self.track_vars = track_vars # NOQA + self.dyn_vars = dyn_vars # NOQA + if millRule is not None: # To prevent overwriting of method-based millRules + self.millRule = millRule + if calcDynamics is not None: # Ditto for calcDynamics + self.calcDynamics = calcDynamics + self.act_T = act_T # NOQA + self.tolerance = tolerance # NOQA + self.max_loops = 1000 # NOQA + self.assignParameters(**kwds) + + self.print_parallel_error_once = True + # Print the error associated with calling the parallel method + # "solveAgents" one time. If set to false, the error will never + # print. See "solveAgents" for why this prints once or never. + + def solveAgents(self): + ''' + Solves the microeconomic problem for all AgentTypes in this market. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + # for this_type in self.agents: + # this_type.solve() + try: + multiThreadCommands(self.agents, ['solve()']) + except Exception as err: + if self.print_parallel_error_once: + # Set flag to False so this is only printed once. + self.print_parallel_error_once = False + print("**** WARNING: could not execute multiThreadCommands in HARK.core.Market.solveAgents() ", + "so using the serial version instead. This will likely be slower. " + "The multiTreadCommands() functions failed with the following error:", '\n', + sys.exc_info()[0], ':', err) # sys.exc_info()[0]) + multiThreadCommandsFake(self.agents, ['solve()']) + + def solve(self): + ''' + "Solves" the market by finding a "dynamic rule" that governs the aggregate + market state such that when agents believe in these dynamics, their actions + collectively generate the same dynamic rule. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + go = True + max_loops = self.max_loops # Failsafe against infinite solution loop + completed_loops = 0 + old_dynamics = None + + while go: # Loop until the dynamic process converges or we hit the loop cap + self.solveAgents() # Solve each AgentType's micro problem + self.makeHistory() # "Run" the model while tracking aggregate variables + new_dynamics = self.updateDynamics() # Find a new aggregate dynamic rule + + # Check to see if the dynamic rule has converged (if this is not the first loop) + if completed_loops > 0: + distance = new_dynamics.distance(old_dynamics) + else: + distance = 1000000.0 + + # Move to the next loop if the terminal conditions are not met + old_dynamics = new_dynamics + completed_loops += 1 + go = distance >= self.tolerance and completed_loops < max_loops + + self.dynamics = new_dynamics # Store the final dynamic rule in self + + def reap(self): + ''' + Collects attributes named in reap_vars from each AgentType in the market, + storing them in respectively named attributes of self. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + for var_name in self.reap_vars: + harvest = [] + for this_type in self.agents: + harvest.append(getattr(this_type, var_name)) + setattr(self, var_name, harvest) + + def sow(self): + ''' + Distributes attrributes named in sow_vars from self to each AgentType + in the market, storing them in respectively named attributes. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + for var_name in self.sow_vars: + this_seed = getattr(self, var_name) + for this_type in self.agents: + setattr(this_type, var_name, this_seed) + + def mill(self): + ''' + Processes the variables collected from agents using the function millRule, + storing the results in attributes named in aggr_sow. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + # Make a dictionary of inputs for the millRule + reap_vars_string = '' + for name in self.reap_vars: + reap_vars_string += ' \'' + name + '\' : self.' + name + ',' + const_vars_string = '' + for name in self.const_vars: + const_vars_string += ' \'' + name + '\' : self.' + name + ',' + mill_dict = eval('{' + reap_vars_string + const_vars_string + '}') + + # Run the millRule and store its output in self + product = self.millRule(**mill_dict) + for j in range(len(self.sow_vars)): + this_var = self.sow_vars[j] + this_product = getattr(product, this_var) + setattr(self, this_var, this_product) + + def cultivate(self): + ''' + Has each AgentType in agents perform their marketAction method, using + variables sown from the market (and maybe also "private" variables). + The marketAction method should store new results in attributes named in + reap_vars to be reaped later. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + for this_type in self.agents: + this_type.marketAction() + + def reset(self): + ''' + Reset the state of the market (attributes in sow_vars, etc) to some + user-defined initial state, and erase the histories of tracked variables. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + for var_name in self.track_vars: # Reset the history of tracked variables + setattr(self, var_name + '_hist', []) + for var_name in self.sow_vars: # Set the sow variables to their initial levels + initial_val = getattr(self, var_name + '_init') + setattr(self, var_name, initial_val) + for this_type in self.agents: # Reset each AgentType in the market + this_type.reset() + + def store(self): + ''' + Record the current value of each variable X named in track_vars in an + attribute named X_hist. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + for var_name in self.track_vars: + value_now = getattr(self, var_name) + getattr(self, var_name + '_hist').append(value_now) + + def makeHistory(self): + ''' + Runs a loop of sow-->cultivate-->reap-->mill act_T times, tracking the + evolution of variables X named in track_vars in attributes named X_hist. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + self.reset() # Initialize the state of the market + for t in range(self.act_T): + self.sow() # Distribute aggregated information/state to agents + self.cultivate() # Agents take action + self.reap() # Collect individual data from agents + self.mill() # Process individual data into aggregate data + self.store() # Record variables of interest + + def updateDynamics(self): + ''' + Calculates a new "aggregate dynamic rule" using the history of variables + named in track_vars, and distributes this rule to AgentTypes in agents. + + Parameters + ---------- + none + + Returns + ------- + dynamics : instance + The new "aggregate dynamic rule" that agents believe in and act on. + Should have attributes named in dyn_vars. + ''' + # Make a dictionary of inputs for the dynamics calculator + history_vars_string = '' + arg_names = list(getArgNames(self.calcDynamics)) + if 'self' in arg_names: + arg_names.remove('self') + for name in arg_names: + history_vars_string += ' \'' + name + '\' : self.' + name + '_hist,' + update_dict = eval('{' + history_vars_string + '}') + + # Calculate a new dynamic rule and distribute it to the agents in agent_list + dynamics = self.calcDynamics(**update_dict) # User-defined dynamics calculator + for var_name in self.dyn_vars: + this_obj = getattr(dynamics, var_name) + for this_type in self.agents: + setattr(this_type, var_name, this_obj) + return dynamics + + +# ------------------------------------------------------------------------------ +# Code to copy entire modules to a local directory +# ------------------------------------------------------------------------------ + +# Define a function to run the copying: +def copy_module(target_path, my_directory_full_path, my_module): + ''' + Helper function for copy_module_to_local(). Provides the actual copy + functionality, with highly cautious safeguards against copying over + important things. + + Parameters + ---------- + target_path : string + String, file path to target location + + my_directory_full_path: string + String, full pathname to this file's directory + + my_module : string + String, name of the module to copy + + Returns + ------- + none + ''' + + if target_path == 'q' or target_path == 'Q': + print("Goodbye!") + return + elif target_path == os.path.expanduser("~") or os.path.normpath(target_path) == os.path.expanduser("~"): + print("You have indicated that the target location is " + target_path + + " -- that is, you want to wipe out your home directory with the contents of " + my_module + + ". My programming does not allow me to do that.\n\nGoodbye!") + return + elif os.path.exists(target_path): + print("There is already a file or directory at the location " + target_path + + ". For safety reasons this code does not overwrite existing files.\n Please remove the file at " + + target_path + + " and try again.") + return + else: + user_input = input("""You have indicated you want to copy module:\n """ + my_module + + """\nto:\n """ + target_path + """\nIs that correct? Please indicate: y / [n]\n\n""") + if user_input == 'y' or user_input == 'Y': + # print("copy_tree(",my_directory_full_path,",", target_path,")") + copy_tree(my_directory_full_path, target_path) + else: + print("Goodbye!") + return + + +def print_helper(): + + my_directory_full_path = os.path.dirname(os.path.realpath(__file__)) + + print(my_directory_full_path) + + +def copy_module_to_local(full_module_name): + ''' + This function contains simple code to copy a submodule to a location on + your hard drive, as specified by you. The purpose of this code is to provide + users with a simple way to access a *copy* of code that usually sits deep in + the Econ-ARK package structure, for purposes of tinkering and experimenting + directly. This is meant to be a simple way to explore HARK code. To interact + with the codebase under active development, please refer to the documentation + under github.com/econ-ark/HARK/ + + To execute, do the following on the Python command line: + + from HARK.core import copy_module_to_local + copy_module_to_local("FULL-HARK-MODULE-NAME-HERE") + + For example, if you want SolvingMicroDSOPs you would enter + + from HARK.core import copy_module_to_local + copy_module_to_local("HARK.SolvingMicroDSOPs") + + ''' + + # Find a default directory -- user home directory: + home_directory_RAW = os.path.expanduser("~") + # Thanks to https://stackoverflow.com/a/4028943 + + # Find the directory of the HARK.core module: + # my_directory_full_path = os.path.dirname(os.path.realpath(__file__)) + hark_core_directory_full_path = os.path.dirname(os.path.realpath(__file__)) + # From https://stackoverflow.com/a/5137509 + # Important note from that answer: + # (Note that the incantation above won't work if you've already used os.chdir() + # to change your current working directory, + # since the value of the __file__ constant is relative to the current working directory and is not changed by an + # os.chdir() call.) + # + # NOTE: for this specific file that I am testing, the path should be: + # '/home/npalmer/anaconda3/envs/py3fresh/lib/python3.6/site-packages/HARK/SolvingMicroDSOPs/---example-file--- + + # Split out the name of the module. Break if proper format is not followed: + all_module_names_list = full_module_name.split('.') # Assume put in at correct format + if all_module_names_list[0] != "HARK": + print("\nWarning: the module name does not start with 'HARK'. Instead it is: '" + + all_module_names_list[0]+"' --please format the full namespace of the module you want. \n" + "For example, 'HARK.SolvingMicroDSOPs'") + print("\nGoodbye!") + return + + # Construct the pathname to the module to copy: + my_directory_full_path = hark_core_directory_full_path + for a_directory_name in all_module_names_list[1:]: + my_directory_full_path = os.path.join(my_directory_full_path, a_directory_name) + + head_path, my_module = os.path.split(my_directory_full_path) + + home_directory_with_module = os.path.join(home_directory_RAW, my_module) + + print("\n\n\nmy_directory_full_path:", my_directory_full_path, '\n\n\n') + + # Interact with the user: + # - Ask the user for the target place to copy the directory + # - Offer use "q/y/other" option + # - Check if there is something there already + # - If so, ask if should replace + # - If not, just copy there + # - Quit + + target_path = input("""You have invoked the 'replicate' process for the current module:\n """ + + my_module + """\nThe default copy location is your home directory:\n """ + + home_directory_with_module + """\nPlease enter one of the three options in single quotes below, excluding the quotes: + + 'q' or return/enter to quit the process + 'y' to accept the default home directory: """+home_directory_with_module+""" + 'n' to specify your own pathname\n\n""") + + if target_path == 'n' or target_path == 'N': + target_path = input("""Please enter the full pathname to your target directory location: """) + + # Clean up: + target_path = os.path.expanduser(target_path) + target_path = os.path.expandvars(target_path) + target_path = os.path.normpath(target_path) + + # Check to see if they included the module name; if not add it here: + temp_head, temp_tail = os.path.split(target_path) + if temp_tail != my_module: + target_path = os.path.join(target_path, my_module) + + elif target_path == 'y' or target_path == 'Y': + # Just using the default path: + target_path = home_directory_with_module + else: + # Assume "quit" + return + + if target_path != 'q' and target_path != 'Q' or target_path == '': + # Run the copy command: + copy_module(target_path, my_directory_full_path, my_module) + + return + + if target_path != 'q' and target_path != 'Q' or target_path == '': + # Run the copy command: + copy_module(target_path, my_directory_full_path, my_module) + + return + + +def main(): + print("Sorry, HARK.core doesn't actually do anything on its own.") + print("To see some examples of its frameworks in action, try running a model module.") + print("Several interesting model modules can be found in /ConsumptionSavingModel.") + print('For an extraordinarily simple model that demonstrates the "microeconomic" and') + print('"macroeconomic" frameworks, see /FashionVictim/FashionVictimModel.') + + +if __name__ == '__main__': + main() diff --git a/HARK/#scratch# b/HARK/#scratch# new file mode 100644 index 000000000..e69de29bb diff --git a/HARK/ConsumptionSaving/ConsAggShockModel.py b/HARK/ConsumptionSaving/ConsAggShockModel.py index d5b99ea65..5a0ca601b 100644 --- a/HARK/ConsumptionSaving/ConsAggShockModel.py +++ b/HARK/ConsumptionSaving/ConsAggShockModel.py @@ -10,7 +10,7 @@ from builtins import range import numpy as np import scipy.stats as stats -from HARK.distribution import DiscreteDistribution, combineIndepDstns, approxMeanOneLognormal +from HARK.distribution import DiscreteDistribution, combineIndepDstns, MeanOneLogNormal from HARK.interpolation import LinearInterp, LinearInterpOnInterp1D, ConstantFunction, IdentityFunction,\ VariableLowerBoundFunc2D, BilinearInterp, LowerEnvelope2D, UpperEnvelope from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv,\ @@ -1095,8 +1095,12 @@ def makeAggShkDstn(self): ------- None ''' - self.TranShkAggDstn = approxMeanOneLognormal(sigma=self.TranShkAggStd, N=self.TranShkAggCount) - self.PermShkAggDstn = approxMeanOneLognormal(sigma=self.PermShkAggStd, N=self.PermShkAggCount) + self.TranShkAggDstn = MeanOneLogNormal( + sigma=self.TranShkAggStd + ).approx(N=self.TranShkAggCount) + self.PermShkAggDstn = MeanOneLognormal( + sigma=self.PermShkAggStd + ).approx(N=self.PermShkAggCount) self.AggShkDstn = combineIndepDstns(self.PermShkAggDstn, self.TranShkAggDstn) def reset(self): @@ -1312,8 +1316,12 @@ def makeAggShkDstn(self): ------- None ''' - self.TranShkAggDstn = approxMeanOneLognormal(sigma=self.TranShkAggStd, N=self.TranShkAggCount) - self.PermShkAggDstn = approxMeanOneLognormal(sigma=self.PermShkAggStd, N=self.PermShkAggCount) + self.TranShkAggDstn = MeanOneLogNormal( + sigma=self.TranShkAggStd + ).approx(N=self.TranShkAggCount) + self.PermShkAggDstn = MeanOneLogNormal( + sigma=self.PermShkAggStd + ).approx(N=self.PermShkAggCount) self.AggShkDstn = combineIndepDstns(self.PermShkAggDstn, self.TranShkAggDstn) def millRule(self): @@ -1526,8 +1534,8 @@ def makeAggShkDstn(self): StateCount = self.MrkvArray.shape[0] for i in range(StateCount): - TranShkAggDstn.append(approxMeanOneLognormal(sigma=self.TranShkAggStd[i], N=self.TranShkAggCount)) - PermShkAggDstn.append(approxMeanOneLognormal(sigma=self.PermShkAggStd[i], N=self.PermShkAggCount)) + TranShkAggDstn.append(MeanOneLogNormal(sigma=self.TranShkAggStd[i]).approx(N=self.TranShkAggCount)) + PermShkAggDstn.append(MeanOneLogNormal(sigma=self.PermShkAggStd[i]).approx(N=self.PermShkAggCount)) AggShkDstn.append(combineIndepDstns(PermShkAggDstn[-1], TranShkAggDstn[-1])) self.TranShkAggDstn = TranShkAggDstn diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 8a474f720..7c4e17ca7 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -24,8 +24,8 @@ from HARK import AgentType, Solution, NullFunc, HARKobject from HARK.utilities import warnings # Because of "patch" to warnings modules from HARK.interpolation import CubicInterp, LowerEnvelope, LinearInterp -from HARK.distribution import Lognormal, Uniform -from HARK.distribution import DiscreteDistribution, approxMeanOneLognormal, addDiscreteOutcomeConstantMean, combineIndepDstns +from HARK.distribution import Lognormal, MeanOneLogNormal, Uniform +from HARK.distribution import DiscreteDistribution, addDiscreteOutcomeConstantMean, combineIndepDstns from HARK.utilities import makeGridExpMult, CRRAutility, CRRAutilityP, \ CRRAutilityPP, CRRAutilityP_inv, CRRAutility_invP, CRRAutility_inv, \ CRRAutilityP_invP @@ -2257,11 +2257,19 @@ def makeEulerErrorFunc(self,mMax=100,approx_inc_dstn=True): if approx_inc_dstn: IncomeDstn = self.IncomeDstn[0] else: - TranShkDstn = approxMeanOneLognormal(N=200,sigma=self.TranShkStd[0], - tail_N=50,tail_order=1.3, tail_bound=[0.05,0.95]) + TranShkDstn = MeanOneLognormal( + sigma=self.TranShkStd[0] + ).approx(N=200, + tail_N=50, + tail_order=1.3, + tail_bound=[0.05,0.95]) TranShkDstn = addDiscreteOutcomeConstantMean(TranShkDstn,self.UnempPrb,self.IncUnemp) - PermShkDstn = approxMeanOneLognormal(N=200,sigma=self.PermShkStd[0], - tail_N=50,tail_order=1.3, tail_bound=[0.05,0.95]) + PermShkDstn = MeanOneLogNormal( + sigma=self.PermShkStd[0] + ).approx(N=200, + tail_N=50, + tail_order=1.3, + tail_bound=[0.05,0.95]) IncomeDstn = combineIndepDstns(PermShkDstn,TranShkDstn) # Make a grid of market resources @@ -2796,10 +2804,14 @@ def constructLognormalIncomeProcessUnemployment(parameters): TranShkDstn.append([ShkPrbsRet,TranShkValsRet]) else: # We are in the "working life" periods. - TranShkDstn_t = approxMeanOneLognormal(N=TranShkCount, sigma=TranShkStd[t], tail_N=0) + TranShkDstn_t = MeanOneLogNormal( + sigma=TranShkStd[t] + ).approx(TranShkCount, tail_N=0) if UnempPrb > 0: TranShkDstn_t = addDiscreteOutcomeConstantMean(TranShkDstn_t, p=UnempPrb, x=IncUnemp) - PermShkDstn_t = approxMeanOneLognormal(N=PermShkCount, sigma=PermShkStd[t], tail_N=0) + PermShkDstn_t = MeanOneLogNormal( + sigma=PermShkStd[t] + ).approx(PermShkCount, tail_N=0) IncomeDstn.append(combineIndepDstns(PermShkDstn_t,TranShkDstn_t)) # mix the independent distributions PermShkDstn.append(PermShkDstn_t) TranShkDstn.append(TranShkDstn_t) diff --git a/HARK/ConsumptionSaving/ConsMarkovModel.py.orig b/HARK/ConsumptionSaving/ConsMarkovModel.py.orig new file mode 100644 index 000000000..f482770d7 --- /dev/null +++ b/HARK/ConsumptionSaving/ConsMarkovModel.py.orig @@ -0,0 +1,1023 @@ +''' +Classes to solve and simulate consumption-savings model with a discrete, exogenous, +stochastic Markov state. The only solver here extends ConsIndShockModel to +include a Markov state; the interest factor, permanent growth factor, and income +distribution can vary with the discrete state. +''' +from __future__ import division, print_function +from __future__ import absolute_import +from builtins import range +from copy import deepcopy +import numpy as np +from HARK import AgentType +from HARK.ConsumptionSaving.ConsIndShockModel import ConsIndShockSolver, ValueFunc, \ + MargValueFunc, ConsumerSolution, IndShockConsumerType +<<<<<<< HEAD +from HARK.distribution import Uniform +======= + + +from HARK.distribution import DiscreteDistribution +from HARK.simulation import Uniform +>>>>>>> 8fc1cd2b34dbb8b29e0c088182bc1d737c7f07ad +from HARK.interpolation import CubicInterp, LowerEnvelope, LinearInterp +from HARK.utilities import CRRAutility, CRRAutilityP, CRRAutilityPP, CRRAutilityP_inv, \ + CRRAutility_invP, CRRAutility_inv, CRRAutilityP_invP + +__all__ = ['ConsMarkovSolver', 'MarkovConsumerType'] + +utility = CRRAutility +utilityP = CRRAutilityP +utilityPP = CRRAutilityPP +utilityP_inv = CRRAutilityP_inv +utility_invP = CRRAutility_invP +utility_inv = CRRAutility_inv +utilityP_invP = CRRAutilityP_invP + +class ConsMarkovSolver(ConsIndShockSolver): + ''' + A class to solve a single period consumption-saving problem with risky income + and stochastic transitions between discrete states, in a Markov fashion. + Extends ConsIndShockSolver, with identical inputs but for a discrete + Markov state, whose transition rule is summarized in MrkvArray. Markov + states can differ in their interest factor, permanent growth factor, live probability, and + income distribution, so the inputs Rfree, PermGroFac, IncomeDstn, and LivPrb are + now arrays or lists specifying those values in each (succeeding) Markov state. + ''' + def __init__(self,solution_next,IncomeDstn_list,LivPrb,DiscFac, + CRRA,Rfree_list,PermGroFac_list,MrkvArray,BoroCnstArt, + aXtraGrid,vFuncBool,CubicBool): + ''' + Constructor for a new solver for a one period problem with risky income + and transitions between discrete Markov states. In the descriptions below, + N is the number of discrete states. + + Parameters + ---------- + solution_next : ConsumerSolution + The solution to next period's one period problem. + IncomeDstn_list : [[np.array]] + A length N list of income distributions in each succeeding Markov + state. Each income distribution contains three arrays of floats, + representing a discrete approximation to the income process at the + beginning of the succeeding period. Order: event probabilities, + permanent shocks, transitory shocks. + LivPrb : np.array + Survival probability; likelihood of being alive at the beginning of + the succeeding period for each Markov state. + DiscFac : float + Intertemporal discount factor for future utility. + CRRA : float + Coefficient of relative risk aversion. + Rfree_list : np.array + Risk free interest factor on end-of-period assets for each Markov + state in the succeeding period. + PermGroFac_list : np.array + Expected permanent income growth factor at the end of this period + for each Markov state in the succeeding period. + MrkvArray : np.array + An NxN array representing a Markov transition matrix between discrete + states. The i,j-th element of MrkvArray is the probability of + moving from state i in period t to state j in period t+1. + BoroCnstArt: float or None + Borrowing constraint for the minimum allowable assets to end the + period with. If it is less than the natural borrowing constraint, + then it is irrelevant; BoroCnstArt=None indicates no artificial bor- + rowing constraint. + aXtraGrid: np.array + Array of "extra" end-of-period asset values-- assets above the + absolute minimum acceptable level. + vFuncBool: boolean + An indicator for whether the value function should be computed and + included in the reported solution. + CubicBool: boolean + An indicator for whether the solver should use cubic or linear inter- + polation. + + Returns + ------- + None + ''' + # Set basic attributes of the problem + ConsIndShockSolver.assignParameters(self,solution_next,np.nan,LivPrb,DiscFac,CRRA,np.nan, + np.nan,BoroCnstArt,aXtraGrid,vFuncBool,CubicBool) + self.defUtilityFuncs() + + # Set additional attributes specific to the Markov model + self.IncomeDstn_list = IncomeDstn_list + self.Rfree_list = Rfree_list + self.PermGroFac_list = PermGroFac_list + self.MrkvArray = MrkvArray + self.StateCount = MrkvArray.shape[0] + + def solve(self): + ''' + Solve the one period problem of the consumption-saving model with a Markov state. + + Parameters + ---------- + none + + Returns + ------- + solution : ConsumerSolution + The solution to the single period consumption-saving problem. Includes + a consumption function cFunc (using cubic or linear splines), a marg- + inal value function vPfunc, a minimum acceptable level of normalized + market resources mNrmMin, normalized human wealth hNrm, and bounding + MPCs MPCmin and MPCmax. It might also have a value function vFunc + and marginal marginal value function vPPfunc. All of these attributes + are lists or arrays, with elements corresponding to the current + Markov state. E.g. solution.cFunc[0] is the consumption function + when in the i=0 Markov state this period. + ''' + # Find the natural borrowing constraint in each current state + self.defBoundary() + + # Initialize end-of-period (marginal) value functions + self.EndOfPrdvFunc_list = [] + self.EndOfPrdvPfunc_list = [] + self.ExIncNextAll = np.zeros(self.StateCount) + np.nan # expected income conditional on the next state + self.WorstIncPrbAll = np.zeros(self.StateCount) + np.nan # probability of getting the worst income shock in each next period state + + # Loop through each next-period-state and calculate the end-of-period + # (marginal) value function + for j in range(self.StateCount): + # Condition values on next period's state (and record a couple for later use) + self.conditionOnState(j) + self.ExIncNextAll[j] = np.dot(self.ShkPrbsNext,self.PermShkValsNext*self.TranShkValsNext) + self.WorstIncPrbAll[j] = self.WorstIncPrb + + # Construct the end-of-period marginal value function conditional + # on next period's state and add it to the list of value functions + EndOfPrdvPfunc_cond = self.makeEndOfPrdvPfuncCond() + self.EndOfPrdvPfunc_list.append(EndOfPrdvPfunc_cond) + + # Construct the end-of-period value functional conditional on next + # period's state and add it to the list of value functions + if self.vFuncBool: + EndOfPrdvFunc_cond = self.makeEndOfPrdvFuncCond() + self.EndOfPrdvFunc_list.append(EndOfPrdvFunc_cond) + + # EndOfPrdvP_cond is EndOfPrdvP conditional on *next* period's state. + # Take expectations to get EndOfPrdvP conditional on *this* period's state. + self.calcEndOfPrdvP() + + # Calculate the bounding MPCs and PDV of human wealth for each state + self.calcHumWealthAndBoundingMPCs() + + # Find consumption and market resources corresponding to each end-of-period + # assets point for each state (and add an additional point at the lower bound) + aNrm = np.asarray(self.aXtraGrid)[np.newaxis,:] + np.array(self.BoroCnstNat_list)[:,np.newaxis] + self.getPointsForInterpolation(self.EndOfPrdvP,aNrm) + cNrm = np.hstack((np.zeros((self.StateCount,1)),self.cNrmNow)) + mNrm = np.hstack((np.reshape(self.mNrmMin_list,(self.StateCount,1)),self.mNrmNow)) + + # Package and return the solution for this period + self.BoroCnstNat = self.BoroCnstNat_list + solution = self.makeSolution(cNrm,mNrm) + return solution + + def defBoundary(self): + ''' + Find the borrowing constraint for each current state and save it as an + attribute of self for use by other methods. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + self.BoroCnstNatAll = np.zeros(self.StateCount) + np.nan + # Find the natural borrowing constraint conditional on next period's state + for j in range(self.StateCount): + PermShkMinNext = np.min(self.IncomeDstn_list[j].X[0]) + TranShkMinNext = np.min(self.IncomeDstn_list[j].X[1]) + self.BoroCnstNatAll[j] = (self.solution_next.mNrmMin[j] - TranShkMinNext)*\ + (self.PermGroFac_list[j]*PermShkMinNext)/self.Rfree_list[j] + + self.BoroCnstNat_list = np.zeros(self.StateCount) + np.nan + self.mNrmMin_list = np.zeros(self.StateCount) + np.nan + self.BoroCnstDependency = np.zeros((self.StateCount,self.StateCount)) + np.nan + # The natural borrowing constraint in each current state is the *highest* + # among next-state-conditional natural borrowing constraints that could + # occur from this current state. + for i in range(self.StateCount): + possible_next_states = self.MrkvArray[i,:] > 0 + self.BoroCnstNat_list[i] = np.max(self.BoroCnstNatAll[possible_next_states]) + + # Explicitly handle the "None" case: + if self.BoroCnstArt is None: + self.mNrmMin_list[i] = self.BoroCnstNat_list[i] + else: + self.mNrmMin_list[i] = np.max([self.BoroCnstNat_list[i],self.BoroCnstArt]) + self.BoroCnstDependency[i,:] = self.BoroCnstNat_list[i] == self.BoroCnstNatAll + # Also creates a Boolean array indicating whether the natural borrowing + # constraint *could* be hit when transitioning from i to j. + + def conditionOnState(self,state_index): + ''' + Temporarily assume that a particular Markov state will occur in the + succeeding period, and condition solver attributes on this assumption. + Allows the solver to construct the future-state-conditional marginal + value function (etc) for that future state. + + Parameters + ---------- + state_index : int + Index of the future Markov state to condition on. + + Returns + ------- + none + ''' + # Set future-state-conditional values as attributes of self + self.IncomeDstn = self.IncomeDstn_list[state_index] + self.Rfree = self.Rfree_list[state_index] + self.PermGroFac = self.PermGroFac_list[state_index] + self.vPfuncNext = self.solution_next.vPfunc[state_index] + self.mNrmMinNow = self.mNrmMin_list[state_index] + self.BoroCnstNat = self.BoroCnstNatAll[state_index] + self.setAndUpdateValues(self.solution_next,self.IncomeDstn,self.LivPrb,self.DiscFac) + self.DiscFacEff = self.DiscFac # survival probability LivPrb represents probability from + # *current* state, so DiscFacEff is just DiscFac for now + + # These lines have to come after setAndUpdateValues to override the definitions there + self.vPfuncNext = self.solution_next.vPfunc[state_index] + if self.CubicBool: + self.vPPfuncNext= self.solution_next.vPPfunc[state_index] + if self.vFuncBool: + self.vFuncNext = self.solution_next.vFunc[state_index] + + def calcEndOfPrdvPP(self): + ''' + Calculates end-of-period marginal marginal value using a pre-defined + array of next period market resources in self.mNrmNext. + + Parameters + ---------- + none + + Returns + ------- + EndOfPrdvPP : np.array + End-of-period marginal marginal value of assets at each value in + the grid of assets. + ''' + EndOfPrdvPP = self.DiscFacEff*self.Rfree*self.Rfree*self.PermGroFac**(-self.CRRA-1.0)*\ + np.sum(self.PermShkVals_temp**(-self.CRRA-1.0)*self.vPPfuncNext(self.mNrmNext) + *self.ShkPrbs_temp,axis=0) + return EndOfPrdvPP + + def makeEndOfPrdvFuncCond(self): + ''' + Construct the end-of-period value function conditional on next period's + state. NOTE: It might be possible to eliminate this method and replace + it with ConsIndShockSolver.makeEndOfPrdvFunc, but the self.X_cond + variables must be renamed. + + Parameters + ---------- + none + + Returns + ------- + EndofPrdvFunc_cond : ValueFunc + The end-of-period value function conditional on a particular state + occuring in the next period. + ''' + VLvlNext = (self.PermShkVals_temp**(1.0-self.CRRA)* + self.PermGroFac**(1.0-self.CRRA))*self.vFuncNext(self.mNrmNext) + EndOfPrdv_cond = self.DiscFacEff*np.sum(VLvlNext*self.ShkPrbs_temp,axis=0) + EndOfPrdvNvrs_cond = self.uinv(EndOfPrdv_cond) + EndOfPrdvNvrsP_cond = self.EndOfPrdvP_cond*self.uinvP(EndOfPrdv_cond) + EndOfPrdvNvrs_cond = np.insert(EndOfPrdvNvrs_cond,0,0.0) + EndOfPrdvNvrsP_cond = np.insert(EndOfPrdvNvrsP_cond,0,EndOfPrdvNvrsP_cond[0]) + aNrm_temp = np.insert(self.aNrm_cond,0,self.BoroCnstNat) + EndOfPrdvNvrsFunc_cond = CubicInterp(aNrm_temp,EndOfPrdvNvrs_cond,EndOfPrdvNvrsP_cond) + EndofPrdvFunc_cond = ValueFunc(EndOfPrdvNvrsFunc_cond,self.CRRA) + return EndofPrdvFunc_cond + + + def calcEndOfPrdvPcond(self): + ''' + Calculate end-of-period marginal value of assets at each point in aNrmNow + conditional on a particular state occuring in the next period. + + Parameters + ---------- + None + + Returns + ------- + EndOfPrdvP : np.array + A 1D array of end-of-period marginal value of assets. + ''' + EndOfPrdvPcond = ConsIndShockSolver.calcEndOfPrdvP(self) + return EndOfPrdvPcond + + + def makeEndOfPrdvPfuncCond(self): + ''' + Construct the end-of-period marginal value function conditional on next + period's state. + + Parameters + ---------- + None + + Returns + ------- + EndofPrdvPfunc_cond : MargValueFunc + The end-of-period marginal value function conditional on a particular + state occuring in the succeeding period. + ''' + # Get data to construct the end-of-period marginal value function (conditional on next state) + self.aNrm_cond = self.prepareToCalcEndOfPrdvP() + self.EndOfPrdvP_cond= self.calcEndOfPrdvPcond() + EndOfPrdvPnvrs_cond = self.uPinv(self.EndOfPrdvP_cond) # "decurved" marginal value + if self.CubicBool: + EndOfPrdvPP_cond = self.calcEndOfPrdvPP() + EndOfPrdvPnvrsP_cond = EndOfPrdvPP_cond*self.uPinvP(self.EndOfPrdvP_cond) # "decurved" marginal marginal value + + # Construct the end-of-period marginal value function conditional on the next state. + if self.CubicBool: + EndOfPrdvPnvrsFunc_cond = CubicInterp(self.aNrm_cond,EndOfPrdvPnvrs_cond, + EndOfPrdvPnvrsP_cond,lower_extrap=True) + else: + EndOfPrdvPnvrsFunc_cond = LinearInterp(self.aNrm_cond,EndOfPrdvPnvrs_cond, + lower_extrap=True) + EndofPrdvPfunc_cond = MargValueFunc(EndOfPrdvPnvrsFunc_cond,self.CRRA) # "recurve" the interpolated marginal value function + return EndofPrdvPfunc_cond + + def calcEndOfPrdvP(self): + ''' + Calculates end of period marginal value (and marginal marginal) value + at each aXtra gridpoint for each current state, unconditional on the + future Markov state (i.e. weighting conditional end-of-period marginal + value by transition probabilities). + + Parameters + ---------- + none + + Returns + ------- + none + ''' + # Find unique values of minimum acceptable end-of-period assets (and the + # current period states for which they apply). + aNrmMin_unique, state_inverse = np.unique(self.BoroCnstNat_list,return_inverse=True) + self.possible_transitions = self.MrkvArray > 0 + + # Calculate end-of-period marginal value (and marg marg value) at each + # asset gridpoint for each current period state + EndOfPrdvP = np.zeros((self.StateCount,self.aXtraGrid.size)) + EndOfPrdvPP = np.zeros((self.StateCount,self.aXtraGrid.size)) + for k in range(aNrmMin_unique.size): + aNrmMin = aNrmMin_unique[k] # minimum assets for this pass + which_states = state_inverse == k # the states for which this minimum applies + aGrid = aNrmMin + self.aXtraGrid # assets grid for this pass + EndOfPrdvP_all = np.zeros((self.StateCount,self.aXtraGrid.size)) + EndOfPrdvPP_all = np.zeros((self.StateCount,self.aXtraGrid.size)) + for j in range(self.StateCount): + if np.any(np.logical_and(self.possible_transitions[:,j],which_states)): # only consider a future state if one of the relevant states could transition to it + EndOfPrdvP_all[j,:] = self.EndOfPrdvPfunc_list[j](aGrid) + if self.CubicBool: # Add conditional end-of-period (marginal) marginal value to the arrays + EndOfPrdvPP_all[j,:] = self.EndOfPrdvPfunc_list[j].derivative(aGrid) + # Weight conditional marginal (marginal) values by transition probs + # to get unconditional marginal (marginal) value at each gridpoint. + EndOfPrdvP_temp = np.dot(self.MrkvArray,EndOfPrdvP_all) + EndOfPrdvP[which_states,:] = EndOfPrdvP_temp[which_states,:] # only take the states for which this asset minimum applies + if self.CubicBool: + EndOfPrdvPP_temp = np.dot(self.MrkvArray,EndOfPrdvPP_all) + EndOfPrdvPP[which_states,:] = EndOfPrdvPP_temp[which_states,:] + + # Store the results as attributes of self, scaling end of period marginal value by survival probability from each current state + LivPrb_tiled = np.tile(np.reshape(self.LivPrb,(self.StateCount,1)),(1,self.aXtraGrid.size)) + self.EndOfPrdvP = LivPrb_tiled*EndOfPrdvP + if self.CubicBool: + self.EndOfPrdvPP = LivPrb_tiled*EndOfPrdvPP + + def calcHumWealthAndBoundingMPCs(self): + ''' + Calculates human wealth and the maximum and minimum MPC for each current + period state, then stores them as attributes of self for use by other methods. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + # Upper bound on MPC at lower m-bound + WorstIncPrb_array = self.BoroCnstDependency*np.tile(np.reshape(self.WorstIncPrbAll, + (1,self.StateCount)),(self.StateCount,1)) + temp_array = self.MrkvArray*WorstIncPrb_array + WorstIncPrbNow = np.sum(temp_array,axis=1) # Probability of getting the "worst" income shock and transition from each current state + ExMPCmaxNext = (np.dot(temp_array,self.Rfree_list**(1.0-self.CRRA)* + self.solution_next.MPCmax**(-self.CRRA))/WorstIncPrbNow)**\ + (-1.0/self.CRRA) + DiscFacEff_temp = self.DiscFac*self.LivPrb + self.MPCmaxNow = 1.0/(1.0 + ((DiscFacEff_temp*WorstIncPrbNow)** + (1.0/self.CRRA))/ExMPCmaxNext) + self.MPCmaxEff = self.MPCmaxNow + self.MPCmaxEff[self.BoroCnstNat_list < self.mNrmMin_list] = 1.0 + # State-conditional PDV of human wealth + hNrmPlusIncNext = self.ExIncNextAll + self.solution_next.hNrm + self.hNrmNow = np.dot(self.MrkvArray,(self.PermGroFac_list/self.Rfree_list)* + hNrmPlusIncNext) + # Lower bound on MPC as m gets arbitrarily large + temp = (DiscFacEff_temp*np.dot(self.MrkvArray,self.solution_next.MPCmin** + (-self.CRRA)*self.Rfree_list**(1.0-self.CRRA)))**(1.0/self.CRRA) + self.MPCminNow = 1.0/(1.0 + temp) + + def makeSolution(self,cNrm,mNrm): + ''' + Construct an object representing the solution to this period's problem. + + Parameters + ---------- + cNrm : np.array + Array of normalized consumption values for interpolation. Each row + corresponds to a Markov state for this period. + mNrm : np.array + Array of normalized market resource values for interpolation. Each + row corresponds to a Markov state for this period. + + Returns + ------- + solution : ConsumerSolution + The solution to the single period consumption-saving problem. Includes + a consumption function cFunc (using cubic or linear splines), a marg- + inal value function vPfunc, a minimum acceptable level of normalized + market resources mNrmMin, normalized human wealth hNrm, and bounding + MPCs MPCmin and MPCmax. It might also have a value function vFunc + and marginal marginal value function vPPfunc. All of these attributes + are lists or arrays, with elements corresponding to the current + Markov state. E.g. solution.cFunc[0] is the consumption function + when in the i=0 Markov state this period. + ''' + solution = ConsumerSolution() # An empty solution to which we'll add state-conditional solutions + # Calculate the MPC at each market resource gridpoint in each state (if desired) + if self.CubicBool: + dcda = self.EndOfPrdvPP/self.uPP(np.array(self.cNrmNow)) + MPC = dcda/(dcda+1.0) + self.MPC_temp = np.hstack((np.reshape(self.MPCmaxNow,(self.StateCount,1)),MPC)) + interpfunc = self.makeCubiccFunc + else: + interpfunc = self.makeLinearcFunc + + # Loop through each current period state and add its solution to the overall solution + for i in range(self.StateCount): + # Set current-period-conditional human wealth and MPC bounds + self.hNrmNow_j = self.hNrmNow[i] + self.MPCminNow_j = self.MPCminNow[i] + if self.CubicBool: + self.MPC_temp_j = self.MPC_temp[i,:] + + # Construct the consumption function by combining the constrained and unconstrained portions + self.cFuncNowCnst = LinearInterp([self.mNrmMin_list[i], self.mNrmMin_list[i]+1.0], + [0.0,1.0]) + cFuncNowUnc = interpfunc(mNrm[i,:],cNrm[i,:]) + cFuncNow = LowerEnvelope(cFuncNowUnc,self.cFuncNowCnst) + + # Make the marginal value function and pack up the current-state-conditional solution + vPfuncNow = MargValueFunc(cFuncNow,self.CRRA) + solution_cond = ConsumerSolution(cFunc=cFuncNow, vPfunc=vPfuncNow, + mNrmMin=self.mNrmMinNow) + if self.CubicBool: # Add the state-conditional marginal marginal value function (if desired) + solution_cond = self.addvPPfunc(solution_cond) + + # Add the current-state-conditional solution to the overall period solution + solution.appendSolution(solution_cond) + + # Add the lower bounds of market resources, MPC limits, human resources, + # and the value functions to the overall solution + solution.mNrmMin = self.mNrmMin_list + solution = self.addMPCandHumanWealth(solution) + if self.vFuncBool: + vFuncNow = self.makevFunc(solution) + solution.vFunc = vFuncNow + + # Return the overall solution to this period + return solution + + + def makeLinearcFunc(self,mNrm,cNrm): + ''' + Make a linear interpolation to represent the (unconstrained) consumption + function conditional on the current period state. + + Parameters + ---------- + mNrm : np.array + Array of normalized market resource values for interpolation. + cNrm : np.array + Array of normalized consumption values for interpolation. + + Returns + ------- + cFuncUnc: an instance of HARK.interpolation.LinearInterp + ''' + cFuncUnc = LinearInterp(mNrm,cNrm,self.MPCminNow_j*self.hNrmNow_j,self.MPCminNow_j) + return cFuncUnc + + + def makeCubiccFunc(self,mNrm,cNrm): + ''' + Make a cubic interpolation to represent the (unconstrained) consumption + function conditional on the current period state. + + Parameters + ---------- + mNrm : np.array + Array of normalized market resource values for interpolation. + cNrm : np.array + Array of normalized consumption values for interpolation. + + Returns + ------- + cFuncUnc: an instance of HARK.interpolation.CubicInterp + ''' + cFuncUnc = CubicInterp(mNrm,cNrm,self.MPC_temp_j,self.MPCminNow_j*self.hNrmNow_j, + self.MPCminNow_j) + return cFuncUnc + + def makevFunc(self,solution): + ''' + Construct the value function for each current state. + + Parameters + ---------- + solution : ConsumerSolution + The solution to the single period consumption-saving problem. Must + have a consumption function cFunc (using cubic or linear splines) as + a list with elements corresponding to the current Markov state. E.g. + solution.cFunc[0] is the consumption function when in the i=0 Markov + state this period. + + Returns + ------- + vFuncNow : [ValueFunc] + A list of value functions (defined over normalized market resources + m) for each current period Markov state. + ''' + vFuncNow = [] # Initialize an empty list of value functions + # Loop over each current period state and construct the value function + for i in range(self.StateCount): + # Make state-conditional grids of market resources and consumption + mNrmMin = self.mNrmMin_list[i] + mGrid = mNrmMin + self.aXtraGrid + cGrid = solution.cFunc[i](mGrid) + aGrid = mGrid - cGrid + + # Calculate end-of-period value at each gridpoint + EndOfPrdv_all = np.zeros((self.StateCount,self.aXtraGrid.size)) + for j in range(self.StateCount): + if self.possible_transitions[i,j]: + EndOfPrdv_all[j,:] = self.EndOfPrdvFunc_list[j](aGrid) + EndOfPrdv = np.dot(self.MrkvArray[i,:],EndOfPrdv_all) + + # Calculate (normalized) value and marginal value at each gridpoint + vNrmNow = self.u(cGrid) + EndOfPrdv + vPnow = self.uP(cGrid) + + # Make a "decurved" value function with the inverse utility function + vNvrs = self.uinv(vNrmNow) # value transformed through inverse utility + vNvrsP = vPnow*self.uinvP(vNrmNow) + mNrm_temp = np.insert(mGrid,0,mNrmMin) # add the lower bound + vNvrs = np.insert(vNvrs,0,0.0) + vNvrsP = np.insert(vNvrsP,0,self.MPCmaxEff[i]**(-self.CRRA/(1.0-self.CRRA))) + MPCminNvrs = self.MPCminNow[i]**(-self.CRRA/(1.0-self.CRRA)) + vNvrsFunc_i = CubicInterp(mNrm_temp,vNvrs,vNvrsP,MPCminNvrs*self.hNrmNow[i],MPCminNvrs) + + # "Recurve" the decurved value function and add it to the list + vFunc_i = ValueFunc(vNvrsFunc_i,self.CRRA) + vFuncNow.append(vFunc_i) + return vFuncNow + + +def _solveConsMarkov(solution_next,IncomeDstn,LivPrb,DiscFac,CRRA,Rfree,PermGroFac, + MrkvArray,BoroCnstArt,aXtraGrid,vFuncBool,CubicBool): + ''' + Solves a single period consumption-saving problem with risky income and + stochastic transitions between discrete states, in a Markov fashion. Has + identical inputs as solveConsIndShock, except for a discrete + Markov transitionrule MrkvArray. Markov states can differ in their interest + factor, permanent growth factor, and income distribution, so the inputs Rfree, + PermGroFac, and IncomeDstn are arrays or lists specifying those values in each + (succeeding) Markov state. + + Parameters + ---------- + solution_next : ConsumerSolution + The solution to next period's one period problem. + IncomeDstn_list : [[np.array]] + A length N list of income distributions in each succeeding Markov + state. Each income distribution contains three arrays of floats, + representing a discrete approximation to the income process at the + beginning of the succeeding period. Order: event probabilities, + permanent shocks, transitory shocks. + LivPrb : float + Survival probability; likelihood of being alive at the beginning of + the succeeding period. + DiscFac : float + Intertemporal discount factor for future utility. + CRRA : float + Coefficient of relative risk aversion. + Rfree_list : np.array + Risk free interest factor on end-of-period assets for each Markov + state in the succeeding period. + PermGroGac_list : float + Expected permanent income growth factor at the end of this period + for each Markov state in the succeeding period. + MrkvArray : numpy.array + An NxN array representing a Markov transition matrix between discrete + states. The i,j-th element of MrkvArray is the probability of + moving from state i in period t to state j in period t+1. + BoroCnstArt: float or None + Borrowing constraint for the minimum allowable assets to end the + period with. If it is less than the natural borrowing constraint, + then it is irrelevant; BoroCnstArt=None indicates no artificial bor- + rowing constraint. + aXtraGrid: np.array + Array of "extra" end-of-period asset values-- assets above the + absolute minimum acceptable level. + vFuncBool: boolean + An indicator for whether the value function should be computed and + included in the reported solution. + CubicBool: boolean + An indicator for whether the solver should use cubic or linear inter- + polation. + + Returns + ------- + solution : ConsumerSolution + The solution to the single period consumption-saving problem. Includes + a consumption function cFunc (using cubic or linear splines), a marg- + inal value function vPfunc, a minimum acceptable level of normalized + market resources mNrmMin, normalized human wealth hNrm, and bounding + MPCs MPCmin and MPCmax. It might also have a value function vFunc + and marginal marginal value function vPPfunc. All of these attributes + are lists or arrays, with elements corresponding to the current + Markov state. E.g. solution.cFunc[0] is the consumption function + when in the i=0 Markov state this period. + ''' + solver = ConsMarkovSolver(solution_next,IncomeDstn,LivPrb,DiscFac,CRRA,Rfree, + PermGroFac,MrkvArray,BoroCnstArt,aXtraGrid,vFuncBool,CubicBool) + solution_now = solver.solve() + return solution_now + + + +#################################################################################################### +#################################################################################################### + +class MarkovConsumerType(IndShockConsumerType): + ''' + An agent in the Markov consumption-saving model. His problem is defined by a sequence + of income distributions, survival probabilities, discount factors, and permanent + income growth rates, as well as time invariant values for risk aversion, the + interest rate, the grid of end-of-period assets, and how he is borrowing constrained. + ''' + time_vary_ = IndShockConsumerType.time_vary_ + ['MrkvArray'] + shock_vars_ = IndShockConsumerType.shock_vars_ + ['MrkvNow'] + + def __init__(self,cycles=1,time_flow=True,**kwds): + IndShockConsumerType.__init__(self,cycles=1,time_flow=True,**kwds) + self.solveOnePeriod = _solveConsMarkov + self.poststate_vars += ['MrkvNow'] + if not hasattr(self, 'global_markov'): + self.global_markov = False + + def checkMarkovInputs(self): + ''' + Many parameters used by MarkovConsumerType are arrays. Make sure those arrays are the + right shape. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + StateCount = self.MrkvArray[0].shape[0] + + # Check that arrays are the right shape + if not isinstance(self.Rfree, np.ndarray) or self.Rfree.shape != (StateCount, ): + raise ValueError('Rfree not the right shape, it should an array of Rfree of all the states.') + + # Check that arrays in lists are the right shape + for MrkvArray_t in self.MrkvArray: + if not isinstance(MrkvArray_t, np.ndarray) or MrkvArray_t.shape != (StateCount, StateCount): + raise ValueError('MrkvArray not the right shape, it should be of the size states*statres.') + for LivPrb_t in self.LivPrb: + if not isinstance(LivPrb_t, np.ndarray) or LivPrb_t.shape != (StateCount, ): + raise ValueError('Array in LivPrb is not the right shape, it should be an array of length equal to number of states') + for PermGroFac_t in self.PermGroFac: + if not isinstance(PermGroFac_t, np.ndarray) or PermGroFac_t.shape != (StateCount, ): + raise ValueError('Array in PermGroFac is not the right shape, it should be an array of length equal to number of states') + + # Now check the income distribution. + # Note IncomeDstn is (potentially) time-varying, so it is in time_vary. + # Therefore it is a list, and each element of that list responds to the income distribution + # at a particular point in time. Each income distribution at a point in time should itself + # be a list, with each element corresponding to the income distribution + # conditional on a particular Markov state. + # TODO: should this be a numpy array too? + for IncomeDstn_t in self.IncomeDstn: + if not isinstance(IncomeDstn_t, list) or len(IncomeDstn_t) != StateCount: + raise ValueError('List in IncomeDstn is not the right length, it should be length equal to number of states') + + def preSolve(self): + """ + Check to make sure that the inputs that are specific to MarkovConsumerType + are of the right shape (if arrays) or length (if lists). + + Parameters + ---------- + None + + Returns + ------- + None + """ + AgentType.preSolve(self) + self.checkMarkovInputs() + + def updateSolutionTerminal(self): + ''' + Update the terminal period solution. This method should be run when a + new AgentType is created or when CRRA changes. + + Parameters + ---------- + none + + Returns + ------- + none + ''' + IndShockConsumerType.updateSolutionTerminal(self) + + # Make replicated terminal period solution: consume all resources, no human wealth, minimum m is 0 + StateCount = self.MrkvArray[0].shape[0] + self.solution_terminal.cFunc = StateCount*[self.cFunc_terminal_] + self.solution_terminal.vFunc = StateCount*[self.solution_terminal.vFunc] + self.solution_terminal.vPfunc = StateCount*[self.solution_terminal.vPfunc] + self.solution_terminal.vPPfunc = StateCount*[self.solution_terminal.vPPfunc] + self.solution_terminal.mNrmMin = np.zeros(StateCount) + self.solution_terminal.hRto = np.zeros(StateCount) + self.solution_terminal.MPCmax = np.ones(StateCount) + self.solution_terminal.MPCmin = np.ones(StateCount) + + def initializeSim(self): + IndShockConsumerType.initializeSim(self) + if self.global_markov: #Need to initialize markov state to be the same for all agents + base_draw = Uniform().draw(1,seed=self.RNG.randint(0,2**31-1)) + Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit)) + self.MrkvNow = np.ones(self.AgentCount)*np.searchsorted(Cutoffs,base_draw).astype(int) + self.MrkvNow = self.MrkvNow.astype(int) + + def simDeath(self): + ''' + Determines which agents die this period and must be replaced. Uses the sequence in LivPrb + to determine survival probabilities for each agent. + + Parameters + ---------- + None + + Returns + ------- + which_agents : np.array(bool) + Boolean array of size AgentCount indicating which agents die. + ''' + # Determine who dies + LivPrb = np.array(self.LivPrb)[self.t_cycle-1,self.MrkvNow] # Time has already advanced, so look back one + DiePrb = 1.0 - LivPrb + DeathShks = Uniform().draw(N=self.AgentCount,seed=self.RNG.randint(0,2**31-1)) + which_agents = DeathShks < DiePrb + if self.T_age is not None: # Kill agents that have lived for too many periods + too_old = self.t_age >= self.T_age + which_agents = np.logical_or(which_agents,too_old) + return which_agents + + def simBirth(self,which_agents): + ''' + Makes new Markov consumer by drawing initial normalized assets, permanent income levels, and + discrete states. Calls IndShockConsumerType.simBirth, then draws from initial Markov distribution. + + Parameters + ---------- + which_agents : np.array(Bool) + Boolean array of size self.AgentCount indicating which agents should be "born". + + Returns + ------- + None + ''' + IndShockConsumerType.simBirth(self,which_agents) # Get initial assets and permanent income + if not self.global_markov: #Markov state is not changed if it is set at the global level + N = np.sum(which_agents) + base_draws = Uniform().draw(N,seed=self.RNG.randint(0,2**31-1)) + Cutoffs = np.cumsum(np.array(self.MrkvPrbsInit)) + self.MrkvNow[which_agents] = np.searchsorted(Cutoffs,base_draws).astype(int) + + + def getMarkovStates(self): + ''' + Draw new Markov states for each agent in the simulated population, using + the attribute MrkvArray to determine transition probabilities. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + # Draw random numbers that will be used to determine the next Markov state + if self.global_markov: + base_draws = np.ones(self.AgentCount)*Uniform().draw(1,seed=self.RNG.randint(0,2**31-1)) + else: + base_draws = Uniform().draw(self.AgentCount,seed=self.RNG.randint(0,2**31-1)) + dont_change = self.t_age == 0 # Don't change Markov state for those who were just born (unless global_markov) + if self.t_sim == 0: # Respect initial distribution of Markov states + dont_change[:] = True + + # Determine which agents are in which states right now + J = self.MrkvArray[0].shape[0] + MrkvPrev = self.MrkvNow + MrkvNow = np.zeros(self.AgentCount,dtype=int) + MrkvBoolArray = np.zeros((J,self.AgentCount)) + for j in range(J): + MrkvBoolArray[j,:] = MrkvPrev == j + + # Draw new Markov states for each agent + for t in range(self.T_cycle): + Cutoffs = np.cumsum(self.MrkvArray[t],axis=1) + right_age = self.t_cycle == t + for j in range(J): + these = np.logical_and(right_age, MrkvBoolArray[j,:]) + MrkvNow[these] = np.searchsorted(Cutoffs[j,:],base_draws[these]).astype(int) + if not self.global_markov: + MrkvNow[dont_change] = MrkvPrev[dont_change] + self.MrkvNow = MrkvNow.astype(int) + + + def getShocks(self): + ''' + Gets new Markov states and permanent and transitory income shocks for this period. Samples + from IncomeDstn for each period-state in the cycle. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + self.getMarkovStates() + MrkvNow = self.MrkvNow + + # Now get income shocks for each consumer, by cycle-time and discrete state + PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays + TranShkNow = np.zeros(self.AgentCount) + for t in range(self.T_cycle): + for j in range(self.MrkvArray[t].shape[0]): + these = np.logical_and(t == self.t_cycle, j == MrkvNow) + N = np.sum(these) + if N > 0: + IncomeDstnNow = self.IncomeDstn[t-1][j] # set current income distribution + PermGroFacNow = self.PermGroFac[t-1][j] # and permanent growth factor + Indices = np.arange(IncomeDstnNow.pmf.size) # just a list of integers + # Get random draws of income shocks from the discrete distribution + EventDraws = IncomeDstnNow.draw_events( + N, + seed=self.RNG.randint(0,2**31-1)) + PermShkNow[these] = IncomeDstnNow.X[0][EventDraws]*PermGroFacNow # permanent "shock" includes expected growth + TranShkNow[these] = IncomeDstnNow.X[1][EventDraws] + newborn = self.t_age == 0 + PermShkNow[newborn] = 1.0 + TranShkNow[newborn] = 1.0 + self.PermShkNow = PermShkNow + self.TranShkNow = TranShkNow + + def readShocks(self): + ''' + A slight modification of AgentType.readShocks that makes sure that MrkvNow is int, not float. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + IndShockConsumerType.readShocks(self) + self.MrkvNow = self.MrkvNow.astype(int) + + def getRfree(self): + ''' + Returns an array of size self.AgentCount with interest factor that varies with discrete state. + + Parameters + ---------- + None + + Returns + ------- + RfreeNow : np.array + Array of size self.AgentCount with risk free interest rate for each agent. + ''' + RfreeNow = self.Rfree[self.MrkvNow] + return RfreeNow + + def getControls(self): + ''' + Calculates consumption for each consumer of this type using the consumption functions. + + Parameters + ---------- + None + + Returns + ------- + None + ''' + cNrmNow = np.zeros(self.AgentCount) + np.nan + MPCnow = np.zeros(self.AgentCount) + np.nan + J = self.MrkvArray[0].shape[0] + + MrkvBoolArray = np.zeros((J,self.AgentCount), dtype=bool) + for j in range(J): + MrkvBoolArray[j,:] = j == self.MrkvNow + + for t in range(self.T_cycle): + right_t = t == self.t_cycle + for j in range(J): + these = np.logical_and(right_t, MrkvBoolArray[j,:]) + cNrmNow[these], MPCnow[these] = self.solution[t].cFunc[j].eval_with_derivative(self.mNrmNow[these]) + self.cNrmNow = cNrmNow + self.MPCnow = MPCnow + + + def calcBoundingValues(self): + ''' + Calculate human wealth plus minimum and maximum MPC in an infinite + horizon model with only one period repeated indefinitely. Store results + as attributes of self. Human wealth is the present discounted value of + expected future income after receiving income this period, ignoring mort- + ality. The maximum MPC is the limit of the MPC as m --> mNrmMin. The + minimum MPC is the limit of the MPC as m --> infty. Results are all + np.array with elements corresponding to each Markov state. + + NOT YET IMPLEMENTED FOR THIS CLASS + + Parameters + ---------- + None + + Returns + ------- + None + ''' + raise NotImplementedError() + + def makeEulerErrorFunc(self,mMax=100,approx_inc_dstn=True): + ''' + Creates a "normalized Euler error" function for this instance, mapping + from market resources to "consumption error per dollar of consumption." + Stores result in attribute eulerErrorFunc as an interpolated function. + Has option to use approximate income distribution stored in self.IncomeDstn + or to use a (temporary) very dense approximation. + + NOT YET IMPLEMENTED FOR THIS CLASS + + Parameters + ---------- + mMax : float + Maximum normalized market resources for the Euler error function. + approx_inc_dstn : Boolean + Indicator for whether to use the approximate discrete income distri- + bution stored in self.IncomeDstn[0], or to use a very accurate + discrete approximation instead. When True, uses approximation in + IncomeDstn; when False, makes and uses a very dense approximation. + + Returns + ------- + None + ''' + raise NotImplementedError() diff --git a/HARK/ConsumptionSaving/ConsPortfolioModel.py b/HARK/ConsumptionSaving/ConsPortfolioModel.py index 7bdbf42da..60b926b91 100644 --- a/HARK/ConsumptionSaving/ConsPortfolioModel.py +++ b/HARK/ConsumptionSaving/ConsPortfolioModel.py @@ -22,7 +22,7 @@ ValueFunc2D, # For representing 2D value function MargValueFunc2D # For representing 2D marginal value function ) -from HARK.distribution import approxLognormal, combineIndepDstns +from HARK.distribution import combineIndepDstns from HARK.distribution import Lognormal, Bernoulli # Random draws for simulating agents from HARK.interpolation import( LinearInterp, # Piecewise linear interpolation @@ -242,7 +242,7 @@ def updateRiskyDstn(self): RiskyVar = self.RiskyStd[t] ** 2 mu = np.log(self.RiskyAvg[t] / (np.sqrt(1. + RiskyVar / RiskyAvgSqrd))) sigma = np.sqrt(np.log(1. + RiskyVar / RiskyAvgSqrd)) - RiskyDstn.append(approxLognormal(self.RiskyCount, mu=mu, sigma=sigma)) + RiskyDstn.append(Lognormal(mu=mu, sigma=sigma).approx(self.RiskyCount)) self.RiskyDstn = RiskyDstn self.addToTimeVary('RiskyDstn') @@ -253,7 +253,7 @@ def updateRiskyDstn(self): RiskyVar = self.RiskyStd ** 2 mu = np.log(self.RiskyAvg / (np.sqrt(1. + RiskyVar / RiskyAvgSqrd))) sigma = np.sqrt(np.log(1. + RiskyVar / RiskyAvgSqrd)) - self.RiskyDstn = approxLognormal(self.RiskyCount, mu=mu, sigma=sigma) + self.RiskyDstn = Lognormal(mu=mu, sigma=sigma).approx(self.RiskyCount) self.addToTimeInv('RiskyDstn') diff --git a/HARK/ConsumptionSaving/ConsPrefShockModel.py b/HARK/ConsumptionSaving/ConsPrefShockModel.py index a8e1b3925..62c67c3f2 100644 --- a/HARK/ConsumptionSaving/ConsPrefShockModel.py +++ b/HARK/ConsumptionSaving/ConsPrefShockModel.py @@ -11,7 +11,7 @@ from builtins import str from builtins import range import numpy as np -from HARK.distribution import approxMeanOneLognormal +from HARK.distribution import MeanOneLogNormal from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType, ConsumerSolution, ConsIndShockSolver, \ ValueFunc, MargValueFunc, KinkedRconsumerType, ConsKinkedRsolver, \ init_idiosyncratic_shocks, init_kinked_R @@ -111,8 +111,11 @@ def updatePrefShockProcess(self): PrefShkDstn = [] # discrete distributions of preference shocks for t in range(len(self.PrefShkStd)): PrefShkStd = self.PrefShkStd[t] - PrefShkDstn.append(approxMeanOneLognormal(N=self.PrefShkCount, - sigma=PrefShkStd,tail_N=self.PrefShk_tail_N)) + PrefShkDstn.append( + MeanOneLogNormal( + sigma=PrefShkStd + ).approx(N=self.PrefShkCount, + tail_N=self.PrefShk_tail_N)) # Store the preference shocks in self (time-varying) and restore time flow self.PrefShkDstn = PrefShkDstn @@ -137,8 +140,9 @@ def getShocks(self): N = np.sum(these) if N > 0: PrefShkNow[these] = self.RNG.permutation( - approxMeanOneLognormal(N, - sigma=self.PrefShkStd[t]).X) + MeanOneLogNormal( + sigma=self.PrefShkStd[t] + ).approx(N).X) self.PrefShkNow = PrefShkNow def getControls(self): diff --git a/HARK/distribution.py b/HARK/distribution.py index ca722bbc9..cf694cb4d 100644 --- a/HARK/distribution.py +++ b/HARK/distribution.py @@ -72,6 +72,90 @@ def draw(self, N, seed=0): size=N)) return draws + def approx(self, N, tail_N=0, tail_bound=[0.02,0.98], tail_order=np.e): + ''' + Construct a discrete approximation to a lognormal distribution with underlying + normal distribution N(mu,sigma). Makes an equiprobable distribution by + default, but user can optionally request augmented tails with exponentially + sized point masses. This can improve solution accuracy in some models. + + Parameters + ---------- + N: int + Number of discrete points in the "main part" of the approximation. + tail_N: int + Number of points in each "tail part" of the approximation; 0 = no tail. + tail_bound: [float] + CDF boundaries of the tails vs main portion; tail_bound[0] is the lower + tail bound, tail_bound[1] is the upper tail bound. Inoperative when + tail_N = 0. Can make "one tailed" approximations with 0.0 or 1.0. + tail_order: float + Factor by which consecutive point masses in a "tail part" differ in + probability. Should be >= 1 for sensible spacing. + + Returns + ------- + d : DiscreteDistribution + Probability associated with each point in array of discrete + points for discrete probability mass function. + ''' + # Find the CDF boundaries of each segment + if self.sigma > 0.0: + if tail_N > 0: + lo_cut = tail_bound[0] + hi_cut = tail_bound[1] + else: + lo_cut = 0.0 + hi_cut = 1.0 + inner_size = hi_cut - lo_cut + inner_CDF_vals = [lo_cut + x*N**(-1.0)*inner_size for x in range(1, N)] + if inner_size < 1.0: + scale = 1.0/tail_order + mag = (1.0-scale**tail_N)/(1.0-scale) + lower_CDF_vals = [0.0] + if lo_cut > 0.0: + for x in range(tail_N-1,-1,-1): + lower_CDF_vals.append(lower_CDF_vals[-1] + lo_cut*scale**x/mag) + upper_CDF_vals = [hi_cut] + if hi_cut < 1.0: + for x in range(tail_N): + upper_CDF_vals.append(upper_CDF_vals[-1] + (1.0-hi_cut)*scale**x/mag) + CDF_vals = lower_CDF_vals + inner_CDF_vals + upper_CDF_vals + temp_cutoffs = list(stats.lognorm.ppf(CDF_vals[1:-1], s=self.sigma, loc=0, scale=np.exp(self.mu))) + cutoffs = [0] + temp_cutoffs + [np.inf] + CDF_vals = np.array(CDF_vals) + + # Construct the discrete approximation by finding the average value within each segment. + # This codeblock ignores warnings because it throws a "divide by zero encountered in log" + # warning due to computing erf(infty) at the tail boundary. This is irrelevant and + # apparently freaks new users out. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + K = CDF_vals.size-1 # number of points in approximation + pmf = CDF_vals[1:(K+1)] - CDF_vals[0:K] + X = np.zeros(K) + for i in range(K): + zBot = cutoffs[i] + zTop = cutoffs[i+1] + tempBot = (self.mu+self.sigma**2-np.log(zBot))/(np.sqrt(2)*self.sigma) + tempTop = (self.mu+self.sigma**2-np.log(zTop))/(np.sqrt(2)*self.sigma) + if tempBot <= 4: + X[i] = -0.5*np.exp(self.mu+(self.sigma**2)*0.5)*(erf(tempTop) - erf(tempBot))/pmf[i] + else: + X[i] = -0.5*np.exp(self.mu+(self.sigma**2)*0.5)*(erfc(tempBot) - erfc(tempTop))/pmf[i] + + else: + pmf = np.ones(N)/N + X = np.exp(self.mu)*np.ones(N) + return DiscreteDistribution(pmf, X) + +class MeanOneLogNormal(Lognormal): + + def __init__(self, sigma=1.0): + mu = -0.5*sigma**2 + super().__init__(mu=mu, sigma=sigma) + + class Normal(): """ A Normal distribution. @@ -277,34 +361,6 @@ def approx(self, N): X = center + width*np.linspace(-(N-1.0)/2.0,(N-1.0)/2.0,N)/(N/2.0) return DiscreteDistribution(pmf,X) -def drawMeanOneLognormal(N, sigma=1.0, seed=0): - ''' - Generate arrays of mean one lognormal draws. The sigma input can be a number - or list-like. If a number, output is a length N array of draws from the - lognormal distribution with standard deviation sigma. If a list, output is - a length T list whose t-th entry is a length N array of draws from the - lognormal with standard deviation sigma[t]. - - Parameters - ---------- - N : int - Number of draws in each row. - sigma : float or [float] - One or more standard deviations. Number of elements T in sigma - determines number of rows of output. - seed : int - Seed for random number generator. - - Returns: - ------------ - draws : np.array or [np.array] - T-length list of arrays of mean one lognormal draws each of size N, or - a single array of size N (if sigma is a scalar). - ''' - mu = -0.5*sigma**2 - - return Lognormal(mu,sigma).draw(N,seed=seed) - ### DISCRETE DISTRIBUTIONS @@ -463,112 +519,6 @@ def drawDiscrete(self, N,X=None,exact_match=False,seed=0): return draws -def approxLognormal(N, mu=0.0, sigma=1.0, tail_N=0, tail_bound=[0.02,0.98], tail_order=np.e): - ''' - Construct a discrete approximation to a lognormal distribution with underlying - normal distribution N(mu,sigma). Makes an equiprobable distribution by - default, but user can optionally request augmented tails with exponentially - sized point masses. This can improve solution accuracy in some models. - - Parameters - ---------- - N: int - Number of discrete points in the "main part" of the approximation. - mu: float - Mean of underlying normal distribution. - sigma: float - Standard deviation of underlying normal distribution. - tail_N: int - Number of points in each "tail part" of the approximation; 0 = no tail. - tail_bound: [float] - CDF boundaries of the tails vs main portion; tail_bound[0] is the lower - tail bound, tail_bound[1] is the upper tail bound. Inoperative when - tail_N = 0. Can make "one tailed" approximations with 0.0 or 1.0. - tail_order: float - Factor by which consecutive point masses in a "tail part" differ in - probability. Should be >= 1 for sensible spacing. - - Returns - ------- - d : DiscreteDistribution - Probability associated with each point in array of discrete - points for discrete probability mass function. - ''' - # Find the CDF boundaries of each segment - if sigma > 0.0: - if tail_N > 0: - lo_cut = tail_bound[0] - hi_cut = tail_bound[1] - else: - lo_cut = 0.0 - hi_cut = 1.0 - inner_size = hi_cut - lo_cut - inner_CDF_vals = [lo_cut + x*N**(-1.0)*inner_size for x in range(1, N)] - if inner_size < 1.0: - scale = 1.0/tail_order - mag = (1.0-scale**tail_N)/(1.0-scale) - lower_CDF_vals = [0.0] - if lo_cut > 0.0: - for x in range(tail_N-1,-1,-1): - lower_CDF_vals.append(lower_CDF_vals[-1] + lo_cut*scale**x/mag) - upper_CDF_vals = [hi_cut] - if hi_cut < 1.0: - for x in range(tail_N): - upper_CDF_vals.append(upper_CDF_vals[-1] + (1.0-hi_cut)*scale**x/mag) - CDF_vals = lower_CDF_vals + inner_CDF_vals + upper_CDF_vals - temp_cutoffs = list(stats.lognorm.ppf(CDF_vals[1:-1], s=sigma, loc=0, scale=np.exp(mu))) - cutoffs = [0] + temp_cutoffs + [np.inf] - CDF_vals = np.array(CDF_vals) - - # Construct the discrete approximation by finding the average value within each segment. - # This codeblock ignores warnings because it throws a "divide by zero encountered in log" - # warning due to computing erf(infty) at the tail boundary. This is irrelevant and - # apparently freaks new users out. - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - K = CDF_vals.size-1 # number of points in approximation - pmf = CDF_vals[1:(K+1)] - CDF_vals[0:K] - X = np.zeros(K) - for i in range(K): - zBot = cutoffs[i] - zTop = cutoffs[i+1] - tempBot = (mu+sigma**2-np.log(zBot))/(np.sqrt(2)*sigma) - tempTop = (mu+sigma**2-np.log(zTop))/(np.sqrt(2)*sigma) - if tempBot <= 4: - X[i] = -0.5*np.exp(mu+(sigma**2)*0.5)*(erf(tempTop) - erf(tempBot))/pmf[i] - else: - X[i] = -0.5*np.exp(mu+(sigma**2)*0.5)*(erfc(tempBot) - erfc(tempTop))/pmf[i] - - else: - pmf = np.ones(N)/N - X = np.exp(mu)*np.ones(N) - return DiscreteDistribution(pmf, X) - -@memoize -def approxMeanOneLognormal(N, sigma=1.0, **kwargs): - ''' - Calculate a discrete approximation to a mean one lognormal distribution. - Based on function approxLognormal; see that function's documentation for - further notes. - - Parameters - ---------- - N : int - Size of discrete space vector to be returned. - sigma : float - standard deviation associated with underlying normal probability distribution. - - Returns - ------- - d : DiscreteDistribution - Probability associated with each point in array of discrete - points for discrete probability mass function. - ''' - mu_adj = - 0.5*sigma**2; - dist = approxLognormal(N=N, mu=mu_adj, sigma=sigma, **kwargs) - return dist - - def approxLognormalGaussHermite(N, mu=0.0, sigma=1.0): d = Normal(mu, sigma).approx(N) return DiscreteDistribution(d.pmf, np.exp(d.X)) diff --git a/HARK/lifecycle.py b/HARK/lifecycle.py new file mode 100644 index 000000000..73016ad09 --- /dev/null +++ b/HARK/lifecycle.py @@ -0,0 +1,29 @@ +from HARK.ConsumptionSaving.ConsIndShockModel import \ + IndShockConsumerType, ConsIndShockSolverBasic +import HARK.ConsumptionSaving.ConsumerParameters as Params + +LifecycleExample = IndShockConsumerType( + **Params.init_lifecycle) +LifecycleExample.cycles = 1 +LifecycleExample.solve() + +# test the solution_terminal +assert( + LifecycleExample.solution[10].cFunc(2).tolist() == 2) + +print(LifecycleExample.solution[9].cFunc(1)) +assert(LifecycleExample.solution[9].cFunc(1) == 0.97769632) +self.assertAlmostEqual(LifecycleExample.solution[8].cFunc(1), + 0.96624445) +self.assertAlmostEqual(LifecycleExample.solution[7].cFunc(1), + 0.95691449) + +self.assertAlmostEqual( + LifecycleExample.solution[0].cFunc(1).tolist(), + 0.87362789) +self.assertAlmostEqual( + LifecycleExample.solution[1].cFunc(1).tolist(), + 0.9081621) +self.assertAlmostEqual( + LifecycleExample.solution[2].cFunc(1).tolist(), + 0.9563899) diff --git a/HARK/policy.py b/HARK/policy.py new file mode 100644 index 000000000..2c2fe720e --- /dev/null +++ b/HARK/policy.py @@ -0,0 +1,11 @@ +from core import HARKobject + + +class Solution(HARKobject): + ''' + A superclass for representing the "solution" to a single period problem in a + dynamic microeconomic model. + + NOTE: This can be deprecated now that HARKobject exists, but this requires + replacing each instance of Solution with HARKobject in the other modules. + ''' diff --git a/HARK/solver.py b/HARK/solver.py new file mode 100644 index 000000000..1f8223f91 --- /dev/null +++ b/HARK/solver.py @@ -0,0 +1,10 @@ +from core import HARKobject + +class Solution(HARKobject): + ''' + A superclass for representing the "solution" to a single period problem in a + dynamic microeconomic model. + + NOTE: This can be deprecated now that HARKobject exists, but this requires + replacing each instance of Solution with HARKobject in the other modules. + ''' diff --git a/HARK/tests/test_distribution.py b/HARK/tests/test_distribution.py index c8b27aa54..b0fecef20 100644 --- a/HARK/tests/test_distribution.py +++ b/HARK/tests/test_distribution.py @@ -26,7 +26,7 @@ class DistributionClassTests(unittest.TestCase): def test_drawMeanOneLognormal(self): self.assertEqual( - distribution.drawMeanOneLognormal(1)[0], + distribution.MeanOneLogNormal().draw(1)[0], 3.5397367004222002) def test_Lognormal(self): From b549fb8bca8758ee5afd5c0fc353c39e4eaf19d5 Mon Sep 17 00:00:00 2001 From: sb Date: Tue, 28 Apr 2020 13:38:33 -0400 Subject: [PATCH 3/3] remove emacs hash files, add class of files to gitignore --- .gitignore | 3 +- HARK/#core.py# | 1409 ------------------------------------------------ HARK/#scratch# | 0 3 files changed, 2 insertions(+), 1410 deletions(-) delete mode 100644 HARK/#core.py# delete mode 100644 HARK/#scratch# diff --git a/.gitignore b/.gitignore index f7720865b..8d6f08d5f 100644 --- a/.gitignore +++ b/.gitignore @@ -263,8 +263,9 @@ Web # Mac system files .DS_Store -# Emacs automatic backup files +# Emacs automatic backup and hash files *~ +\#*# # Vim swap files *.swp diff --git a/HARK/#core.py# b/HARK/#core.py# deleted file mode 100644 index b4b8c2a03..000000000 --- a/HARK/#core.py# +++ /dev/null @@ -1,1409 +0,0 @@ -''' -High-level functions and classes for solving a wide variety of economic models. -The "core" of HARK is a framework for "microeconomic" and "macroeconomic" -models. A micro model concerns the dynamic optimization problem for some type -of agents, where agents take the inputs to their problem as exogenous. A macro -model adds an additional layer, endogenizing some of the inputs to the micro -problem by finding a general equilibrium dynamic rule. -''' -from __future__ import print_function, division -from __future__ import absolute_import - -from builtins import str -from builtins import range -from builtins import object -import sys -import os -from distutils.dir_util import copy_tree -from .utilities import getArgNames, NullFunc -from copy import copy, deepcopy -import numpy as np -from time import time -from .parallel import multiThreadCommands, multiThreadCommandsFake - - -def distanceMetric(thing_A, thing_B): - ''' - A "universal distance" metric that can be used as a default in many settings. - - Parameters - ---------- - thing_A : object - A generic object. - thing_B : object - Another generic object. - - Returns: - ------------ - distance : float - The "distance" between thing_A and thing_B. - ''' - # Get the types of the two inputs - typeA = type(thing_A) - typeB = type(thing_B) - - if typeA is list and typeB is list: - lenA = len(thing_A) # If both inputs are lists, then the distance between - lenB = len(thing_B) # them is the maximum distance between corresponding - if lenA == lenB: # elements in the lists. If they differ in length, - distance_temp = [] # the distance is the difference in lengths. - for n in range(lenA): - distance_temp.append(distanceMetric(thing_A[n], thing_B[n])) - distance = max(distance_temp) - else: - distance = float(abs(lenA - lenB)) - # If both inputs are numbers, return their difference - elif isinstance(thing_A, (int, float)) and isinstance(thing_B, (int, float)): - distance = float(abs(thing_A - thing_B)) - # If both inputs are array-like, return the maximum absolute difference b/w - # corresponding elements (if same shape); return largest difference in dimensions - # if shapes do not align. - elif hasattr(thing_A, 'shape') and hasattr(thing_B, 'shape'): - if thing_A.shape == thing_B.shape: - distance = np.max(abs(thing_A - thing_B)) - else: - # Flatten arrays so they have the same dimensions - distance = np.max(abs(thing_A.flatten().shape[0] - thing_B.flatten().shape[0])) - # If none of the above cases, but the objects are of the same class, call - # the distance method of one on the other - elif thing_A.__class__.__name__ == thing_B.__class__.__name__: - if thing_A.__class__.__name__ == 'function': - distance = 0.0 - else: - distance = thing_A.distance(thing_B) - else: # Failsafe: the inputs are very far apart - distance = 1000.0 - return distance - - -class HARKobject(object): - ''' - A superclass for object classes in HARK. Comes with two useful methods: - a generic/universal distance method and an attribute assignment method. - ''' - def distance(self, other): - ''' - A generic distance method, which requires the existence of an attribute - called distance_criteria, giving a list of strings naming the attributes - to be considered by the distance metric. - - Parameters - ---------- - other : object - Another object to compare this instance to. - - Returns - ------- - (unnamed) : float - The distance between this object and another, using the "universal - distance" metric. - ''' - distance_list = [0.0] - for attr_name in self.distance_criteria: - try: - obj_A = getattr(self, attr_name) - obj_B = getattr(other, attr_name) - distance_list.append(distanceMetric(obj_A, obj_B)) - except AttributeError: - distance_list.append(1000.0) # if either object lacks attribute, they are not the same - return max(distance_list) - - def assignParameters(self, **kwds): - ''' - Assign an arbitrary number of attributes to this agent. - - Parameters - ---------- - **kwds : keyword arguments - Any number of keyword arguments of the form key=value. Each value - will be assigned to the attribute named in self. - - Returns - ------- - none - ''' - for key in kwds: - setattr(self, key, kwds[key]) - - def __call__(self, **kwds): - ''' - Assign an arbitrary number of attributes to this agent, as a convenience. - See assignParameters. - ''' - self.assignParameters(**kwds) - - def getAvg(self, varname, **kwds): - ''' - Calculates the average of an attribute of this instance. Returns NaN if no such attribute. - - Parameters - ---------- - varname : string - The name of the attribute whose average is to be calculated. This attribute must be an - np.array or other class compatible with np.mean. - - Returns - ------- - avg : float or np.array - The average of this attribute. Might be an array if the axis keyword is passed. - ''' - if hasattr(self, varname): - return np.mean(getattr(self, varname), **kwds) - else: - return np.nan - - -class Solution(HARKobject): - ''' - A superclass for representing the "solution" to a single period problem in a - dynamic microeconomic model. - - NOTE: This can be deprecated now that HARKobject exists, but this requires - replacing each instance of Solution with HARKobject in the other modules. - ''' - - -class AgentType(HARKobject): - ''' - A superclass for economic agents in the HARK framework. Each model should - specify its own subclass of AgentType, inheriting its methods and overwriting - as necessary. Critically, every subclass of AgentType should define class- - specific static values of the attributes time_vary and time_inv as lists of - strings. Each element of time_vary is the name of a field in AgentSubType - that varies over time in the model. Each element of time_inv is the name of - a field in AgentSubType that is constant over time in the model. - ''' - def __init__(self, solution_terminal=None, cycles=1, time_flow=True, pseudo_terminal=True, - tolerance=0.000001, seed=0, **kwds): - ''' - Initialize an instance of AgentType by setting attributes. - - Parameters - ---------- - solution_terminal : Solution - A representation of the solution to the terminal period problem of - this AgentType instance, or an initial guess of the solution if this - is an infinite horizon problem. - cycles : int - The number of times the sequence of periods is experienced by this - AgentType in their "lifetime". cycles=1 corresponds to a lifecycle - model, with a certain sequence of one period problems experienced - once before terminating. cycles=0 corresponds to an infinite horizon - model, with a sequence of one period problems repeating indefinitely. - time_flow : boolean - Whether time is currently "flowing" forward(True) or backward(False) for this - instance. Used to flip between solving (using backward iteration) - and simulating (etc). - pseudo_terminal : boolean - Indicates whether solution_terminal isn't actually part of the - solution to the problem (as a known solution to the terminal period - problem), but instead represents a "scrap value"-style termination. - When True, solution_terminal is not included in the solution; when - false, solution_terminal is the last element of the solution. - tolerance : float - Maximum acceptable "distance" between successive solutions to the - one period problem in an infinite horizon (cycles=0) model in order - for the solution to be considered as having "converged". Inoperative - when cycles>0. - seed : int - A seed for this instance's random number generator. - - Returns - ------- - None - ''' - if solution_terminal is None: - solution_terminal = NullFunc() - self.solution_terminal = solution_terminal # NOQA - self.cycles = cycles # NOQA - self.time_flow = time_flow # NOQA - self.pseudo_terminal = pseudo_terminal # NOQA - self.solveOnePeriod = NullFunc() # NOQA - self.tolerance = tolerance # NOQA - self.seed = seed # NOQA - self.track_vars = [] # NOQA - self.poststate_vars = [] # NOQA - self.read_shocks = False # NOQA - self.assignParameters(**kwds) # NOQA - self.resetRNG() # NOQA - - def timeReport(self): - ''' - Report to the user the direction that time is currently "flowing" for - this instance. Only exists as a reminder of how time_flow works. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - if self.time_flow: - print('Time varying objects are listed in ordinary chronological order.') - else: - print('Time varying objects are listed in reverse chronological order.') - - def timeFlip(self): - ''' - Reverse the flow of time for this instance. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - for name in self.time_vary: - self.__dict__[name].reverse() - self.time_flow = not self.time_flow - - def timeFwd(self): - ''' - Make time flow forward for this instance. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - if not self.time_flow: - self.timeFlip() - - def timeRev(self): - ''' - Make time flow backward for this instance. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - if self.time_flow: - self.timeFlip() - - def addToTimeVary(self, *params): - ''' - Adds any number of parameters to time_vary for this instance. - - Parameters - ---------- - params : string - Any number of strings naming attributes to be added to time_vary - - Returns - ------- - None - ''' - for param in params: - if param not in self.time_vary: - self.time_vary.append(param) - - def addToTimeInv(self, *params): - ''' - Adds any number of parameters to time_inv for this instance. - - Parameters - ---------- - params : string - Any number of strings naming attributes to be added to time_inv - - Returns - ------- - None - ''' - for param in params: - if param not in self.time_inv: - self.time_inv.append(param) - - def delFromTimeVary(self, *params): - ''' - Removes any number of parameters from time_vary for this instance. - - Parameters - ---------- - params : string - Any number of strings naming attributes to be removed from time_vary - - Returns - ------- - None - ''' - for param in params: - if param in self.time_vary: - self.time_vary.remove(param) - - def delFromTimeInv(self, *params): - ''' - Removes any number of parameters from time_inv for this instance. - - Parameters - ---------- - params : string - Any number of strings naming attributes to be removed from time_inv - - Returns - ------- - None - ''' - for param in params: - if param in self.time_inv: - self.time_inv.remove(param) - - def solve(self, verbose=False): - ''' - Solve the model for this instance of an agent type by backward induction. - Loops through the sequence of one period problems, passing the solution - from period t+1 to the problem for period t. - - Parameters - ---------- - verbose : boolean - If True, solution progress is printed to screen. - - Returns - ------- - none - ''' - - # Ignore floating point "errors". Numpy calls it "errors", but really it's excep- - # tions with well-defined answers such as 1.0/0.0 that is np.inf, -1.0/0.0 that is - # -np.inf, np.inf/np.inf is np.nan and so on. - with np.errstate(divide='ignore', over='ignore', under='ignore', invalid='ignore'): - self.preSolve() # Do pre-solution stuff - self.solution = solveAgent(self, verbose) # Solve the model by backward induction - if self.time_flow: # Put the solution in chronological order if this instance's time flow runs that way - self.solution.reverse() - self.addToTimeVary('solution') # Add solution to the list of time-varying attributes - self.postSolve() # Do post-solution stuff - - def resetRNG(self): - ''' - Reset the random number generator for this type. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - self.RNG = np.random.RandomState(self.seed) - - def checkElementsOfTimeVaryAreLists(self): - """ - A method to check that elements of time_vary are lists. - """ - for param in self.time_vary: - assert type(getattr(self, param)) == list, param + ' is not a list, but should be' + \ - ' because it is in time_vary' - - def checkRestrictions(self): - """ - A method to check that various restrictions are met for the model class. - """ - return - - def preSolve(self): - ''' - A method that is run immediately before the model is solved, to check inputs or to prepare - the terminal solution, perhaps. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - self.checkRestrictions() - self.checkElementsOfTimeVaryAreLists() - return None - - def postSolve(self): - ''' - A method that is run immediately after the model is solved, to finalize - the solution in some way. Does nothing here. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - return None - - def initializeSim(self): - ''' - Prepares this AgentType for a new simulation. Resets the internal random number generator, - makes initial states for all agents (using simBirth), clears histories of tracked variables. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - if not hasattr(self, 'T_sim'): - raise Exception('To initialize simulation variables it is necessary to first ' + - 'set the attribute T_sim to the largest number of observations ' + - 'you plan to simulate for each agent including re-births.') - elif self.T_sim <= 0: - raise Exception('T_sim represents the largest number of observations ' + - 'that can be simulated for an agent, and must be a positive number.') - - self.resetRNG() - self.t_sim = 0 - all_agents = np.ones(self.AgentCount, dtype=bool) - blank_array = np.zeros(self.AgentCount) - for var_name in self.poststate_vars: - setattr(self, var_name, copy(blank_array)) - self.t_age = np.zeros(self.AgentCount, dtype=int) # Number of periods since agent entry - self.t_cycle = np.zeros(self.AgentCount, dtype=int) # Which cycle period each agent is on - self.simBirth(all_agents) - self.clearHistory() - return None - - def simOnePeriod(self): - ''' - Simulates one period for this type. Calls the methods getMortality(), getShocks() or - readShocks, getStates(), getControls(), and getPostStates(). These should be defined for - AgentType subclasses, except getMortality (define its components simDeath and simBirth - instead) and readShocks. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - if not hasattr(self, 'solution'): - raise Exception('Model instance does not have a solution stored. To simulate, it is necessary' - ' to run the `solve()` method of the class first.') - - self.getMortality() # Replace some agents with "newborns" - if self.read_shocks: # If shock histories have been pre-specified, use those - self.readShocks() - else: # Otherwise, draw shocks as usual according to subclass-specific method - self.getShocks() - self.getStates() # Determine each agent's state at decision time - self.getControls() # Determine each agent's choice or control variables based on states - self.getPostStates() # Determine each agent's post-decision / end-of-period states using states and controls - - # Advance time for all agents - self.t_age = self.t_age + 1 # Age all consumers by one period - self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle - self.t_cycle[self.t_cycle == self.T_cycle] = 0 # Resetting to zero for those who have reached the end - - def makeShockHistory(self): - ''' - Makes a pre-specified history of shocks for the simulation. Shock variables should be named - in self.shock_vars, a list of strings that is subclass-specific. This method runs a subset - of the standard simulation loop by simulating only mortality and shocks; each variable named - in shock_vars is stored in a T_sim x AgentCount array in an attribute of self named X_hist. - Automatically sets self.read_shocks to True so that these pre-specified shocks are used for - all subsequent calls to simulate(). - - Parameters - ---------- - None - - Returns - ------- - None - ''' - # Make sure time is flowing forward and re-initialize the simulation - orig_time = self.time_flow - self.timeFwd() - self.initializeSim() - - # Make blank history arrays for each shock variable (and mortality) - for var_name in self.shock_vars: - setattr(self, var_name+'_hist', np.zeros((self.T_sim, self.AgentCount)) + np.nan) - setattr(self, 'who_dies_hist', np.zeros((self.T_sim, self.AgentCount), dtype=bool)) - - # Make and store the history of shocks for each period - for t in range(self.T_sim): - self.getMortality() - self.who_dies_hist[t,:] = self.who_dies - self.getShocks() - for var_name in self.shock_vars: - getattr(self, var_name + '_hist')[self.t_sim,:] = getattr(self, var_name) - self.t_sim += 1 - self.t_age = self.t_age + 1 # Age all consumers by one period - self.t_cycle = self.t_cycle + 1 # Age all consumers within their cycle - self.t_cycle[self.t_cycle == self.T_cycle] = 0 # Resetting to zero for those who have reached the end - - # Restore the flow of time and flag that shocks can be read rather than simulated - self.read_shocks = True - if not orig_time: - self.timeRev() - - def getMortality(self): - ''' - Simulates mortality or agent turnover according to some model-specific rules named simDeath - and simBirth (methods of an AgentType subclass). simDeath takes no arguments and returns - a Boolean array of size AgentCount, indicating which agents of this type have "died" and - must be replaced. simBirth takes such a Boolean array as an argument and generates initial - post-decision states for those agent indices. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - if self.read_shocks: - who_dies = self.who_dies_hist[self.t_sim,:] - else: - who_dies = self.simDeath() - self.simBirth(who_dies) - self.who_dies = who_dies - return None - - def simDeath(self): - ''' - Determines which agents in the current population "die" or should be replaced. Takes no - inputs, returns a Boolean array of size self.AgentCount, which has True for agents who die - and False for those that survive. Returns all False by default, must be overwritten by a - subclass to have replacement events. - - Parameters - ---------- - None - - Returns - ------- - who_dies : np.array - Boolean array of size self.AgentCount indicating which agents die and are replaced. - ''' - who_dies = np.zeros(self.AgentCount, dtype=bool) - return who_dies - - def simBirth(self, which_agents): - ''' - Makes new agents for the simulation. Takes a boolean array as an input, indicating which - agent indices are to be "born". Does nothing by default, must be overwritten by a subclass. - - Parameters - ---------- - which_agents : np.array(Bool) - Boolean array of size self.AgentCount indicating which agents should be "born". - - Returns - ------- - None - ''' - print('AgentType subclass must define method simBirth!') - return None - - def getShocks(self): - ''' - Gets values of shock variables for the current period. Does nothing by default, but can - be overwritten by subclasses of AgentType. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - return None - - def readShocks(self): - ''' - Reads values of shock variables for the current period from history arrays. For each var- - iable X named in self.shock_vars, this attribute of self is set to self.X_hist[self.t_sim,:]. - - This method is only ever called if self.read_shocks is True. This can be achieved by using - the method makeShockHistory() (or manually after storing a "handcrafted" shock history). - - Parameters - ---------- - None - - Returns - ------- - None - ''' - for var_name in self.shock_vars: - setattr(self, var_name, getattr(self, var_name + '_hist')[self.t_sim, :]) - - def getStates(self): - ''' - Gets values of state variables for the current period, probably by using post-decision states - from last period, current period shocks, and maybe market-level events. Does nothing by - default, but can be overwritten by subclasses of AgentType. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - return None - - def getControls(self): - ''' - Gets values of control variables for the current period, probably by using current states. - Does nothing by default, but can be overwritten by subclasses of AgentType. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - return None - - def getPostStates(self): - ''' - Gets values of post-decision state variables for the current period, probably by current - states and controls and maybe market-level events or shock variables. Does nothing by - default, but can be overwritten by subclasses of AgentType. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - return None - - def simulate(self, sim_periods=None): - ''' - Simulates this agent type for a given number of periods. Defaults to self.T_sim if no input. - Records histories of attributes named in self.track_vars in attributes named varname_hist. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - if not hasattr(self, 't_sim'): - raise Exception('It seems that the simulation variables were not initialize before calling ' + - 'simulate(). Call initializeSim() to initialize the variables before calling simulate() again.') - - if not hasattr(self, 'T_sim'): - raise Exception('This agent type instance must have the attribute T_sim set to a positive integer.' + - 'Set T_sim to match the largest dataset you might simulate, and run this agent\'s' + - 'initalizeSim() method before running simulate() again.') - - if sim_periods is not None and self.T_sim < sim_periods: - raise Exception('To simulate, sim_periods has to be larger than the maximum data set size ' + - 'T_sim. Either increase the attribute T_sim of this agent type instance ' + - 'and call the initializeSim() method again, or set sim_periods <= T_sim.') - - - # Ignore floating point "errors". Numpy calls it "errors", but really it's excep- - # tions with well-defined answers such as 1.0/0.0 that is np.inf, -1.0/0.0 that is - # -np.inf, np.inf/np.inf is np.nan and so on. - with np.errstate(divide='ignore', over='ignore', under='ignore', invalid='ignore'): - orig_time = self.time_flow - self.timeFwd() - if sim_periods is None: - sim_periods = self.T_sim - - for t in range(sim_periods): - self.simOnePeriod() - for var_name in self.track_vars: - getattr(self, var_name + '_hist')[self.t_sim,:] = getattr(self,var_name) - self.t_sim += 1 - - if not orig_time: - self.timeRev() - - def clearHistory(self): - ''' - Clears the histories of the attributes named in self.track_vars. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - for var_name in self.track_vars: - setattr(self, var_name + '_hist', np.zeros((self.T_sim,self.AgentCount)) + np.nan) - - -def solveAgent(agent, verbose): - ''' - Solve the dynamic model for one agent type. This function iterates on "cycles" - of an agent's model either a given number of times or until solution convergence - if an infinite horizon model is used (with agent.cycles = 0). - - Parameters - ---------- - agent : AgentType - The microeconomic AgentType whose dynamic problem is to be solved. - verbose : boolean - If True, solution progress is printed to screen (when cycles != 1). - - Returns - ------- - solution : [Solution] - A list of solutions to the one period problems that the agent will - encounter in his "lifetime". Returns in reverse chronological order. - ''' - # Record the flow of time when the Agent began the process, and make sure time is flowing backwards - original_time_flow = agent.time_flow - agent.timeRev() - - # Check to see whether this is an (in)finite horizon problem - cycles_left = agent.cycles # NOQA - infinite_horizon = cycles_left == 0 # NOQA - # Initialize the solution, which includes the terminal solution if it's not a pseudo-terminal period - solution = [] - if not agent.pseudo_terminal: - solution.append(deepcopy(agent.solution_terminal)) - - - # Initialize the process, then loop over cycles - solution_last = agent.solution_terminal # NOQA - go = True # NOQA - completed_cycles = 0 # NOQA - max_cycles = 5000 # NOQA - escape clause - if verbose: - t_last = time() - while go: - # Solve a cycle of the model, recording it if horizon is finite - solution_cycle = solveOneCycle(agent, solution_last) - if not infinite_horizon: - solution += solution_cycle - - # Check for termination: identical solutions across cycle iterations or run out of cycles - solution_now = solution_cycle[-1] - if infinite_horizon: - if completed_cycles > 0: - solution_distance = solution_now.distance(solution_last) - agent.solution_distance = solution_distance # Add these attributes so users can - agent.completed_cycles = completed_cycles # query them to see if solution is ready - go = (solution_distance > agent.tolerance and completed_cycles < max_cycles) - else: # Assume solution does not converge after only one cycle - solution_distance = 100.0 - go = True - else: - cycles_left += -1 - go = cycles_left > 0 - - # Update the "last period solution" - solution_last = solution_now - completed_cycles += 1 - - # Display progress if requested - if verbose: - t_now = time() - if infinite_horizon: - print('Finished cycle #' + str(completed_cycles) + ' in ' + str(t_now-t_last) + - ' seconds, solution distance = ' + str(solution_distance)) - else: - print('Finished cycle #' + str(completed_cycles) + ' of ' + str(agent.cycles) + - ' in ' + str(t_now-t_last) + ' seconds.') - t_last = t_now - - # Record the last cycle if horizon is infinite (solution is still empty!) - if infinite_horizon: - solution = solution_cycle # PseudoTerminal=False impossible for infinite horizon - - # Restore the direction of time to its original orientation, then return the solution - if original_time_flow: - agent.timeFwd() - return solution - - -def solveOneCycle(agent, solution_last): - ''' - Solve one "cycle" of the dynamic model for one agent type. This function - iterates over the periods within an agent's cycle, updating the time-varying - parameters and passing them to the single period solver(s). - - Parameters - ---------- - agent : AgentType - The microeconomic AgentType whose dynamic problem is to be solved. - solution_last : Solution - A representation of the solution of the period that comes after the - end of the sequence of one period problems. This might be the term- - inal period solution, a "pseudo terminal" solution, or simply the - solution to the earliest period from the succeeding cycle. - - Returns - ------- - solution_cycle : [Solution] - A list of one period solutions for one "cycle" of the AgentType's - microeconomic model. Returns in reverse chronological order. - ''' - # Calculate number of periods per cycle, defaults to 1 if all variables are time invariant - if len(agent.time_vary) > 0: - # name = agent.time_vary[0] - # T = len(eval('agent.' + name)) - T = len(agent.__dict__[agent.time_vary[0]]) - else: - T = 1 - - # Construct a dictionary to be passed to the solver - # time_inv_string = '' - # for name in agent.time_inv: - # time_inv_string += ' \'' + name + '\' : agent.' + name + ',' - # time_vary_string = '' - # for name in agent.time_vary: - # time_vary_string += ' \'' + name + '\' : None,' - # solve_dict = eval('{' + time_inv_string + time_vary_string + '}') - solve_dict = {parameter: agent.__dict__[parameter] for parameter in agent.time_inv} - solve_dict.update({parameter: None for parameter in agent.time_vary}) - - # Initialize the solution for this cycle, then iterate on periods - solution_cycle = [] - solution_next = solution_last - for t in range(T): - # Update which single period solver to use (if it depends on time) - if hasattr(agent.solveOnePeriod, "__getitem__"): - solveOnePeriod = agent.solveOnePeriod[t] - else: - solveOnePeriod = agent.solveOnePeriod - - these_args = getArgNames(solveOnePeriod) - - # Update time-varying single period inputs - for name in agent.time_vary: - if name in these_args: - # solve_dict[name] = eval('agent.' + name + '[t]') - solve_dict[name] = agent.__dict__[name][t] - solve_dict['solution_next'] = solution_next - - # Make a temporary dictionary for this period - temp_dict = {name: solve_dict[name] for name in these_args} - - # Solve one period, add it to the solution, and move to the next period - solution_t = solveOnePeriod(**temp_dict) - solution_cycle.append(solution_t) - solution_next = solution_t - - # Return the list of per-period solutions - return solution_cycle - - -# ======================================================================== -# ======================================================================== - -class Market(HARKobject): - ''' - A superclass to represent a central clearinghouse of information. Used for - dynamic general equilibrium models to solve the "macroeconomic" model as a - layer on top of the "microeconomic" models of one or more AgentTypes. - ''' - def __init__(self, agents=[], sow_vars=[], reap_vars=[], const_vars=[], track_vars=[], dyn_vars=[], - millRule=None, calcDynamics=None, act_T=1000, tolerance=0.000001,**kwds): - ''' - Make a new instance of the Market class. - - Parameters - ---------- - agents : [AgentType] - A list of all the AgentTypes in this market. - sow_vars : [string] - Names of variables generated by the "aggregate market process" that should - be "sown" to the agents in the market. Aggregate state, etc. - reap_vars : [string] - Names of variables to be collected ("reaped") from agents in the market - to be used in the "aggregate market process". - const_vars : [string] - Names of attributes of the Market instance that are used in the "aggregate - market process" but do not come from agents-- they are constant or simply - parameters inherent to the process. - track_vars : [string] - Names of variables generated by the "aggregate market process" that should - be tracked as a "history" so that a new dynamic rule can be calculated. - This is often a subset of sow_vars. - dyn_vars : [string] - Names of variables that constitute a "dynamic rule". - millRule : function - A function that takes inputs named in reap_vars and returns an object - with attributes named in sow_vars. The "aggregate market process" that - transforms individual agent actions/states/data into aggregate data to - be sent back to agents. - calcDynamics : function - A function that takes inputs named in track_vars and returns an object - with attributes named in dyn_vars. Looks at histories of aggregate - variables and generates a new "dynamic rule" for agents to believe and - act on. - act_T : int - The number of times that the "aggregate market process" should be run - in order to generate a history of aggregate variables. - tolerance: float - Minimum acceptable distance between "dynamic rules" to consider the - Market solution process converged. Distance is a user-defined metric. - - Returns - ------- - None - ''' - self.agents = agents # NOQA - self.reap_vars = reap_vars # NOQA - self.sow_vars = sow_vars # NOQA - self.const_vars = const_vars # NOQA - self.track_vars = track_vars # NOQA - self.dyn_vars = dyn_vars # NOQA - if millRule is not None: # To prevent overwriting of method-based millRules - self.millRule = millRule - if calcDynamics is not None: # Ditto for calcDynamics - self.calcDynamics = calcDynamics - self.act_T = act_T # NOQA - self.tolerance = tolerance # NOQA - self.max_loops = 1000 # NOQA - self.assignParameters(**kwds) - - self.print_parallel_error_once = True - # Print the error associated with calling the parallel method - # "solveAgents" one time. If set to false, the error will never - # print. See "solveAgents" for why this prints once or never. - - def solveAgents(self): - ''' - Solves the microeconomic problem for all AgentTypes in this market. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - # for this_type in self.agents: - # this_type.solve() - try: - multiThreadCommands(self.agents, ['solve()']) - except Exception as err: - if self.print_parallel_error_once: - # Set flag to False so this is only printed once. - self.print_parallel_error_once = False - print("**** WARNING: could not execute multiThreadCommands in HARK.core.Market.solveAgents() ", - "so using the serial version instead. This will likely be slower. " - "The multiTreadCommands() functions failed with the following error:", '\n', - sys.exc_info()[0], ':', err) # sys.exc_info()[0]) - multiThreadCommandsFake(self.agents, ['solve()']) - - def solve(self): - ''' - "Solves" the market by finding a "dynamic rule" that governs the aggregate - market state such that when agents believe in these dynamics, their actions - collectively generate the same dynamic rule. - - Parameters - ---------- - None - - Returns - ------- - None - ''' - go = True - max_loops = self.max_loops # Failsafe against infinite solution loop - completed_loops = 0 - old_dynamics = None - - while go: # Loop until the dynamic process converges or we hit the loop cap - self.solveAgents() # Solve each AgentType's micro problem - self.makeHistory() # "Run" the model while tracking aggregate variables - new_dynamics = self.updateDynamics() # Find a new aggregate dynamic rule - - # Check to see if the dynamic rule has converged (if this is not the first loop) - if completed_loops > 0: - distance = new_dynamics.distance(old_dynamics) - else: - distance = 1000000.0 - - # Move to the next loop if the terminal conditions are not met - old_dynamics = new_dynamics - completed_loops += 1 - go = distance >= self.tolerance and completed_loops < max_loops - - self.dynamics = new_dynamics # Store the final dynamic rule in self - - def reap(self): - ''' - Collects attributes named in reap_vars from each AgentType in the market, - storing them in respectively named attributes of self. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - for var_name in self.reap_vars: - harvest = [] - for this_type in self.agents: - harvest.append(getattr(this_type, var_name)) - setattr(self, var_name, harvest) - - def sow(self): - ''' - Distributes attrributes named in sow_vars from self to each AgentType - in the market, storing them in respectively named attributes. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - for var_name in self.sow_vars: - this_seed = getattr(self, var_name) - for this_type in self.agents: - setattr(this_type, var_name, this_seed) - - def mill(self): - ''' - Processes the variables collected from agents using the function millRule, - storing the results in attributes named in aggr_sow. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - # Make a dictionary of inputs for the millRule - reap_vars_string = '' - for name in self.reap_vars: - reap_vars_string += ' \'' + name + '\' : self.' + name + ',' - const_vars_string = '' - for name in self.const_vars: - const_vars_string += ' \'' + name + '\' : self.' + name + ',' - mill_dict = eval('{' + reap_vars_string + const_vars_string + '}') - - # Run the millRule and store its output in self - product = self.millRule(**mill_dict) - for j in range(len(self.sow_vars)): - this_var = self.sow_vars[j] - this_product = getattr(product, this_var) - setattr(self, this_var, this_product) - - def cultivate(self): - ''' - Has each AgentType in agents perform their marketAction method, using - variables sown from the market (and maybe also "private" variables). - The marketAction method should store new results in attributes named in - reap_vars to be reaped later. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - for this_type in self.agents: - this_type.marketAction() - - def reset(self): - ''' - Reset the state of the market (attributes in sow_vars, etc) to some - user-defined initial state, and erase the histories of tracked variables. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - for var_name in self.track_vars: # Reset the history of tracked variables - setattr(self, var_name + '_hist', []) - for var_name in self.sow_vars: # Set the sow variables to their initial levels - initial_val = getattr(self, var_name + '_init') - setattr(self, var_name, initial_val) - for this_type in self.agents: # Reset each AgentType in the market - this_type.reset() - - def store(self): - ''' - Record the current value of each variable X named in track_vars in an - attribute named X_hist. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - for var_name in self.track_vars: - value_now = getattr(self, var_name) - getattr(self, var_name + '_hist').append(value_now) - - def makeHistory(self): - ''' - Runs a loop of sow-->cultivate-->reap-->mill act_T times, tracking the - evolution of variables X named in track_vars in attributes named X_hist. - - Parameters - ---------- - none - - Returns - ------- - none - ''' - self.reset() # Initialize the state of the market - for t in range(self.act_T): - self.sow() # Distribute aggregated information/state to agents - self.cultivate() # Agents take action - self.reap() # Collect individual data from agents - self.mill() # Process individual data into aggregate data - self.store() # Record variables of interest - - def updateDynamics(self): - ''' - Calculates a new "aggregate dynamic rule" using the history of variables - named in track_vars, and distributes this rule to AgentTypes in agents. - - Parameters - ---------- - none - - Returns - ------- - dynamics : instance - The new "aggregate dynamic rule" that agents believe in and act on. - Should have attributes named in dyn_vars. - ''' - # Make a dictionary of inputs for the dynamics calculator - history_vars_string = '' - arg_names = list(getArgNames(self.calcDynamics)) - if 'self' in arg_names: - arg_names.remove('self') - for name in arg_names: - history_vars_string += ' \'' + name + '\' : self.' + name + '_hist,' - update_dict = eval('{' + history_vars_string + '}') - - # Calculate a new dynamic rule and distribute it to the agents in agent_list - dynamics = self.calcDynamics(**update_dict) # User-defined dynamics calculator - for var_name in self.dyn_vars: - this_obj = getattr(dynamics, var_name) - for this_type in self.agents: - setattr(this_type, var_name, this_obj) - return dynamics - - -# ------------------------------------------------------------------------------ -# Code to copy entire modules to a local directory -# ------------------------------------------------------------------------------ - -# Define a function to run the copying: -def copy_module(target_path, my_directory_full_path, my_module): - ''' - Helper function for copy_module_to_local(). Provides the actual copy - functionality, with highly cautious safeguards against copying over - important things. - - Parameters - ---------- - target_path : string - String, file path to target location - - my_directory_full_path: string - String, full pathname to this file's directory - - my_module : string - String, name of the module to copy - - Returns - ------- - none - ''' - - if target_path == 'q' or target_path == 'Q': - print("Goodbye!") - return - elif target_path == os.path.expanduser("~") or os.path.normpath(target_path) == os.path.expanduser("~"): - print("You have indicated that the target location is " + target_path + - " -- that is, you want to wipe out your home directory with the contents of " + my_module + - ". My programming does not allow me to do that.\n\nGoodbye!") - return - elif os.path.exists(target_path): - print("There is already a file or directory at the location " + target_path + - ". For safety reasons this code does not overwrite existing files.\n Please remove the file at " - + target_path + - " and try again.") - return - else: - user_input = input("""You have indicated you want to copy module:\n """ + my_module - + """\nto:\n """ + target_path + """\nIs that correct? Please indicate: y / [n]\n\n""") - if user_input == 'y' or user_input == 'Y': - # print("copy_tree(",my_directory_full_path,",", target_path,")") - copy_tree(my_directory_full_path, target_path) - else: - print("Goodbye!") - return - - -def print_helper(): - - my_directory_full_path = os.path.dirname(os.path.realpath(__file__)) - - print(my_directory_full_path) - - -def copy_module_to_local(full_module_name): - ''' - This function contains simple code to copy a submodule to a location on - your hard drive, as specified by you. The purpose of this code is to provide - users with a simple way to access a *copy* of code that usually sits deep in - the Econ-ARK package structure, for purposes of tinkering and experimenting - directly. This is meant to be a simple way to explore HARK code. To interact - with the codebase under active development, please refer to the documentation - under github.com/econ-ark/HARK/ - - To execute, do the following on the Python command line: - - from HARK.core import copy_module_to_local - copy_module_to_local("FULL-HARK-MODULE-NAME-HERE") - - For example, if you want SolvingMicroDSOPs you would enter - - from HARK.core import copy_module_to_local - copy_module_to_local("HARK.SolvingMicroDSOPs") - - ''' - - # Find a default directory -- user home directory: - home_directory_RAW = os.path.expanduser("~") - # Thanks to https://stackoverflow.com/a/4028943 - - # Find the directory of the HARK.core module: - # my_directory_full_path = os.path.dirname(os.path.realpath(__file__)) - hark_core_directory_full_path = os.path.dirname(os.path.realpath(__file__)) - # From https://stackoverflow.com/a/5137509 - # Important note from that answer: - # (Note that the incantation above won't work if you've already used os.chdir() - # to change your current working directory, - # since the value of the __file__ constant is relative to the current working directory and is not changed by an - # os.chdir() call.) - # - # NOTE: for this specific file that I am testing, the path should be: - # '/home/npalmer/anaconda3/envs/py3fresh/lib/python3.6/site-packages/HARK/SolvingMicroDSOPs/---example-file--- - - # Split out the name of the module. Break if proper format is not followed: - all_module_names_list = full_module_name.split('.') # Assume put in at correct format - if all_module_names_list[0] != "HARK": - print("\nWarning: the module name does not start with 'HARK'. Instead it is: '" - + all_module_names_list[0]+"' --please format the full namespace of the module you want. \n" - "For example, 'HARK.SolvingMicroDSOPs'") - print("\nGoodbye!") - return - - # Construct the pathname to the module to copy: - my_directory_full_path = hark_core_directory_full_path - for a_directory_name in all_module_names_list[1:]: - my_directory_full_path = os.path.join(my_directory_full_path, a_directory_name) - - head_path, my_module = os.path.split(my_directory_full_path) - - home_directory_with_module = os.path.join(home_directory_RAW, my_module) - - print("\n\n\nmy_directory_full_path:", my_directory_full_path, '\n\n\n') - - # Interact with the user: - # - Ask the user for the target place to copy the directory - # - Offer use "q/y/other" option - # - Check if there is something there already - # - If so, ask if should replace - # - If not, just copy there - # - Quit - - target_path = input("""You have invoked the 'replicate' process for the current module:\n """ + - my_module + """\nThe default copy location is your home directory:\n """ + - home_directory_with_module + """\nPlease enter one of the three options in single quotes below, excluding the quotes: - - 'q' or return/enter to quit the process - 'y' to accept the default home directory: """+home_directory_with_module+""" - 'n' to specify your own pathname\n\n""") - - if target_path == 'n' or target_path == 'N': - target_path = input("""Please enter the full pathname to your target directory location: """) - - # Clean up: - target_path = os.path.expanduser(target_path) - target_path = os.path.expandvars(target_path) - target_path = os.path.normpath(target_path) - - # Check to see if they included the module name; if not add it here: - temp_head, temp_tail = os.path.split(target_path) - if temp_tail != my_module: - target_path = os.path.join(target_path, my_module) - - elif target_path == 'y' or target_path == 'Y': - # Just using the default path: - target_path = home_directory_with_module - else: - # Assume "quit" - return - - if target_path != 'q' and target_path != 'Q' or target_path == '': - # Run the copy command: - copy_module(target_path, my_directory_full_path, my_module) - - return - - if target_path != 'q' and target_path != 'Q' or target_path == '': - # Run the copy command: - copy_module(target_path, my_directory_full_path, my_module) - - return - - -def main(): - print("Sorry, HARK.core doesn't actually do anything on its own.") - print("To see some examples of its frameworks in action, try running a model module.") - print("Several interesting model modules can be found in /ConsumptionSavingModel.") - print('For an extraordinarily simple model that demonstrates the "microeconomic" and') - print('"macroeconomic" frameworks, see /FashionVictim/FashionVictimModel.') - - -if __name__ == '__main__': - main() diff --git a/HARK/#scratch# b/HARK/#scratch# deleted file mode 100644 index e69de29bb..000000000