@@ -26,7 +26,7 @@ kernelspec:
2626
2727In addition to what's in Anaconda, this lecture will need the following libraries:
2828
29- ``` {code-cell} python
29+ ``` {code-cell} python
3030:tags: ["hide-output"]
3131!pip install --upgrade quantecon
3232```
@@ -74,7 +74,7 @@ Such dynamics are consistent with experiences of many countries.
7474
7575Let's start with some imports:
7676
77- ``` {code-cell} python
77+ ``` {code-cell} python
7878import matplotlib.pyplot as plt
7979import numpy as np
8080import quantecon as qe
@@ -346,21 +346,21 @@ As we have in other places, we accelerate our code using Numba.
346346We define a class that will store parameters, grids and transition
347347probabilities.
348348
349- ``` {code-cell} python
349+ ``` {code-cell} python
350350class Arellano_Economy:
351351 " Stores data and creates primitives for the Arellano economy. "
352352
353- def __init__(self,
353+ def __init__(self,
354354 B_grid_size= 251, # Grid size for bonds
355- B_grid_min=-0.45, # Smallest B value
355+ B_grid_min=-0.45, # Smallest B value
356356 B_grid_max=0.45, # Largest B value
357- y_grid_size=51, # Grid size for income
357+ y_grid_size=51, # Grid size for income
358358 β=0.953, # Time discount parameter
359359 γ=2.0, # Utility parameter
360360 r=0.017, # Lending rate
361361 ρ=0.945, # Persistence in the income process
362362 η=0.025, # Standard deviation of the income process
363- θ=0.282, # Prob of re-entering financial markets
363+ θ=0.282, # Prob of re-entering financial markets
364364 def_y_param=0.969): # Parameter governing income in default
365365
366366 # Save parameters
@@ -370,7 +370,7 @@ class Arellano_Economy:
370370 self.y_grid_size = y_grid_size
371371 self.B_grid_size = B_grid_size
372372 self.B_grid = np.linspace(B_grid_min, B_grid_max, B_grid_size)
373- mc = qe.markov.tauchen(ρ, η, 0, 3, y_grid_size )
373+ mc = qe.markov.tauchen(y_grid_size, ρ, η, 0, 3)
374374 self.y_grid, self.P = np.exp(mc.state_values), mc.P
375375
376376 # The index at which B_grid is (close to) zero
@@ -380,7 +380,7 @@ class Arellano_Economy:
380380 self.def_y = np.minimum(def_y_param * np.mean(self.y_grid), self.y_grid)
381381
382382 def params(self):
383- return self.β, self.γ, self.r, self.ρ, self.η, self.θ
383+ return self.β, self.γ, self.r, self.ρ, self.η, self.θ
384384
385385 def arrays(self):
386386 return self.P, self.y_grid, self.B_grid, self.def_y, self.B0_idx
@@ -389,14 +389,14 @@ class Arellano_Economy:
389389Notice how the class returns the data it stores as simple numerical values and
390390arrays via the methods ` params ` and ` arrays ` .
391391
392- We will use this data in the Numba-jitted functions defined below.
392+ We will use this data in the Numba-jitted functions defined below.
393393
394394Jitted functions prefer simple arguments, since type inference is easier.
395395
396396Here is the utility function.
397397
398398
399- ``` {code-cell} python
399+ ``` {code-cell} python
400400@njit
401401def u(c, γ):
402402 return c**(1-γ)/(1-γ)
@@ -405,7 +405,7 @@ def u(c, γ):
405405Here is a function to compute the bond price at each state, given $v_c$ and
406406$v_d$.
407407
408- ``` {code-cell} python
408+ ``` {code-cell} python
409409@njit
410410def compute_q(v_c, v_d, q, params, arrays):
411411 """
@@ -414,29 +414,29 @@ def compute_q(v_c, v_d, q, params, arrays):
414414 This function writes to the array q that is passed in as an argument.
415415 """
416416
417- # Unpack
418- β, γ, r, ρ, η, θ = params
417+ # Unpack
418+ β, γ, r, ρ, η, θ = params
419419 P, y_grid, B_grid, def_y, B0_idx = arrays
420420
421421 for B_idx in range(len(B_grid)):
422422 for y_idx in range(len(y_grid)):
423423 # Compute default probability and corresponding bond price
424424 default_states = 1.0 * (v_c[B_idx, :] < v_d)
425- delta = np.dot(default_states, P[y_idx, :])
425+ delta = np.dot(default_states, P[y_idx, :])
426426 q[B_idx, y_idx] = (1 - delta ) / (1 + r)
427427```
428428
429429Next we introduce Bellman operators that updated $v_d$ and $v_c$.
430430
431- ``` {code-cell} python
431+ ``` {code-cell} python
432432@njit
433433def T_d(y_idx, v_c, v_d, params, arrays):
434434 """
435435 The RHS of the Bellman equation when income is at index y_idx and
436436 the country has chosen to default. Returns an update of v_d.
437437 """
438- # Unpack
439- β, γ, r, ρ, η, θ = params
438+ # Unpack
439+ β, γ, r, ρ, η, θ = params
440440 P, y_grid, B_grid, def_y, B0_idx = arrays
441441
442442 current_utility = u(def_y[y_idx], γ)
@@ -453,8 +453,8 @@ def T_c(B_idx, y_idx, v_c, v_d, q, params, arrays):
453453 defaulted state on their debt. Returns a value that corresponds to
454454 v_c[B_idx, y_idx], as well as the optimal level of bond sales B'.
455455 """
456- # Unpack
457- β, γ, r, ρ, η, θ = params
456+ # Unpack
457+ β, γ, r, ρ, η, θ = params
458458 P, y_grid, B_grid, def_y, B0_idx = arrays
459459 B = B_grid[B_idx]
460460 y = y_grid[y_idx]
@@ -475,14 +475,14 @@ def T_c(B_idx, y_idx, v_c, v_d, q, params, arrays):
475475
476476Here is a fast function that calls these operators in the right sequence.
477477
478- ``` {code-cell} python
478+ ``` {code-cell} python
479479@njit(parallel=True)
480480def update_values_and_prices(v_c, v_d, # Current guess of value functions
481481 B_star, q, # Arrays to be written to
482482 params, arrays):
483483
484- # Unpack
485- β, γ, r, ρ, η, θ = params
484+ # Unpack
485+ β, γ, r, ρ, η, θ = params
486486 P, y_grid, B_grid, def_y, B0_idx = arrays
487487 y_grid_size = len(y_grid)
488488 B_grid_size = len(B_grid)
@@ -498,7 +498,7 @@ def update_values_and_prices(v_c, v_d, # Current guess of value functions
498498 for y_idx in prange(y_grid_size):
499499 new_v_d[y_idx] = T_d(y_idx, v_c, v_d, params, arrays)
500500 for B_idx in range(B_grid_size):
501- new_v_c[B_idx, y_idx], Bp_idx = T_c(B_idx, y_idx,
501+ new_v_c[B_idx, y_idx], Bp_idx = T_c(B_idx, y_idx,
502502 v_c, v_d, q, params, arrays)
503503 B_star[B_idx, y_idx] = Bp_idx
504504
@@ -515,7 +515,7 @@ In fact, one of the jobs of this function is to take an instance of
515515` Arellano_Economy ` , which is hard for the JIT compiler to handle, and strip it
516516down to more basic objects, which are then passed out to jitted functions.
517517
518- ``` {code-cell} python
518+ ``` {code-cell} python
519519def solve(model, tol=1e-8, max_iter=10_000):
520520 """
521521 Given an instance of Arellano_Economy, this function computes the optimal
@@ -524,7 +524,7 @@ def solve(model, tol=1e-8, max_iter=10_000):
524524 # Unpack
525525 params = model.params()
526526 arrays = model.arrays()
527- y_grid_size, B_grid_size = model.y_grid_size, model.B_grid_size
527+ y_grid_size, B_grid_size = model.y_grid_size, model.B_grid_size
528528
529529 # Initial conditions for v_c and v_d
530530 v_c = np.zeros((B_grid_size, y_grid_size))
@@ -549,13 +549,13 @@ def solve(model, tol=1e-8, max_iter=10_000):
549549 current_iter += 1
550550
551551 print(f"Terminating at iteration {current_iter}.")
552- return v_c, v_d, q, B_star
552+ return v_c, v_d, q, B_star
553553```
554554
555555Finally, we write a function that will allow us to simulate the economy once
556556we have the policy functions
557557
558- ``` {code-cell} python
558+ ``` {code-cell} python
559559def simulate(model, T, v_c, v_d, q, B_star, y_idx=None, B_idx=None):
560560 """
561561 Simulates the Arellano 2008 model of sovereign debt
@@ -707,17 +707,17 @@ To the extent that you can, replicate the figures shown above
707707
708708Compute the value function, policy and equilibrium prices
709709
710- ``` {code-cell} python
710+ ``` {code-cell} python
711711ae = Arellano_Economy()
712712```
713713
714- ``` {code-cell} python
714+ ``` {code-cell} python
715715v_c, v_d, q, B_star = solve(ae)
716716```
717717
718718Compute the bond price schedule as seen in figure 3 of Arellano (2008)
719719
720- ``` {code-cell} python
720+ ``` {code-cell} python
721721# Unpack some useful names
722722B_grid, y_grid, P = ae.B_grid, ae.y_grid, ae.P
723723B_grid_size, y_grid_size = len(B_grid), len(y_grid)
@@ -748,7 +748,7 @@ plt.show()
748748
749749Draw a plot of the value functions
750750
751- ``` {code-cell} python
751+ ``` {code-cell} python
752752v = np.maximum(v_c, np.reshape(v_d, (1, y_grid_size)))
753753
754754fig, ax = plt.subplots(figsize=(10, 6.5))
@@ -763,14 +763,14 @@ plt.show()
763763
764764Draw a heat map for default probability
765765
766- ``` {code-cell} python
766+ ``` {code-cell} python
767767xx, yy = B_grid, y_grid
768768zz = np.empty_like(v_c)
769769
770770for B_idx in range(B_grid_size):
771771 for y_idx in range(y_grid_size):
772772 default_states = 1.0 * (v_c[B_idx, :] < v_d)
773- zz[B_idx, y_idx] = np.dot(default_states, P[y_idx, :])
773+ zz[B_idx, y_idx] = np.dot(default_states, P[y_idx, :])
774774
775775# Create figure
776776fig, ax = plt.subplots(figsize=(10, 6.5))
@@ -784,7 +784,7 @@ plt.show()
784784
785785Plot a time series of major variables simulated from the model
786786
787- ``` {code-cell} python
787+ ``` {code-cell} python
788788T = 250
789789np.random.seed(42)
790790y_sim, y_a_sim, B_sim, q_sim, d_sim = simulate(ae, T, v_c, v_d, q, B_star)
@@ -826,4 +826,4 @@ for ax, series, title in zip(axes, plot_series, titles):
826826plt.show()
827827```
828828``` {solution-end}
829- ```
829+ ```
0 commit comments