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