@@ -353,10 +353,11 @@ u_prime_inv = lambda c, γ: c**(-1/γ)
353353``` {code-cell} python3
354354def K(σ: jnp.ndarray, ifp: IFP) -> jnp.ndarray:
355355 """
356- The Coleman-Reffett operator for the IFP model using the Endogenous Grid Method.
356+ The Coleman-Reffett operator for the IFP model using the
357+ Endogenous Grid Method.
357358
358- This operator implements one iteration of the EGM algorithm to update the
359- consumption policy function.
359+ This operator implements one iteration of the EGM algorithm to
360+ update the consumption policy function.
360361
361362 Parameters
362363 ----------
@@ -374,8 +375,10 @@ def K(σ: jnp.ndarray, ifp: IFP) -> jnp.ndarray:
374375 Algorithm
375376 ---------
376377 The EGM works backwards from next period:
377- 1. Given σ(a', z'), compute current consumption c that satisfies Euler equation
378- 2. Compute the endogenous current asset level a that leads to (c, a')
378+ 1. Given σ(a', z'), compute current consumption c that
379+ satisfies Euler equation
380+ 2. Compute the endogenous current asset level a that leads
381+ to (c, a')
379382 3. Interpolate back to exogenous grid to get σ_new(a, z)
380383 """
381384 R, β, γ, Π, z_grid, asset_grid = ifp
@@ -386,7 +389,8 @@ def K(σ: jnp.ndarray, ifp: IFP) -> jnp.ndarray:
386389 """
387390 Compute updated consumption policy for income state z_j.
388391
389- The asset_grid here represents a' (next period assets), not current assets.
392+ The asset_grid here represents a' (next period assets),
393+ not current assets.
390394 """
391395
392396 # Step 1: Compute expected marginal utility of consumption tomorrow
@@ -396,12 +400,14 @@ def K(σ: jnp.ndarray, ifp: IFP) -> jnp.ndarray:
396400 # where the expectation is over tomorrow's income state z'
397401 # conditional on today's income state z_j
398402
399- u_prime_vals = u_prime(σ, γ) # u'(σ(a', z')) for all (a', z')
400- # Shape: (n_a, n_z) where n_a is # of a' values
403+ # u'(σ(a', z')) for all (a', z')
404+ # Shape: (n_a, n_z) where n_a is # of a' values
405+ u_prime_vals = u_prime(σ, γ)
401406
402- expected_marginal = u_prime_vals @ Π[j, :] # Matrix multiply to get expectation
403- # Π[j, :] are transition probs from z_j
404- # Result shape: (n_a,) - one value per a'
407+ # Matrix multiply to get expectation
408+ # Π[j, :] are transition probs from z_j
409+ # Result shape: (n_a,) - one value per a'
410+ expected_marginal = u_prime_vals @ Π[j, :]
405411
406412 # Step 2: Use Euler equation to find today's consumption
407413 # -------------------------------------------------------
@@ -418,13 +424,14 @@ def K(σ: jnp.ndarray, ifp: IFP) -> jnp.ndarray:
418424 # --------------------------------------------------
419425 # The budget constraint is: a_{t+1} + c_t = R * a_t + Y_t
420426 # Rearranging: a_t = (a_{t+1} + c_t - Y_t) / R
421- # For each (a'_i, c_i) pair, find the current asset level a^e_i that
422- # makes this budget constraint hold
427+ # For each (a'_i, c_i) pair, find the current asset
428+ # level a^e_i that makes this budget constraint hold
423429
430+ # asset_grid[i] is a'_i, c_vals[i] is c_i
431+ # y(z_grid[j]) is income today
432+ # a_endogenous[i] is the current asset level that
433+ # leads to this (c_i, a'_i) pair. Shape: (n_a,)
424434 a_endogenous = (1/R) * (asset_grid + c_vals - y(z_grid[j]))
425- # asset_grid[i] is a'_i, c_vals[i] is c_i, y(z_grid[j]) is income today
426- # a_endogenous[i] is the current asset level that leads to this (c_i, a'_i) pair
427- # Shape: (n_a,)
428435
429436 # Step 4: Interpolate back to exogenous grid
430437 # -------------------------------------------
@@ -547,7 +554,9 @@ R, β, γ, Π, z_grid, asset_grid = ifp_cake_eating
547554
548555fig, ax = plt.subplots()
549556ax.plot(asset_grid, σ_star[:, 0], label='numerical')
550- ax.plot(asset_grid, c_star(asset_grid, ifp_cake_eating.β, ifp_cake_eating.γ), '--', label='analytical')
557+ ax.plot(asset_grid,
558+ c_star(asset_grid, ifp_cake_eating.β, ifp_cake_eating.γ),
559+ '--', label='analytical')
551560ax.set(xlabel='assets', ylabel='consumption')
552561ax.legend()
553562plt.show()
@@ -581,13 +590,15 @@ Your figure should show that higher interest rates boost savings and suppress co
581590Here's one solution:
582591
583592``` {code-cell} python3
584- r_vals = np.linspace(0, 0.04, 2.0 )
593+ r_vals = np.linspace(0, 0.04, 4 )
585594
586595fig, ax = plt.subplots()
587596for r_val in r_vals:
588- ifp = IFP(r=r_val)
597+ ifp = create_ifp(r=r_val)
598+ R, β, γ, Π, z_grid, asset_grid = ifp
599+ σ_init = R * asset_grid[:, None] + y(z_grid)
589600 σ_star = solve_model(ifp, σ_init)
590- ax.plot(ifp. asset_grid, σ_star[:, 0], label=f'$r = {r_val:.3f}$')
601+ ax.plot(asset_grid, σ_star[:, 0], label=f'$r = {r_val:.3f}$')
591602
592603ax.set(xlabel='asset level', ylabel='consumption (low income)')
593604ax.legend()
@@ -608,11 +619,11 @@ default parameters.
608619The following figure is a 45 degree diagram showing the law of motion for assets when consumption is optimal
609620
610621``` {code-cell} python3
611- ifp = IFP()
612-
622+ ifp = create_ifp()
623+ R, β, γ, Π, z_grid, asset_grid = ifp
624+ σ_init = R * asset_grid[:, None] + y(z_grid)
613625σ_star = solve_model(ifp, σ_init)
614- a = ifp.asset_grid
615- R = ifp.R
626+ a = asset_grid
616627
617628fig, ax = plt.subplots()
618629for z, lb in zip((0, 1), ('low income', 'high income')):
@@ -663,36 +674,39 @@ Your task is to generate such a histogram.
663674First we write a function to compute a long asset series.
664675
665676``` {code-cell} python3
666- def compute_asset_series(ifp, T=500_000, seed=1234):
677+ def compute_asset_series(ifp, σ_init, T=500_000, seed=1234):
667678 """
668679 Simulates a time series of length T for assets, given optimal
669680 savings behavior.
670681
671682 ifp is an instance of IFP
672683 """
673- P, y, R = ifp.P, ifp.y, ifp.R # Simplify names
684+ R, β, γ, Π, z_grid, asset_grid = ifp
674685
675686 # Solve for the optimal policy
676- σ_star = solve_model_time_iter (ifp, σ_init, verbose=False )
677- σ = lambda a, z: np.interp(a, ifp. asset_grid, σ_star[:, z])
687+ σ_star = solve_model (ifp, σ_init)
688+ σ = lambda a, z: np.interp(a, asset_grid, σ_star[:, z])
678689
679690 # Simulate the exogeneous state process
680- mc = MarkovChain(P )
691+ mc = MarkovChain(Π )
681692 z_seq = mc.simulate(T, random_state=seed)
682693
683694 # Simulate the asset path
684695 a = np.zeros(T+1)
685696 for t in range(T):
686- z = z_seq[t]
687- a[t+1] = R * (a[t] - σ(a[t], z)) + y[z]
697+ z_idx = z_seq[t]
698+ z_val = z_grid[z_idx]
699+ a[t+1] = R * a[t] + y(z_val) - σ(a[t], z_idx)
688700 return a
689701```
690702
691703Now we call the function, generate the series and then histogram it:
692704
693705``` {code-cell} python3
694- ifp = IFP()
695- a = compute_asset_series(ifp)
706+ ifp = create_ifp()
707+ R, β, γ, Π, z_grid, asset_grid = ifp
708+ σ_init = R * asset_grid[:, None] + y(z_grid)
709+ a = compute_asset_series(ifp, σ_init)
696710
697711fig, ax = plt.subplots()
698712ax.hist(a, bins=20, alpha=0.5, density=True)
@@ -750,8 +764,10 @@ fig, ax = plt.subplots()
750764asset_mean = []
751765for r in r_vals:
752766 print(f'Solving model at r = {r}')
753- ifp = IFP(r=r)
754- mean = np.mean(compute_asset_series(ifp, T=250_000))
767+ ifp = create_ifp(r=r)
768+ R, β, γ, Π, z_grid, asset_grid = ifp
769+ σ_init = R * asset_grid[:, None] + y(z_grid)
770+ mean = np.mean(compute_asset_series(ifp, σ_init, T=250_000))
755771 asset_mean.append(mean)
756772ax.plot(asset_mean, r_vals)
757773
0 commit comments