Coverage for HARK / ConsumptionSaving / ConsWealthPortfolioModel.py: 96%
159 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-25 05:22 +0000
1from copy import deepcopy
3import numpy as np
4from HARK.ConsumptionSaving.ConsPortfolioModel import (
5 PortfolioConsumerType,
6 PortfolioSolution,
7 make_portfolio_solution_terminal,
8)
9from HARK.distributions import expected
10from HARK.interpolation import (
11 BilinearInterp,
12 ConstantFunction,
13 CubicInterp,
14 LinearInterp,
15 MargValueFuncCRRA,
16 ValueFuncCRRA,
17)
18from HARK.Calibration.Assets.AssetProcesses import (
19 make_lognormal_RiskyDstn,
20 combine_IncShkDstn_and_RiskyDstn,
21 calc_ShareLimit_for_CRRA,
22)
23from HARK.Calibration.Income.IncomeProcesses import (
24 construct_lognormal_income_process_unemployment,
25 get_PermShkDstn_from_IncShkDstn,
26 get_TranShkDstn_from_IncShkDstn,
27)
28from HARK.ConsumptionSaving.ConsRiskyAssetModel import (
29 make_simple_ShareGrid,
30 make_AdjustDstn,
31)
32from HARK.ConsumptionSaving.ConsWealthUtilityModel import (
33 make_ChiFromOmega_function,
34)
35from HARK.ConsumptionSaving.ConsIndShockModel import (
36 make_lognormal_kNrm_init_dstn,
37 make_lognormal_pLvl_init_dstn,
38)
39from HARK.rewards import UtilityFuncCRRA
40from HARK.utilities import NullFunc, make_assets_grid
43def utility(c, a, CRRA, share=0.0, intercept=0.0):
44 w = a + intercept
45 return (c ** (1 - share) * w**share) ** (1 - CRRA) / (1 - CRRA)
48def dudc(c, a, CRRA, share=0.0, intercept=0.0):
49 u = utility(c, a, CRRA, share, intercept)
50 return u * (1 - CRRA) * (1 - share) / c
53def calc_m_nrm_next(shocks, b_nrm, perm_gro_fac):
54 """
55 Calculate future realizations of market resources mNrm from the income
56 shock distribution S and normalized bank balances b.
57 """
58 return b_nrm / (shocks["PermShk"] * perm_gro_fac) + shocks["TranShk"]
61def calc_dvdm_next(shocks, b_nrm, perm_gro_fac, crra, vp_func):
62 """
63 Evaluate realizations of marginal value of market resources next period,
64 based on the income distribution S and values of bank balances bNrm
65 """
66 m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac)
67 perm_shk_fac = shocks["PermShk"] * perm_gro_fac
68 return perm_shk_fac ** (-crra) * vp_func(m_nrm)
71def calc_end_dvdx(shocks, a_nrm, share, rfree, dvdb_func):
72 ex_ret = shocks - rfree # Excess returns
73 rport = rfree + share * ex_ret # Portfolio return
74 b_nrm = rport * a_nrm
76 # Calculate and return dvds (second term is all zeros)
77 dvdb = dvdb_func(b_nrm)
78 dvda = rport * dvdb
79 dvds = ex_ret * a_nrm * dvdb
80 return dvda, dvds
83def calc_med_v(shocks, b_nrm, perm_gro_fac, crra, v_func):
84 """
85 Calculate "intermediate" value from next period's bank balances, the
86 income shocks S, and the risky asset share.
87 """
88 m_nrm = calc_m_nrm_next(shocks, b_nrm, perm_gro_fac)
89 v_next = v_func(m_nrm)
90 return (shocks["PermShk"] * perm_gro_fac) ** (1.0 - crra) * v_next
93def calc_end_v(shocks, a_nrm, share, rfree, v_func):
94 # Calculate future realizations of bank balances bNrm
95 ex_ret = shocks - rfree
96 rport = rfree + share * ex_ret
97 b_nrm = rport * a_nrm
98 return v_func(b_nrm)
101###############################################################################
104def solve_one_period_WealthPortfolio(
105 solution_next,
106 IncShkDstn,
107 RiskyDstn,
108 LivPrb,
109 DiscFac,
110 CRRA,
111 Rfree,
112 PermGroFac,
113 BoroCnstArt,
114 aXtraGrid,
115 ShareGrid,
116 ShareLimit,
117 vFuncBool,
118 WealthShare,
119 WealthShift,
120 ChiFunc,
121):
122 """
123 TODO: Fill in this missing docstring.
125 Parameters
126 ----------
127 solution_next : TYPE
128 DESCRIPTION.
129 IncShkDstn : TYPE
130 DESCRIPTION.
131 RiskyDstn : TYPE
132 DESCRIPTION.
133 LivPrb : TYPE
134 DESCRIPTION.
135 DiscFac : TYPE
136 DESCRIPTION.
137 CRRA : TYPE
138 DESCRIPTION.
139 Rfree : TYPE
140 DESCRIPTION.
141 PermGroFac : TYPE
142 DESCRIPTION.
143 BoroCnstArt : TYPE
144 DESCRIPTION.
145 aXtraGrid : TYPE
146 DESCRIPTION.
147 ShareGrid : TYPE
148 DESCRIPTION.
149 ShareLimit : TYPE
150 DESCRIPTION.
151 vFuncBool : TYPE
152 DESCRIPTION.
153 WealthShare : TYPE
154 DESCRIPTION.
155 WealthShift : TYPE
156 DESCRIPTION.
157 ChiFunc : TYPE
158 DESCRIPTION.
160 Returns
161 -------
162 solution_now : TYPE
163 DESCRIPTION.
165 """
166 # Make sure the individual is liquidity constrained. Allowing a consumer to
167 # borrow *and* invest in an asset with unbounded (negative) returns is a bad mix.
168 if BoroCnstArt != 0.0:
169 raise ValueError("PortfolioConsumerType must have BoroCnstArt=0.0!")
171 # Define the current period utility function and effective discount factor
172 uFunc = UtilityFuncCRRA(CRRA)
173 DiscFacEff = DiscFac * LivPrb # "effective" discount factor
175 # Unpack next period's solution for easier access
176 vp_func_next = solution_next.vPfuncAdj
177 v_func_next = solution_next.vFuncAdj
179 # Set a flag for whether the natural borrowing constraint is zero, which
180 # depends on whether the smallest transitory income shock is zero
181 BoroCnstNat_iszero = (np.min(IncShkDstn.atoms[1]) == 0.0) or (
182 WealthShare != 0.0 and WealthShift == 0.0
183 )
185 # Prepare to calculate end-of-period marginal values by creating an array
186 # of market resources that the agent could have next period, considering
187 # the grid of end-of-period assets and the distribution of shocks he might
188 # experience next period.
190 # Unpack the risky return shock distribution
191 Risky_next = RiskyDstn.atoms
192 RiskyMax = np.max(Risky_next)
193 RiskyMin = np.min(Risky_next)
195 # bNrm represents R*a, balances after asset return shocks but before income.
196 # This just uses the highest risky return as a rough shifter for the aXtraGrid.
197 if BoroCnstNat_iszero:
198 aNrmGrid = aXtraGrid
199 bNrmGrid = np.insert(RiskyMax * aXtraGrid, 0, RiskyMin * aXtraGrid[0])
200 else:
201 # Add an asset point at exactly zero
202 aNrmGrid = np.insert(aXtraGrid, 0, 0.0)
203 bNrmGrid = RiskyMax * np.insert(aXtraGrid, 0, 0.0)
205 # Get grid and shock sizes, for easier indexing
206 aNrmCount = aNrmGrid.size
208 # Make tiled arrays to calculate future realizations of mNrm and Share when integrating over IncShkDstn
209 bNrmNext = bNrmGrid
211 # Calculate end-of-period marginal value of assets and shares at each point
212 # in aNrm and ShareGrid. Does so by taking expectation of next period marginal
213 # values across income and risky return shocks.
215 # Calculate intermediate marginal value of bank balances by taking expectations over income shocks
216 med_dvdb = expected(
217 calc_dvdm_next,
218 IncShkDstn,
219 args=(bNrmNext, PermGroFac, CRRA, vp_func_next),
220 )
221 med_dvdb_nvrs = uFunc.derinv(med_dvdb, order=(1, 0))
222 med_dvdb_nvrs_func = LinearInterp(bNrmGrid, med_dvdb_nvrs)
223 med_dvdb_func = MargValueFuncCRRA(med_dvdb_nvrs_func, CRRA)
225 # Make tiled arrays to calculate future realizations of bNrm and Share when integrating over RiskyDstn
226 aNrmNow, ShareNext = np.meshgrid(aNrmGrid, ShareGrid, indexing="ij")
228 # Evaluate realizations of value and marginal value after asset returns are realized
229 end_dvda, end_dvds = DiscFacEff * expected(
230 calc_end_dvdx,
231 RiskyDstn,
232 args=(aNrmNow, ShareNext, Rfree, med_dvdb_func),
233 )
234 end_dvda_nvrs = uFunc.derinv(end_dvda)
236 # Now find the optimal (continuous) risky share on [0,1] by solving the first
237 # order condition end_dvds == 0.
238 focs = end_dvds # Relabel for convenient typing
240 # For each value of aNrm, find the value of Share such that focs == 0
241 crossing = np.logical_and(focs[:, 1:] <= 0.0, focs[:, :-1] >= 0.0)
242 share_idx = np.argmax(crossing, axis=1)
243 # This represents the index of the segment of the share grid where dvds flips
244 # from positive to negative, indicating that there's a zero *on* the segment
246 # Calculate the fractional distance between those share gridpoints where the
247 # zero should be found, assuming a linear function; call it alpha
248 a_idx = np.arange(aNrmCount)
249 bot_s = ShareGrid[share_idx]
250 top_s = ShareGrid[share_idx + 1]
251 bot_f = focs[a_idx, share_idx]
252 top_f = focs[a_idx, share_idx + 1]
253 bot_c = end_dvda_nvrs[a_idx, share_idx]
254 top_c = end_dvda_nvrs[a_idx, share_idx + 1]
255 bot_dvda = end_dvda[a_idx, share_idx]
256 top_dvda = end_dvda[a_idx, share_idx + 1]
257 alpha = 1.0 - top_f / (top_f - bot_f)
259 # Calculate the continuous optimal risky share and optimal consumption
260 Share_now = (1.0 - alpha) * bot_s + alpha * top_s
261 end_dvda_nvrs_now = (1.0 - alpha) * bot_c + alpha * top_c
262 end_dvda_now = (1.0 - alpha) * bot_dvda + alpha * top_dvda
264 # If agent wants to put more than 100% into risky asset, he is constrained.
265 # Likewise if he wants to put less than 0% into risky asset, he is constrained.
266 constrained_top = focs[:, -1] > 0.0
267 constrained_bot = focs[:, 0] < 0.0
269 # Apply those constraints to both risky share and consumption (but lower
270 # constraint should never be relevant)
271 Share_now[constrained_top] = 1.0
272 Share_now[constrained_bot] = 0.0
273 end_dvda_nvrs_now[constrained_top] = end_dvda_nvrs[constrained_top, -1]
274 end_dvda_nvrs_now[constrained_bot] = end_dvda_nvrs[constrained_bot, 0]
275 end_dvda_now[constrained_top] = end_dvda[constrained_top, -1]
276 end_dvda_now[constrained_bot] = end_dvda[constrained_bot, 0]
278 # When the natural borrowing constraint is *not* zero, then aNrm=0 is in the
279 # grid, but there's no way to "optimize" the portfolio if a=0, and consumption
280 # can't depend on the risky share if it doesn't meaningfully exist. Apply
281 # a small fix to the bottom gridpoint (aNrm=0) when this happens.
282 if not BoroCnstNat_iszero:
283 Share_now[0] = 1.0
284 end_dvda_nvrs_now[0] = end_dvda_nvrs[0, -1]
285 end_dvda_now[0] = end_dvda[0, -1]
287 # Now this is where we look for optimal C
288 # for each a in the agrid find corresponding c that satisfies the euler equation
290 if WealthShare == 0.0:
291 cNrm_now = end_dvda_nvrs_now
292 else:
293 omega = end_dvda_nvrs_now / (aNrmGrid + WealthShift)
294 cNrm_now = ChiFunc(omega) * (aNrmGrid + WealthShift)
296 # Calculate the endogenous mNrm gridpoints when the agent adjusts his portfolio,
297 # then construct the consumption function when the agent can adjust his share
298 mNrm_now = np.insert(aNrmGrid + cNrm_now, 0, 0.0)
299 cNrm_now = np.insert(cNrm_now, 0, 0.0)
300 cFuncNow = LinearInterp(mNrm_now, cNrm_now)
302 dudc_now = dudc(cNrm_now, mNrm_now - cNrm_now, CRRA, WealthShare, WealthShift)
303 dudc_nvrs_now = uFunc.derinv(dudc_now, order=(1, 0))
304 dudc_nvrs_func_now = LinearInterp(mNrm_now, dudc_nvrs_now)
306 # Construct the marginal value (of mNrm) function
307 vPfuncNow = MargValueFuncCRRA(dudc_nvrs_func_now, CRRA)
309 # If the share choice is continuous, just make an ordinary interpolating function
310 if BoroCnstNat_iszero:
311 Share_lower_bound = ShareLimit
312 else:
313 Share_lower_bound = 1.0
314 Share_now = np.insert(Share_now, 0, Share_lower_bound)
315 ShareFuncNow = LinearInterp(mNrm_now, Share_now, ShareLimit, 0.0)
317 # Add the value function if requested
318 if vFuncBool:
319 # Calculate intermediate value by taking expectations over income shocks
320 med_v = expected(
321 calc_med_v, IncShkDstn, args=(bNrmNext, PermGroFac, CRRA, v_func_next)
322 )
324 # Construct the "intermediate value function" for this period
325 med_v_nvrs = uFunc.inv(med_v)
326 med_v_nvrs_func = LinearInterp(bNrmGrid, med_v_nvrs)
327 med_v_func = ValueFuncCRRA(med_v_nvrs_func, CRRA)
329 # Calculate end-of-period value by taking expectations
330 end_v = DiscFacEff * expected(
331 calc_end_v,
332 RiskyDstn,
333 args=(aNrmNow, ShareNext, Rfree, med_v_func),
334 )
335 end_v_nvrs = uFunc.inv(end_v)
337 # Now make an end-of-period value function over aNrm and Share
338 end_v_nvrs_func = BilinearInterp(end_v_nvrs, aNrmGrid, ShareGrid)
339 end_v_func = ValueFuncCRRA(end_v_nvrs_func, CRRA)
340 # This will be used later to make the value function for this period
342 # Create the value functions for this period, defined over market resources
343 # mNrm when agent can adjust his portfolio, and over market resources and
344 # fixed share when agent can not adjust his portfolio.
346 # Construct the value function
347 mNrm_temp = aXtraGrid # Just use aXtraGrid as our grid of mNrm values
348 cNrm_temp = cFuncNow(mNrm_temp)
349 aNrm_temp = np.maximum(mNrm_temp - cNrm_temp, 0.0) # Fix tiny violations
350 Share_temp = ShareFuncNow(mNrm_temp)
351 v_temp = uFunc(cNrm_temp) + end_v_func(aNrm_temp, Share_temp)
352 vNvrs_temp = uFunc.inv(v_temp)
353 vNvrsP_temp = uFunc.der(cNrm_temp) * uFunc.inverse(v_temp, order=(0, 1))
354 vNvrsFunc = CubicInterp(
355 np.insert(mNrm_temp, 0, 0.0), # x_list
356 np.insert(vNvrs_temp, 0, 0.0), # f_list
357 np.insert(vNvrsP_temp, 0, vNvrsP_temp[0]), # dfdx_list
358 )
359 # Re-curve the pseudo-inverse value function
360 vFuncNow = ValueFuncCRRA(vNvrsFunc, CRRA)
362 else: # If vFuncBool is False, fill in dummy values
363 vFuncNow = NullFunc()
365 # Package and return the solution
366 solution_now = PortfolioSolution(
367 cFuncAdj=cFuncNow,
368 ShareFuncAdj=ShareFuncNow,
369 vPfuncAdj=vPfuncNow,
370 vFuncAdj=vFuncNow,
371 )
372 return solution_now
375###############################################################################
377# Make a dictionary of constructors for the wealth-in-utility portfolio choice consumer type
378WealthPortfolioConsumerType_constructors_default = {
379 "IncShkDstn": construct_lognormal_income_process_unemployment,
380 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
381 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
382 "aXtraGrid": make_assets_grid,
383 "RiskyDstn": make_lognormal_RiskyDstn,
384 "ShockDstn": combine_IncShkDstn_and_RiskyDstn,
385 "ShareLimit": calc_ShareLimit_for_CRRA,
386 "ShareGrid": make_simple_ShareGrid,
387 "ChiFunc": make_ChiFromOmega_function,
388 "AdjustDstn": make_AdjustDstn,
389 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
390 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
391 "solution_terminal": make_portfolio_solution_terminal,
392}
394# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
395WealthPortfolioConsumerType_IncShkDstn_default = {
396 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
397 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
398 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
399 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
400 "UnempPrb": 0.05, # Probability of unemployment while working
401 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
402 "T_retire": 0, # Period of retirement (0 --> no retirement)
403 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
404 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
405}
407# Default parameters to make aXtraGrid using make_assets_grid
408WealthPortfolioConsumerType_aXtraGrid_default = {
409 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
410 "aXtraMax": 100, # Maximum end-of-period "assets above minimum" value
411 "aXtraNestFac": 1, # Exponential nesting factor for aXtraGrid
412 "aXtraCount": 200, # Number of points in the grid of "assets above minimum"
413 "aXtraExtra": None, # Additional other values to add in grid (optional)
414}
416# Default parameters to make RiskyDstn with make_lognormal_RiskyDstn (and uniform ShareGrid)
417WealthPortfolioConsumerType_RiskyDstn_default = {
418 "RiskyAvg": 1.08, # Mean return factor of risky asset
419 "RiskyStd": 0.18362634887, # Stdev of log returns on risky asset
420 "RiskyCount": 5, # Number of integration nodes to use in approximation of risky returns
421}
423WealthPortfolioConsumerType_ShareGrid_default = {
424 "ShareCount": 25 # Number of discrete points in the risky share approximation
425}
427# Default parameters to make ChiFunc with make_ChiFromOmega_function
428WealthPortfolioConsumerType_ChiFunc_default = {
429 "ChiFromOmega_N": 501, # Number of gridpoints in chi-from-omega function
430 "ChiFromOmega_bound": 15, # Highest gridpoint to use for it
431}
433# Make a dictionary with parameters for the default constructor for kNrmInitDstn
434WealthPortfolioConsumerType_kNrmInitDstn_default = {
435 "kLogInitMean": -12.0, # Mean of log initial capital
436 "kLogInitStd": 0.0, # Stdev of log initial capital
437 "kNrmInitCount": 15, # Number of points in initial capital discretization
438}
440# Make a dictionary with parameters for the default constructor for pLvlInitDstn
441WealthPortfolioConsumerType_pLvlInitDstn_default = {
442 "pLogInitMean": 0.0, # Mean of log permanent income
443 "pLogInitStd": 0.0, # Stdev of log permanent income
444 "pLvlInitCount": 15, # Number of points in initial capital discretization
445}
447# Make a dictionary to specify a risky asset consumer type
448WealthPortfolioConsumerType_solving_default = {
449 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
450 "cycles": 1, # Finite, non-cyclic model
451 "T_cycle": 1, # Number of periods in the cycle for this agent type
452 "constructors": WealthPortfolioConsumerType_constructors_default, # See dictionary above
453 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
454 "CRRA": 5.0, # Coefficient of relative risk aversion
455 "Rfree": [1.03], # Return factor on risk free asset
456 "DiscFac": 0.90, # Intertemporal discount factor
457 "LivPrb": [0.98], # Survival probability after each period
458 "PermGroFac": [1.01], # Permanent income growth factor
459 "BoroCnstArt": 0.0, # Artificial borrowing constraint
460 "WealthShare": 0.5, # Share of wealth in Cobb-Douglas aggregator in utility function
461 "WealthShift": 0.1, # Shifter for wealth in utility function
462 "DiscreteShareBool": False, # Whether risky asset share is restricted to discrete values
463 "PortfolioBool": True, # Whether there is portfolio choice
464 "PortfolioBisect": False, # This is a mystery parameter
465 "IndepDstnBool": True, # Whether income and return shocks are independent
466 "vFuncBool": False, # Whether to calculate the value function during solution
467 "CubicBool": False, # Whether to use cubic spline interpolation when True
468 # (Uses linear spline interpolation for cFunc when False)
469 "AdjustPrb": 1.0, # Probability that the agent can update their risky portfolio share each period
470 "RiskyShareFixed": None, # This just needs to exist because of inheritance, does nothing
471 "sim_common_Rrisky": True, # Whether risky returns have a shared/common value across agents
472}
473WealthPortfolioConsumerType_simulation_default = {
474 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
475 "AgentCount": 10000, # Number of agents of this type
476 "T_age": None, # Age after which simulated agents are automatically killed
477 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
478 # (The portion of PermGroFac attributable to aggregate productivity growth)
479 "NewbornTransShk": False, # Whether Newborns have transitory shock
480 # ADDITIONAL OPTIONAL PARAMETERS
481 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
482 # (Forces Newborns to follow solution path of the agent they replaced if True)
483 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
484}
486# Assemble the default dictionary
487WealthPortfolioConsumerType_default = {}
488WealthPortfolioConsumerType_default.update(WealthPortfolioConsumerType_solving_default)
489WealthPortfolioConsumerType_default.update(
490 WealthPortfolioConsumerType_simulation_default
491)
492WealthPortfolioConsumerType_default.update(
493 WealthPortfolioConsumerType_aXtraGrid_default
494)
495WealthPortfolioConsumerType_default.update(
496 WealthPortfolioConsumerType_ShareGrid_default
497)
498WealthPortfolioConsumerType_default.update(
499 WealthPortfolioConsumerType_IncShkDstn_default
500)
501WealthPortfolioConsumerType_default.update(
502 WealthPortfolioConsumerType_RiskyDstn_default
503)
504WealthPortfolioConsumerType_default.update(WealthPortfolioConsumerType_ChiFunc_default)
505WealthPortfolioConsumerType_default.update(
506 WealthPortfolioConsumerType_kNrmInitDstn_default
507)
508WealthPortfolioConsumerType_default.update(
509 WealthPortfolioConsumerType_pLvlInitDstn_default
510)
511init_wealth_portfolio = WealthPortfolioConsumerType_default
513###############################################################################
516class WealthPortfolioConsumerType(PortfolioConsumerType):
517 """
518 TODO: This docstring is missing and needs to be written.
519 """
521 time_inv_ = deepcopy(PortfolioConsumerType.time_inv_)
522 time_inv_ = time_inv_ + [
523 "WealthShare",
524 "WealthShift",
525 "ChiFunc",
526 "RiskyDstn",
527 ]
528 default_ = {
529 "params": init_wealth_portfolio,
530 "solver": solve_one_period_WealthPortfolio,
531 "model": "ConsRiskyAsset.yaml",
532 "track_vars": ["aNrm", "cNrm", "mNrm", "Share", "pLvl"],
533 }
535 def pre_solve(self):
536 self.construct("solution_terminal")
537 self.solution_terminal.ShareFunc = ConstantFunction(1.0)