From 48ef8d22dde0cdcc735f0a850446151752e695fa Mon Sep 17 00:00:00 2001 From: Alan Lujan Date: Thu, 2 Jun 2022 12:34:15 -0500 Subject: [PATCH 1/4] Update population.py initial PortfolioSharkFinAgentType --- sharkfin/population.py | 85 ++++++++++++++++++++++++++---------------- 1 file changed, 52 insertions(+), 33 deletions(-) diff --git a/sharkfin/population.py b/sharkfin/population.py index cbf88f0..0b15537 100644 --- a/sharkfin/population.py +++ b/sharkfin/population.py @@ -1,10 +1,29 @@ -import HARK.ConsumptionSaving.ConsPortfolioModel as cpm -import HARK.ConsumptionSaving.ConsIndShockModel as cism -from sharkfin.utilities import * import math -import pandas as pd import random +import HARK.ConsumptionSaving.ConsIndShockModel as cism +import HARK.ConsumptionSaving.ConsPortfolioModel as cpm +import pandas as pd + +from sharkfin.utilities import * + + +class PortfolioSharkFinAgentType(cpm.SequentialPortfolioConsumerType): + """ + SHARKFin agent class extends HARK's SequentialPortfolioConsumerType class + with SHARK features. In particular, it takes an external solution in the + simulation step. This external solution can come from interpolated + population solutions. + """ + + def simulate(self, solution, sim_periods=None): + """ + Provide solution to SHARKFin Agent before internal HARK simulation step. + """ + self.solution = solution + super().simulate(sim_periods) + + class AgentPopulation: """ A class encapsulating the population of 'macroeconomy' agents. @@ -45,21 +64,21 @@ def __init__(self, base_parameters, dist_params, n_per_class): ) def agent_df(self): - ''' + """ Output a dataframe for agent attributes returns agent_df from class_stats - ''' + """ records = [] for agent in self.agents: - for i, aLvl in enumerate(agent.state_now['aLvl']): + for i, aLvl in enumerate(agent.state_now["aLvl"]): record = { - 'aLvl': aLvl, - 'mNrm': agent.state_now['mNrm'][i], - 'cNrm': agent.controls['cNrm'][i] - if 'cNrm' in agent.controls + "aLvl": aLvl, + "mNrm": agent.state_now["mNrm"][i], + "cNrm": agent.controls["cNrm"][i] + if "cNrm" in agent.controls else None, } @@ -80,12 +99,12 @@ def class_stats(self, store=False): records = [] for agent in self.agents: - for i, aLvl in enumerate(agent.state_now['aLvl']): + for i, aLvl in enumerate(agent.state_now["aLvl"]): record = { - 'aLvl': aLvl, - 'mNrm': agent.state_now['mNrm'][i], + "aLvl": aLvl, + "mNrm": agent.state_now["mNrm"][i], # difference between mNrm and the equilibrium mNrm from BST - 'mNrm_ratio_StE': agent.state_now['mNrm'][i] / agent.mNrmStE, + "mNrm_ratio_StE": agent.state_now["mNrm"][i] / agent.mNrmStE, } for dp in self.dist_params: @@ -97,20 +116,20 @@ def class_stats(self, store=False): class_stats = ( agent_df.groupby(list(self.dist_params.keys())) - .aggregate(['mean', 'std']) + .aggregate(["mean", "std"]) .reset_index() ) cs = class_stats - cs['label'] = round(cs['CRRA'], 2).apply(lambda x: f'CRRA: {x}, ') + round( - cs['DiscFac'], 2 + cs["label"] = round(cs["CRRA"], 2).apply(lambda x: f"CRRA: {x}, ") + round( + cs["DiscFac"], 2 ).apply(lambda x: f"DiscFac: {x}") - cs['aLvl_mean'] = cs['aLvl']['mean'] - cs['aLvl_std'] = cs['aLvl']['std'] - cs['mNrm_mean'] = cs['mNrm']['mean'] - cs['mNrm_std'] = cs['mNrm']['std'] - cs['mNrm_ratio_StE_mean'] = cs['mNrm_ratio_StE']['mean'] - cs['mNrm_ratio_StE_std'] = cs['mNrm_ratio_StE']['std'] + cs["aLvl_mean"] = cs["aLvl"]["mean"] + cs["aLvl_std"] = cs["aLvl"]["std"] + cs["mNrm_mean"] = cs["mNrm"]["mean"] + cs["mNrm_std"] = cs["mNrm"]["std"] + cs["mNrm_ratio_StE_mean"] = cs["mNrm_ratio_StE"]["mean"] + cs["mNrm_ratio_StE_std"] = cs["mNrm_ratio_StE"]["std"] if store: self.stored_class_stats = class_stats @@ -133,8 +152,8 @@ def create_distributed_agents(self, agent_parameters, dist_params, n_per_class): n_per_class: int number of agents to instantiate per class """ - num_classes = math.prod([dist_params[dp]['n'] for dp in dist_params]) - agent_batches = [{'AgentCount': num_classes}] * n_per_class + num_classes = math.prod([dist_params[dp]["n"] for dp in dist_params]) + agent_batches = [{"AgentCount": num_classes}] * n_per_class agents = [ cpm.PortfolioConsumerType( @@ -153,7 +172,7 @@ def init_simulation(self): Sets up the agents with their state for the state of the simulation """ for agent in self.agents: - agent.track_vars += ['pLvl', 'mNrm', 'cNrm', 'Share', 'Risky'] + agent.track_vars += ["pLvl", "mNrm", "cNrm", "Share", "Risky"] agent.assign_parameters(AdjustPrb=1.0) agent.T_sim = 1000 # arbitrary! @@ -171,7 +190,7 @@ def init_simulation(self): ind_shock_double.solve() mNrmStE = ind_shock_double.solution[0].mNrmStE - agent.state_now['mNrm'][:] = mNrmStE + agent.state_now["mNrm"][:] = mNrmStE agent.mNrmStE = ( mNrmStE # saving this for later, in case we do the analysis. ) @@ -180,11 +199,11 @@ def init_simulation(self): mNrm = ( self.stored_class_stats.copy() .set_index([dp for dp in self.dist_params]) - .xs((idx))['mNrm']['mean'] + .xs((idx))["mNrm"]["mean"] ) - agent.state_now['mNrm'][:] = mNrm + agent.state_now["mNrm"][:] = mNrm - agent.state_now['aNrm'] = agent.state_now['mNrm'] - agent.solution[ + agent.state_now["aNrm"] = agent.state_now["mNrm"] - agent.solution[ 0 - ].cFuncAdj(agent.state_now['mNrm']) - agent.state_now['aLvl'] = agent.state_now['aNrm'] * agent.state_now['pLvl'] + ].cFuncAdj(agent.state_now["mNrm"]) + agent.state_now["aLvl"] = agent.state_now["aNrm"] * agent.state_now["pLvl"] From 5f6c60b042c224ea5585d79cd173ac1fc41b8063 Mon Sep 17 00:00:00 2001 From: Alan Lujan Date: Thu, 3 Nov 2022 09:27:33 -0400 Subject: [PATCH 2/4] update rfree --- sharkfin/population.py | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/sharkfin/population.py b/sharkfin/population.py index f6aabb3..9f43f8c 100644 --- a/sharkfin/population.py +++ b/sharkfin/population.py @@ -688,3 +688,37 @@ def _merge_solutions_3d(self, continuous_states): self.solution_database = self.solution_database.set_index(discrete_params) return self.solution_database + + +class PortfolioSharkFinAgentType(cpm.SequentialPortfolioConsumerType): + """ + SHARKFin agent class extends HARK's SequentialPortfolioConsumerType class + with SHARK features. In particular, it takes an external solution in the + simulation step. This external solution can come from interpolated + population solutions. + """ + + def get_Rfree(self): + """ + Calculates realized return factor for each agent, using the attributes Rfree, + RiskyNow, and ShareNow. This method is a bit of a misnomer, as the return + factor is not riskless, but would more accurately be labeled as Rport. However, + this method makes the portfolio model compatible with its parent class. + + Parameters + ---------- + None + + Returns + ------- + Rport : np.array + Array of size AgentCount with each simulated agent's realized portfolio + return factor. Will be used by get_states() to calculate mNrmNow, where it + will be mislabeled as "Rfree". + """ + + Rfree = np.asarray(self.Rfree) + RfreeNow = Rfree[self.t_cycle - 1] + + Rport = self.RiskyAvg * np.ones_like(RfreeNow) + return Rport From 3582a5300f459d5ba081ce2515da8d46246e6320 Mon Sep 17 00:00:00 2001 From: Alan Lujan Date: Mon, 16 Jan 2023 16:45:30 -0500 Subject: [PATCH 3/4] organize important simulation methods simulate method for SHARKFIN needs to do 2 kinds of simulations, macro day and intraday --- sharkfin/population.py | 249 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 236 insertions(+), 13 deletions(-) diff --git a/sharkfin/population.py b/sharkfin/population.py index 9f43f8c..37d2a6c 100644 --- a/sharkfin/population.py +++ b/sharkfin/population.py @@ -3,12 +3,19 @@ from pprint import pprint from typing import NewType -import HARK.ConsumptionSaving.ConsIndShockModel as cism -import HARK.ConsumptionSaving.ConsPortfolioModel as cpm import numpy as np import pandas as pd +from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType +from HARK.ConsumptionSaving.ConsPortfolioModel import SequentialPortfolioConsumerType +from HARK.ConsumptionSaving.ConsRiskyAssetModel import RiskyAssetConsumerType from HARK.core import AgentType -from HARK.distribution import Distribution, IndexDistribution, combine_indep_dstns +from HARK.distribution import ( + Bernoulli, + Distribution, + IndexDistribution, + Lognormal, + combine_indep_dstns, +) from HARK.interpolation import BilinearInterpOnInterp1D, TrilinearInterpOnInterp1D from xarray import DataArray @@ -263,17 +270,18 @@ def init_simulation(self, T_sim=1000): if self.stored_class_stats is None: - ## build an IndShockConsumerType "double" of this agent, with the same parameters - ind_shock_double = cism.IndShockConsumerType(**agent.parameters) + # build an IndShockConsumerType "double" of this agent, with the same parameters + ind_shock_double = IndShockConsumerType(**agent.parameters) - ## solve to get the mNrmStE value - ## that is, the Steady-state Equilibrium value of mNrm, for the IndShockModel + # solve to get the mNrmStE value + # that is, the Steady-state Equilibrium value of mNrm, for the IndShockModel ind_shock_double.solve() mNrmStE = ind_shock_double.solution[0].mNrmStE agent.state_now["mNrm"][:] = mNrmStE agent.mNrmStE = ( - mNrmStE # saving this for later, in case we do the analysis. + # saving this for later, in case we do the analysis. + mNrmStE ) else: idx = [agent.parameters[dp] for dp in self.dist_params] @@ -459,7 +467,7 @@ def macro_update(self, agent, price): print("ERROR: Agent has share > 1 after macro update.") print(true_risky_expectations) - ## put back the expectations that include capital gains now + # put back the expectations that include capital gains now agent.assign_parameters(**true_risky_expectations) # Selling off shares if necessary to @@ -519,7 +527,7 @@ def update_agent_wealth_capital_gains(self, new_share_price, pror, dividend): ) print("Setting normalize assets and shares to 0.") agent.state_now["aNrm"][(agent.state_now["aNrm"] < 0)] = 0.0 - ## TODO: This change in shares needs to be registered with the Broker. + # TODO: This change in shares needs to be registered with the Broker. agent.shares[(agent.state_now["aNrm"] == 0)] = 0 # update non-normalized market assets @@ -690,7 +698,7 @@ def _merge_solutions_3d(self, continuous_states): return self.solution_database -class PortfolioSharkFinAgentType(cpm.SequentialPortfolioConsumerType): +class PortfolioSharkFinAgentType(SequentialPortfolioConsumerType): """ SHARKFin agent class extends HARK's SequentialPortfolioConsumerType class with SHARK features. In particular, it takes an external solution in the @@ -698,6 +706,101 @@ class PortfolioSharkFinAgentType(cpm.SequentialPortfolioConsumerType): population solutions. """ + def initialize_sim(self): + """ + Initialize the state of simulation attributes. Simply calls the same method + for IndShockConsumerType, then sets the type of AdjustNow to bool. + + Parameters + ---------- + None + + Returns + ------- + None + """ + # these need to be set because "post states", + # but are a control variable and shock, respectively + self.controls["Share"] = np.zeros(self.AgentCount) + RiskyAssetConsumerType.initialize_sim(self) + + def transition(self): + pLvlPrev = self.state_prev["pLvl"] + aNrmPrev = self.state_prev["aNrm"] + RfreeNow = self.get_Rfree() + + # Calculate new states: normalized market resources and permanent income level + # Updated permanent income level + pLvlNow = pLvlPrev * self.shocks["PermShk"] + # Updated aggregate permanent productivity level + PlvlAggNow = self.state_prev["PlvlAgg"] * self.PermShkAggNow + # "Effective" interest factor on normalized assets + ReffNow = RfreeNow / self.shocks["PermShk"] + bNrmNow = ReffNow * aNrmPrev # Bank balances before labor income + # Market resources after income + mNrmNow = bNrmNow + self.shocks["TranShk"] + + return pLvlNow, PlvlAggNow, bNrmNow, mNrmNow, None + + def get_controls(self): + """ + Calculates consumption cNrmNow and risky portfolio share ShareNow using + the policy functions in the attribute solution. These are stored as attributes. + + Parameters + ---------- + None + + Returns + ------- + None + """ + cNrmNow = np.zeros(self.AgentCount) + np.nan + ShareNow = np.zeros(self.AgentCount) + np.nan + + # Loop over each period of the cycle, getting controls separately depending on "age" + for t in range(self.T_cycle): + these = t == self.t_cycle + + # Get controls for agents who *can* adjust their portfolio share + those = np.logical_and(these, self.shocks["Adjust"]) + cNrmNow[those] = self.solution[t].cFuncAdj(self.state_now["mNrm"][those]) + ShareNow[those] = self.solution[t].ShareFuncAdj( + self.state_now["mNrm"][those] + ) + + # Get Controls for agents who *can't* adjust their portfolio share + those = np.logical_and(these, np.logical_not(self.shocks["Adjust"])) + cNrmNow[those] = self.solution[t].cFuncFxd( + self.state_now["mNrm"][those], ShareNow[those] + ) + ShareNow[those] = self.solution[t].ShareFuncFxd( + self.state_now["mNrm"][those], ShareNow[those] + ) + + # Store controls as attributes of self + self.controls["cNrm"] = cNrmNow + self.controls["Share"] = ShareNow + + def get_poststates(self): + """ + Calculates end-of-period assets for each consumer of this type. + + Parameters + ---------- + None + + Returns + ------- + None + """ + # should this be "Now", or "Prev"?!? + self.state_now["aNrm"] = self.state_now["mNrm"] - self.controls["cNrm"] + # Useful in some cases to precalculate asset level + self.state_now["aLvl"] = self.state_now["aNrm"] * self.state_now["pLvl"] + + return None + def get_Rfree(self): """ Calculates realized return factor for each agent, using the attributes Rfree, @@ -717,8 +820,128 @@ def get_Rfree(self): will be mislabeled as "Rfree". """ - Rfree = np.asarray(self.Rfree) + Rfree = np.array(self.Rfree) RfreeNow = Rfree[self.t_cycle - 1] - Rport = self.RiskyAvg * np.ones_like(RfreeNow) + Rport = ( + self.controls["Share"] * self.shocks["Risky"] + + (1.0 - self.controls["Share"]) * RfreeNow + ) + self.Rport = Rport return Rport + + def get_Risky(self): + """ + Sets the attribute Risky as a single draw from a lognormal distribution. + Uses the attributes RiskyAvgTrue and RiskyStdTrue if RiskyAvg is time-varying, + else just uses the single values from RiskyAvg and RiskyStd. + + Parameters + ---------- + None + + Returns + ------- + None + """ + if "RiskyDstn" in self.time_vary: + RiskyAvg = self.RiskyAvgTrue + RiskyStd = self.RiskyStdTrue + else: + RiskyAvg = self.RiskyAvg + RiskyStd = self.RiskyStd + + self.shocks["Risky"] = Lognormal.from_mean_std( + RiskyAvg, RiskyStd, seed=self.RNG.randint(0, 2**31 - 1) + ).draw(1) + + def get_Adjust(self): + """ + Sets the attribute Adjust as a boolean array of size AgentCount, indicating + whether each agent is able to adjust their risky portfolio share this period. + Uses the attribute AdjustPrb to draw from a Bernoulli distribution. + + Parameters + ---------- + None + + Returns + ------- + None + """ + self.shocks["Adjust"] = IndexDistribution( + Bernoulli, {"p": self.AdjustPrb}, seed=self.RNG.randint(0, 2**31 - 1) + ).draw(self.t_cycle) + + def get_shocks(self): + """ + Draw idiosyncratic income shocks, just as for IndShockConsumerType, then draw + a single common value for the risky asset return. Also draws whether each + agent is able to adjust their portfolio this period. + + Parameters + ---------- + None + + Returns + ------- + None + """ + IndShockConsumerType.get_shocks(self) + self.get_Risky() + self.get_Adjust() + + def sim_birth(self, which_agents): + """ + Create new agents to replace ones who have recently died; takes draws of + initial aNrm and pLvl, as in ConsIndShockModel, then sets Share and Adjust + to zero as initial values. + Parameters + ---------- + which_agents : np.array + Boolean array of size AgentCount indicating which agents should be "born". + + Returns + ------- + None + """ + IndShockConsumerType.sim_birth(self, which_agents) + + self.controls["Share"][which_agents] = 0 + # here a shock is being used as a 'post state' + self.shocks["Adjust"][which_agents] = False + + def sim_death(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 + DiePrb_by_t_cycle = 1.0 - np.asarray(self.LivPrb) + DiePrb = DiePrb_by_t_cycle[ + self.t_cycle - 1 if self.cycles == 1 else self.t_cycle + ] # Time has already advanced, so look back one + + # In finite-horizon problems the previous line gives newborns the + # survival probability of the last non-terminal period. This is okay, + # however, since they will be instantly replaced by new newborns if + # they die. + # See: https://github.com/econ-ark/HARK/pull/981 + + DeathShks = Uniform(seed=self.RNG.randint(0, 2**31 - 1)).draw( + N=self.AgentCount + ) + 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 From c0ba986092e2d2a30f91839e5605a43610205d89 Mon Sep 17 00:00:00 2001 From: Alan Lujan Date: Mon, 16 Jan 2023 16:50:33 -0500 Subject: [PATCH 4/4] update discretize --- sharkfin/population.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sharkfin/population.py b/sharkfin/population.py index 37d2a6c..be4843d 100644 --- a/sharkfin/population.py +++ b/sharkfin/population.py @@ -96,7 +96,7 @@ def approx_distributions(self, approx_params: dict): for key in approx_params: if key in param_dict and isinstance(param_dict[key], Distribution): discrete_points = approx_params[key] - discrete_distribution = param_dict[key].approx(discrete_points) + discrete_distribution = param_dict[key].discretize(discrete_points) self.continuous_distributions[key] = param_dict[key] self.discrete_distributions[key] = discrete_distribution else: