@@ -4,7 +4,7 @@ jupytext:
4
4
extension : .md
5
5
format_name : myst
6
6
format_version : 0.13
7
- jupytext_version : 1.14.5
7
+ jupytext_version : 1.16.1
8
8
kernelspec :
9
9
display_name : Python 3 (ipykernel)
10
10
language : python
@@ -31,12 +31,10 @@ you can jump to [](eq:ntecx2).
31
31
32
32
The code outputs below are generated by machine connected to the following GPU
33
33
34
-
35
34
``` {code-cell} ipython3
36
35
!nvidia-smi
37
36
```
38
37
39
-
40
38
In addition to what's in Anaconda, this lecture will need the following libraries:
41
39
42
40
``` {code-cell} ipython3
@@ -410,7 +408,7 @@ def create_model(N=100, # size of state space for Markov chain
410
408
σ=0.01, # persistence parameter for Markov chain
411
409
β=0.98, # discount factor
412
410
γ=2.5, # coefficient of risk aversion
413
- μ_c=0.01, # mean growth of consumtion
411
+ μ_c=0.01, # mean growth of consumption
414
412
μ_d=0.01, # mean growth of dividends
415
413
σ_c=0.02, # consumption volatility
416
414
σ_d=0.04): # dividend volatility
@@ -431,14 +429,14 @@ Here's a function that does this using loops.
431
429
432
430
```{code-cell} ipython3
433
431
def compute_K_loop(model):
434
- # Setp up
432
+ # unpack
435
433
P, S, β, γ, μ_c, μ_d, σ_c, σ_d = model
436
434
N = len(S)
437
435
K = np.empty((N, N))
438
436
a = μ_d - γ * μ_c
439
437
for i, x in enumerate(S):
440
438
for j, y in enumerate(S):
441
- e = jnp .exp(a + (1 - γ) * x + (σ_d**2 + γ**2 * σ_c**2) / 2)
439
+ e = np .exp(a + (1 - γ) * x + (σ_d**2 + γ**2 * σ_c**2) / 2)
442
440
K[i, j] = β * e * P[i, j]
443
441
return K
444
442
```
@@ -447,24 +445,24 @@ To exploit the parallelization capabilities of JAX, let's also write a vectorize
447
445
448
446
```{code-cell} ipython3
449
447
def compute_K(model):
450
- # Setp up
448
+ # unpack
451
449
P, S, β, γ, μ_c, μ_d, σ_c, σ_d = model
452
450
N = len(S)
453
451
# Reshape and multiply pointwise using broadcasting
454
- x = jnp .reshape(S, (N, 1))
452
+ x = np .reshape(S, (N, 1))
455
453
a = μ_d - γ * μ_c
456
- e = jnp .exp(a + (1 - γ) * x + (σ_d**2 + γ**2 * σ_c**2) / 2)
454
+ e = np .exp(a + (1 - γ) * x + (σ_d**2 + γ**2 * σ_c**2) / 2)
457
455
K = β * e * P
458
456
return K
459
457
```
460
458
461
- These two functions produce the saem output:
459
+ These two functions produce the same output:
462
460
463
461
```{code-cell} ipython3
464
462
model = create_model(N=10)
465
463
K1 = compute_K(model)
466
464
K2 = compute_K_loop(model)
467
- jnp .allclose(K1, K2)
465
+ np .allclose(K1, K2)
468
466
```
469
467
470
468
Now we can compute the price-dividend ratio:
@@ -492,9 +490,9 @@ def price_dividend_ratio(model, test_stable=True):
492
490
test_stability(K)
493
491
494
492
# Compute v
495
- I = jnp .identity(N)
496
- ones_vec = jnp .ones(N)
497
- v = jnp .linalg.solve(I - K, K @ ones_vec)
493
+ I = np .identity(N)
494
+ ones_vec = np .ones(N)
495
+ v = np .linalg.solve(I - K, K @ ones_vec)
498
496
499
497
return v
500
498
```
@@ -504,7 +502,7 @@ Here's a plot of $v$ as a function of the state for several values of $\gamma$.
504
502
```{code-cell} ipython3
505
503
model = create_model()
506
504
S = model.S
507
- γs = jnp .linspace(2.0, 3.0, 5)
505
+ γs = np .linspace(2.0, 3.0, 5)
508
506
509
507
fig, ax = plt.subplots()
510
508
@@ -701,7 +699,7 @@ def create_sv_model(β=0.98, # discount factor
701
699
bar_σ=0.01, # volatility scaling parameter
702
700
ρ_z=0.9, # persistence parameter for z
703
701
σ_z=0.01, # persistence parameter for z
704
- μ_c=0.001, # mean growth of consumtion
702
+ μ_c=0.001, # mean growth of consumption
705
703
μ_d=0.005): # mean growth of dividends
706
704
707
705
mc = qe.tauchen(I, ρ_c, σ_c)
@@ -763,7 +761,7 @@ def sv_pd_ratio(sv_model, test_stable=True):
763
761
price-dividend ratio
764
762
765
763
"""
766
- # Setp up
764
+ # unpack
767
765
P, hc_grid, Q, hd_grid, R, z_grid, β, γ, bar_σ, μ_c, μ_d = sv_model
768
766
I, J, K = len(hc_grid), len(hd_grid), len(z_grid)
769
767
N = I * J * K
@@ -862,12 +860,12 @@ def compute_A_jax(sv_model, shapes):
862
860
I, J, K = shapes
863
861
N = I * J * K
864
862
# Reshape and broadcast over (i, j, k, i', j', k')
865
- hc = np .reshape(hc_grid, (I, 1, 1, 1, 1, 1))
866
- hd = np .reshape(hd_grid, (1, J, 1, 1, 1, 1))
867
- z = np .reshape(z_grid, (1, 1, K, 1, 1, 1))
868
- P = np .reshape(P, (I, 1, 1, I, 1, 1))
869
- Q = np .reshape(Q, (1, J, 1, 1, J, 1))
870
- R = np .reshape(R, (1, 1, K, 1, 1, K))
863
+ hc = jnp .reshape(hc_grid, (I, 1, 1, 1, 1, 1))
864
+ hd = jnp .reshape(hd_grid, (1, J, 1, 1, 1, 1))
865
+ z = jnp .reshape(z_grid, (1, 1, K, 1, 1, 1))
866
+ P = jnp .reshape(P, (I, 1, 1, I, 1, 1))
867
+ Q = jnp .reshape(Q, (1, J, 1, 1, J, 1))
868
+ R = jnp .reshape(R, (1, 1, K, 1, 1, K))
871
869
# Compute A and then reshape to create a matrix
872
870
a = μ_d - γ * μ_c
873
871
b = bar_σ**2 * (jnp.exp(2 * hd) + γ**2 * jnp.exp(2 * hc)) / 2
@@ -896,7 +894,7 @@ def sv_pd_ratio_jax(sv_model, shapes):
896
894
price-dividend ratio
897
895
898
896
"""
899
- # Setp up
897
+ # unpack
900
898
P, hc_grid, Q, hd_grid, R, z_grid, β, γ, bar_σ, μ_c, μ_d = sv_model
901
899
I, J, K = len(hc_grid), len(hd_grid), len(z_grid)
902
900
shapes = I, J, K
@@ -971,13 +969,13 @@ def A(g, sv_model, shapes):
971
969
P, hc_grid, Q, hd_grid, R, z_grid, β, γ, bar_σ, μ_c, μ_d = sv_model
972
970
I, J, K = shapes
973
971
# Reshape and broadcast over (i, j, k, i', j', k')
974
- hc = np .reshape(hc_grid, (I, 1, 1, 1, 1, 1))
975
- hd = np .reshape(hd_grid, (1, J, 1, 1, 1, 1))
976
- z = np .reshape(z_grid, (1, 1, K, 1, 1, 1))
977
- P = np .reshape(P, (I, 1, 1, I, 1, 1))
978
- Q = np .reshape(Q, (1, J, 1, 1, J, 1))
979
- R = np .reshape(R, (1, 1, K, 1, 1, K))
980
- g = np .reshape(g, (1, 1, 1, I, J, K))
972
+ hc = jnp .reshape(hc_grid, (I, 1, 1, 1, 1, 1))
973
+ hd = jnp .reshape(hd_grid, (1, J, 1, 1, 1, 1))
974
+ z = jnp .reshape(z_grid, (1, 1, K, 1, 1, 1))
975
+ P = jnp .reshape(P, (I, 1, 1, I, 1, 1))
976
+ Q = jnp .reshape(Q, (1, J, 1, 1, J, 1))
977
+ R = jnp .reshape(R, (1, 1, K, 1, 1, K))
978
+ g = jnp .reshape(g, (1, 1, 1, I, J, K))
981
979
a = μ_d - γ * μ_c
982
980
b = bar_σ**2 * (jnp.exp(2 * hd) + γ**2 * jnp.exp(2 * hc)) / 2
983
981
κ = jnp.exp(a + (1 - γ) * z + b)
@@ -991,12 +989,10 @@ that acts directly on the linear operator `A`.
991
989
992
990
```{code-cell} ipython3
993
991
def sv_pd_ratio_linop(sv_model, shapes):
994
-
995
- # Setp up
996
992
P, hc_grid, Q, hd_grid, R, z_grid, β, γ, bar_σ, μ_c, μ_d = sv_model
997
993
I, J, K = shapes
998
994
999
- ones_array = np .ones((I, J, K))
995
+ ones_array = jnp .ones((I, J, K))
1000
996
# Set up the operator g -> (I - A) g
1001
997
J = lambda g: g - A(g, sv_model, shapes)
1002
998
# Solve v = (I - A)^{-1} A 1
0 commit comments