Coverage for HARK / ConsumptionSaving / ConsAggShockModel.py: 97%
925 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-25 05:22 +0000
1"""
2Consumption-saving models with aggregate productivity shocks as well as idiosyn-
3cratic income shocks. Currently only contains one microeconomic model with a
4basic solver. Also includes a subclass of Market called CobbDouglas economy,
5used for solving "macroeconomic" models with aggregate shocks.
6"""
8from copy import deepcopy
10import numpy as np
11import scipy.stats as stats
13from HARK import AgentType, Market
14from HARK.Calibration.Income.IncomeProcesses import (
15 construct_lognormal_income_process_unemployment,
16 get_PermShkDstn_from_IncShkDstn,
17 get_TranShkDstn_from_IncShkDstn,
18)
19from HARK.ConsumptionSaving.ConsIndShockModel import (
20 ConsumerSolution,
21 IndShockConsumerType,
22 make_lognormal_kNrm_init_dstn,
23 make_lognormal_pLvl_init_dstn,
24)
25from HARK.distributions import (
26 MarkovProcess,
27 MeanOneLogNormal,
28 Uniform,
29 calc_expectation,
30 combine_indep_dstns,
31)
32from HARK.interpolation import (
33 BilinearInterp,
34 ConstantFunction,
35 IdentityFunction,
36 LinearInterp,
37 LinearInterpOnInterp1D,
38 LowerEnvelope2D,
39 MargValueFuncCRRA,
40 UpperEnvelope,
41 VariableLowerBoundFunc2D,
42)
43from HARK.metric import MetricObject
44from HARK.rewards import (
45 CRRAutility,
46 CRRAutility_inv,
47 CRRAutility_invP,
48 CRRAutilityP,
49 CRRAutilityP_inv,
50 CRRAutilityPP,
51)
52from HARK.utilities import make_assets_grid, get_it_from, NullFunc
54__all__ = [
55 "AggShockConsumerType",
56 "AggShockMarkovConsumerType",
57 "CobbDouglasEconomy",
58 "SmallOpenEconomy",
59 "CobbDouglasMarkovEconomy",
60 "SmallOpenMarkovEconomy",
61 "AggregateSavingRule",
62 "AggShocksDynamicRule",
63 "init_agg_shocks",
64 "init_agg_mrkv_shocks",
65 "init_cobb_douglas",
66 "init_mrkv_cobb_douglas",
67]
69utility = CRRAutility
70utilityP = CRRAutilityP
71utilityPP = CRRAutilityPP
72utilityP_inv = CRRAutilityP_inv
73utility_invP = CRRAutility_invP
74utility_inv = CRRAutility_inv
77def make_aggshock_solution_terminal(CRRA):
78 """
79 Creates the terminal period solution for an aggregate shock consumer.
80 Only fills in the consumption function and marginal value function.
82 Parameters
83 ----------
84 CRRA : float
85 Coefficient of relative risk aversion.
87 Returns
88 -------
89 solution_terminal : ConsumerSolution
90 Solution to the terminal period problem.
91 """
92 cFunc_terminal = IdentityFunction(i_dim=0, n_dims=2)
93 vPfunc_terminal = MargValueFuncCRRA(cFunc_terminal, CRRA)
94 mNrmMin_terminal = ConstantFunction(0)
95 solution_terminal = ConsumerSolution(
96 cFunc=cFunc_terminal, vPfunc=vPfunc_terminal, mNrmMin=mNrmMin_terminal
97 )
98 return solution_terminal
101def make_aggmrkv_solution_terminal(CRRA, MrkvArray):
102 """
103 Creates the terminal period solution for an aggregate shock consumer with
104 discrete Markov state. Only fills in the consumption function and marginal
105 value function.
107 Parameters
108 ----------
109 CRRA : float
110 Coefficient of relative risk aversion.
111 MrkvArray : np.array
112 Transition probability array.
114 Returns
115 -------
116 solution_terminal : ConsumerSolution
117 Solution to the terminal period problem.
118 """
119 solution_terminal = make_aggshock_solution_terminal(CRRA)
121 # Make replicated terminal period solution
122 StateCount = MrkvArray.shape[0]
123 solution_terminal.cFunc = StateCount * [solution_terminal.cFunc]
124 solution_terminal.vPfunc = StateCount * [solution_terminal.vPfunc]
125 solution_terminal.mNrmMin = StateCount * [solution_terminal.mNrmMin]
127 return solution_terminal
130def make_exponential_MgridBase(MaggCount, MaggPerturb, MaggExpFac):
131 """
132 Constructor function for MgridBase, the grid of aggregate market resources
133 relative to the steady state. This grid is always centered around 1.0.
135 Parameters
136 ----------
137 MaggCount : int
138 Number of gridpoints for aggregate market resources. Should be odd.
139 MaggPerturb : float
140 Small perturbation around the steady state; the grid will always include
141 1+perturb and 1-perturb.
142 MaggExpFac : float
143 Log growth factor for gridpoints beyond the two adjacent to the steady state.
145 Returns
146 -------
147 MgridBase : np.array
148 Grid of aggregate market resources relative to the steady state.
149 """
150 N = int((MaggCount - 1) / 2)
151 gridpoints = [1.0 - MaggPerturb, 1.0, 1.0 + MaggPerturb]
152 fac = np.exp(MaggExpFac)
153 for n in range(N - 1):
154 new_hi = gridpoints[-1] * fac
155 new_lo = gridpoints[0] / fac
156 gridpoints.append(new_hi)
157 gridpoints.insert(0, new_lo)
158 MgridBase = np.array(gridpoints)
159 return MgridBase
162###############################################################################
165def solveConsAggShock(
166 solution_next,
167 IncShkDstn,
168 LivPrb,
169 DiscFac,
170 CRRA,
171 PermGroFac,
172 PermGroFacAgg,
173 aXtraGrid,
174 BoroCnstArt,
175 Mgrid,
176 AFunc,
177 Rfunc,
178 wFunc,
179 DeprRte,
180):
181 """
182 Solve one period of a consumption-saving problem with idiosyncratic and
183 aggregate shocks (transitory and permanent). This is a basic solver that
184 can't handle cubic splines, nor can it calculate a value function. This
185 version uses calc_expectation to reduce code clutter.
187 Parameters
188 ----------
189 solution_next : ConsumerSolution
190 The solution to the succeeding one period problem.
191 IncShkDstn : distribution.Distribution
192 A discrete
193 approximation to the income process between the period being solved
194 and the one immediately following (in solution_next). Order:
195 idiosyncratic permanent shocks, idiosyncratic transitory
196 shocks, aggregate permanent shocks, aggregate transitory shocks.
197 LivPrb : float
198 Survival probability; likelihood of being alive at the beginning of
199 the succeeding period.
200 DiscFac : float
201 Intertemporal discount factor for future utility.
202 CRRA : float
203 Coefficient of relative risk aversion.
204 PermGroFac : float
205 Expected permanent income growth factor at the end of this period.
206 PermGroFacAgg : float
207 Expected aggregate productivity growth factor.
208 aXtraGrid : np.array
209 Array of "extra" end-of-period asset values-- assets above the
210 absolute minimum acceptable level.
211 BoroCnstArt : float
212 Artificial borrowing constraint; minimum allowable end-of-period asset-to-
213 permanent-income ratio. Unlike other models, this *can't* be None.
214 Mgrid : np.array
215 A grid of aggregate market resourses to permanent income in the economy.
216 AFunc : function
217 Aggregate savings as a function of aggregate market resources.
218 Rfunc : function
219 The net interest factor on assets as a function of capital ratio k.
220 wFunc : function
221 The wage rate for labor as a function of capital-to-labor ratio k.
222 DeprRte : float
223 Capital Depreciation Rate
225 Returns
226 -------
227 solution_now : ConsumerSolution
228 The solution to the single period consumption-saving problem. Includes
229 a consumption function cFunc (linear interpolation over linear interpola-
230 tions) and marginal value function vPfunc.
231 """
232 # Unpack the income shocks and get grid sizes
233 PermShkValsNext = IncShkDstn.atoms[0]
234 TranShkValsNext = IncShkDstn.atoms[1]
235 PermShkAggValsNext = IncShkDstn.atoms[2]
236 TranShkAggValsNext = IncShkDstn.atoms[3]
237 aCount = aXtraGrid.size
238 Mcount = Mgrid.size
240 # Define a function that calculates M_{t+1} from M_t and the aggregate shocks;
241 # the function also returns the wage rate and effective interest factor
242 def calcAggObjects(M, Psi, Theta):
243 A = AFunc(M) # End-of-period aggregate assets (normalized)
244 # Next period's aggregate capital/labor ratio
245 kNext = A / (PermGroFacAgg * Psi)
246 kNextEff = kNext / Theta # Same thing, but account for *transitory* shock
247 R = Rfunc(kNextEff) # Interest factor on aggregate assets
248 wEff = (
249 wFunc(kNextEff) * Theta
250 ) # Effective wage rate (accounts for labor supply)
251 Reff = R / LivPrb # Account for redistribution of decedents' wealth
252 Mnext = kNext * R + wEff # Next period's aggregate market resources
253 return Mnext, Reff, wEff
255 # Define a function that evaluates R*v'(m_{t+1},M_{t+1}) from a_t, M_t, and the income shocks
256 def vPnextFunc(S, a, M):
257 psi = S[0]
258 theta = S[1]
259 Psi = S[2]
260 Theta = S[3]
262 Mnext, Reff, wEff = calcAggObjects(M, Psi, Theta)
263 PermShkTotal = (
264 PermGroFac * PermGroFacAgg * psi * Psi
265 ) # Total / combined permanent shock
266 mNext = Reff * a / PermShkTotal + theta * wEff # Idiosyncratic market resources
267 vPnext = Reff * PermShkTotal ** (-CRRA) * solution_next.vPfunc(mNext, Mnext)
268 return vPnext
270 # Make an array of a_t values at which to calculate end-of-period marginal value of assets
271 # Natural borrowing constraint at each M_t
272 BoroCnstNat_vec = np.zeros(Mcount)
273 aNrmNow = np.zeros((aCount, Mcount))
274 for j in range(Mcount):
275 Mnext, Reff, wEff = calcAggObjects(
276 Mgrid[j], PermShkAggValsNext, TranShkAggValsNext
277 )
278 aNrmMin_cand = (
279 PermGroFac * PermGroFacAgg * PermShkValsNext * PermShkAggValsNext / Reff
280 ) * (solution_next.mNrmMin(Mnext) - wEff * TranShkValsNext)
281 aNrmMin = np.max(aNrmMin_cand) # Lowest valid a_t value for this M_t
282 aNrmNow[:, j] = aNrmMin + aXtraGrid
283 BoroCnstNat_vec[j] = aNrmMin
285 # Compute end-of-period marginal value of assets
286 MaggNow = np.tile(np.reshape(Mgrid, (1, Mcount)), (aCount, 1)) # Tiled Mgrid
287 EndOfPrdvP = (
288 DiscFac * LivPrb * calc_expectation(IncShkDstn, vPnextFunc, *(aNrmNow, MaggNow))
289 )
291 # Calculate optimal consumption from each asset gridpoint and endogenous m_t gridpoint
292 cNrmNow = EndOfPrdvP ** (-1.0 / CRRA)
293 mNrmNow = aNrmNow + cNrmNow
295 # Loop through the values in Mgrid and make a linear consumption function for each
296 cFuncBaseByM_list = []
297 for j in range(Mcount):
298 c_temp = np.insert(cNrmNow[:, j], 0, 0.0) # Add point at bottom
299 m_temp = np.insert(mNrmNow[:, j] - BoroCnstNat_vec[j], 0, 0.0)
300 cFuncBaseByM_list.append(LinearInterp(m_temp, c_temp))
302 # Construct the overall unconstrained consumption function by combining the M-specific functions
303 BoroCnstNat = LinearInterp(
304 np.insert(Mgrid, 0, 0.0), np.insert(BoroCnstNat_vec, 0, 0.0)
305 )
306 cFuncBase = LinearInterpOnInterp1D(cFuncBaseByM_list, Mgrid)
307 cFuncUnc = VariableLowerBoundFunc2D(cFuncBase, BoroCnstNat)
309 # Make the constrained consumption function and combine it with the unconstrained component
310 cFuncCnst = BilinearInterp(
311 np.array([[0.0, 0.0], [1.0, 1.0]]),
312 np.array([BoroCnstArt, BoroCnstArt + 1.0]),
313 np.array([0.0, 1.0]),
314 )
315 cFuncNow = LowerEnvelope2D(cFuncUnc, cFuncCnst)
317 # Make the minimum m function as the greater of the natural and artificial constraints
318 mNrmMinNow = UpperEnvelope(BoroCnstNat, ConstantFunction(BoroCnstArt))
320 # Construct the marginal value function using the envelope condition
321 vPfuncNow = MargValueFuncCRRA(cFuncNow, CRRA)
323 # Pack up and return the solution
324 solution_now = ConsumerSolution(
325 cFunc=cFuncNow, vPfunc=vPfuncNow, mNrmMin=mNrmMinNow
326 )
327 return solution_now
330###############################################################################
333def solve_ConsAggMarkov(
334 solution_next,
335 IncShkDstn,
336 LivPrb,
337 DiscFac,
338 CRRA,
339 MrkvArray,
340 PermGroFac,
341 PermGroFacAgg,
342 aXtraGrid,
343 BoroCnstArt,
344 Mgrid,
345 AFunc,
346 Rfunc,
347 wFunc,
348):
349 """
350 Solve one period of a consumption-saving problem with idiosyncratic and
351 aggregate shocks (transitory and permanent). Moreover, the macroeconomic
352 state follows a Markov process that determines the income distribution and
353 aggregate permanent growth factor. This is a basic solver that can't handle
354 cubic splines, nor can it calculate a value function.
356 Parameters
357 ----------
358 solution_next : ConsumerSolution
359 The solution to the succeeding one period problem.
360 IncShkDstn : [distribution.Distribution]
361 A list of
362 discrete approximations to the income process between the period being
363 solved and the one immediately following (in solution_next). Order:
364 idisyncratic permanent shocks, idiosyncratic transitory
365 shocks, aggregate permanent shocks, aggregate transitory shocks.
366 LivPrb : float
367 Survival probability; likelihood of being alive at the beginning of
368 the succeeding period.
369 DiscFac : float
370 Intertemporal discount factor for future utility.
371 CRRA : float
372 Coefficient of relative risk aversion.
373 MrkvArray : np.array
374 Markov transition matrix between discrete macroeconomic states.
375 MrkvArray[i,j] is probability of being in state j next period conditional
376 on being in state i this period.
377 PermGroFac : float
378 Expected permanent income growth factor at the end of this period,
379 for the *individual*'s productivity.
380 PermGroFacAgg : [float]
381 Expected aggregate productivity growth in each Markov macro state.
382 aXtraGrid : np.array
383 Array of "extra" end-of-period asset values-- assets above the
384 absolute minimum acceptable level.
385 BoroCnstArt : float
386 Artificial borrowing constraint; minimum allowable end-of-period asset-to-
387 permanent-income ratio. Unlike other models, this *can't* be None.
388 Mgrid : np.array
389 A grid of aggregate market resourses to permanent income in the economy.
390 AFunc : [function]
391 Aggregate savings as a function of aggregate market resources, for each
392 Markov macro state.
393 Rfunc : function
394 The net interest factor on assets as a function of capital ratio k.
395 wFunc : function
396 The wage rate for labor as a function of capital-to-labor ratio k.
397 DeprRte : float
398 Capital Depreciation Rate
400 Returns
401 -------
402 solution_now : ConsumerSolution
403 The solution to the single period consumption-saving problem. Includes
404 a consumption function cFunc (linear interpolation over linear interpola-
405 tions) and marginal value function vPfunc.
406 """
407 # Get sizes of grids
408 aCount = aXtraGrid.size
409 Mcount = Mgrid.size
410 StateCount = MrkvArray.shape[0]
412 # Loop through next period's states, assuming we reach each one at a time.
413 # Construct EndOfPrdvP_cond functions for each state.
414 EndOfPrdvPfunc_cond = []
415 BoroCnstNat_cond = []
416 for j in range(StateCount):
417 # Unpack next period's solution
418 vPfuncNext = solution_next.vPfunc[j]
419 mNrmMinNext = solution_next.mNrmMin[j]
421 # Unpack the income shocks
422 ShkPrbsNext = IncShkDstn[j].pmv
423 PermShkValsNext = IncShkDstn[j].atoms[0]
424 TranShkValsNext = IncShkDstn[j].atoms[1]
425 PermShkAggValsNext = IncShkDstn[j].atoms[2]
426 TranShkAggValsNext = IncShkDstn[j].atoms[3]
427 ShkCount = ShkPrbsNext.size
428 aXtra_tiled = np.tile(
429 np.reshape(aXtraGrid, (1, aCount, 1)), (Mcount, 1, ShkCount)
430 )
432 # Make tiled versions of the income shocks
433 # Dimension order: Mnow, aNow, Shk
434 ShkPrbsNext_tiled = np.tile(
435 np.reshape(ShkPrbsNext, (1, 1, ShkCount)), (Mcount, aCount, 1)
436 )
437 PermShkValsNext_tiled = np.tile(
438 np.reshape(PermShkValsNext, (1, 1, ShkCount)), (Mcount, aCount, 1)
439 )
440 TranShkValsNext_tiled = np.tile(
441 np.reshape(TranShkValsNext, (1, 1, ShkCount)), (Mcount, aCount, 1)
442 )
443 PermShkAggValsNext_tiled = np.tile(
444 np.reshape(PermShkAggValsNext, (1, 1, ShkCount)), (Mcount, aCount, 1)
445 )
446 TranShkAggValsNext_tiled = np.tile(
447 np.reshape(TranShkAggValsNext, (1, 1, ShkCount)), (Mcount, aCount, 1)
448 )
450 # Make a tiled grid of end-of-period aggregate assets. These lines use
451 # next prd state j's aggregate saving rule to get a relevant set of Aagg,
452 # which will be used to make an interpolated EndOfPrdvP_cond function.
453 # After constructing these functions, we will use the aggregate saving
454 # rule for *current* state i to get values of Aagg at which to evaluate
455 # these conditional marginal value functions. In the strange, maybe even
456 # impossible case where the aggregate saving rules differ wildly across
457 # macro states *and* there is "anti-persistence", so that the macro state
458 # is very likely to change each period, then this procedure will lead to
459 # an inaccurate solution because the grid of Aagg values on which the
460 # conditional marginal value functions are constructed is not relevant
461 # to the values at which it will actually be evaluated.
462 AaggGrid = AFunc[j](Mgrid)
463 AaggNow_tiled = np.tile(
464 np.reshape(AaggGrid, (Mcount, 1, 1)), (1, aCount, ShkCount)
465 )
467 # Calculate returns to capital and labor in the next period
468 kNext_array = AaggNow_tiled / (
469 PermGroFacAgg[j] * PermShkAggValsNext_tiled
470 ) # Next period's aggregate capital to labor ratio
471 kNextEff_array = (
472 kNext_array / TranShkAggValsNext_tiled
473 ) # Same thing, but account for *transitory* shock
474 R_array = Rfunc(kNextEff_array) # Interest factor on aggregate assets
475 Reff_array = (
476 R_array / LivPrb
477 ) # Effective interest factor on individual assets *for survivors*
478 wEff_array = (
479 wFunc(kNextEff_array) * TranShkAggValsNext_tiled
480 ) # Effective wage rate (accounts for labor supply)
481 PermShkTotal_array = (
482 PermGroFac
483 * PermGroFacAgg[j]
484 * PermShkValsNext_tiled
485 * PermShkAggValsNext_tiled
486 ) # total / combined permanent shock
487 Mnext_array = (
488 kNext_array * R_array + wEff_array
489 ) # next period's aggregate market resources
491 # Find the natural borrowing constraint for each value of M in the Mgrid.
492 # There is likely a faster way to do this, but someone needs to do the math:
493 # is aNrmMin determined by getting the worst shock of all four types?
494 aNrmMin_candidates = (
495 PermGroFac
496 * PermGroFacAgg[j]
497 * PermShkValsNext_tiled[:, 0, :]
498 * PermShkAggValsNext_tiled[:, 0, :]
499 / Reff_array[:, 0, :]
500 * (
501 mNrmMinNext(Mnext_array[:, 0, :])
502 - wEff_array[:, 0, :] * TranShkValsNext_tiled[:, 0, :]
503 )
504 )
505 aNrmMin_vec = np.max(aNrmMin_candidates, axis=1)
506 BoroCnstNat_vec = aNrmMin_vec
507 aNrmMin_tiled = np.tile(
508 np.reshape(aNrmMin_vec, (Mcount, 1, 1)), (1, aCount, ShkCount)
509 )
510 aNrmNow_tiled = aNrmMin_tiled + aXtra_tiled
512 # Calculate market resources next period (and a constant array of capital-to-labor ratio)
513 mNrmNext_array = (
514 Reff_array * aNrmNow_tiled / PermShkTotal_array
515 + TranShkValsNext_tiled * wEff_array
516 )
518 # Find marginal value next period at every income shock
519 # realization and every aggregate market resource gridpoint
520 vPnext_array = (
521 Reff_array
522 * PermShkTotal_array ** (-CRRA)
523 * vPfuncNext(mNrmNext_array, Mnext_array)
524 )
526 # Calculate expectated marginal value at the end of the period at every asset gridpoint
527 EndOfPrdvP = DiscFac * LivPrb * np.sum(vPnext_array * ShkPrbsNext_tiled, axis=2)
529 # Make the conditional end-of-period marginal value function
530 BoroCnstNat = LinearInterp(
531 np.insert(AaggGrid, 0, 0.0), np.insert(BoroCnstNat_vec, 0, 0.0)
532 )
533 EndOfPrdvPnvrs = np.concatenate(
534 (np.zeros((Mcount, 1)), EndOfPrdvP ** (-1.0 / CRRA)), axis=1
535 )
536 EndOfPrdvPnvrsFunc_base = BilinearInterp(
537 np.transpose(EndOfPrdvPnvrs), np.insert(aXtraGrid, 0, 0.0), AaggGrid
538 )
539 EndOfPrdvPnvrsFunc = VariableLowerBoundFunc2D(
540 EndOfPrdvPnvrsFunc_base, BoroCnstNat
541 )
542 EndOfPrdvPfunc_cond.append(MargValueFuncCRRA(EndOfPrdvPnvrsFunc, CRRA))
543 BoroCnstNat_cond.append(BoroCnstNat)
545 # Prepare some objects that are the same across all current states
546 aXtra_tiled = np.tile(np.reshape(aXtraGrid, (1, aCount)), (Mcount, 1))
547 cFuncCnst = BilinearInterp(
548 np.array([[0.0, 0.0], [1.0, 1.0]]),
549 np.array([BoroCnstArt, BoroCnstArt + 1.0]),
550 np.array([0.0, 1.0]),
551 )
553 # Now loop through *this* period's discrete states, calculating end-of-period
554 # marginal value (weighting across state transitions), then construct consumption
555 # and marginal value function for each state.
556 cFuncNow = []
557 vPfuncNow = []
558 mNrmMinNow = []
559 for i in range(StateCount):
560 # Find natural borrowing constraint for this state by Aagg
561 AaggNow = AFunc[i](Mgrid)
562 aNrmMin_candidates = np.zeros((StateCount, Mcount)) + np.nan
563 for j in range(StateCount):
564 if MrkvArray[i, j] > 0.0: # Irrelevant if transition is impossible
565 aNrmMin_candidates[j, :] = BoroCnstNat_cond[j](AaggNow)
566 aNrmMin_vec = np.nanmax(aNrmMin_candidates, axis=0)
567 BoroCnstNat_vec = aNrmMin_vec
569 # Make tiled grids of aNrm and Aagg
570 aNrmMin_tiled = np.tile(np.reshape(aNrmMin_vec, (Mcount, 1)), (1, aCount))
571 aNrmNow_tiled = aNrmMin_tiled + aXtra_tiled
572 AaggNow_tiled = np.tile(np.reshape(AaggNow, (Mcount, 1)), (1, aCount))
574 # Loop through feasible transitions and calculate end-of-period marginal value
575 EndOfPrdvP = np.zeros((Mcount, aCount))
576 for j in range(StateCount):
577 if MrkvArray[i, j] > 0.0:
578 temp = EndOfPrdvPfunc_cond[j](aNrmNow_tiled, AaggNow_tiled)
579 EndOfPrdvP += MrkvArray[i, j] * temp
581 # Calculate consumption and the endogenous mNrm gridpoints for this state
582 cNrmNow = EndOfPrdvP ** (-1.0 / CRRA)
583 mNrmNow = aNrmNow_tiled + cNrmNow
585 # Loop through the values in Mgrid and make a piecewise linear consumption function for each
586 cFuncBaseByM_list = []
587 for n in range(Mcount):
588 c_temp = np.insert(cNrmNow[n, :], 0, 0.0) # Add point at bottom
589 m_temp = np.insert(mNrmNow[n, :] - BoroCnstNat_vec[n], 0, 0.0)
590 cFuncBaseByM_list.append(LinearInterp(m_temp, c_temp))
591 # Add the M-specific consumption function to the list
593 # Construct the unconstrained consumption function by combining the M-specific functions
594 BoroCnstNat = LinearInterp(
595 np.insert(Mgrid, 0, 0.0), np.insert(BoroCnstNat_vec, 0, 0.0)
596 )
597 cFuncBase = LinearInterpOnInterp1D(cFuncBaseByM_list, Mgrid)
598 cFuncUnc = VariableLowerBoundFunc2D(cFuncBase, BoroCnstNat)
600 # Combine the constrained consumption function with unconstrained component
601 cFuncNow.append(LowerEnvelope2D(cFuncUnc, cFuncCnst))
603 # Make the minimum m function as the greater of the natural and artificial constraints
604 mNrmMinNow.append(UpperEnvelope(BoroCnstNat, ConstantFunction(BoroCnstArt)))
606 # Construct the marginal value function using the envelope condition
607 vPfuncNow.append(MargValueFuncCRRA(cFuncNow[-1], CRRA))
609 # Pack up and return the solution
610 solution_now = ConsumerSolution(
611 cFunc=cFuncNow, vPfunc=vPfuncNow, mNrmMin=mNrmMinNow
612 )
613 return solution_now
616###############################################################################
619def solve_KrusellSmith(
620 solution_next,
621 DiscFac,
622 CRRA,
623 aGrid,
624 Mgrid,
625 mNextArray,
626 MnextArray,
627 ProbArray,
628 RnextArray,
629):
630 """
631 Solve the one period problem of an agent in Krusell & Smith's canonical 1998 model.
632 Because this model is so specialized and only intended to be used with a very narrow
633 case, many arrays can be precomputed, making the code here very short. See the
634 default constructor functions for details.
636 Parameters
637 ----------
638 solution_next : ConsumerSolution
639 Representation of the solution to next period's problem, including the
640 discrete-state-conditional consumption function and marginal value function.
641 DiscFac : float
642 Intertemporal discount factor.
643 CRRA : float
644 Coefficient of relative risk aversion.
645 aGrid : np.array
646 Array of end-of-period asset values.
647 Mgrid : np.array
648 A grid of aggregate market resources in the economy.
649 mNextArray : np.array
650 Precomputed array of next period's market resources attained from every
651 end-of-period state in the exogenous grid crossed with every shock that
652 might attain. Has shape [aCount, Mcount, 4, 4] ~ [a, M, s, s'].
653 MnextArray : np.array
654 Precomputed array of next period's aggregate market resources attained
655 from every end-of-period state in the exogenous grid crossed with every
656 shock that might attain. Corresponds to mNextArray.
657 ProbArray : np.array
658 Tiled array of transition probabilities among discrete states. Every
659 slice [i,j,:,:] is identical and translated from MrkvIndArray.
660 RnextArray : np.array
661 Tiled array of net interest factors next period, attained from every
662 end-of-period state crossed with every shock that might attain.
664 Returns
665 -------
666 solution_now : ConsumerSolution
667 Representation of this period's solution to the Krusell-Smith model.
668 """
669 # Loop over next period's state realizations, computing marginal value of market resources
670 vPnext = np.zeros_like(mNextArray)
671 for j in range(4):
672 vPnext[:, :, :, j] = solution_next.vPfunc[j](
673 mNextArray[:, :, :, j], MnextArray[:, :, :, j]
674 )
676 # Compute end-of-period marginal value of assets
677 EndOfPrdvP = DiscFac * np.sum(RnextArray * vPnext * ProbArray, axis=3)
679 # Invert the first order condition to find optimal consumption
680 cNow = EndOfPrdvP ** (-1.0 / CRRA)
682 # Find the endogenous gridpoints
683 aCount = aGrid.size
684 Mcount = Mgrid.size
685 aNow = np.tile(np.reshape(aGrid, [aCount, 1, 1]), [1, Mcount, 4])
686 mNow = aNow + cNow
688 # Insert zeros at the bottom of both cNow and mNow arrays (consume nothing)
689 cNow = np.concatenate([np.zeros([1, Mcount, 4]), cNow], axis=0)
690 mNow = np.concatenate([np.zeros([1, Mcount, 4]), mNow], axis=0)
692 # Construct the consumption and marginal value function for each discrete state
693 cFunc_by_state = []
694 vPfunc_by_state = []
695 for j in range(4):
696 cFunc_by_M = [LinearInterp(mNow[:, k, j], cNow[:, k, j]) for k in range(Mcount)]
697 cFunc_j = LinearInterpOnInterp1D(cFunc_by_M, Mgrid)
698 vPfunc_j = MargValueFuncCRRA(cFunc_j, CRRA)
699 cFunc_by_state.append(cFunc_j)
700 vPfunc_by_state.append(vPfunc_j)
702 # Package and return the solution
703 solution_now = ConsumerSolution(cFunc=cFunc_by_state, vPfunc=vPfunc_by_state)
704 return solution_now
707###############################################################################
709# Make a dictionary of constructors for the aggregate income shocks model
710aggshock_constructor_dict = {
711 "IncShkDstn": construct_lognormal_income_process_unemployment,
712 "PermShkDstn": get_PermShkDstn_from_IncShkDstn,
713 "TranShkDstn": get_TranShkDstn_from_IncShkDstn,
714 "aXtraGrid": make_assets_grid,
715 "MgridBase": make_exponential_MgridBase,
716 "kNrmInitDstn": make_lognormal_kNrm_init_dstn,
717 "pLvlInitDstn": make_lognormal_pLvl_init_dstn,
718 "solution_terminal": make_aggshock_solution_terminal,
719}
721# Make a dictionary with parameters for the default constructor for kNrmInitDstn
722default_kNrmInitDstn_params = {
723 "kLogInitMean": 0.0, # Mean of log initial capital
724 "kLogInitStd": 0.0, # Stdev of log initial capital
725 "kNrmInitCount": 15, # Number of points in initial capital discretization
726}
728# Make a dictionary with parameters for the default constructor for pLvlInitDstn
729default_pLvlInitDstn_params = {
730 "pLogInitMean": 0.0, # Mean of log permanent income
731 "pLogInitStd": 0.0, # Stdev of log permanent income
732 "pLvlInitCount": 15, # Number of points in initial capital discretization
733}
736# Default parameters to make IncShkDstn using construct_lognormal_income_process_unemployment
737default_IncShkDstn_params = {
738 "PermShkStd": [0.1], # Standard deviation of log permanent income shocks
739 "PermShkCount": 7, # Number of points in discrete approximation to permanent income shocks
740 "TranShkStd": [0.1], # Standard deviation of log transitory income shocks
741 "TranShkCount": 7, # Number of points in discrete approximation to transitory income shocks
742 "UnempPrb": 0.05, # Probability of unemployment while working
743 "IncUnemp": 0.3, # Unemployment benefits replacement rate while working
744 "T_retire": 0, # Period of retirement (0 --> no retirement)
745 "UnempPrbRet": 0.005, # Probability of "unemployment" while retired
746 "IncUnempRet": 0.0, # "Unemployment" benefits when retired
747}
749# Default parameters to make aXtraGrid using make_assets_grid
750default_aXtraGrid_params = {
751 "aXtraMin": 0.001, # Minimum end-of-period "assets above minimum" value
752 "aXtraMax": 20, # Maximum end-of-period "assets above minimum" value
753 "aXtraNestFac": 3, # Exponential nesting factor for aXtraGrid
754 "aXtraCount": 24, # Number of points in the grid of "assets above minimum"
755 "aXtraExtra": None, # Additional other values to add in grid (optional)
756}
758# Default parameters to make MgridBase using make_exponential_MgridBase
759default_MgridBase_params = {
760 "MaggCount": 17,
761 "MaggPerturb": 0.01,
762 "MaggExpFac": 0.15,
763}
765# Make a dictionary to specify an aggregate income shocks consumer type
766init_agg_shocks = {
767 # BASIC HARK PARAMETERS REQUIRED TO SOLVE THE MODEL
768 "cycles": 1, # Finite, non-cyclic model
769 "T_cycle": 1, # Number of periods in the cycle for this agent type
770 "constructors": aggshock_constructor_dict, # See dictionary above
771 "pseudo_terminal": False, # Terminal period is real
772 # PRIMITIVE RAW PARAMETERS REQUIRED TO SOLVE THE MODEL
773 "CRRA": 2.0, # Coefficient of relative risk aversion
774 "DiscFac": 0.96, # Intertemporal discount factor
775 "LivPrb": [0.98], # Survival probability after each period
776 "PermGroFac": [1.00], # Permanent income growth factor
777 "BoroCnstArt": 0.0, # Artificial borrowing constraint
778 # PARAMETERS REQUIRED TO SIMULATE THE MODEL
779 "AgentCount": 10000, # Number of agents of this type
780 "T_age": None, # Age after which simulated agents are automatically killed
781 "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor
782 # (The portion of PermGroFac attributable to aggregate productivity growth)
783 "NewbornTransShk": False, # Whether Newborns have transitory shock
784 # ADDITIONAL OPTIONAL PARAMETERS
785 "PerfMITShk": False, # Do Perfect Foresight MIT Shock
786 # (Forces Newborns to follow solution path of the agent they replaced if True)
787 "neutral_measure": False, # Whether to use permanent income neutral measure (see Harmenberg 2021)
788}
789init_agg_shocks.update(default_kNrmInitDstn_params)
790init_agg_shocks.update(default_pLvlInitDstn_params)
791init_agg_shocks.update(default_IncShkDstn_params)
792init_agg_shocks.update(default_aXtraGrid_params)
793init_agg_shocks.update(default_MgridBase_params)
796class AggShockConsumerType(IndShockConsumerType):
797 """
798 A class to represent consumers who face idiosyncratic (transitory and per-
799 manent) shocks to their income and live in an economy that has aggregate
800 (transitory and permanent) shocks to labor productivity. As the capital-
801 to-labor ratio varies in the economy, so does the wage rate and interest
802 rate. "Aggregate shock consumers" have beliefs about how the capital ratio
803 evolves over time and take aggregate shocks into account when making their
804 decision about how much to consume.
805 """
807 default_ = {
808 "params": init_agg_shocks,
809 "solver": solveConsAggShock,
810 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
811 }
812 time_inv_ = IndShockConsumerType.time_inv_.copy()
813 try:
814 time_inv_.remove("vFuncBool")
815 time_inv_.remove("CubicBool")
816 except: # pragma: nocover
817 pass
819 def reset(self):
820 """
821 Initialize this type for a new simulated history of K/L ratio.
823 Parameters
824 ----------
825 None
827 Returns
828 -------
829 None
830 """
831 self.initialize_sim()
832 self.state_now["aLvlNow"] = self.kInit * np.ones(
833 self.AgentCount
834 ) # Start simulation near SS
835 self.state_now["aNrm"] = self.state_now["aLvlNow"] / self.state_now["pLvl"]
837 def pre_solve(self):
838 self.construct("solution_terminal")
840 def get_economy_data(self, economy):
841 """
842 Imports economy-determined objects into self from a Market.
843 Instances of AggShockConsumerType "live" in some macroeconomy that has
844 attributes relevant to their microeconomic model, like the relationship
845 between the capital-to-labor ratio and the interest and wage rates; this
846 method imports those attributes from an "economy" object and makes them
847 attributes of the ConsumerType.
849 Parameters
850 ----------
851 economy : Market
852 The "macroeconomy" in which this instance "lives". Might be of the
853 subclass CobbDouglasEconomy, which has methods to generate the
854 relevant attributes.
856 Returns
857 -------
858 None
859 """
860 self.T_sim = (
861 economy.act_T
862 ) # Need to be able to track as many periods as economy runs
863 self.kInit = economy.kSS # Initialize simulation assets to steady state
864 self.aNrmInitMean = np.log(
865 0.00000001
866 ) # Initialize newborn assets to nearly zero
867 self.Mgrid = (
868 economy.MSS * self.MgridBase
869 ) # Aggregate market resources grid adjusted around SS capital ratio
870 self.AFunc = economy.AFunc # Next period's aggregate savings function
871 self.Rfunc = economy.Rfunc # Interest factor as function of capital ratio
872 self.wFunc = economy.wFunc # Wage rate as function of capital ratio
873 self.DeprRte = economy.DeprRte # Rate of capital depreciation
874 self.PermGroFacAgg = (
875 economy.PermGroFacAgg
876 ) # Aggregate permanent productivity growth
877 self.add_AggShkDstn(
878 economy.AggShkDstn
879 ) # Combine idiosyncratic and aggregate shocks into one dstn
880 self.add_to_time_inv(
881 "Mgrid", "AFunc", "Rfunc", "wFunc", "DeprRte", "PermGroFacAgg"
882 )
884 def add_AggShkDstn(self, AggShkDstn):
885 """
886 Updates attribute IncShkDstn by combining idiosyncratic shocks with aggregate shocks.
888 Parameters
889 ----------
890 AggShkDstn : [np.array]
891 Aggregate productivity shock distribution. First element is proba-
892 bilities, second element is agg permanent shocks, third element is
893 agg transitory shocks.
895 Returns
896 -------
897 None
898 """
899 if len(self.IncShkDstn[0].atoms) > 2:
900 self.IncShkDstn = self.IncShkDstnWithoutAggShocks
901 else:
902 self.IncShkDstnWithoutAggShocks = self.IncShkDstn
903 self.IncShkDstn = [
904 combine_indep_dstns(
905 self.IncShkDstn[t], AggShkDstn, seed=self.RNG.integers(0, 2**31 - 1)
906 )
907 for t in range(self.T_cycle)
908 ]
910 def sim_birth(self, which_agents):
911 """
912 Makes new consumers for the given indices. Initialized variables include
913 aNrm and pLvl, as well as time variables t_age and t_cycle.
915 Parameters
916 ----------
917 which_agents : np.array(Bool)
918 Boolean array of size self.AgentCount indicating which agents should be "born".
920 Returns
921 -------
922 None
923 """
924 IndShockConsumerType.sim_birth(self, which_agents)
925 if "aLvl" in self.state_now and self.state_now["aLvl"] is not None:
926 self.state_now["aLvl"][which_agents] = (
927 self.state_now["aNrm"][which_agents]
928 * self.state_now["pLvl"][which_agents]
929 )
930 else:
931 self.state_now["aLvl"] = self.state_now["aNrm"] * self.state_now["pLvl"]
933 def sim_death(self):
934 """
935 Randomly determine which consumers die, and distribute their wealth among the survivors.
936 This method only works if there is only one period in the cycle.
938 Parameters
939 ----------
940 None
942 Returns
943 -------
944 who_dies : np.array(bool)
945 Boolean array of size AgentCount indicating which agents die.
946 """
947 # Just select a random set of agents to die
948 how_many_die = int(round(self.AgentCount * (1.0 - self.LivPrb[0])))
949 base_bool = np.zeros(self.AgentCount, dtype=bool)
950 base_bool[0:how_many_die] = True
951 who_dies = self.RNG.permutation(base_bool)
952 if self.T_age is not None:
953 who_dies[self.t_age >= self.T_age] = True
955 # Divide up the wealth of those who die, giving it to those who survive
956 who_lives = np.logical_not(who_dies)
957 wealth_living = np.sum(self.state_now["aLvl"][who_lives])
958 wealth_dead = np.sum(self.state_now["aLvl"][who_dies])
959 Ractuarial = 1.0 + wealth_dead / wealth_living
960 self.state_now["aNrm"][who_lives] = (
961 self.state_now["aNrm"][who_lives] * Ractuarial
962 )
963 self.state_now["aLvl"][who_lives] = (
964 self.state_now["aLvl"][who_lives] * Ractuarial
965 )
966 return who_dies
968 def get_Rport(self):
969 """
970 Returns an array of size self.AgentCount with self.RfreeNow in every entry.
971 This is the risk-free portfolio return in this model.
973 Parameters
974 ----------
975 None
977 Returns
978 -------
979 RfreeNow : np.array
980 Array of size self.AgentCount with risk free interest rate for each agent.
981 """
982 RfreeNow = self.RfreeNow * np.ones(self.AgentCount)
983 return RfreeNow
985 def get_shocks(self):
986 """
987 Finds the effective permanent and transitory shocks this period by combining the aggregate
988 and idiosyncratic shocks of each type.
990 Parameters
991 ----------
992 None
994 Returns
995 -------
996 None
997 """
998 IndShockConsumerType.get_shocks(self) # Update idiosyncratic shocks
999 self.shocks["TranShk"] *= self.TranShkAggNow * self.wRteNow
1000 self.shocks["PermShk"] *= self.PermShkAggNow
1002 def get_controls(self):
1003 """
1004 Calculates consumption for each consumer of this type using the consumption functions.
1006 Parameters
1007 ----------
1008 None
1010 Returns
1011 -------
1012 None
1013 """
1014 cNrmNow = np.zeros(self.AgentCount) + np.nan
1015 MPCnow = np.zeros(self.AgentCount) + np.nan
1016 MaggNow = self.get_MaggNow()
1017 for t in range(self.T_cycle):
1018 these = t == self.t_cycle
1019 cNrmNow[these] = self.solution[t].cFunc(
1020 self.state_now["mNrm"][these], MaggNow[these]
1021 )
1022 MPCnow[these] = self.solution[t].cFunc.derivativeX(
1023 self.state_now["mNrm"][these], MaggNow[these]
1024 ) # Marginal propensity to consume
1026 self.controls["cNrm"] = cNrmNow
1027 self.MPCnow = MPCnow
1029 def get_MaggNow(self): # This function exists to be overwritten in StickyE model
1030 return self.MaggNow * np.ones(self.AgentCount)
1032 def market_action(self):
1033 """
1034 In the aggregate shocks model, the "market action" is to simulate one
1035 period of receiving income and choosing how much to consume.
1037 Parameters
1038 ----------
1039 None
1041 Returns
1042 -------
1043 None
1044 """
1045 self.simulate(1)
1047 def calc_bounding_values(self): # pragma: nocover
1048 """
1049 NOT YET IMPLEMENTED FOR THIS CLASS
1050 """
1051 raise NotImplementedError()
1053 def make_euler_error_func(self, mMax=100, approx_inc_dstn=True): # pragma: nocover
1054 """
1055 NOT YET IMPLEMENTED FOR THIS CLASS
1056 """
1057 raise NotImplementedError()
1059 def check_conditions(self, verbose=None): # pragma: nocover
1060 raise NotImplementedError()
1062 def calc_limiting_values(self): # pragma: nocover
1063 raise NotImplementedError()
1066###############################################################################
1069# This example makes a high risk, low growth state and a low risk, high growth state
1070MrkvArray = np.array([[0.90, 0.10], [0.04, 0.96]])
1072# Make a dictionary to specify a Markov aggregate shocks consumer
1073init_agg_mrkv_shocks = init_agg_shocks.copy()
1074init_agg_mrkv_shocks["MrkvArray"] = MrkvArray
1075aggmrkv_constructor_dict = aggshock_constructor_dict.copy()
1076aggmrkv_constructor_dict["solution_terminal"] = make_aggmrkv_solution_terminal
1077init_agg_mrkv_shocks["constructors"] = aggmrkv_constructor_dict
1080class AggShockMarkovConsumerType(AggShockConsumerType):
1081 """
1082 A class for representing ex ante heterogeneous "types" of consumers who
1083 experience both aggregate and idiosyncratic shocks to productivity (both
1084 permanent and transitory), who lives in an environment where the macroeconomic
1085 state is subject to Markov-style discrete state evolution.
1086 """
1088 time_inv_ = AggShockConsumerType.time_inv_ + ["MrkvArray"]
1089 shock_vars_ = AggShockConsumerType.shock_vars_ + ["Mrkv"]
1090 default_ = {
1091 "params": init_agg_mrkv_shocks,
1092 "solver": solve_ConsAggMarkov,
1093 "track_vars": ["aNrm", "cNrm", "mNrm", "pLvl"],
1094 }
1096 def add_AggShkDstn(self, AggShkDstn):
1097 """
1098 Variation on AggShockConsumerType.add_AggShkDstn that handles the Markov
1099 state. AggShkDstn is a list of aggregate productivity shock distributions
1100 for each Markov state.
1101 """
1102 if len(self.IncShkDstn[0][0].atoms) > 2:
1103 self.IncShkDstn = self.IncShkDstnWithoutAggShocks
1104 else:
1105 self.IncShkDstnWithoutAggShocks = self.IncShkDstn
1107 IncShkDstnOut = []
1108 N = self.MrkvArray.shape[0]
1109 for t in range(self.T_cycle):
1110 IncShkDstnOut.append(
1111 [
1112 combine_indep_dstns(
1113 self.IncShkDstn[t][n],
1114 AggShkDstn[n],
1115 seed=self.RNG.integers(0, 2**31 - 1),
1116 )
1117 for n in range(N)
1118 ]
1119 )
1120 self.IncShkDstn = IncShkDstnOut
1122 def initialize_sim(self):
1123 self.shocks["Mrkv"] = 0
1124 AggShockConsumerType.initialize_sim(self)
1126 def get_shocks(self):
1127 """
1128 Gets permanent and transitory income shocks for this period. Samples from IncShkDstn for
1129 each period in the cycle. This is a copy-paste from IndShockConsumerType, with the
1130 addition of the Markov macroeconomic state. Unfortunately, the get_shocks method for
1131 MarkovConsumerType cannot be used, as that method assumes that MrkvNow is a vector
1132 with a value for each agent, not just a single int.
1134 Parameters
1135 ----------
1136 None
1138 Returns
1139 -------
1140 None
1141 """
1142 PermShkNow = np.zeros(self.AgentCount) # Initialize shock arrays
1143 TranShkNow = np.zeros(self.AgentCount)
1144 newborn = self.t_age == 0
1145 for t in range(self.T_cycle):
1146 these = t == self.t_cycle
1147 N = np.sum(these)
1148 if N > 0:
1149 IncShkDstnNow = self.IncShkDstn[t - 1][
1150 self.shocks["Mrkv"]
1151 ] # set current income distribution
1152 # and permanent growth factor
1153 PermGroFacNow = self.PermGroFac[t - 1]
1155 # Get random draws of income shocks from the discrete distribution
1156 ShockDraws = IncShkDstnNow.draw(N, shuffle=True)
1157 # Permanent "shock" includes expected growth
1158 PermShkNow[these] = ShockDraws[0] * PermGroFacNow
1159 TranShkNow[these] = ShockDraws[1]
1161 # That procedure used the *last* period in the sequence for newborns, but that's not right
1162 # Redraw shocks for newborns, using the *first* period in the sequence. Approximation.
1163 N = np.sum(newborn)
1164 if N > 0:
1165 these = newborn
1166 IncShkDstnNow = self.IncShkDstn[0][
1167 self.shocks["Mrkv"]
1168 ] # set current income distribution
1169 PermGroFacNow = self.PermGroFac[0] # and permanent growth factor
1171 # Get random draws of income shocks from the discrete distribution
1172 ShockDraws = IncShkDstnNow.draw(N, shuffle=True)
1174 # Permanent "shock" includes expected growth
1175 PermShkNow[these] = ShockDraws[0] * PermGroFacNow
1176 TranShkNow[these] = ShockDraws[1]
1178 # Store the shocks in self
1179 self.EmpNow = np.ones(self.AgentCount, dtype=bool)
1180 self.EmpNow[TranShkNow == self.IncUnemp] = False
1181 self.shocks["TranShk"] = TranShkNow * self.TranShkAggNow * self.wRteNow
1182 self.shocks["PermShk"] = PermShkNow * self.PermShkAggNow
1184 def get_controls(self):
1185 """
1186 Calculates consumption for each consumer of this type using the consumption functions.
1187 For this AgentType class, MrkvNow is the same for all consumers. However, in an
1188 extension with "macroeconomic inattention", consumers might misperceive the state
1189 and thus act as if they are in different states.
1191 Parameters
1192 ----------
1193 None
1195 Returns
1196 -------
1197 None
1198 """
1199 cNrmNow = np.zeros(self.AgentCount) + np.nan
1200 MPCnow = np.zeros(self.AgentCount) + np.nan
1201 MaggNow = self.get_MaggNow()
1202 MrkvNow = self.getMrkvNow()
1204 StateCount = self.MrkvArray.shape[0]
1205 MrkvBoolArray = np.zeros((StateCount, self.AgentCount), dtype=bool)
1206 for i in range(StateCount):
1207 MrkvBoolArray[i, :] = i == MrkvNow
1209 for t in range(self.T_cycle):
1210 these = t == self.t_cycle
1211 for i in range(StateCount):
1212 those = np.logical_and(these, MrkvBoolArray[i, :])
1213 cNrmNow[those] = self.solution[t].cFunc[i](
1214 self.state_now["mNrm"][those], MaggNow[those]
1215 )
1216 # Marginal propensity to consume
1217 MPCnow[those] = (
1218 self.solution[t]
1219 .cFunc[i]
1220 .derivativeX(self.state_now["mNrm"][those], MaggNow[those])
1221 )
1222 self.controls["cNrm"] = cNrmNow
1223 self.MPCnow = MPCnow
1224 return None
1226 def getMrkvNow(self): # This function exists to be overwritten in StickyE model
1227 return self.shocks["Mrkv"] * np.ones(self.AgentCount, dtype=int)
1230###############################################################################
1232# Define some constructor functions for the basic Krusell-Smith model
1235def make_solution_terminal_KS(CRRA):
1236 cFunc_terminal = 4 * [IdentityFunction(n_dims=2)]
1237 vPfunc_terminal = [MargValueFuncCRRA(cFunc_terminal[j], CRRA) for j in range(4)]
1238 solution_terminal = ConsumerSolution(cFunc=cFunc_terminal, vPfunc=vPfunc_terminal)
1239 return solution_terminal
1242def make_assets_grid_KS(aMin, aMax, aCount, aNestFac):
1243 return make_assets_grid(aMin, aMax, aCount, None, aNestFac)
1246def make_KS_transition_arrays(
1247 aGrid, Mgrid, AFunc, LbrInd, UrateB, UrateG, ProdB, ProdG, MrkvIndArray
1248):
1249 """
1250 Construct the attributes ProbArray, mNextArray, MnextArray, and RnextArray,
1251 which will be used by the one period solver. The information for this method
1252 is usually obtained by the get_economy_data method. Output is returned as a
1253 *list* of four arrays, which are later assigned to their appropriate attributes.
1255 Parameters
1256 ----------
1257 aGrid : np.array
1258 Grid of end-of-period individual assets.
1259 MGrid : np.array
1260 Grid of aggregate market resources.
1261 AFunc : function
1262 End-of-period aggregate assets as a function of aggregate market resources.
1263 LbrInd : float
1264 Individual labor supply measure.
1265 UrateB : float
1266 Unemployment rate in the "bad" aggregate state.
1267 UrateG : float
1268 Unemployment rate in the "good" aggregate state.
1269 ProdB : float
1270 TFP in the "bad" aggregate state.
1271 ProdG : float
1272 TFP in the "good" aggregate state.
1273 MrkvIndArray : np.array
1274 Markov transition probabilities from the perspective of the individual.
1276 Returns
1277 -------
1278 ProbArray : np.array
1279 Array of discrete future outcome probabilities.
1280 mNextArray : np.array
1281 Array of discrete realizations of next-period idiosyncratic market resources.
1282 MnextArray : np.array
1283 Array of discrete realizations of next-period aggregate market resources.
1284 RnextArray : np.array
1285 Array of discrete realizations of next-period rate of return.
1286 """
1287 # Get array sizes
1288 aCount = aGrid.size
1289 Mcount = Mgrid.size
1291 # Make tiled array of end-of-period idiosyncratic assets (order: a, M, s, s')
1292 aNow_tiled = np.tile(np.reshape(aGrid, [aCount, 1, 1, 1]), [1, Mcount, 4, 4])
1294 # Make arrays of end-of-period aggregate assets (capital next period)
1295 AnowB = AFunc[0](Mgrid)
1296 AnowG = AFunc[1](Mgrid)
1297 KnextB = np.tile(np.reshape(AnowB, [1, Mcount, 1, 1]), [1, 1, 1, 4])
1298 KnextG = np.tile(np.reshape(AnowG, [1, Mcount, 1, 1]), [1, 1, 1, 4])
1299 Knext = np.concatenate((KnextB, KnextB, KnextG, KnextG), axis=2)
1301 # Make arrays of aggregate labor and TFP next period
1302 Lnext = np.zeros((1, Mcount, 4, 4)) # shape (1,Mcount,4,4)
1303 Lnext[0, :, :, 0:2] = (1.0 - UrateB) * LbrInd
1304 Lnext[0, :, :, 2:4] = (1.0 - UrateG) * LbrInd
1305 Znext = np.zeros((1, Mcount, 4, 4))
1306 Znext[0, :, :, 0:2] = ProdB
1307 Znext[0, :, :, 2:4] = ProdG
1309 # Calculate (net) interest factor and wage rate next period
1310 KtoLnext = Knext / Lnext
1311 Rnext = 1.0 + Znext * CapShare * KtoLnext ** (CapShare - 1.0) - DeprRte
1312 Wnext = Znext * (1.0 - CapShare) * KtoLnext**CapShare
1314 # Calculate aggregate market resources next period
1315 Ynext = Znext * Knext**CapShare * Lnext ** (1.0 - CapShare)
1316 Mnext = (1.0 - DeprRte) * Knext + Ynext
1318 # Tile the interest, wage, and aggregate market resources arrays
1319 Rnext_tiled = np.tile(Rnext, [aCount, 1, 1, 1])
1320 Wnext_tiled = np.tile(Wnext, [aCount, 1, 1, 1])
1321 Mnext_tiled = np.tile(Mnext, [aCount, 1, 1, 1])
1323 # Make an array of idiosyncratic labor supply next period
1324 lNext_tiled = np.zeros([aCount, Mcount, 4, 4])
1325 lNext_tiled[:, :, :, 1] = LbrInd
1326 lNext_tiled[:, :, :, 3] = LbrInd
1328 # Calculate idiosyncratic market resources next period
1329 mNext = Rnext_tiled * aNow_tiled + Wnext_tiled * lNext_tiled
1331 # Make a tiled array of transition probabilities
1332 Probs_tiled = np.tile(
1333 np.reshape(MrkvIndArray, [1, 1, 4, 4]), [aCount, Mcount, 1, 1]
1334 )
1336 # Return the attributes that will be used by the solver
1337 transition_arrays = {
1338 "ProbArray": Probs_tiled,
1339 "mNextArray": mNext,
1340 "MnextArray": Mnext_tiled,
1341 "RnextArray": Rnext_tiled,
1342 }
1343 return transition_arrays
1346###############################################################################
1348# Make a dictionary for Krusell-Smith agents
1349KS_constructor_dict = {
1350 "solution_terminal": make_solution_terminal_KS,
1351 "aGrid": make_assets_grid_KS,
1352 "transition_arrays": make_KS_transition_arrays,
1353 "ProbArray": get_it_from("transition_arrays"),
1354 "mNextArray": get_it_from("transition_arrays"),
1355 "MnextArray": get_it_from("transition_arrays"),
1356 "RnextArray": get_it_from("transition_arrays"),
1357 "MgridBase": make_exponential_MgridBase,
1358}
1360init_KS_agents = {
1361 "T_cycle": 1,
1362 "pseudo_terminal": False,
1363 "constructors": KS_constructor_dict,
1364 "DiscFac": 0.99,
1365 "CRRA": 1.0,
1366 "LbrInd": 1.0,
1367 "aMin": 0.001,
1368 "aMax": 50.0,
1369 "aCount": 32,
1370 "aNestFac": 2,
1371 "MaggCount": 25,
1372 "MaggPerturb": 0.01,
1373 "MaggExpFac": 0.12,
1374 "MgridBase": np.array([0.99, 1.0, 1.01]), ## dummy, this will be overwritten
1375 "AgentCount": 5000,
1376}
1379class KrusellSmithType(AgentType):
1380 """
1381 A class for representing agents in the seminal Krusell-Smith (1998) model from
1382 the paper "Income and Wealth Heterogeneity in the Macroeconomy". All default
1383 parameters have been set to match those in the paper, but the equilibrium object
1384 is perceptions of aggregate assets as a function of aggregate market resources
1385 in each macroeconomic state (bad=0, good=1), rather than aggregate capital as
1386 a function of previous aggregate capital. This choice was made so that some
1387 of the code from HARK's other HA-macro models can be used.
1389 To make this class work properly, instantiate both this class and an instance
1390 of KrusellSmithEconomy, then use this class' get_economy_data method with the
1391 economy object.
1392 """
1394 time_inv_ = [
1395 "DiscFac",
1396 "CRRA",
1397 "aGrid",
1398 "ProbArray",
1399 "mNextArray",
1400 "MnextArray",
1401 "RnextArray",
1402 ]
1403 time_vary_ = []
1404 shock_vars_ = ["Mrkv"]
1405 state_vars = ["aNow", "mNow", "EmpNow"]
1406 default_ = {
1407 "params": init_KS_agents,
1408 "solver": solve_KrusellSmith,
1409 "track_vars": ["aNow", "cNow", "mNow", "EmpNow"],
1410 }
1412 def __init__(self, **kwds):
1413 temp = kwds.copy()
1414 temp["construct"] = False
1415 AgentType.__init__(self, **temp)
1416 self.construct("MgridBase")
1418 # Special case: this type *must* be initialized with construct=False
1419 # because the data required to make its solution attributes is obtained
1420 # from the associated economy, not passed as part of its parameters.
1421 # To make it work properly, instantiate both this class and an instance
1422 # of KrusellSmithEconomy, then use this class' get_economy_data method.
1423 # Exception: MgridBase must exist
1425 def pre_solve(self):
1426 self.construct("solution_terminal")
1428 def get_economy_data(self, Economy):
1429 """
1430 Imports economy-determined objects into self from a Market.
1432 Parameters
1433 ----------
1434 Economy : KrusellSmithEconomy
1435 The "macroeconomy" in which this instance "lives".
1437 Returns
1438 -------
1439 None
1440 """
1441 self.T_sim = (
1442 Economy.act_T
1443 ) # Need to be able to track as many periods as economy runs
1444 self.kInit = Economy.KSS # Initialize simulation assets to steady state
1445 self.MrkvInit = Economy.sow_init[
1446 "Mrkv"
1447 ] # Starting Markov state for the macroeconomy
1448 self.Mgrid = (
1449 Economy.MSS * self.MgridBase
1450 ) # Aggregate market resources grid adjusted around SS capital ratio
1451 self.AFunc = Economy.AFunc # Next period's aggregate savings function
1452 self.DeprRte = Economy.DeprRte # Rate of capital depreciation
1453 self.CapShare = Economy.CapShare # Capital's share of production
1454 # Idiosyncratic labor supply (when employed)
1455 self.LbrInd = Economy.LbrInd
1456 self.UrateB = Economy.UrateB # Unemployment rate in bad state
1457 self.UrateG = Economy.UrateG # Unemployment rate in good state
1458 self.ProdB = Economy.ProdB # Total factor productivity in bad state
1459 self.ProdG = Economy.ProdG # Total factor productivity in good state
1460 self.MrkvIndArray = (
1461 Economy.MrkvIndArray
1462 ) # Transition probabilities among discrete states
1463 self.MrkvAggArray = (
1464 Economy.MrkvArray
1465 ) # Transition probabilities among aggregate discrete states
1466 self.add_to_time_inv(
1467 "Mgrid",
1468 "AFunc",
1469 "DeprRte",
1470 "CapShare",
1471 "UrateB",
1472 "LbrInd",
1473 "UrateG",
1474 "ProdB",
1475 "ProdG",
1476 "MrkvIndArray",
1477 "MrkvAggArray",
1478 )
1480 def make_emp_idx_arrays(self):
1481 """
1482 Construct the attributes emp_permute and unemp_permute, each of which is
1483 a 2x2 nested list of boolean arrays. The j,k-th element of emp_permute
1484 represents the employment states this period for agents who were employed
1485 last period when the macroeconomy is transitioning from state j to state k.
1486 Likewise, j,k-th element of unemp_permute represents the employment states
1487 this period for agents who were unemployed last period when the macro-
1488 economy is transitioning from state j to state k. These attributes are
1489 referenced during simulation, when they are randomly permuted in order to
1490 maintain exact unemployment rates in each period.
1491 """
1492 # Get counts of employed and unemployed agents in each macroeconomic state
1493 B_unemp_N = int(np.round(self.UrateB * self.AgentCount))
1494 B_emp_N = self.AgentCount - B_unemp_N
1495 G_unemp_N = int(np.round(self.UrateG * self.AgentCount))
1496 G_emp_N = self.AgentCount - G_unemp_N
1498 # Bad-bad transition indices
1499 BB_stay_unemp_N = int(
1500 np.round(B_unemp_N * self.MrkvIndArray[0, 0] / self.MrkvAggArray[0, 0])
1501 )
1502 BB_become_unemp_N = B_unemp_N - BB_stay_unemp_N
1503 BB_stay_emp_N = int(
1504 np.round(B_emp_N * self.MrkvIndArray[1, 1] / self.MrkvAggArray[0, 0])
1505 )
1506 BB_become_emp_N = B_emp_N - BB_stay_emp_N
1507 BB_unemp_permute = np.concatenate(
1508 [
1509 np.ones(BB_become_emp_N, dtype=bool),
1510 np.zeros(BB_stay_unemp_N, dtype=bool),
1511 ]
1512 )
1513 BB_emp_permute = np.concatenate(
1514 [
1515 np.ones(BB_stay_emp_N, dtype=bool),
1516 np.zeros(BB_become_unemp_N, dtype=bool),
1517 ]
1518 )
1520 # Bad-good transition indices
1521 BG_stay_unemp_N = int(
1522 np.round(B_unemp_N * self.MrkvIndArray[0, 2] / self.MrkvAggArray[0, 1])
1523 )
1524 BG_become_unemp_N = G_unemp_N - BG_stay_unemp_N
1525 BG_stay_emp_N = int(
1526 np.round(B_emp_N * self.MrkvIndArray[1, 3] / self.MrkvAggArray[0, 1])
1527 )
1528 BG_become_emp_N = G_emp_N - BG_stay_emp_N
1529 BG_unemp_permute = np.concatenate(
1530 [
1531 np.ones(BG_become_emp_N, dtype=bool),
1532 np.zeros(BG_stay_unemp_N, dtype=bool),
1533 ]
1534 )
1535 BG_emp_permute = np.concatenate(
1536 [
1537 np.ones(BG_stay_emp_N, dtype=bool),
1538 np.zeros(BG_become_unemp_N, dtype=bool),
1539 ]
1540 )
1542 # Good-bad transition indices
1543 GB_stay_unemp_N = int(
1544 np.round(G_unemp_N * self.MrkvIndArray[2, 0] / self.MrkvAggArray[1, 0])
1545 )
1546 GB_become_unemp_N = B_unemp_N - GB_stay_unemp_N
1547 GB_stay_emp_N = int(
1548 np.round(G_emp_N * self.MrkvIndArray[3, 1] / self.MrkvAggArray[1, 0])
1549 )
1550 GB_become_emp_N = B_emp_N - GB_stay_emp_N
1551 GB_unemp_permute = np.concatenate(
1552 [
1553 np.ones(GB_become_emp_N, dtype=bool),
1554 np.zeros(GB_stay_unemp_N, dtype=bool),
1555 ]
1556 )
1557 GB_emp_permute = np.concatenate(
1558 [
1559 np.ones(GB_stay_emp_N, dtype=bool),
1560 np.zeros(GB_become_unemp_N, dtype=bool),
1561 ]
1562 )
1564 # Good-good transition indices
1565 GG_stay_unemp_N = int(
1566 np.round(G_unemp_N * self.MrkvIndArray[2, 2] / self.MrkvAggArray[1, 1])
1567 )
1568 GG_become_unemp_N = G_unemp_N - GG_stay_unemp_N
1569 GG_stay_emp_N = int(
1570 np.round(G_emp_N * self.MrkvIndArray[3, 3] / self.MrkvAggArray[1, 1])
1571 )
1572 GG_become_emp_N = G_emp_N - GG_stay_emp_N
1573 GG_unemp_permute = np.concatenate(
1574 [
1575 np.ones(GG_become_emp_N, dtype=bool),
1576 np.zeros(GG_stay_unemp_N, dtype=bool),
1577 ]
1578 )
1579 GG_emp_permute = np.concatenate(
1580 [
1581 np.ones(GG_stay_emp_N, dtype=bool),
1582 np.zeros(GG_become_unemp_N, dtype=bool),
1583 ]
1584 )
1586 # Store transition matrices as attributes of self
1587 self.unemp_permute = [
1588 [BB_unemp_permute, BG_unemp_permute],
1589 [GB_unemp_permute, GG_unemp_permute],
1590 ]
1591 self.emp_permute = [
1592 [BB_emp_permute, BG_emp_permute],
1593 [GB_emp_permute, GG_emp_permute],
1594 ]
1596 def reset(self):
1597 self.initialize_sim()
1599 def market_action(self):
1600 self.simulate(1)
1602 def initialize_sim(self):
1603 self.shocks["Mrkv"] = self.MrkvInit
1604 AgentType.initialize_sim(self)
1605 self.state_now["EmpNow"] = self.state_now["EmpNow"].astype(bool)
1606 self.make_emp_idx_arrays()
1608 def sim_birth(self, which):
1609 """
1610 Create newborn agents with randomly drawn employment states. This will
1611 only ever be called by initialize_sim() at the start of a new simulation
1612 history, as the Krusell-Smith model does not have death and replacement.
1613 The sim_death() method does not exist, as AgentType's default of "no death"
1614 is the correct behavior for the model.
1615 """
1616 N = np.sum(which)
1617 if N == 0:
1618 return
1620 if self.shocks["Mrkv"] == 0:
1621 unemp_N = int(np.round(self.UrateB * N))
1622 emp_N = self.AgentCount - unemp_N
1623 elif self.shocks["Mrkv"] == 1:
1624 unemp_N = int(np.round(self.UrateG * N))
1625 emp_N = self.AgentCount - unemp_N
1626 else:
1627 assert False, "Illegal macroeconomic state: MrkvNow must be 0 or 1"
1628 EmpNew = np.concatenate(
1629 [np.zeros(unemp_N, dtype=bool), np.ones(emp_N, dtype=bool)]
1630 )
1632 self.state_now["EmpNow"][which] = self.RNG.permutation(EmpNew)
1633 self.state_now["aNow"][which] = self.kInit
1635 def get_shocks(self):
1636 """
1637 Get new idiosyncratic employment states based on the macroeconomic state.
1638 """
1639 # Get boolean arrays for current employment states
1640 employed = self.state_prev["EmpNow"].copy().astype(bool)
1641 unemployed = np.logical_not(employed)
1643 # derive from past employment rate rather than store previous value
1644 mrkv_prev = int((unemployed.sum() / float(self.AgentCount)) != self.UrateB)
1646 # Transition some agents between unemployment and employment
1647 emp_permute = self.emp_permute[mrkv_prev][self.shocks["Mrkv"]]
1648 unemp_permute = self.unemp_permute[mrkv_prev][self.shocks["Mrkv"]]
1649 # TODO: replace poststate_vars functionality with shocks here
1650 EmpNow = self.state_now["EmpNow"]
1652 # It's really this permutation that is the shock...
1653 # This apparatus is trying to 'exact match' the 'internal' Markov process.
1654 EmpNow[employed] = self.RNG.permutation(emp_permute)
1655 EmpNow[unemployed] = self.RNG.permutation(unemp_permute)
1657 def get_states(self):
1658 """
1659 Get each agent's idiosyncratic state, their household market resources.
1660 """
1661 self.state_now["mNow"] = (
1662 self.Rnow * self.state_prev["aNow"]
1663 + self.Wnow * self.LbrInd * self.state_now["EmpNow"]
1664 )
1666 def get_controls(self):
1667 """
1668 Get each agent's consumption given their current state.'
1669 """
1670 employed = self.state_now["EmpNow"].copy().astype(bool)
1671 unemployed = np.logical_not(employed)
1673 # Get the discrete index for (un)employed agents
1674 if self.shocks["Mrkv"] == 0: # Bad macroeconomic conditions
1675 unemp_idx = 0
1676 emp_idx = 1
1677 elif self.shocks["Mrkv"] == 1: # Good macroeconomic conditions
1678 unemp_idx = 2
1679 emp_idx = 3
1680 else:
1681 assert False, "Illegal macroeconomic state: MrkvNow must be 0 or 1"
1683 # Get consumption for each agent using the appropriate consumption function
1684 cNow = np.zeros(self.AgentCount)
1685 Mnow = self.Mnow * np.ones(self.AgentCount)
1686 cNow[unemployed] = self.solution[0].cFunc[unemp_idx](
1687 self.state_now["mNow"][unemployed], Mnow[unemployed]
1688 )
1689 cNow[employed] = self.solution[0].cFunc[emp_idx](
1690 self.state_now["mNow"][employed], Mnow[employed]
1691 )
1692 self.controls["cNow"] = cNow
1694 def get_poststates(self):
1695 """
1696 Gets each agent's retained assets after consumption.
1697 """
1698 self.state_now["aNow"] = self.state_now["mNow"] - self.controls["cNow"]
1701CRRA = 2.0
1702DiscFac = 0.96
1704# Parameters for a Cobb-Douglas economy
1705PermGroFacAgg = 1.00 # Aggregate permanent income growth factor
1706PermShkAggCount = (
1707 3 # Number of points in discrete approximation to aggregate permanent shock dist
1708)
1709TranShkAggCount = (
1710 3 # Number of points in discrete approximation to aggregate transitory shock dist
1711)
1712PermShkAggStd = 0.0063 # Standard deviation of log aggregate permanent shocks
1713TranShkAggStd = 0.0031 # Standard deviation of log aggregate transitory shocks
1714DeprRte = 0.025 # Capital depreciation rate
1715CapShare = 0.36 # Capital's share of income
1716DiscFacPF = DiscFac # Discount factor of perfect foresight calibration
1717CRRAPF = CRRA # Coefficient of relative risk aversion of perfect foresight calibration
1718intercept_prev = 0.0 # Intercept of aggregate savings function
1719slope_prev = 1.0 # Slope of aggregate savings function
1720verbose_cobb_douglas = (
1721 True # Whether to print solution progress to screen while solving
1722)
1723T_discard = 200 # Number of simulated "burn in" periods to discard when updating AFunc
1724# Damping factor when updating AFunc; puts DampingFac weight on old params, rest on new
1725DampingFac = 0.5
1726max_loops = 20 # Maximum number of AFunc updating loops to allow
1729# Make a dictionary to specify a Cobb-Douglas economy
1730init_cobb_douglas = {
1731 "PermShkAggCount": PermShkAggCount,
1732 "TranShkAggCount": TranShkAggCount,
1733 "PermShkAggStd": PermShkAggStd,
1734 "TranShkAggStd": TranShkAggStd,
1735 "DeprRte": DeprRte,
1736 "CapShare": CapShare,
1737 "DiscFac": DiscFacPF,
1738 "CRRA": CRRAPF,
1739 "PermGroFacAgg": PermGroFacAgg,
1740 "AggregateL": 1.0,
1741 "intercept_prev": intercept_prev,
1742 "slope_prev": slope_prev,
1743 "verbose": verbose_cobb_douglas,
1744 "T_discard": T_discard,
1745 "DampingFac": DampingFac,
1746 "max_loops": max_loops,
1747}
1750class CobbDouglasEconomy(Market):
1751 """
1752 A class to represent an economy with a Cobb-Douglas aggregate production
1753 function over labor and capital, extending HARK.Market. The "aggregate
1754 market process" for this market combines all individuals' asset holdings
1755 into aggregate capital, yielding the interest factor on assets and the wage
1756 rate for the upcoming period.
1758 Note: The current implementation assumes a constant labor supply, but
1759 this will be generalized in the future.
1761 Parameters
1762 ----------
1763 agents : [ConsumerType]
1764 List of types of consumers that live in this economy.
1765 tolerance: float
1766 Minimum acceptable distance between "dynamic rules" to consider the
1767 solution process converged. Distance depends on intercept and slope
1768 of the log-linear "next capital ratio" function.
1769 act_T : int
1770 Number of periods to simulate when making a history of of the market.
1771 """
1773 def __init__(self, agents=None, tolerance=0.0001, act_T=1200, **kwds):
1774 agents = agents if agents is not None else list()
1775 params = init_cobb_douglas.copy()
1776 params["sow_vars"] = [
1777 "MaggNow",
1778 "AaggNow",
1779 "RfreeNow",
1780 "wRteNow",
1781 "PermShkAggNow",
1782 "TranShkAggNow",
1783 "KtoLnow",
1784 ]
1785 params["reap_vars"] = ["aLvl", "pLvl"]
1786 params["track_vars"] = ["MaggNow", "AaggNow"]
1787 params["dyn_vars"] = ["AFunc"]
1788 params["distributions"] = ["PermShkAggDstn", "TranShkAggDstn", "AggShkDstn"]
1789 params.update(kwds)
1791 Market.__init__(self, agents=agents, tolerance=tolerance, act_T=act_T, **params)
1792 self.update()
1794 def mill_rule(self, aLvl, pLvl):
1795 """
1796 Function to calculate the capital to labor ratio, interest factor, and
1797 wage rate based on each agent's current state. Just calls calc_R_and_W().
1799 See documentation for calc_R_and_W for more information.
1800 """
1801 return self.calc_R_and_W(aLvl, pLvl)
1803 def calc_dynamics(self, MaggNow, AaggNow):
1804 """
1805 Calculates a new dynamic rule for the economy: end of period savings as
1806 a function of aggregate market resources. Just calls calc_AFunc().
1808 See documentation for calc_AFunc for more information.
1809 """
1810 return self.calc_AFunc(MaggNow, AaggNow)
1812 def update(self):
1813 """
1814 Use primitive parameters (and perfect foresight calibrations) to make
1815 interest factor and wage rate functions (of capital to labor ratio),
1816 as well as discrete approximations to the aggregate shock distributions.
1818 Parameters
1819 ----------
1820 None
1822 Returns
1823 -------
1824 None
1825 """
1826 self.kSS = (
1827 (
1828 self.get_PermGroFacAggLR() ** (self.CRRA) / self.DiscFac
1829 - (1.0 - self.DeprRte)
1830 )
1831 / self.CapShare
1832 ) ** (1.0 / (self.CapShare - 1.0))
1833 self.KtoYSS = self.kSS ** (1.0 - self.CapShare)
1834 self.wRteSS = (1.0 - self.CapShare) * self.kSS ** (self.CapShare)
1835 self.RfreeSS = (
1836 1.0 + self.CapShare * self.kSS ** (self.CapShare - 1.0) - self.DeprRte
1837 )
1838 self.MSS = self.kSS * self.RfreeSS + self.wRteSS
1839 self.convertKtoY = lambda KtoY: KtoY ** (
1840 1.0 / (1.0 - self.CapShare)
1841 ) # converts K/Y to K/L
1842 self.Rfunc = lambda k: (
1843 1.0 + self.CapShare * k ** (self.CapShare - 1.0) - self.DeprRte
1844 )
1845 self.wFunc = lambda k: ((1.0 - self.CapShare) * k ** (self.CapShare))
1847 self.sow_init["KtoLnow"] = self.kSS
1848 self.sow_init["MaggNow"] = self.kSS
1849 self.sow_init["AaggNow"] = self.kSS
1850 self.sow_init["RfreeNow"] = self.Rfunc(self.kSS)
1851 self.sow_init["wRteNow"] = self.wFunc(self.kSS)
1852 self.sow_init["PermShkAggNow"] = 1.0
1853 self.sow_init["TranShkAggNow"] = 1.0
1854 self.make_AggShkDstn()
1855 self.AFunc = AggregateSavingRule(self.intercept_prev, self.slope_prev)
1857 def get_PermGroFacAggLR(self):
1858 """
1859 A trivial function that returns self.PermGroFacAgg. Exists to be overwritten
1860 and extended by ConsAggShockMarkov model.
1862 Parameters
1863 ----------
1864 None
1866 Returns
1867 -------
1868 PermGroFacAggLR : float
1869 Long run aggregate permanent income growth, which is the same thing
1870 as aggregate permanent income growth.
1871 """
1872 return self.PermGroFacAgg
1874 def make_AggShkDstn(self):
1875 """
1876 Creates the attributes TranShkAggDstn, PermShkAggDstn, and AggShkDstn.
1877 Draws on attributes TranShkAggStd, PermShkAddStd, TranShkAggCount, PermShkAggCount.
1879 Parameters
1880 ----------
1881 None
1883 Returns
1884 -------
1885 None
1886 """
1887 RNG = self.RNG
1888 TranShkAggDstn = MeanOneLogNormal(
1889 sigma=self.TranShkAggStd, seed=RNG.integers(0, 2**31 - 1)
1890 )
1891 self.TranShkAggDstn = TranShkAggDstn.discretize(N=self.TranShkAggCount)
1892 PermShkAggDstn = MeanOneLogNormal(
1893 sigma=self.PermShkAggStd, seed=RNG.integers(0, 2**31 - 1)
1894 )
1895 self.PermShkAggDstn = PermShkAggDstn.discretize(N=self.PermShkAggCount)
1896 self.AggShkDstn = combine_indep_dstns(
1897 self.PermShkAggDstn, self.TranShkAggDstn, seed=RNG.integers(0, 2**31 - 1)
1898 )
1900 def reset(self):
1901 """
1902 Reset the economy to prepare for a new simulation. Sets the time index
1903 of aggregate shocks to zero and runs Market.reset().
1905 Parameters
1906 ----------
1907 None
1909 Returns
1910 -------
1911 None
1912 """
1913 self.Shk_idx = 0
1914 Market.reset(self)
1916 def make_AggShkHist(self):
1917 """
1918 Make simulated histories of aggregate transitory and permanent shocks.
1919 Histories are of length self.act_T, for use in the general equilibrium
1920 simulation.
1922 Parameters
1923 ----------
1924 None
1926 Returns
1927 -------
1928 None
1929 """
1930 sim_periods = self.act_T
1931 Events = np.arange(self.AggShkDstn.pmv.size) # just a list of integers
1932 EventDraws = self.AggShkDstn.draw(N=sim_periods, atoms=Events)
1933 PermShkAggHist = self.AggShkDstn.atoms[0][EventDraws]
1934 TranShkAggHist = self.AggShkDstn.atoms[1][EventDraws]
1936 # Store the histories
1937 self.PermShkAggHist = PermShkAggHist * self.PermGroFacAgg
1938 self.TranShkAggHist = TranShkAggHist
1940 def calc_R_and_W(self, aLvlNow, pLvlNow):
1941 """
1942 Calculates the interest factor and wage rate this period using each agent's
1943 capital stock to get the aggregate capital ratio.
1945 Parameters
1946 ----------
1947 aLvlNow : [np.array]
1948 Agents' current end-of-period assets. Elements of the list correspond
1949 to types in the economy, entries within arrays to agents of that type.
1951 Returns
1952 -------
1953 MaggNow : float
1954 Aggregate market resources for this period normalized by mean permanent income
1955 AaggNow : float
1956 Aggregate savings for this period normalized by mean permanent income
1957 RfreeNow : float
1958 Interest factor on assets in the economy this period.
1959 wRteNow : float
1960 Wage rate for labor in the economy this period.
1961 PermShkAggNow : float
1962 Permanent shock to aggregate labor productivity this period.
1963 TranShkAggNow : float
1964 Transitory shock to aggregate labor productivity this period.
1965 KtoLnow : float
1966 Capital-to-labor ratio in the economy this period.
1967 """
1968 # Calculate aggregate savings
1969 AaggPrev = np.mean(np.array(aLvlNow)) / np.mean(
1970 pLvlNow
1971 ) # End-of-period savings from last period
1972 # Calculate aggregate capital this period
1973 AggregateK = np.mean(np.array(aLvlNow)) # ...becomes capital today
1974 # This version uses end-of-period assets and
1975 # permanent income to calculate aggregate capital, unlike the Mathematica
1976 # version, which first applies the idiosyncratic permanent income shocks
1977 # and then aggregates. Obviously this is mathematically equivalent.
1979 # Get this period's aggregate shocks
1980 PermShkAggNow = self.PermShkAggHist[self.Shk_idx]
1981 TranShkAggNow = self.TranShkAggHist[self.Shk_idx]
1982 self.Shk_idx += 1
1984 AggregateL = np.mean(pLvlNow) * PermShkAggNow
1986 # Calculate the interest factor and wage rate this period
1987 KtoLnow = AggregateK / AggregateL
1988 self.KtoYnow = KtoLnow ** (1.0 - self.CapShare)
1989 RfreeNow = self.Rfunc(KtoLnow / TranShkAggNow)
1990 wRteNow = self.wFunc(KtoLnow / TranShkAggNow)
1991 MaggNow = KtoLnow * RfreeNow + wRteNow * TranShkAggNow
1992 self.KtoLnow = KtoLnow # Need to store this as it is a sow variable
1994 # Package the results into an object and return it
1995 return (
1996 MaggNow,
1997 AaggPrev,
1998 RfreeNow,
1999 wRteNow,
2000 PermShkAggNow,
2001 TranShkAggNow,
2002 KtoLnow,
2003 )
2005 def calc_AFunc(self, MaggNow, AaggNow):
2006 """
2007 Calculate a new aggregate savings rule based on the history
2008 of the aggregate savings and aggregate market resources from a simulation.
2010 Parameters
2011 ----------
2012 MaggNow : [float]
2013 List of the history of the simulated aggregate market resources for an economy.
2014 AaggNow : [float]
2015 List of the history of the simulated aggregate savings for an economy.
2017 Returns
2018 -------
2019 (unnamed) : CapDynamicRule
2020 Object containing a new savings rule
2021 """
2022 verbose = self.verbose
2023 discard_periods = (
2024 self.T_discard
2025 ) # Throw out the first T periods to allow the simulation to approach the SS
2026 update_weight = (
2027 1.0 - self.DampingFac
2028 ) # Proportional weight to put on new function vs old function parameters
2029 total_periods = len(MaggNow)
2031 # Regress the log savings against log market resources
2032 logAagg = np.log(AaggNow[discard_periods:total_periods])
2033 logMagg = np.log(MaggNow[discard_periods - 1 : total_periods - 1])
2034 slope, intercept, r_value, p_value, std_err = stats.linregress(logMagg, logAagg)
2036 # Make a new aggregate savings rule by combining the new regression parameters
2037 # with the previous guess
2038 intercept = (
2039 update_weight * intercept + (1.0 - update_weight) * self.intercept_prev
2040 )
2041 slope = update_weight * slope + (1.0 - update_weight) * self.slope_prev
2042 AFunc = AggregateSavingRule(
2043 intercept, slope
2044 ) # Make a new next-period capital function
2046 # Save the new values as "previous" values for the next iteration
2047 self.intercept_prev = intercept
2048 self.slope_prev = slope
2050 # Print the new parameters
2051 if verbose:
2052 print(
2053 "intercept="
2054 + str(intercept)
2055 + ", slope="
2056 + str(slope)
2057 + ", r-sq="
2058 + str(r_value**2)
2059 )
2061 return AggShocksDynamicRule(AFunc)
2064class SmallOpenEconomy(Market):
2065 """
2066 A class for representing a small open economy, where the wage rate and interest rate are
2067 exogenously determined by some "global" rate. However, the economy is still subject to
2068 aggregate productivity shocks.
2070 Parameters
2071 ----------
2072 agents : [ConsumerType]
2073 List of types of consumers that live in this economy.
2074 tolerance: float
2075 Minimum acceptable distance between "dynamic rules" to consider the
2076 solution process converged. Distance depends on intercept and slope
2077 of the log-linear "next capital ratio" function.
2078 act_T : int
2079 Number of periods to simulate when making a history of of the market.
2080 """
2082 def __init__(self, agents=None, tolerance=0.0001, act_T=1000, **kwds):
2083 agents = agents if agents is not None else list()
2084 Market.__init__(
2085 self,
2086 agents=agents,
2087 sow_vars=[
2088 "MaggNow",
2089 "AaggNow",
2090 "RfreeNow",
2091 "wRteNow",
2092 "PermShkAggNow",
2093 "TranShkAggNow",
2094 "KtoLnow",
2095 ],
2096 reap_vars=[],
2097 track_vars=["MaggNow", "AaggNow", "KtoLnow"],
2098 dyn_vars=[],
2099 distributions=["PermShkAggDstn", "TranShkAggDstn", "AggShkDstn"],
2100 tolerance=tolerance,
2101 act_T=act_T,
2102 )
2103 self.assign_parameters(**kwds)
2104 self.update()
2106 def update(self):
2107 """
2108 Use primitive parameters to set basic objects.
2109 This is an extremely stripped-down version
2110 of update for CobbDouglasEconomy.
2112 Parameters
2113 ----------
2114 none
2116 Returns
2117 -------
2118 none
2119 """
2120 self.kSS = 1.0
2121 self.MSS = 1.0
2122 self.sow_init["KtoLnow_init"] = self.kSS
2123 self.Rfunc = ConstantFunction(self.Rfree)
2124 self.wFunc = ConstantFunction(self.wRte)
2125 self.sow_init["RfreeNow"] = self.Rfunc(self.kSS)
2126 self.sow_init["wRteNow"] = self.wFunc(self.kSS)
2127 self.sow_init["MaggNow"] = self.kSS
2128 self.sow_init["AaggNow"] = self.kSS
2129 self.sow_init["PermShkAggNow"] = 1.0
2130 self.sow_init["TranShkAggNow"] = 1.0
2131 self.make_AggShkDstn()
2132 self.AFunc = ConstantFunction(1.0)
2134 def make_AggShkDstn(self):
2135 """
2136 Creates the attributes TranShkAggDstn, PermShkAggDstn, and AggShkDstn.
2137 Draws on attributes TranShkAggStd, PermShkAddStd, TranShkAggCount, PermShkAggCount.
2139 Parameters
2140 ----------
2141 None
2143 Returns
2144 -------
2145 None
2146 """
2147 RNG = self.RNG
2148 TranShkAggDstn = MeanOneLogNormal(
2149 sigma=self.TranShkAggStd, seed=RNG.integers(0, 2**31 - 1)
2150 )
2151 self.TranShkAggDstn = TranShkAggDstn.discretize(N=self.TranShkAggCount)
2152 PermShkAggDstn = MeanOneLogNormal(
2153 sigma=self.PermShkAggStd, seed=RNG.integers(0, 2**31 - 1)
2154 )
2155 self.PermShkAggDstn = PermShkAggDstn.discretize(N=self.PermShkAggCount)
2156 self.AggShkDstn = combine_indep_dstns(
2157 self.PermShkAggDstn, self.TranShkAggDstn, seed=RNG.integers(0, 2**31 - 1)
2158 )
2160 def mill_rule(self):
2161 """
2162 No aggregation occurs for a small open economy, because the wage and interest rates are
2163 exogenously determined. However, aggregate shocks may occur.
2165 See documentation for get_AggShocks() for more information.
2166 """
2167 return self.get_AggShocks()
2169 def calc_dynamics(self, KtoLnow):
2170 """
2171 Calculates a new dynamic rule for the economy, which is just an empty object.
2172 There is no "dynamic rule" for a small open economy, because K/L does not generate w and R.
2173 """
2174 return MetricObject()
2176 def reset(self):
2177 """
2178 Reset the economy to prepare for a new simulation. Sets the time index of aggregate shocks
2179 to zero and runs Market.reset(). This replicates the reset method for CobbDouglasEconomy;
2180 future version should create parent class of that class and this one.
2182 Parameters
2183 ----------
2184 None
2186 Returns
2187 -------
2188 None
2189 """
2190 self.Shk_idx = 0
2191 Market.reset(self)
2193 def make_AggShkHist(self):
2194 """
2195 Make simulated histories of aggregate transitory and permanent shocks. Histories are of
2196 length self.act_T, for use in the general equilibrium simulation. This replicates the same
2197 method for CobbDouglasEconomy; future version should create parent class.
2199 Parameters
2200 ----------
2201 None
2203 Returns
2204 -------
2205 None
2206 """
2207 sim_periods = self.act_T
2208 Events = np.arange(self.AggShkDstn.pmv.size) # just a list of integers
2209 EventDraws = self.AggShkDstn.draw(N=sim_periods, atoms=Events)
2210 PermShkAggHist = self.AggShkDstn.atoms[0][EventDraws]
2211 TranShkAggHist = self.AggShkDstn.atoms[1][EventDraws]
2213 # Store the histories
2214 self.PermShkAggHist = PermShkAggHist
2215 self.TranShkAggHist = TranShkAggHist
2217 def get_AggShocks(self):
2218 """
2219 Returns aggregate state variables and shocks for this period. The capital-to-labor ratio
2220 is irrelevant and thus treated as constant, and the wage and interest rates are also
2221 constant. However, aggregate shocks are assigned from a prespecified history.
2223 Parameters
2224 ----------
2225 None
2227 Returns
2228 -------
2229 MaggNow : float
2230 Aggregate market resources for this period normalized by mean permanent income
2231 AaggNow : float
2232 Aggregate savings for this period normalized by mean permanent income
2233 RfreeNow : float
2234 Interest factor on assets in the economy this period.
2235 wRteNow : float
2236 Wage rate for labor in the economy this period.
2237 PermShkAggNow : float
2238 Permanent shock to aggregate labor productivity this period.
2239 TranShkAggNow : float
2240 Transitory shock to aggregate labor productivity this period.
2241 KtoLnow : float
2242 Capital-to-labor ratio in the economy this period.
2244 """
2245 # Get this period's aggregate shocks
2246 PermShkAggNow = self.PermShkAggHist[self.Shk_idx]
2247 TranShkAggNow = self.TranShkAggHist[self.Shk_idx]
2248 self.Shk_idx += 1
2250 # Factor prices are constant
2251 RfreeNow = self.Rfunc(1.0 / TranShkAggNow)
2252 wRteNow = self.wFunc(1.0 / TranShkAggNow)
2254 # Aggregates are irrelavent
2255 AaggNow = 1.0
2256 MaggNow = 1.0
2257 KtoLnow = 1.0 / PermShkAggNow
2259 return (
2260 MaggNow,
2261 AaggNow,
2262 RfreeNow,
2263 wRteNow,
2264 PermShkAggNow,
2265 TranShkAggNow,
2266 KtoLnow,
2267 )
2270# Make a dictionary to specify a Markov Cobb-Douglas economy
2271init_mrkv_cobb_douglas = init_cobb_douglas.copy()
2272init_mrkv_cobb_douglas["PermShkAggStd"] = [0.012, 0.006]
2273init_mrkv_cobb_douglas["TranShkAggStd"] = [0.006, 0.003]
2274init_mrkv_cobb_douglas["PermGroFacAgg"] = [0.98, 1.02]
2275init_mrkv_cobb_douglas["MrkvArray"] = MrkvArray
2276init_mrkv_cobb_douglas["MrkvNow_init"] = 0
2277init_mrkv_cobb_douglas["slope_prev"] = 2 * [slope_prev]
2278init_mrkv_cobb_douglas["intercept_prev"] = 2 * [intercept_prev]
2281class CobbDouglasMarkovEconomy(CobbDouglasEconomy):
2282 """
2283 A class to represent an economy with a Cobb-Douglas aggregate production
2284 function over labor and capital, extending HARK.Market. The "aggregate
2285 market process" for this market combines all individuals' asset holdings
2286 into aggregate capital, yielding the interest factor on assets and the wage
2287 rate for the upcoming period. This small extension incorporates a Markov
2288 state for the "macroeconomy", so that the shock distribution and aggregate
2289 productivity growth factor can vary over time.
2291 Parameters
2292 ----------
2293 agents : [ConsumerType]
2294 List of types of consumers that live in this economy.
2295 tolerance: float
2296 Minimum acceptable distance between "dynamic rules" to consider the
2297 solution process converged. Distance depends on intercept and slope
2298 of the log-linear "next capital ratio" function.
2299 act_T : int
2300 Number of periods to simulate when making a history of of the market.
2301 """
2303 def __init__(
2304 self,
2305 agents=None,
2306 tolerance=0.0001,
2307 act_T=1200,
2308 **kwds,
2309 ):
2310 agents = agents if agents is not None else list()
2311 params = init_mrkv_cobb_douglas.copy()
2312 params.update(kwds)
2314 CobbDouglasEconomy.__init__(
2315 self,
2316 agents=agents,
2317 tolerance=tolerance,
2318 act_T=act_T,
2319 sow_vars=[
2320 "MaggNow",
2321 "AaggNow",
2322 "RfreeNow",
2323 "wRteNow",
2324 "PermShkAggNow",
2325 "TranShkAggNow",
2326 "KtoLnow",
2327 "Mrkv", # This one is new
2328 ],
2329 **params,
2330 )
2332 self.sow_init["Mrkv"] = params["MrkvNow_init"]
2334 def update(self):
2335 """
2336 Use primitive parameters (and perfect foresight calibrations) to make
2337 interest factor and wage rate functions (of capital to labor ratio),
2338 as well as discrete approximations to the aggregate shock distributions.
2340 Parameters
2341 ----------
2342 None
2344 Returns
2345 -------
2346 None
2347 """
2348 CobbDouglasEconomy.update(self)
2349 StateCount = self.MrkvArray.shape[0]
2350 AFunc_all = []
2351 for i in range(StateCount):
2352 AFunc_all.append(
2353 AggregateSavingRule(self.intercept_prev[i], self.slope_prev[i])
2354 )
2355 self.AFunc = AFunc_all
2357 def get_PermGroFacAggLR(self):
2358 """
2359 Calculates and returns the long run permanent income growth factor. This
2360 is the average growth factor in self.PermGroFacAgg, weighted by the long
2361 run distribution of Markov states (as determined by self.MrkvArray).
2363 Parameters
2364 ----------
2365 None
2367 Returns
2368 -------
2369 PermGroFacAggLR : float
2370 Long run aggregate permanent income growth factor
2371 """
2372 # Find the long run distribution of Markov states
2373 w, v = np.linalg.eig(np.transpose(self.MrkvArray))
2374 idx = (np.abs(w - 1.0)).argmin()
2375 x = v[:, idx].astype(float)
2376 LR_dstn = x / np.sum(x)
2378 # Return the weighted average of aggregate permanent income growth factors
2379 PermGroFacAggLR = np.dot(LR_dstn, np.array(self.PermGroFacAgg))
2380 return PermGroFacAggLR
2382 def make_AggShkDstn(self):
2383 """
2384 Creates the attributes TranShkAggDstn, PermShkAggDstn, and AggShkDstn.
2385 Draws on attributes TranShkAggStd, PermShkAddStd, TranShkAggCount, PermShkAggCount.
2386 This version accounts for the Markov macroeconomic state.
2388 Parameters
2389 ----------
2390 None
2392 Returns
2393 -------
2394 None
2395 """
2396 TranShkAggDstn = []
2397 PermShkAggDstn = []
2398 AggShkDstn = []
2399 StateCount = self.MrkvArray.shape[0]
2400 RNG = self.RNG
2402 for i in range(StateCount):
2403 TranShkAggDstn.append(
2404 MeanOneLogNormal(
2405 sigma=self.TranShkAggStd[i], seed=RNG.integers(0, 2**31 - 1)
2406 ).discretize(
2407 N=self.TranShkAggCount,
2408 )
2409 )
2410 PermShkAggDstn.append(
2411 MeanOneLogNormal(
2412 sigma=self.PermShkAggStd[i], seed=RNG.integers(0, 2**31 - 1)
2413 ).discretize(
2414 N=self.PermShkAggCount,
2415 )
2416 )
2417 AggShkDstn.append(
2418 combine_indep_dstns(
2419 PermShkAggDstn[-1],
2420 TranShkAggDstn[-1],
2421 seed=RNG.integers(0, 2**31 - 1),
2422 )
2423 )
2425 self.TranShkAggDstn = TranShkAggDstn
2426 self.PermShkAggDstn = PermShkAggDstn
2427 self.AggShkDstn = AggShkDstn
2429 def make_AggShkHist(self):
2430 """
2431 Make simulated histories of aggregate transitory and permanent shocks.
2432 Histories are of length self.act_T, for use in the general equilibrium
2433 simulation. Draws on history of aggregate Markov states generated by
2434 internal call to make_Mrkv_history().
2436 Parameters
2437 ----------
2438 None
2440 Returns
2441 -------
2442 None
2443 """
2444 self.make_Mrkv_history() # Make a (pseudo)random sequence of Markov states
2445 sim_periods = self.act_T
2447 # For each Markov state in each simulated period, draw the aggregate shocks
2448 # that would occur in that state in that period
2449 StateCount = self.MrkvArray.shape[0]
2450 PermShkAggHistAll = np.zeros((StateCount, sim_periods))
2451 TranShkAggHistAll = np.zeros((StateCount, sim_periods))
2452 for i in range(StateCount):
2453 AggShockDraws = self.AggShkDstn[i].draw(N=sim_periods)
2454 PermShkAggHistAll[i, :] = AggShockDraws[0, :]
2455 TranShkAggHistAll[i, :] = AggShockDraws[1, :]
2457 # Select the actual history of aggregate shocks based on the sequence
2458 # of Markov states that the economy experiences
2459 PermShkAggHist = np.zeros(sim_periods)
2460 TranShkAggHist = np.zeros(sim_periods)
2461 for i in range(StateCount):
2462 these = i == self.MrkvNow_hist
2463 PermShkAggHist[these] = PermShkAggHistAll[i, these] * self.PermGroFacAgg[i]
2464 TranShkAggHist[these] = TranShkAggHistAll[i, these]
2466 # Store the histories
2467 self.PermShkAggHist = PermShkAggHist
2468 self.TranShkAggHist = TranShkAggHist
2470 def make_Mrkv_history(self):
2471 """
2472 Makes a history of macroeconomic Markov states, stored in the attribute
2473 MrkvNow_hist. This version ensures that each state is reached a sufficient
2474 number of times to have a valid sample for calc_dynamics to produce a good
2475 dynamic rule. It will sometimes cause act_T to be increased beyond its
2476 initially specified level.
2478 Parameters
2479 ----------
2480 None
2482 Returns
2483 -------
2484 None
2485 """
2486 loops_max = getattr(self, "loops_max", 10)
2487 state_T_min = 50 # Choose minimum number of periods in each state for a valid Markov sequence
2488 logit_scale = (
2489 0.2 # Scaling factor on logit choice shocks when jumping to a new state
2490 )
2491 # Values close to zero make the most underrepresented states very likely to visit, while
2492 # large values of logit_scale make any state very likely to be jumped to.
2494 # Reset act_T to the level actually specified by the user
2495 if hasattr(self, "act_T_orig"):
2496 act_T = self.act_T_orig
2497 else: # Or store it for the first time
2498 self.act_T_orig = self.act_T
2499 act_T = self.act_T
2501 # Find the long run distribution of Markov states
2502 w, v = np.linalg.eig(np.transpose(self.MrkvArray))
2503 idx = (np.abs(w - 1.0)).argmin()
2504 x = v[:, idx].astype(float)
2505 LR_dstn = x / np.sum(x)
2507 # Initialize the Markov history and set up transitions
2508 MrkvNow_hist = np.zeros(self.act_T_orig, dtype=int)
2509 loops = 0
2510 go = True
2511 MrkvNow = self.sow_init["Mrkv"]
2512 t = 0
2513 StateCount = self.MrkvArray.shape[0]
2515 # Add histories until each state has been visited at least state_T_min times
2516 while go:
2517 draws = Uniform(seed=loops).draw(N=self.act_T_orig)
2518 markov_process = MarkovProcess(self.MrkvArray, seed=loops)
2519 for s in range(self.act_T_orig): # Add act_T_orig more periods
2520 MrkvNow_hist[t] = MrkvNow
2521 MrkvNow = markov_process.draw(MrkvNow)
2522 t += 1
2524 # Calculate the empirical distribution
2525 state_T = np.zeros(StateCount)
2526 for i in range(StateCount):
2527 state_T[i] = np.sum(MrkvNow_hist == i)
2529 # Check whether each state has been visited state_T_min times
2530 if np.all(state_T >= state_T_min):
2531 go = False # If so, terminate the loop
2532 continue
2534 # Choose an underrepresented state to "jump" to
2535 if np.any(
2536 state_T == 0
2537 ): # If any states have *never* been visited, randomly choose one of those
2538 never_visited = np.where(np.array(state_T == 0))[0]
2539 MrkvNow = np.random.choice(never_visited)
2540 else: # Otherwise, use logit choice probabilities to visit an underrepresented state
2541 emp_dstn = state_T / act_T
2542 ratios = LR_dstn / emp_dstn
2543 ratios_adj = ratios - np.max(ratios)
2544 ratios_exp = np.exp(ratios_adj / logit_scale)
2545 ratios_sum = np.sum(ratios_exp)
2546 jump_probs = ratios_exp / ratios_sum
2547 cum_probs = np.cumsum(jump_probs)
2548 MrkvNow = np.searchsorted(cum_probs, draws[-1])
2550 loops += 1
2551 # Make the Markov state history longer by act_T_orig periods
2552 if loops >= loops_max:
2553 go = False
2554 print(
2555 "make_Mrkv_history reached maximum number of loops without generating a valid sequence!"
2556 )
2557 else:
2558 MrkvNow_new = np.zeros(self.act_T_orig, dtype=int)
2559 MrkvNow_hist = np.concatenate((MrkvNow_hist, MrkvNow_new))
2560 act_T += self.act_T_orig
2562 # Store the results as attributes of self
2563 self.MrkvNow_hist = MrkvNow_hist
2564 self.act_T = act_T
2566 def mill_rule(self, aLvl, pLvl):
2567 """
2568 Function to calculate the capital to labor ratio, interest factor, and
2569 wage rate based on each agent's current state. Just calls calc_R_and_W()
2570 and adds the Markov state index.
2572 See documentation for calc_R_and_W for more information.
2574 Params
2575 -------
2576 aLvl : float
2577 pLvl : float
2579 Returns
2580 -------
2581 Mnow : float
2582 Aggregate market resources for this period.
2583 Aprev : float
2584 Aggregate savings for the prior period.
2585 KtoLnow : float
2586 Capital-to-labor ratio in the economy this period.
2587 Rnow : float
2588 Interest factor on assets in the economy this period.
2589 Wnow : float
2590 Wage rate for labor in the economy this period.
2591 MrkvNow : int
2592 Binary indicator for bad (0) or good (1) macroeconomic state.
2593 """
2594 MrkvNow = self.MrkvNow_hist[self.Shk_idx]
2595 temp = self.calc_R_and_W(aLvl, pLvl)
2597 return temp + (MrkvNow,)
2599 def calc_AFunc(self, MaggNow, AaggNow):
2600 """
2601 Calculate a new aggregate savings rule based on the history of the
2602 aggregate savings and aggregate market resources from a simulation.
2603 Calculates an aggregate saving rule for each macroeconomic Markov state.
2605 Parameters
2606 ----------
2607 MaggNow : [float]
2608 List of the history of the simulated aggregate market resources for an economy.
2609 AaggNow : [float]
2610 List of the history of the simulated aggregate savings for an economy.
2612 Returns
2613 -------
2614 (unnamed) : CapDynamicRule
2615 Object containing new saving rules for each Markov state.
2616 """
2617 verbose = self.verbose
2618 discard_periods = (
2619 self.T_discard
2620 ) # Throw out the first T periods to allow the simulation to approach the SS
2621 update_weight = (
2622 1.0 - self.DampingFac
2623 ) # Proportional weight to put on new function vs old function parameters
2624 total_periods = len(MaggNow)
2626 # Trim the histories of M_t and A_t and convert them to logs
2627 logAagg = np.log(AaggNow[discard_periods:total_periods])
2628 logMagg = np.log(MaggNow[discard_periods - 1 : total_periods - 1])
2629 MrkvHist = self.MrkvNow_hist[discard_periods - 1 : total_periods - 1]
2631 # For each Markov state, regress A_t on M_t and update the saving rule
2632 AFunc_list = []
2633 rSq_list = []
2634 for i in range(self.MrkvArray.shape[0]):
2635 these = i == MrkvHist
2636 slope, intercept, r_value, p_value, std_err = stats.linregress(
2637 logMagg[these], logAagg[these]
2638 )
2640 # Make a new aggregate savings rule by combining the new regression parameters
2641 # with the previous guess
2642 intercept = (
2643 update_weight * intercept
2644 + (1.0 - update_weight) * self.intercept_prev[i]
2645 )
2646 slope = update_weight * slope + (1.0 - update_weight) * self.slope_prev[i]
2647 AFunc_list.append(
2648 AggregateSavingRule(intercept, slope)
2649 ) # Make a new next-period capital function
2650 rSq_list.append(r_value**2)
2652 # Save the new values as "previous" values for the next iteration
2653 self.intercept_prev[i] = intercept
2654 self.slope_prev[i] = slope
2656 # Print the new parameters
2657 if verbose: # pragma: nocover
2658 print(
2659 "intercept="
2660 + str(self.intercept_prev)
2661 + ", slope="
2662 + str(self.slope_prev)
2663 + ", r-sq="
2664 + str(rSq_list)
2665 )
2667 return AggShocksDynamicRule(AFunc_list)
2670class SmallOpenMarkovEconomy(CobbDouglasMarkovEconomy, SmallOpenEconomy):
2671 """
2672 A class for representing a small open economy, where the wage rate and interest rate are
2673 exogenously determined by some "global" rate. However, the economy is still subject to
2674 aggregate productivity shocks. This version supports a discrete Markov state. All
2675 methods in this class inherit from the two parent classes.
2676 """
2678 def __init__(self, agents=None, tolerance=0.0001, act_T=1000, **kwds):
2679 agents = agents if agents is not None else list()
2680 temp_dict = init_mrkv_cobb_douglas.copy()
2681 temp_dict.update(kwds)
2682 CobbDouglasMarkovEconomy.__init__(
2683 self,
2684 agents=agents,
2685 tolerance=tolerance,
2686 act_T=act_T,
2687 reap_vars=[],
2688 dyn_vars=[],
2689 **temp_dict,
2690 )
2692 def update(self):
2693 SmallOpenEconomy.update(self)
2694 StateCount = self.MrkvArray.shape[0]
2695 self.AFunc = StateCount * [IdentityFunction()]
2697 def make_AggShkDstn(self):
2698 CobbDouglasMarkovEconomy.make_AggShkDstn(self)
2700 def mill_rule(self):
2701 MrkvNow = self.MrkvNow_hist[self.Shk_idx]
2702 temp = SmallOpenEconomy.get_AggShocks(self)
2703 temp += (MrkvNow,)
2704 return temp
2706 def calc_dynamics(self):
2707 return NullFunc()
2709 def make_AggShkHist(self):
2710 CobbDouglasMarkovEconomy.make_AggShkHist(self)
2713init_KS_economy = {
2714 "verbose": True,
2715 "act_T": 11000,
2716 "T_discard": 1000,
2717 "DampingFac": 0.5,
2718 "intercept_prev": [0.0, 0.0],
2719 "slope_prev": [1.0, 1.0],
2720 "DiscFac": 0.99,
2721 "CRRA": 1.0,
2722 # Not listed in KS (1998), but Alan Lujan got this number indirectly from KS
2723 "LbrInd": 0.3271,
2724 "ProdB": 0.99,
2725 "ProdG": 1.01,
2726 "CapShare": 0.36,
2727 "DeprRte": 0.025,
2728 "DurMeanB": 8.0,
2729 "DurMeanG": 8.0,
2730 "SpellMeanB": 2.5,
2731 "SpellMeanG": 1.5,
2732 "UrateB": 0.10,
2733 "UrateG": 0.04,
2734 "RelProbBG": 0.75,
2735 "RelProbGB": 1.25,
2736 "MrkvNow_init": 0,
2737}
2740class KrusellSmithEconomy(Market):
2741 """
2742 A class to represent an economy in the special Krusell-Smith (1998) model.
2743 This model replicates the one presented in the JPE article "Income and Wealth
2744 Heterogeneity in the Macroeconomy", with its default parameters set to match
2745 those in the paper.
2747 Parameters
2748 ----------
2749 agents : [ConsumerType]
2750 List of types of consumers that live in this economy.
2751 tolerance: float
2752 Minimum acceptable distance between "dynamic rules" to consider the
2753 solution process converged. Distance depends on intercept and slope
2754 of the log-linear "next capital ratio" function.
2755 act_T : int
2756 Number of periods to simulate when making a history of of the market.
2757 """
2759 def __init__(self, agents=None, tolerance=0.0001, **kwds):
2760 agents = agents if agents is not None else list()
2761 params = deepcopy(init_KS_economy)
2762 params.update(kwds)
2764 Market.__init__(
2765 self,
2766 agents=agents,
2767 tolerance=tolerance,
2768 sow_vars=["Mnow", "Aprev", "Mrkv", "Rnow", "Wnow"],
2769 reap_vars=["aNow", "EmpNow"],
2770 track_vars=["Mrkv", "Aprev", "Mnow", "Urate"],
2771 dyn_vars=["AFunc"],
2772 **params,
2773 )
2774 self.update()
2776 def update(self):
2777 """
2778 Construct trivial initial guesses of the aggregate saving rules, as well
2779 as the perfect foresight steady state and associated objects.
2780 """
2781 StateCount = 2
2782 AFunc_all = [
2783 AggregateSavingRule(self.intercept_prev[j], self.slope_prev[j])
2784 for j in range(StateCount)
2785 ]
2786 self.AFunc = AFunc_all
2787 self.KtoLSS = (
2788 (1.0**self.CRRA / self.DiscFac - (1.0 - self.DeprRte)) / self.CapShare
2789 ) ** (1.0 / (self.CapShare - 1.0))
2790 self.KSS = self.KtoLSS * self.LbrInd
2791 self.KtoYSS = self.KtoLSS ** (1.0 - self.CapShare)
2792 self.WSS = (1.0 - self.CapShare) * self.KtoLSS ** (self.CapShare)
2793 self.RSS = (
2794 1.0 + self.CapShare * self.KtoLSS ** (self.CapShare - 1.0) - self.DeprRte
2795 )
2796 self.MSS = self.KSS * self.RSS + self.WSS * self.LbrInd
2797 self.convertKtoY = lambda KtoY: KtoY ** (
2798 1.0 / (1.0 - self.CapShare)
2799 ) # converts K/Y to K/L
2800 self.rFunc = lambda k: self.CapShare * k ** (self.CapShare - 1.0)
2801 self.Wfunc = lambda k: ((1.0 - self.CapShare) * k ** (self.CapShare))
2802 self.sow_init["KtoLnow"] = self.KtoLSS
2803 self.sow_init["Mnow"] = self.MSS
2804 self.sow_init["Aprev"] = self.KSS
2805 self.sow_init["Rnow"] = self.RSS
2806 self.sow_init["Wnow"] = self.WSS
2807 self.PermShkAggNow_init = 1.0
2808 self.TranShkAggNow_init = 1.0
2809 self.sow_init["Mrkv"] = 0
2810 self.make_MrkvArray()
2812 def reset(self):
2813 """
2814 Reset the economy to prepare for a new simulation. Sets the time index
2815 of aggregate shocks to zero and runs Market.reset().
2816 """
2817 self.Shk_idx = 0
2818 Market.reset(self)
2820 def make_MrkvArray(self):
2821 """
2822 Construct the attributes MrkvAggArray and MrkvIndArray from the primitive
2823 attributes DurMeanB, DurMeanG, SpellMeanB, SpellMeanG, UrateB, UrateG,
2824 RelProbGB, and RelProbBG.
2825 """
2826 # Construct aggregate Markov transition probabilities
2827 ProbBG = 1.0 / self.DurMeanB
2828 ProbGB = 1.0 / self.DurMeanG
2829 ProbBB = 1.0 - ProbBG
2830 ProbGG = 1.0 - ProbGB
2831 MrkvAggArray = np.array([[ProbBB, ProbBG], [ProbGB, ProbGG]])
2833 # Construct idiosyncratic Markov transition probabilities
2834 # ORDER: BU, BE, GU, GE
2835 MrkvIndArray = np.zeros((4, 4))
2837 # BAD-BAD QUADRANT
2838 MrkvIndArray[0, 1] = ProbBB * 1.0 / self.SpellMeanB
2839 MrkvIndArray[0, 0] = ProbBB * (1 - 1.0 / self.SpellMeanB)
2840 MrkvIndArray[1, 0] = self.UrateB / (1.0 - self.UrateB) * MrkvIndArray[0, 1]
2841 MrkvIndArray[1, 1] = ProbBB - MrkvIndArray[1, 0]
2843 # GOOD-GOOD QUADRANT
2844 MrkvIndArray[2, 3] = ProbGG * 1.0 / self.SpellMeanG
2845 MrkvIndArray[2, 2] = ProbGG * (1 - 1.0 / self.SpellMeanG)
2846 MrkvIndArray[3, 2] = self.UrateG / (1.0 - self.UrateG) * MrkvIndArray[2, 3]
2847 MrkvIndArray[3, 3] = ProbGG - MrkvIndArray[3, 2]
2849 # BAD-GOOD QUADRANT
2850 MrkvIndArray[0, 2] = self.RelProbBG * MrkvIndArray[2, 2] / ProbGG * ProbBG
2851 MrkvIndArray[0, 3] = ProbBG - MrkvIndArray[0, 2]
2852 MrkvIndArray[1, 2] = (
2853 ProbBG * self.UrateG - self.UrateB * MrkvIndArray[0, 2]
2854 ) / (1.0 - self.UrateB)
2855 MrkvIndArray[1, 3] = ProbBG - MrkvIndArray[1, 2]
2857 # GOOD-BAD QUADRANT
2858 MrkvIndArray[2, 0] = self.RelProbGB * MrkvIndArray[0, 0] / ProbBB * ProbGB
2859 MrkvIndArray[2, 1] = ProbGB - MrkvIndArray[2, 0]
2860 MrkvIndArray[3, 0] = (
2861 ProbGB * self.UrateB - self.UrateG * MrkvIndArray[2, 0]
2862 ) / (1.0 - self.UrateG)
2863 MrkvIndArray[3, 1] = ProbGB - MrkvIndArray[3, 0]
2865 # Test for valid idiosyncratic transition probabilities
2866 assert np.all(MrkvIndArray >= 0.0), (
2867 "Invalid idiosyncratic transition probabilities!"
2868 )
2869 self.MrkvArray = MrkvAggArray
2870 self.MrkvIndArray = MrkvIndArray
2872 def make_Mrkv_history(self):
2873 """
2874 Makes a history of macroeconomic Markov states, stored in the attribute
2875 MrkvNow_hist. This variable is binary (0 bad, 1 good) in the KS model.
2876 """
2877 # Initialize the Markov history and set up transitions
2878 self.MrkvNow_hist = np.zeros(self.act_T, dtype=int)
2879 MrkvNow = self.MrkvNow_init
2881 markov_process = MarkovProcess(self.MrkvArray, seed=0)
2882 for s in range(self.act_T): # Add act_T_orig more periods
2883 self.MrkvNow_hist[s] = MrkvNow
2884 MrkvNow = markov_process.draw(MrkvNow)
2886 def mill_rule(self, aNow, EmpNow):
2887 """
2888 Method to calculate the capital to labor ratio, interest factor, and
2889 wage rate based on each agent's current state. Just calls calc_R_and_W().
2891 See documentation for calc_R_and_W for more information.
2893 Returns
2894 -------
2895 Mnow : float
2896 Aggregate market resources for this period.
2897 Aprev : float
2898 Aggregate savings for the prior period.
2899 MrkvNow : int
2900 Binary indicator for bad (0) or good (1) macroeconomic state.
2901 Rnow : float
2902 Interest factor on assets in the economy this period.
2903 Wnow : float
2904 Wage rate for labor in the economy this period.
2905 """
2907 return self.calc_R_and_W(aNow, EmpNow)
2909 def calc_dynamics(self, Mnow, Aprev):
2910 """
2911 Method to update perceptions of the aggregate saving rule in each
2912 macroeconomic state; just calls calc_AFunc.
2913 """
2914 return self.calc_AFunc(Mnow, Aprev)
2916 def calc_R_and_W(self, aNow, EmpNow):
2917 """
2918 Calculates the interest factor and wage rate this period using each agent's
2919 capital stock to get the aggregate capital ratio.
2921 Parameters
2922 ----------
2923 aNow : [np.array]
2924 Agents' current end-of-period assets. Elements of the list correspond
2925 to types in the economy, entries within arrays to agents of that type.
2926 EmpNow [np.array]
2927 Agents' binary employment states. Not actually used in computation of
2928 interest and wage rates, but stored in the history to verify that the
2929 idiosyncratic unemployment probabilities are behaving as expected.
2931 Returns
2932 -------
2933 Mnow : float
2934 Aggregate market resources for this period.
2935 Aprev : float
2936 Aggregate savings for the prior period.
2937 MrkvNow : int
2938 Binary indicator for bad (0) or good (1) macroeconomic state.
2939 Rnow : float
2940 Interest factor on assets in the economy this period.
2941 Wnow : float
2942 Wage rate for labor in the economy this period.
2943 """
2944 # Calculate aggregate savings
2945 # End-of-period savings from last period
2946 Aprev = np.mean(np.array(aNow))
2947 # Calculate aggregate capital this period
2948 AggK = Aprev # ...becomes capital today
2950 # Calculate unemployment rate
2951 Urate = 1.0 - np.mean(np.array(EmpNow))
2952 self.Urate = Urate # This is the unemployment rate for the *prior* period
2954 # Get this period's TFP and labor supply
2955 MrkvNow = self.MrkvNow_hist[self.Shk_idx]
2956 if MrkvNow == 0:
2957 Prod = self.ProdB
2958 AggL = (1.0 - self.UrateB) * self.LbrInd
2959 elif MrkvNow == 1:
2960 Prod = self.ProdG
2961 AggL = (1.0 - self.UrateG) * self.LbrInd
2962 self.Shk_idx += 1
2964 # Calculate the interest factor and wage rate this period
2965 KtoLnow = AggK / AggL
2966 Rnow = 1.0 + Prod * self.rFunc(KtoLnow) - self.DeprRte
2967 Wnow = Prod * self.Wfunc(KtoLnow)
2968 Mnow = Rnow * AggK + Wnow * AggL
2969 self.KtoLnow = KtoLnow # Need to store this as it is a sow variable
2971 # Returns a tuple of these values
2972 return Mnow, Aprev, MrkvNow, Rnow, Wnow
2974 def calc_AFunc(self, Mnow, Aprev):
2975 """
2976 Calculate a new aggregate savings rule based on the history of the
2977 aggregate savings and aggregate market resources from a simulation.
2978 Calculates an aggregate saving rule for each macroeconomic Markov state.
2980 Parameters
2981 ----------
2982 Mnow : [float]
2983 List of the history of the simulated aggregate market resources for an economy.
2984 Anow : [float]
2985 List of the history of the simulated aggregate savings for an economy.
2987 Returns
2988 -------
2989 (unnamed) : CapDynamicRule
2990 Object containing new saving rules for each Markov state.
2991 """
2992 verbose = self.verbose
2993 discard_periods = (
2994 self.T_discard
2995 ) # Throw out the first T periods to allow the simulation to approach the SS
2996 update_weight = (
2997 1.0 - self.DampingFac
2998 ) # Proportional weight to put on new function vs old function parameters
2999 total_periods = len(Mnow)
3001 # Trim the histories of M_t and A_t and convert them to logs
3002 logAagg = np.log(Aprev[discard_periods:total_periods])
3003 logMagg = np.log(Mnow[discard_periods - 1 : total_periods - 1])
3004 MrkvHist = self.MrkvNow_hist[discard_periods - 1 : total_periods - 1]
3006 # For each Markov state, regress A_t on M_t and update the saving rule
3007 AFunc_list = []
3008 rSq_list = []
3009 for i in range(self.MrkvArray.shape[0]):
3010 these = i == MrkvHist
3011 slope, intercept, r_value, p_value, std_err = stats.linregress(
3012 logMagg[these], logAagg[these]
3013 )
3015 # Make a new aggregate savings rule by combining the new regression parameters
3016 # with the previous guess
3017 intercept = (
3018 update_weight * intercept
3019 + (1.0 - update_weight) * self.intercept_prev[i]
3020 )
3021 slope = update_weight * slope + (1.0 - update_weight) * self.slope_prev[i]
3022 AFunc_list.append(
3023 AggregateSavingRule(intercept, slope)
3024 ) # Make a new next-period capital function
3025 rSq_list.append(r_value**2)
3027 # Save the new values as "previous" values for the next iteration
3028 self.intercept_prev[i] = intercept
3029 self.slope_prev[i] = slope
3031 # Print the new parameters
3032 if verbose: # pragma: nocover
3033 print(
3034 "intercept="
3035 + str(self.intercept_prev)
3036 + ", slope="
3037 + str(self.slope_prev)
3038 + ", r-sq="
3039 + str(rSq_list)
3040 )
3042 return AggShocksDynamicRule(AFunc_list)
3045class AggregateSavingRule(MetricObject):
3046 """
3047 A class to represent agent beliefs about aggregate saving at the end of this period (AaggNow) as
3048 a function of (normalized) aggregate market resources at the beginning of the period (MaggNow).
3050 Parameters
3051 ----------
3052 intercept : float
3053 Intercept of the log-linear capital evolution rule.
3054 slope : float
3055 Slope of the log-linear capital evolution rule.
3056 """
3058 def __init__(self, intercept, slope):
3059 self.intercept = intercept
3060 self.slope = slope
3061 self.distance_criteria = ["slope", "intercept"]
3063 def __call__(self, Mnow):
3064 """
3065 Evaluates aggregate savings as a function of the aggregate market resources this period.
3067 Parameters
3068 ----------
3069 Mnow : float
3070 Aggregate market resources this period.
3072 Returns
3073 -------
3074 Aagg : Expected aggregate savings this period.
3075 """
3076 Aagg = np.exp(self.intercept + self.slope * np.log(Mnow))
3077 return Aagg
3080class AggShocksDynamicRule(MetricObject):
3081 """
3082 Just a container class for passing the dynamic rule in the aggregate shocks model to agents.
3084 Parameters
3085 ----------
3086 AFunc : CapitalEvoRule
3087 Aggregate savings as a function of aggregate market resources.
3088 """
3090 def __init__(self, AFunc):
3091 self.AFunc = AFunc
3092 self.distance_criteria = ["AFunc"]