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