From c227f9c60304b445bd332eac3abed27ec55856e7 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Mon, 20 Oct 2025 17:27:54 +0900 Subject: [PATCH 01/16] update markov_asset --- lectures/markov_asset.md | 166 +++++++++++++++++++++------------------ 1 file changed, 91 insertions(+), 75 deletions(-) diff --git a/lectures/markov_asset.md b/lectures/markov_asset.md index 1739dde0f..0c4264df0 100644 --- a/lectures/markov_asset.md +++ b/lectures/markov_asset.md @@ -75,7 +75,10 @@ Let's start with some imports: import matplotlib.pyplot as plt import numpy as np import quantecon as qe -from numpy.linalg import eigvals, solve +import jax +import jax.numpy as jnp +from jax.numpy.linalg import eigvals, solve +from typing import NamedTuple ``` ## {index}`Pricing Models ` @@ -91,7 +94,7 @@ Let $\{d_t\}_{t \geq 0}$ be a stream of dividends Let's look at some equations that we expect to hold for prices of assets under ex-dividend contracts (we will consider cum-dividend pricing in the exercises). -### Risk-Neutral Pricing +### Risk-neutral pricing ```{index} single: Pricing Models; Risk-Neutral ``` @@ -116,7 +119,7 @@ Here ${\mathbb E}_t [y]$ denotes the best forecast of $y$, conditioned on inform More precisely, ${\mathbb E}_t [y]$ is the mathematical expectation of $y$ conditional on information available at time $t$. -### Pricing with Random Discount Factor +### Pricing with random discount factor ```{index} single: Pricing Models; Risk Aversion ``` @@ -145,7 +148,7 @@ This is because such assets pay well when funds are more urgently wanted. We give examples of how the stochastic discount factor has been modeled below. -### Asset Pricing and Covariances +### Asset pricing and covariances Recall that, from the definition of a conditional covariance ${\rm cov}_t (x_{t+1}, y_{t+1})$, we have @@ -174,7 +177,7 @@ Equation {eq}`lteeqs102` asserts that the covariance of the stochastic discount We give examples of some models of stochastic discount factors that have been proposed later in this lecture and also in a [later lecture](https://python-advanced.quantecon.org/lucas_model.html). -### The Price-Dividend Ratio +### The price-dividend ratio Aside from prices, another quantity of interest is the **price-dividend ratio** $v_t := p_t / d_t$. @@ -190,7 +193,7 @@ v_t = {\mathbb E}_t \left[ m_{t+1} \frac{d_{t+1}}{d_t} (1 + v_{t+1}) \right] Below we'll discuss the implication of this equation. -## Prices in the Risk-Neutral Case +## Prices in the risk-neutral case What can we say about price dynamics on the basis of the models described above? @@ -203,7 +206,7 @@ For now we'll study the risk-neutral case in which the stochastic discount fac We'll focus on how an asset price depends on a dividend process. -### Example 1: Constant Dividends +### Example 1: constant dividends The simplest case is risk-neutral price of a constant, non-random dividend stream $d_t = d > 0$. @@ -234,7 +237,7 @@ This is the equilibrium price in the constant dividend case. Indeed, simple algebra shows that setting $p_t = \bar p$ for all $t$ satisfies the difference equation $p_t = \beta (d + p_{t+1})$. -### Example 2: Dividends with Deterministic Growth Paths +### Example 2: dividends with deterministic growth paths Consider a growing, non-random dividend process $d_{t+1} = g d_t$ where $0 < g \beta < 1$. @@ -267,7 +270,7 @@ $$ This is called the *Gordon formula*. (mass_mg)= -### Example 3: Markov Growth, Risk-Neutral Pricing +### Example 3: Markov growth, risk-neutral pricing Next, we consider a dividend process @@ -310,14 +313,14 @@ The next figure shows a simulation, where * $\{X_t\}$ evolves as a discretized AR1 process produced using {ref}`Tauchen's method `. * $g_t = \exp(X_t)$, so that $\ln g_t = X_t$ is the growth rate. -```{code-cell} ipython +```{code-cell} ipython3 n = 7 mc = qe.tauchen(n, 0.96, 0.25) sim_length = 80 -x_series = mc.simulate(sim_length, init=np.median(mc.state_values)) -g_series = np.exp(x_series) -d_series = np.cumprod(g_series) # Assumes d_0 = 1 +x_series = mc.simulate(sim_length, init=jnp.median(mc.state_values)) +g_series = jnp.exp(x_series) +d_series = jnp.cumprod(g_series) # Assumes d_0 = 1 series = [x_series, g_series, d_series, np.log(d_series)] labels = ['$X_t$', '$g_t$', '$d_t$', r'$\log \, d_t$'] @@ -330,7 +333,7 @@ plt.tight_layout() plt.show() ``` -#### Pricing Formula +#### Pricing formula To obtain asset prices in this setting, let's adapt our analysis from the case of deterministic growth. @@ -400,18 +403,18 @@ As before, we'll generate $\{X_t\}$ as a {ref}`discretized AR1 process Here's the code, including a test of the spectral radius condition -```{code-cell} python3 +```{code-cell} ipython3 n = 25 # Size of state space β = 0.9 mc = qe.tauchen(n, 0.96, 0.02) -K = mc.P * np.exp(mc.state_values) +K = mc.P * jnp.exp(mc.state_values) warning_message = "Spectral radius condition fails" -assert np.max(np.abs(eigvals(K))) < 1 / β, warning_message +assert jnp.max(jnp.abs(eigvals(K))) < 1 / β, warning_message -I = np.identity(n) -v = solve(I - β * K, β * K @ np.ones(n)) +I = jnp.identity(n) +v = solve(I - β * K, β * K @ jnp.ones(n)) fig, ax = plt.subplots() ax.plot(mc.state_values, v, 'g-o', lw=2, alpha=0.7, label='$v$') @@ -440,7 +443,7 @@ We'll price several distinct assets, including * A consol (a type of bond issued by the UK government in the 19th century) * Call options on a consol -### Pricing a Lucas Tree +### Pricing a Lucas tree ```{index} single: Finite Markov Asset Pricing; Lucas Tree ``` @@ -539,46 +542,51 @@ v = (I - \beta J)^{-1} \beta J {\mathbb 1} We will define a function tree_price to compute $v$ given parameters stored in the class AssetPriceModel -```{code-cell} python3 -class AssetPriceModel: +```{code-cell} ipython3 +class AssetPriceModel(NamedTuple): """ A class that stores the primitives of the asset pricing model. Parameters ---------- - β : scalar, float - Discount factor mc : MarkovChain - Contains the transition matrix and set of state values for the state - process - γ : scalar(float) - Coefficient of risk aversion + Contains the transition matrix and set of state values g : callable The function mapping states to growth rates - + β : float + Discount factor + γ : float + Coefficient of risk aversion + n: int + The number of states + """ + mc: qe.MarkovChain + g: callable + β: float + γ: float + n: int + + +def create_ap_model(mc=None, g=jnp.exp, β=0.96, γ=2.0): + """Create an AssetPriceModel class""" + if mc is None: + n, ρ, σ = 25, 0.9, 0.02 + mc = qe.tauchen(n, ρ, σ) + else: + mc = mc + n = mc.P.shape[0] + + return AssetPriceModel(mc=mc, g=g, β=β, γ=γ, n=n) + + +def test_stability(Q, β): """ - def __init__(self, β=0.96, mc=None, γ=2.0, g=np.exp): - self.β, self.γ = β, γ - self.g = g - - # A default process for the Markov chain - if mc is None: - self.ρ = 0.9 - self.σ = 0.02 - self.mc = qe.tauchen(n, self.ρ, self.σ) - else: - self.mc = mc - - self.n = self.mc.P.shape[0] - - def test_stability(self, Q): - """ - Stability test for a given matrix Q. - """ - sr = np.max(np.abs(eigvals(Q))) - if not sr < 1 / self.β: - msg = f"Spectral radius condition failed with radius = {sr}" - raise ValueError(msg) + Stability test for a given matrix Q. + """ + sr = np.max(np.abs(eigvals(Q))) + if not sr < 1 / β: + msg = f"Spectral radius condition failed with radius = {sr}" + raise ValueError(msg) def tree_price(ap): @@ -601,11 +609,11 @@ def tree_price(ap): J = P * ap.g(y)**(1 - γ) # Make sure that a unique solution exists - ap.test_stability(J) + test_stability(J, β) # Compute v - I = np.identity(ap.n) - Ones = np.ones(ap.n) + I = jnp.identity(ap.n) + Ones = jnp.ones(ap.n) v = solve(I - β * J, β * J @ Ones) return v @@ -614,16 +622,16 @@ def tree_price(ap): Here's a plot of $v$ as a function of the state for several values of $\gamma$, with a positively correlated Markov process and $g(x) = \exp(x)$ -```{code-cell} python3 +```{code-cell} ipython3 γs = [1.2, 1.4, 1.6, 1.8, 2.0] -ap = AssetPriceModel() +ap = create_ap_model() states = ap.mc.state_values fig, ax = plt.subplots() for γ in γs: - ap.γ = γ - v = tree_price(ap) + tem_ap = create_ap_model(mc=ap.mc, g=ap.g, β=ap.β, γ=γ) + v = tree_price(tem_ap) ax.plot(states, v, lw=2, alpha=0.6, label=rf"$\gamma = {γ}$") ax.set_title('Price-dividend ratio as a function of the state') @@ -706,7 +714,7 @@ p = (I - \beta M)^{-1} \beta M \zeta {\mathbb 1} The above is implemented in the function consol_price. -```{code-cell} python3 +```{code-cell} ipython3 def consol_price(ap, ζ): """ Computes price of a consol bond with payoff ζ @@ -715,7 +723,7 @@ def consol_price(ap, ζ): ---------- ap: AssetPriceModel An instance of AssetPriceModel containing primitives - + ζ : scalar(float) Coupon of the console @@ -723,18 +731,17 @@ def consol_price(ap, ζ): ------- p : array_like(float) Console bond prices - """ # Simplify names, set up matrices β, γ, P, y = ap.β, ap.γ, ap.mc.P, ap.mc.state_values M = P * ap.g(y)**(- γ) # Make sure that a unique solution exists - ap.test_stability(M) + test_stability(M, β) # Compute price - I = np.identity(ap.n) - Ones = np.ones(ap.n) + I = jnp.identity(ap.n) + Ones = jnp.ones(ap.n) p = solve(I - β * M, β * ζ * M @ Ones) return p @@ -812,7 +819,7 @@ Start at some initial $w$ and iterate with $T$ to convergence . We can find the solution with the following function call_option -```{code-cell} python3 +```{code-cell} ipython3 def call_option(ap, ζ, p_s, ϵ=1e-7): """ Computes price of a call option on a consol bond. @@ -828,7 +835,7 @@ def call_option(ap, ζ, p_s, ϵ=1e-7): p_s : scalar(float) Strike price - ϵ : scalar(float), optional(default=1e-8) + ϵ : scalar(float), optional(default=1e-7) Tolerance for infinite horizon problem Returns @@ -842,26 +849,35 @@ def call_option(ap, ζ, p_s, ϵ=1e-7): M = P * ap.g(y)**(- γ) # Make sure that a unique consol price exists - ap.test_stability(M) + test_stability(M, β) # Compute option price p = consol_price(ap, ζ) - w = np.zeros(ap.n) + w = jnp.zeros(ap.n) error = ϵ + 1 - while error > ϵ: + + def step(state): + w, error = state # Maximize across columns - w_new = np.maximum(β * M @ w, p - p_s) + w_new = jnp.maximum(β * M @ w, p - p_s) # Find maximal difference of each component and update - error = np.amax(np.abs(w - w_new)) - w = w_new + error_new = jnp.amax(jnp.abs(w - w_new)) + return (w_new, error_new) - return w + # Check whether converged + def cond(state): + _, error = state + return error > ϵ + + final_w, _ = jax.lax.while_loop(cond, step, (w, error)) + + return final_w ``` Here's a plot of $w$ compared to the consol price when $P_S = 40$ -```{code-cell} python3 -ap = AssetPriceModel(β=0.9) +```{code-cell} ipython3 +ap = create_ap_model(β=0.9) ζ = 1.0 strike_price = 40 From 0317e4a9f51861b3102dec1a15759f3f46ec04f3 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Wed, 22 Oct 2025 15:36:53 +0900 Subject: [PATCH 02/16] modified: lectures/markov_asset.md --- lectures/markov_asset.md | 53 ++++++++++++++++++++++++++-------------- 1 file changed, 34 insertions(+), 19 deletions(-) diff --git a/lectures/markov_asset.md b/lectures/markov_asset.md index 0c4264df0..298d43433 100644 --- a/lectures/markov_asset.md +++ b/lectures/markov_asset.md @@ -35,6 +35,16 @@ kernelspec: "Asset pricing is all about covariances" -- Lars Peter Hansen ``` +```{admonition} GPU +:class: warning + +This lecture is accelerated via [hardware](status:machine-details) that has access to a GPU and JAX for GPU programming. + +Free GPUs are available on Google Colab. To use this option, please click on the play icon top right, select Colab, and set the runtime environment to include a GPU. + +Alternatively, if you have your own GPU, you can follow the [instructions](https://github.com/google/jax) for installing JAX with GPU support. If you would like to install JAX running on the `cpu` only you can use `pip install jax[cpu]` +``` + In addition to what's in Anaconda, this lecture will need the following libraries: ```{code-cell} ipython @@ -976,10 +986,12 @@ $$ Consider the following primitives -```{code-cell} python3 +```{code-cell} ipython3 n = 5 # Size of State Space -P = np.full((n, n), 0.0125) -P[range(n), range(n)] += 1 - P.sum(1) +P = jnp.full((n, n), 0.0125) +P = P.at[jnp.arange(n), jnp.arange(n)].set( + P[jnp.arange(n), jnp.arange(n)] + 1 - P.sum(1) + ) # State values of the Markov chain s = np.array([0.95, 0.975, 1.0, 1.025, 1.05]) γ = 2.0 @@ -1004,11 +1016,13 @@ Do the same for First, let's enter the parameters: -```{code-cell} python3 +```{code-cell} ipython3 n = 5 -P = np.full((n, n), 0.0125) -P[range(n), range(n)] += 1 - P.sum(1) -s = np.array([0.95, 0.975, 1.0, 1.025, 1.05]) # State values +P = jnp.full((n, n), 0.0125) +P = P.at[jnp.arange(n), jnp.arange(n)].set( + P[jnp.arange(n), jnp.arange(n)] + 1 - P.sum(1) + ) +s = jnp.array([0.95, 0.975, 1.0, 1.025, 1.05]) # State values mc = qe.MarkovChain(P, state_values=s) γ = 2.0 @@ -1020,27 +1034,27 @@ p_s = 150.0 Next, we'll create an instance of `AssetPriceModel` to feed into the functions -```{code-cell} python3 -apm = AssetPriceModel(β=β, mc=mc, γ=γ, g=lambda x: x) +```{code-cell} ipython3 +apm = create_ap_model(mc=mc, g=lambda x: x, β=β, γ=γ) ``` Now we just need to call the relevant functions on the data: -```{code-cell} python3 +```{code-cell} ipython3 tree_price(apm) ``` -```{code-cell} python3 +```{code-cell} ipython3 consol_price(apm, ζ) ``` -```{code-cell} python3 +```{code-cell} ipython3 call_option(apm, ζ, p_s) ``` Let's show the last two functions as a plot -```{code-cell} python3 +```{code-cell} ipython3 fig, ax = plt.subplots() ax.plot(s, consol_price(apm, ζ), label='consol') ax.plot(s, call_option(apm, ζ, p_s), label='call option') @@ -1101,7 +1115,7 @@ Is one higher than the other? Can you give intuition? Here's a suitable function: -```{code-cell} python3 +```{code-cell} ipython3 def finite_horizon_call_option(ap, ζ, p_s, k): """ Computes k period option value. @@ -1111,15 +1125,16 @@ def finite_horizon_call_option(ap, ζ, p_s, k): M = P * ap.g(y)**(- γ) # Make sure that a unique solution exists - ap.test_stability(M) - + test_stability(M, β) # Compute option price p = consol_price(ap, ζ) - w = np.zeros(ap.n) - for i in range(k): + def step(i, w): # Maximize across columns - w = np.maximum(β * M @ w, p - p_s) + w = jnp.maximum(β * M @ w, p - p_s) + return w + + w = jax.lax.fori_loop(0, k, step, jnp.zeros(ap.n)) return w ``` From 17396978b84468c0e83f50365659e6762ebf90ab Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Wed, 22 Oct 2025 17:09:33 +0900 Subject: [PATCH 03/16] add metadata cell for figures --- lectures/markov_asset.md | 45 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/lectures/markov_asset.md b/lectures/markov_asset.md index 298d43433..90e20c315 100644 --- a/lectures/markov_asset.md +++ b/lectures/markov_asset.md @@ -324,6 +324,13 @@ The next figure shows a simulation, where * $g_t = \exp(X_t)$, so that $\ln g_t = X_t$ is the growth rate. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + State, growth, and dividend simulation + name: fig_markov_sim +--- n = 7 mc = qe.tauchen(n, 0.96, 0.25) sim_length = 80 @@ -414,6 +421,13 @@ As before, we'll generate $\{X_t\}$ as a {ref}`discretized AR1 process Here's the code, including a test of the spectral radius condition ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Price-dividend ratio risk-neutral case + name: fig_pdv_neutral +--- n = 25 # Size of state space β = 0.9 mc = qe.tauchen(n, 0.96, 0.02) @@ -633,6 +647,13 @@ Here's a plot of $v$ as a function of the state for several values of $\gamma$, with a positively correlated Markov process and $g(x) = \exp(x)$ ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Lucas tree prices for varying risk aversion + name: fig_lucas_gamma +--- γs = [1.2, 1.4, 1.6, 1.8, 2.0] ap = create_ap_model() states = ap.mc.state_values @@ -644,7 +665,6 @@ for γ in γs: v = tree_price(tem_ap) ax.plot(states, v, lw=2, alpha=0.6, label=rf"$\gamma = {γ}$") -ax.set_title('Price-dividend ratio as a function of the state') ax.set_ylabel("price-dividend ratio") ax.set_xlabel("state") ax.legend(loc='upper right') @@ -887,6 +907,13 @@ def call_option(ap, ζ, p_s, ϵ=1e-7): Here's a plot of $w$ compared to the consol price when $P_S = 40$ ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Consol price and call option value + name: fig_consol_call +--- ap = create_ap_model(β=0.9) ζ = 1.0 strike_price = 40 @@ -1055,6 +1082,13 @@ call_option(apm, ζ, p_s) Let's show the last two functions as a plot ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Consol and call option exercise two comparison + name: fig_ex2_prices +--- fig, ax = plt.subplots() ax.plot(s, consol_price(apm, ζ), label='consol') ax.plot(s, call_option(apm, ζ, p_s), label='call option') @@ -1141,7 +1175,14 @@ def finite_horizon_call_option(ap, ζ, p_s, k): Now let's compute the option values at `k=5` and `k=25` -```{code-cell} python3 +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Finite horizon call option values + name: fig_ex3_finite +--- fig, ax = plt.subplots() for k in [5, 25]: w = finite_horizon_call_option(apm, ζ, p_s, k) From 76c64b2a050b0036945d58d4e20ad5b441f75fba Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Wed, 22 Oct 2025 17:27:57 +0900 Subject: [PATCH 04/16] delete metadata cell in solutions --- lectures/markov_asset.md | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/lectures/markov_asset.md b/lectures/markov_asset.md index 90e20c315..a7e04a010 100644 --- a/lectures/markov_asset.md +++ b/lectures/markov_asset.md @@ -1082,13 +1082,6 @@ call_option(apm, ζ, p_s) Let's show the last two functions as a plot ```{code-cell} ipython3 ---- -mystnb: - figure: - caption: | - Consol and call option exercise two comparison - name: fig_ex2_prices ---- fig, ax = plt.subplots() ax.plot(s, consol_price(apm, ζ), label='consol') ax.plot(s, call_option(apm, ζ, p_s), label='call option') @@ -1176,13 +1169,6 @@ def finite_horizon_call_option(ap, ζ, p_s, k): Now let's compute the option values at `k=5` and `k=25` ```{code-cell} ipython3 ---- -mystnb: - figure: - caption: | - Finite horizon call option values - name: fig_ex3_finite ---- fig, ax = plt.subplots() for k in [5, 25]: w = finite_horizon_call_option(apm, ζ, p_s, k) From 02dfe6e50889b1fc45843b2e3c6fa91dc1e6ebac Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Mon, 27 Oct 2025 13:31:16 +0900 Subject: [PATCH 05/16] modify document link --- lectures/markov_asset.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lectures/markov_asset.md b/lectures/markov_asset.md index a7e04a010..ef6c8d3c2 100644 --- a/lectures/markov_asset.md +++ b/lectures/markov_asset.md @@ -185,7 +185,7 @@ It is useful to regard equation {eq}`lteeqs102` as a generalization of equatio Equation {eq}`lteeqs102` asserts that the covariance of the stochastic discount factor with the one period payout $d_{t+1} + p_{t+1}$ is an important determinant of the price $p_t$. -We give examples of some models of stochastic discount factors that have been proposed later in this lecture and also in a [later lecture](https://python-advanced.quantecon.org/lucas_model.html). +We give examples of some models of stochastic discount factors that have been proposed later in this lecture and also in a {doc}`later lecture`. ### The price-dividend ratio @@ -496,7 +496,7 @@ m_{t+1} = \beta \frac{u'(c_{t+1})}{u'(c_t)} where $u$ is a concave utility function and $c_t$ is time $t$ consumption of a representative consumer. -(A derivation of this expression is given in a [later lecture](https://python-advanced.quantecon.org/lucas_model.html)) +(A derivation of this expression is given in a {doc}`later lecture`) Assume the existence of an endowment that follows growth process {eq}`mass_fmce`. From bd4c2f1d128f0e7fbd1726516d7ad1ad88cc8689 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Sat, 15 Nov 2025 17:08:20 +0800 Subject: [PATCH 06/16] modified: lectures/ge_arrow.md --- lectures/ge_arrow.md | 296 +++++++++++++++++++++---------------------- 1 file changed, 148 insertions(+), 148 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index 3886d6ad8..bee76ccb5 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -390,11 +390,11 @@ on $s_t$ via the following string of equalities $$ \begin{aligned} -E \left[ E d(s_{t+j}) | s_{t+1} \right] | s_t +E \left[ E d(s_{t+j} | s_{t+1}) | s_t \right] & = \sum_{s_{t+1}} \left[ \sum_{s_{t+j}} d(s_{t+j}) P_{j-1}(s_{t+j}| s_{t+1} ) \right] P(s_{t+1} | s_t) \\ & = \sum_{s_{t+j}} d(s_{t+j}) \left[ \sum_{s_{t+1}} P_{j-1} ( s_{t+j} |s_{t+1}) P(s_{t+1}| s_t) \right] \\ & = \sum_{s_{t+j}} d(s_{t+j}) P_j (s_{t+j} | s_t ) \\ - & = E d(s_{t+j})| s_t + & = E d(s_{t+j}| s_t) \end{aligned} $$ @@ -416,7 +416,7 @@ $$ The **law of iterated values** states $$ -V \left[ V (d(s_{t+j}) | s_{t+1}) \right] | s_t = V(d(s_{t+j}))| s_t +V \left[ V (d(s_{t+j}) | s_{t+1}) | s_t \right] = V(d(s_{t+j})| s_t ) $$ We verify it by pursuing the following a string of inequalities that are counterparts to those we used @@ -424,11 +424,11 @@ to verify the law of iterated expectations: $$ \begin{aligned} -V \left[ V ( d(s_{t+j}) | s_{t+1} ) \right] | s_t +V \left[ V ( d(s_{t+j}) | s_{t+1} ) | s_t \right] & = \sum_{s_{t+1}} \left[ \sum_{s_{t+j}} d(s_{t+j}) Q_{j-1}(s_{t+j}| s_{t+1} ) \right] Q(s_{t+1} | s_t) \\ & = \sum_{s_{t+j}} d(s_{t+j}) \left[ \sum_{s_{t+1}} Q_{j-1} ( s_{t+j} |s_{t+1}) Q(s_{t+1}| s_t) \right] \\ & = \sum_{s_{t+j}} d(s_{t+j}) Q_j (s_{t+j} | s_t ) \\ - & = E V(d(s_{t+j}))| s_t + & = E V(d(s_{t+j})| s_t) \end{aligned} $$ @@ -544,7 +544,7 @@ A^{K}\left(s\right) \end{array}\right], \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] $$ -and an $n \times 1$ vector of continuation endowment values for each individual $k$ as +and an $n \times 1$ vector-form function of continuation endowment values for each individual $k$ as $$ A^{k}=\left[\begin{array}{c} @@ -557,7 +557,7 @@ $$ $A^k$ of consumer $k$ satisfies $$ -A^k = \left[I - Q\right]^{-1} \left[ y^k\right] +A^k = \left[I - Q\right]^{-1} y^k $$ where @@ -571,7 +571,7 @@ y^{k}\left(\bar{s}_{n}\right) $$ -In a competitive equilibrium of an **infinite horizon** economy with sequential trading of one-period Arrow securities, $A^k(s)$ serves as a state-by-state vector of **debt limits** on the quantities of one-period Arrow securities +In a competitive equilibrium of an *infinite horizon* economy with sequential trading of one-period Arrow securities, $A^k(s)$ serves as a state-by-state vector of *debt limits* on the quantities of one-period Arrow securities paying off in state $s$ at time $t+1$ that individual $k$ can issue at time $t$. @@ -580,8 +580,10 @@ These are often called **natural debt limits**. Evidently, they equal the maximum amount that it is feasible for individual $k$ to repay even if he consumes zero goods forevermore. -**Remark:** If we have an Inada condition at zero consumption or just impose that consumption -be nonnegative, then in a **finite horizon** economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. See the section on a Finite Horizon Economy below. +```{prf:remark} +If we have an Inada condition at zero consumption or just impose that consumption +be nonnegative, then in a *finite horizon* economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. See the section on a [finite orizon economy](#finite-horizon) below. +``` @@ -601,14 +603,14 @@ $$ \end{array}\right], \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] $$ -and an $n \times 1$ vector of continuation wealths for each individual $k$ as +and an $n \times 1$ vector-form function of continuation wealths for each individual $k\in \left[1, \ldots, K\right]$ as $$ \psi^{k}=\left[\begin{array}{c} \psi^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ \psi^{k}\left(\bar{s}_{n}\right) -\end{array}\right], \quad k \in \left[1, \ldots, K\right] +\end{array}\right] $$ Continuation wealth $\psi^k$ of consumer $k$ satisfies @@ -633,12 +635,15 @@ $$ Note that $\sum_{k=1}^K \psi^k = {0}_{n \times 1}$. -**Remark:** At the initial state $s_0 \in \begin{bmatrix} \bar s_1, \ldots, \bar s_n \end{bmatrix}$, +```{prf:remark} +At the initial state $s_0 \in \begin{bmatrix} \bar s_1, \ldots, \bar s_n \end{bmatrix}$, the continuation wealth $\psi^k(s_0) = 0$ for all agents $k = 1, \ldots, K$. This indicates that the economy begins with all agents being debt-free and financial-asset-free at time $0$, state $s_0$. +``` - -**Remark:** Note that all agents' continuation wealths recurrently return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. +```{prf:remark} +Note that all agents' continuation wealths recurrently return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. +``` ### Optimal Portfolios @@ -704,7 +709,7 @@ We now describe a finite-horizon version of the economy that operates for $T+1$ $t \in {\bf T} = \{ 0, 1, \ldots, T\}$. Consequently, we'll want $T+1$ counterparts to objects described above, with one important exception: -we won't need **borrowing limits**. +we won't need *borrowing limits*. * borrowing limits aren't required for a finite horizon economy in which a one-period utility function $u(c)$ satisfies an Inada condition that sets the marginal utility of consumption at zero consumption to zero. @@ -764,12 +769,16 @@ $$ Note that $\sum_{k=1}^K \psi_t^k = {0}_{n \times 1}$ for all $t \in {\bf T}$. -**Remark:** At the initial state $s_0 \in \begin{bmatrix} \bar s_1, \ldots, \bar s_n \end{bmatrix}$, +```{prf:remark} +At the initial state $s_0 \in \begin{bmatrix} \bar s_1, \ldots, \bar s_n \end{bmatrix}$, for all agents $k = 1, \ldots, K$, continuation wealth $\psi_0^k(s_0) = 0$. This indicates that the economy begins with all agents being debt-free and financial-asset-free at time $0$, state $s_0$. +``` -**Remark:** Note that all agents' continuation wealths return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. This will recur if the Markov chain makes the initial state $s_0$ recurrent. +```{prf:remark} +Note that all agents' continuation wealths return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. This will recur if the Markov chain makes the initial state $s_0$ recurrent. +``` @@ -833,7 +842,10 @@ We are ready to dive into some Python code. As usual, we start with Python imports. ```{code-cell} ipython3 +import jax +import jax.numpy as jnp import numpy as np +from typing import NamedTuple import matplotlib.pyplot as plt ``` @@ -849,166 +861,149 @@ In addition to infinite-horizon economies, the code is set up to handle finite-h We'll study examples of finite horizon economies after we first look at some infinite-horizon economies. ```{code-cell} ipython3 -class RecurCompetitive: +class RecurCompetitive(NamedTuple): """ A class that represents a recursive competitive economy with one-period Arrow securities. """ + s: jax.Array # state vector + P: jax.Array # transition matrix + ys: jax.Array # endowments ys = [y1, y2, .., yT] + y: jax.Array # total endownment under each state + n: float # number of states + K: float # number of agents + γ: float # risk aversion + β: float # discount rate + T: float # time horizon, 0 if infinite - def __init__(self, - s, # state vector - P, # transition matrix - ys, # endowments ys = [y1, y2, .., yI] - γ=0.5, # risk aversion - β=0.98, # discount rate - T=None): # time horizon, none if infinite - - # preference parameters - self.γ = γ - self.β = β - # variables dependent on state - self.s = s - self.P = P - self.ys = ys - self.y = np.sum(ys, 1) +def create_rc_model(s, P, ys, γ=0.5, β=0.98, T=0): + n, K = ys.shape + y = jnp.sum(ys, 1) + return RecurCompetitive( + s=s, P=P, ys=ys, y=y n=n, K=K, γ=γ, β=β, T=T + ) - # dimensions - self.n, self.K = ys.shape - # compute pricing kernel - self.Q = self.pricing_kernel() +def u(rc, c): + "The CRRA utility" + return c ** (1 - rc.γ) / (1 - rc.γ) - # compute price of risk-free one-period bond - self.PRF = self.price_risk_free_bond() - # compute risk-free rate - self.R = self.risk_free_rate() +def u_prime(rc, c): + "The first derivative of CRRA utility" + return c ** (-rc.γ) - # V = [I - Q]^{-1} (infinite case) - if T is None: - self.T = None - self.V = np.empty((1, n, n)) - self.V[0] = np.linalg.inv(np.eye(n) - self.Q) - # V = [I + Q + Q^2 + ... + Q^T] (finite case) - else: - self.T = T - self.V = np.empty((T+1, n, n)) - self.V[0] = np.eye(n) - Qt = np.eye(n) - for t in range(1, T+1): - Qt = Qt.dot(self.Q) - self.V[t] = self.V[t-1] + Qt +def pricing_kernel(rc): + "Compute the pricing kernel matrix Q" - # natural debt limit - self.A = self.V[-1] @ ys + c, n, β = rc.y, rc.n, rc.β - def u(self, c): - "The CRRA utility" + Q = jnp.empty((n, n)) + for i in range(n): + for j in range(n): + ratio = u_prime(c[j]) / u_prime(c[i]) + Q[i, j] = β * ratio * P[i, j] - return c ** (1 - self.γ) / (1 - self.γ) + return Q - def u_prime(self, c): - "The first derivative of CRRA utility" - return c ** (-self.γ) +def V(rc, Q): + T = rc.T + n = rc.n + if T==0: + V = jnp.empty((1, n, n)) + V.at[0].set(jnp.linalg.inv(jnp.eye(n) - Q)) + # V = [I + Q + Q^2 + ... + Q^T] (finite case) + else: + V = jnp.empty((T+1, n, n)) + V.at[0].set(jnp.eye(n)) - def pricing_kernel(self): - "Compute the pricing kernel matrix Q" + Qt = jnp.eye(n) + for t in range(1, T+1): + Qt = Qt @ Q + V.at[t].set(V[t-1] + Qt) + return V - c = self.y - n = self.n - Q = np.empty((n, n)) - for i in range(n): - for j in range(n): - ratio = self.u_prime(c[j]) / self.u_prime(c[i]) - Q[i, j] = self.β * ratio * P[i, j] +def A(rc, V): + return V[-1] @ rc.ys - self.Q = Q - return Q +def wealth_distribution(rc, s0_idx): + "Solve for wealth distribution α" - def wealth_distribution(self, s0_idx): - "Solve for wealth distribution α" + # simplify notations + n, T = rc.n, rc.T + Q = pricing_kernel(rc) + y, ys = rc.y, rc.ys - # set initial state - self.s0_idx = s0_idx + V = V(rc, Q) + # row of V corresponding to s0 + Vs0 = V[-1, s0_idx, :] + α = Vs0 @ ys / (Vs0 @ y) - # simplify notations - n = self.n - Q = self.Q - y, ys = self.y, self.ys + return α - # row of V corresponding to s0 - Vs0 = self.V[-1, s0_idx, :] - α = Vs0 @ self.ys / (Vs0 @ self.y) +def continuation_wealths(rc, Q, V, α): + "Given α, compute the continuation wealths ψ" + n, K = rc.n, rc.K + y, ys = rc.y, rc.ys + diff = np.empty((n, K)) + for k in range(K): + diff[:, k] = α[k] * y - ys[:, k] - self.α = α + ψ = V @ diff - return α + return ψ - def continuation_wealths(self): - "Given α, compute the continuation wealths ψ" +def price_risk_free_bond(Q): + "Give Q, compute price of one-period risk free bond" - diff = np.empty((n, K)) - for k in range(K): - diff[:, k] = self.α[k] * self.y - self.ys[:, k] + PRF = jnp.sum(Q, axis=1) - ψ = self.V @ diff - self.ψ = ψ + return PRF - return ψ +def risk_free_rate(Q): + "Given Q, compute one-period gross risk-free interest rate R" - def price_risk_free_bond(self): - "Give Q, compute price of one-period risk free bond" + R = jnp.sum(Q, axis=1) + R = jnp.reciprocal(R) - PRF = np.sum(self.Q, axis=1) - self.PRF = PRF + return R - return PRF +def value_functionss(rc): + "Given α, compute the optimal value functions J in equilibrium" - def risk_free_rate(self): - "Given Q, compute one-period gross risk-free interest rate R" + n, K, T = rc.n, rc.K, rc.T + β = rc.β + P = rc.P - R = np.sum(self.Q, axis=1) - R = np.reciprocal(R) - self.R = R + # compute (I - βP)^(-1) in infinite case + if T==0: + P_seq = jnp.empty((1, n, n)) + P_seq.at[0].set( + jnp.linalg.inv(jnp.eye(n) - β * P) + ) + # and (I + βP + ... + β^T P^T) in finite case + else: + P_seq = jnp.empty((T+1, n, n)) + P_seq.at[0].set(np.eye(n)) - return R + Pt = jnp.eye(n) + for t in range(1, T+1): + Pt = Pt @ P + P_seq.at[t].set(P_seq[t-1] + Pt * β ** t) - def value_functionss(self): - "Given α, compute the optimal value functions J in equilibrium" + # compute the matrix [u(α_1 y), ..., u(α_K, y)] + flow = jnp.empty((n, K)) + for k in range(K): + flow[:, k] = u(α[k] * y) - n, T = self.n, self.T - β = self.β - P = self.P + J = P_seq @ flow - # compute (I - βP)^(-1) in infinite case - if T is None: - P_seq = np.empty((1, n, n)) - P_seq[0] = np.linalg.inv(np.eye(n) - β * P) - # and (I + βP + ... + β^T P^T) in finite case - else: - P_seq = np.empty((T+1, n, n)) - P_seq[0] = np.eye(n) - - Pt = np.eye(n) - for t in range(1, T+1): - Pt = Pt.dot(P) - P_seq[t] = P_seq[t-1] + Pt * β ** t - - # compute the matrix [u(α_1 y), ..., u(α_K, y)] - flow = np.empty((n, K)) - for k in range(K): - flow[:, k] = self.u(self.α[k] * self.y) - - J = P_seq @ flow - - self.J = J - - return J + return J ``` ## Examples @@ -1030,13 +1025,13 @@ Here goes. K, n = 2, 2 # states -s = np.array([0, 1]) +s = jnp.array([0, 1]) # transition -P = np.array([[.5, .5], [.5, .5]]) +P = jnp.array([[.5, .5], [.5, .5]]) # endowments -ys = np.empty((n, K)) +ys = jnp.empty((n, K)) ys[:, 0] = 1 - s # y1 ys[:, 1] = s # y2 ``` @@ -1052,31 +1047,36 @@ ex1.ys ```{code-cell} ipython3 # pricing kernal -ex1.Q +Q = pricing_kernel(ex1) ``` ```{code-cell} ipython3 # Risk free rate R -ex1.R +R = risk_free_rate(Q) ``` ```{code-cell} ipython3 # natural debt limit, A = [A1, A2, ..., AI] -ex1.A +V = V(ex1, Q) +A = A(ex1, V) +A ``` ```{code-cell} ipython3 # when the initial state is state 1 -print(f'α = {ex1.wealth_distribution(s0_idx=0)}') -print(f'ψ = \n{ex1.continuation_wealths()}') +print(f'α = {wealth_distribution(ex1, s0_idx=0)}') +print(f'ψ = \n{continuation_wealths(ex1, )}') print(f'J = \n{ex1.value_functionss()}') ``` ```{code-cell} ipython3 # when the initial state is state 2 -print(f'α = {ex1.wealth_distribution(s0_idx=1)}') -print(f'ψ = \n{ex1.continuation_wealths()}') -print(f'J = \n{ex1.value_functionss()}') +α = wealth_distribution(ex1, s0_idx=1) +ψ = continuation_wealths(ex1, Q, V, α) +J = value_functionss(ex1) +print(f'α = {α}') +print(f'ψ = \n{ψ}') +print(f'J = \n{J}') ``` ### Example 2 From 60a684a14639c8f3d234f94ad10d7114e013c722 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Sun, 16 Nov 2025 14:30:47 +0800 Subject: [PATCH 07/16] modified: lectures/ge_arrow.md --- lectures/ge_arrow.md | 402 +++++++++++++++++++++---------------------- 1 file changed, 200 insertions(+), 202 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index bee76ccb5..4e9024a18 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -875,135 +875,132 @@ class RecurCompetitive(NamedTuple): γ: float # risk aversion β: float # discount rate T: float # time horizon, 0 if infinite + Q: jax.Array # pricing kernek + PRF: jax.Array # price of risk-free bond + R: jax.Array # risk-free rate + A: jax.Array # natural debt limit + α: jax.Array # wealth distribution + ψ: jax.Array # continuation value + J: jax.Array # optimal value -def create_rc_model(s, P, ys, γ=0.5, β=0.98, T=0): +def compute_rc_model(s, P, ys, 0_idx=0, γ=0.5, β=0.98, T=0): n, K = ys.shape y = jnp.sum(ys, 1) - return RecurCompetitive( - s=s, P=P, ys=ys, y=y n=n, K=K, γ=γ, β=β, T=T - ) - - -def u(rc, c): - "The CRRA utility" - return c ** (1 - rc.γ) / (1 - rc.γ) - - -def u_prime(rc, c): - "The first derivative of CRRA utility" - return c ** (-rc.γ) - - -def pricing_kernel(rc): - "Compute the pricing kernel matrix Q" - - c, n, β = rc.y, rc.n, rc.β - - Q = jnp.empty((n, n)) - for i in range(n): - for j in range(n): - ratio = u_prime(c[j]) / u_prime(c[i]) - Q[i, j] = β * ratio * P[i, j] - - return Q - - -def V(rc, Q): - T = rc.T - n = rc.n - if T==0: - V = jnp.empty((1, n, n)) - V.at[0].set(jnp.linalg.inv(jnp.eye(n) - Q)) - # V = [I + Q + Q^2 + ... + Q^T] (finite case) - else: - V = jnp.empty((T+1, n, n)) - V.at[0].set(jnp.eye(n)) - - Qt = jnp.eye(n) - for t in range(1, T+1): - Qt = Qt @ Q - V.at[t].set(V[t-1] + Qt) - return V - - -def A(rc, V): - return V[-1] @ rc.ys - - -def wealth_distribution(rc, s0_idx): - "Solve for wealth distribution α" - - # simplify notations - n, T = rc.n, rc.T - Q = pricing_kernel(rc) - y, ys = rc.y, rc.ys - - V = V(rc, Q) - # row of V corresponding to s0 - Vs0 = V[-1, s0_idx, :] - α = Vs0 @ ys / (Vs0 @ y) - - return α - -def continuation_wealths(rc, Q, V, α): - "Given α, compute the continuation wealths ψ" - n, K = rc.n, rc.K - y, ys = rc.y, rc.ys - diff = np.empty((n, K)) - for k in range(K): - diff[:, k] = α[k] * y - ys[:, k] - - ψ = V @ diff - - return ψ - -def price_risk_free_bond(Q): - "Give Q, compute price of one-period risk free bond" - PRF = jnp.sum(Q, axis=1) + def u(c): + "The CRRA utility" + return c ** (1 - γ) / (1 - γ) + + def u_prime(c): + "The first derivative of CRRA utility" + return c ** (-γ) + + def pricing_kernel(c): + "Compute the pricing kernel matrix Q" + + Q = jnp.empty((n, n)) + for i in range(n): + for j in range(n): + ratio = u_prime(c[j]) / u_prime(c[i]) + Q[i, j] = β * ratio * P[i, j] + + return Q + + def V(Q): + if T==0: + V = jnp.empty((1, n, n)) + V.at[0].set(jnp.linalg.inv(jnp.eye(n) - Q)) + # V = [I + Q + Q^2 + ... + Q^T] (finite case) + else: + V = jnp.empty((T+1, n, n)) + V.at[0].set(jnp.eye(n)) + + Qt = jnp.eye(n) + for t in range(1, T+1): + Qt = Qt @ Q + V.at[t].set(V[t-1] + Qt) + return V + + def A(ys, V): + "Give endownments and pricing kernel, compute natural debt limit" + return V[-1] @ ys + + def wealth_distribution(V, ys, y, 0_idx): + "Solve for wealth distribution α" + + # row of V corresponding to s0 + Vs0 = V[-1, s0_idx, :] + α = Vs0 @ ys / (Vs0 @ y) + + return α + + def continuation_wealths(V, α): + "Given α, compute the continuation wealths ψ" + diff = jnp.empty((n, K)) + for k in range(K): + diff.at[:, k].set(α[k] * y - ys[:, k]) + + ψ = V @ diff + + return ψ + + def price_risk_free_bond(Q): + "Give Q, compute price of one-period risk free bond" + + PRF = jnp.sum(Q, axis=1) + + return PRF + + def risk_free_rate(Q): + "Given Q, compute one-period gross risk-free interest rate R" + + R = jnp.sum(Q, axis=1) + R = jnp.reciprocal(R) + + return R + + def value_functions(α, y): + "Given α, compute the optimal value functions J in equilibrium" + + # compute (I - βP)^(-1) in infinite case + if T==0: + P_seq = jnp.empty((1, n, n)) + P_seq.at[0].set( + jnp.linalg.inv(jnp.eye(n) - β * P) + ) + # and (I + βP + ... + β^T P^T) in finite case + else: + P_seq = jnp.empty((T+1, n, n)) + P_seq.at[0].set(np.eye(n)) + + Pt = jnp.eye(n) + for t in range(1, T+1): + Pt = Pt @ P + P_seq.at[t].set(P_seq[t-1] + Pt * β ** t) + + # compute the matrix [u(α_1 y), ..., u(α_K, y)] + flow = jnp.empty((n, K)) + for k in range(K): + flow.at[:, k].set(u(α[k] * y)) + + J = P_seq @ flow + + return J + + Q = pricing_kernek(y) + V = V(Q) + A = A(ys, Q) + α = wealth_distribution(V, ys, y, s0_idx) + ψ = continuation_wealths(V, α) + PRF = price_risk_free_bond(Q) + R = risk_free_rate(Q) + J = value_functions(α, y) - return PRF - -def risk_free_rate(Q): - "Given Q, compute one-period gross risk-free interest rate R" - - R = jnp.sum(Q, axis=1) - R = jnp.reciprocal(R) - - return R - -def value_functionss(rc): - "Given α, compute the optimal value functions J in equilibrium" - - n, K, T = rc.n, rc.K, rc.T - β = rc.β - P = rc.P - - # compute (I - βP)^(-1) in infinite case - if T==0: - P_seq = jnp.empty((1, n, n)) - P_seq.at[0].set( - jnp.linalg.inv(jnp.eye(n) - β * P) - ) - # and (I + βP + ... + β^T P^T) in finite case - else: - P_seq = jnp.empty((T+1, n, n)) - P_seq.at[0].set(np.eye(n)) - - Pt = jnp.eye(n) - for t in range(1, T+1): - Pt = Pt @ P - P_seq.at[t].set(P_seq[t-1] + Pt * β ** t) - - # compute the matrix [u(α_1 y), ..., u(α_K, y)] - flow = jnp.empty((n, K)) - for k in range(K): - flow[:, k] = u(α[k] * y) - - J = P_seq @ flow - - return J + return RecurCompetitive( + s=s, P=P, ys=ys, y=y n=n, K=K, γ=γ, β=β, T=T, + Q=Q, V=V, A=A, α=α, ψ=ψ, PRF=PRF, R=R, J=J + ) ``` ## Examples @@ -1032,12 +1029,12 @@ P = jnp.array([[.5, .5], [.5, .5]]) # endowments ys = jnp.empty((n, K)) -ys[:, 0] = 1 - s # y1 -ys[:, 1] = s # y2 +ys.at[:, 0].set(1 - s) # y1 +ys.at[:, 1].set(s) # y2 ``` ```{code-cell} ipython3 -ex1 = RecurCompetitive(s, P, ys) +ex1 = compute_rc_model(s, P, ys, s0_idx=0) ``` ```{code-cell} ipython3 @@ -1047,36 +1044,32 @@ ex1.ys ```{code-cell} ipython3 # pricing kernal -Q = pricing_kernel(ex1) +ex1.Q ``` ```{code-cell} ipython3 # Risk free rate R -R = risk_free_rate(Q) +ex1.R ``` ```{code-cell} ipython3 # natural debt limit, A = [A1, A2, ..., AI] -V = V(ex1, Q) -A = A(ex1, V) -A +ex1.A ``` ```{code-cell} ipython3 # when the initial state is state 1 -print(f'α = {wealth_distribution(ex1, s0_idx=0)}') -print(f'ψ = \n{continuation_wealths(ex1, )}') -print(f'J = \n{ex1.value_functionss()}') +print(f'α = {ex1.α}') +print(f'ψ = \n{ex1.ψ}') +print(f'J = \n{ex1.J}') ``` ```{code-cell} ipython3 # when the initial state is state 2 -α = wealth_distribution(ex1, s0_idx=1) -ψ = continuation_wealths(ex1, Q, V, α) -J = value_functionss(ex1) -print(f'α = {α}') -print(f'ψ = \n{ψ}') -print(f'J = \n{J}') +ex1 = compute_rc_model(s, P, ys, s0_idx=1) +print(f'α = {ex1.α}') +print(f'ψ = \n{wx1.ψ}') +print(f'J = \n{ex1.J}') ``` ### Example 2 @@ -1086,19 +1079,19 @@ print(f'J = \n{J}') K, n = 2, 2 # states -s = np.array([1, 2]) +s = jnp.array([1, 2]) # transition -P = np.array([[.5, .5], [.5, .5]]) +P = jnp.array([[.5, .5], [.5, .5]]) # endowments -ys = np.empty((n, K)) -ys[:, 0] = 1.5 # y1 -ys[:, 1] = s # y2 +ys = jnp.empty((n, K)) +ys.at[:, 0].set(1.5) # y1 +ys.at[:, 1].set(s) # y2 ``` ```{code-cell} ipython3 -ex2 = RecurCompetitive(s, P, ys) +ex2 = compute_rc_model(s, P, ys) ``` ```{code-cell} ipython3 @@ -1123,11 +1116,15 @@ Note that the pricing kernal in example economies 1 and 2 differ. This comes from differences in the aggregate endowments in state 1 and 2 in example 1. ```{code-cell} ipython3 -ex2.β * ex2.u_prime(3.5) / ex2.u_prime(2.5) * ex2.P[0,1] +def u_prime(c, γ=0.5): + "The first derivative of CRRA utility" + return c ** (-γ) + +ex2.β * u_prime(3.5) / u_prime(2.5) * ex2.P[0,1] ``` ```{code-cell} ipython3 -ex2.β * ex2.u_prime(2.5) / ex2.u_prime(3.5) * ex2.P[1,0] +ex2.β * u_prime(2.5) / u_prime(3.5) * ex2.P[1,0] ``` @@ -1144,16 +1141,17 @@ ex2.A ```{code-cell} ipython3 # when the initial state is state 1 -print(f'α = {ex2.wealth_distribution(s0_idx=0)}') -print(f'ψ = \n{ex2.continuation_wealths()}') -print(f'J = \n{ex2.value_functionss()}') +print(f'α = {ex2.α}') +print(f'ψ = \n{ex2.ψ}') +print(f'J = \n{ex2.J}') ``` ```{code-cell} ipython3 # when the initial state is state 1 -print(f'α = {ex2.wealth_distribution(s0_idx=1)}') -print(f'ψ = \n{ex2.continuation_wealths()}') -print(f'J = \n{ex2.value_functionss()}') +ex2 = compute_rc_model(s, P, ys, s0_idx=1) +print(f'α = {ex2.α}') +print(f'ψ = \n{ex2.ψ}') +print(f'J = \n{ex2.J}') ``` ### Example 3 @@ -1163,20 +1161,20 @@ print(f'J = \n{ex2.value_functionss()}') K, n = 2, 2 # states -s = np.array([1, 2]) +s = jnp.array([1, 2]) # transition λ = 0.9 -P = np.array([[1-λ, λ], [0, 1]]) +P = jnp.array([[1-λ, λ], [0, 1]]) # endowments -ys = np.empty((n, K)) -ys[:, 0] = [1, 0] # y1 -ys[:, 1] = [0, 1] # y2 +ys = jnp.empty((n, K)) +ys.at[:, 0].set([1, 0]) # y1 +ys.at[:, 1].set([0, 1]) # y2 ``` ```{code-cell} ipython3 -ex3 = RecurCompetitive(s, P, ys) +ex3 = compute_rc_model(s, P, ys) ``` ```{code-cell} ipython3 @@ -1205,38 +1203,38 @@ Note that the natural debt limit for agent $1$ in state $2$ is $0$. ```{code-cell} ipython3 # when the initial state is state 1 -print(f'α = {ex3.wealth_distribution(s0_idx=0)}') -print(f'ψ = \n{ex3.continuation_wealths()}') -print(f'J = \n{ex3.value_functionss()}') +print(f'α = {ex3.α}') +print(f'ψ = \n{ex3.ψ}') +print(f'J = \n{ex3.J}') ``` ```{code-cell} ipython3 -# when the initial state is state 1 -print(f'α = {ex3.wealth_distribution(s0_idx=1)}') -print(f'ψ = \n{ex3.continuation_wealths()}') -print(f'J = \n{ex3.value_functionss()}') +# when the initial state is state 2 +ex3 = compute_rc_model(s, P, ys, s0_idx=1) +print(f'α = {ex3.α}') +print(f'ψ = \n{ex3.ψ}') +print(f'J = \n{ex3.J}') ``` For the specification of the Markov chain in example 3, let's take a look at how the equilibrium allocation changes as a function of transition probability $\lambda$. ```{code-cell} ipython3 -λ_seq = np.linspace(0, 0.99, 100) +λ_seq = jnp.linspace(0, 0.99, 100) # prepare containers -αs0_seq = np.empty((len(λ_seq), 2)) -αs1_seq = np.empty((len(λ_seq), 2)) +αs0_seq = jnp.empty((len(λ_seq), 2)) +αs1_seq = jnp.empty((len(λ_seq), 2)) for i, λ in enumerate(λ_seq): - P = np.array([[1-λ, λ], [0, 1]]) - ex3 = RecurCompetitive(s, P, ys) + P = jnp.array([[1-λ, λ], [0, 1]]) + ex3 = compute_rc_model(s, P, ys) # initial state s0 = 1 - α = ex3.wealth_distribution(s0_idx=0) - αs0_seq[i, :] = α + αs0_seq.at[i, :].set(ex3.α) # initial state s0 = 2 - α = ex3.wealth_distribution(s0_idx=1) - αs1_seq[i, :] = α + ex3 = compute_rc_model(s, P, ys, s0_idx=1) + αs1_seq.at[i, :].set(ex3.α) ``` ```{code-cell} ipython3 @@ -1259,7 +1257,7 @@ plt.show() K, n = 2, 3 # states -s = np.array([1, 2, 3]) +s = jnp.array([1, 2, 3]) # transition λ = .9 @@ -1267,16 +1265,16 @@ s = np.array([1, 2, 3]) δ = .05 # prosperous, moderate, and recession states -P = np.array([[1-λ, λ, 0], [μ/2, μ, μ/2], [(1-δ)/2, (1-δ)/2, δ]]) +P = jnp.array([[1-λ, λ, 0], [μ/2, μ, μ/2], [(1-δ)/2, (1-δ)/2, δ]]) # endowments -ys = np.empty((n, K)) -ys[:, 0] = [.25, .75, .2] # y1 -ys[:, 1] = [1.25, .25, .2] # y2 +ys = jnp.empty((n, K)) +ys.at[:, 0].set([.25, .75, .2]) # y1 +ys.at[:, 1].set([1.25, .25, .2]) # y2 ``` ```{code-cell} ipython3 -ex4 = RecurCompetitive(s, P, ys) +ex4 = compute_rc_model(s, P, ys) ``` ```{code-cell} ipython3 @@ -1296,10 +1294,11 @@ print('') for i in range(1, 4): # when the initial state is state i + ex4 = compute_rc_model(s, P, ys, s0_idx=i-1) print(f"when the initial state is state {i}") - print(f'α = {ex4.wealth_distribution(s0_idx=i-1)}') - print(f'ψ = \n{ex4.continuation_wealths()}') - print(f'J = \n{ex4.value_functionss()}\n') + print(f'α = {ex4.α}') + print(f'ψ = \n{ex4.ψ}') + print(f'J = \n{ex4.J}\n') ``` @@ -1312,19 +1311,19 @@ We now revisit the economy defined in example 1, but set the time horizon to be K, n = 2, 2 # states -s = np.array([0, 1]) +s = jnp.array([0, 1]) # transition -P = np.array([[.5, .5], [.5, .5]]) +P = jnp.array([[.5, .5], [.5, .5]]) # endowments -ys = np.empty((n, K)) -ys[:, 0] = 1 - s # y1 -ys[:, 1] = s # y2 +ys = jnp.empty((n, K)) +ys.at[:, 0].set(1 - s) # y1 +ys.at[:, 1].set(s) # y2 ``` ```{code-cell} ipython3 -ex1_finite = RecurCompetitive(s, P, ys, T=10) +ex1_finite = compute_rc_model(s, P, ys, T=10) ``` ```{code-cell} ipython3 @@ -1352,24 +1351,25 @@ In the finite time horizon case, `ψ` and `J` are returned as sequences. Components are ordered from $t=T$ to $t=0$. ```{code-cell} ipython3 -# when the initial state is state 2 -print(f'α = {ex1_finite.wealth_distribution(s0_idx=0)}') -print(f'ψ = \n{ex1_finite.continuation_wealths()}\n') -print(f'J = \n{ex1_finite.value_functionss()}') +# when the initial state is state 1 +print(f'α = {ex1_finite.α}') +print(f'ψ = \n{ex1_finite.ψ}\n') +print(f'J = \n{ex1_finite.J}') ``` ```{code-cell} ipython3 # when the initial state is state 2 -print(f'α = {ex1_finite.wealth_distribution(s0_idx=1)}') -print(f'ψ = \n{ex1_finite.continuation_wealths()}\n') -print(f'J = \n{ex1_finite.value_functionss()}') +ex1_finite = compute_rc_model(s, P, ys, s0_idx=1, T=10) +print(f'α = {ex1_finite.α}') +print(f'ψ = \n{ex1_finite.ψ}\n') +print(f'J = \n{ex1_finite.J}') ``` We can check the results with finite horizon converges to the ones with infinite horizon as $T \rightarrow \infty$. ```{code-cell} ipython3 -ex1_large = RecurCompetitive(s, P, ys, T=10000) -ex1_large.wealth_distribution(s0_idx=1) +ex1_large = compute_rc_model(s, P, ys, s0_idx=1, T=10000) +ex1_large.α ``` ```{code-cell} ipython3 @@ -1377,11 +1377,9 @@ ex1.V, ex1_large.V[-1] ``` ```{code-cell} ipython3 -ex1_large.continuation_wealths() ex1.ψ, ex1_large.ψ[-1] ``` ```{code-cell} ipython3 -ex1_large.value_functionss() ex1.J, ex1_large.J[-1] ``` From 598f541a872f3f5dcbacaf655a8e23d041327e26 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Sun, 16 Nov 2025 17:43:57 +0800 Subject: [PATCH 08/16] modified: lectures/ge_arrow.md --- lectures/ge_arrow.md | 237 +++++++++++++++++++++++++++---------------- 1 file changed, 149 insertions(+), 88 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index 4e9024a18..f7543991a 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -34,7 +34,7 @@ This lecture presents Python code for experimenting with competitive equilibria * Differences in their endowments make individuals want to reallocate consumption goods across time and Markov states -We impose restrictions that allow us to **Bellmanize** competitive equilibrium prices and quantities +We impose restrictions that allow us to *Bellmanize* competitive equilibrium prices and quantities We use Bellman equations to describe @@ -47,15 +47,15 @@ We use Bellman equations to describe In the course of presenting the model we shall encounter these important ideas -* a **resolvent operator** widely used in this class of models +* a *resolvent operator* widely used in this class of models -* absence of **borrowing limits** in finite horizon economies +* absence of *borrowing limits* in finite horizon economies -* state-by-state **borrowing limits** required in infinite horizon economies +* state-by-state *borrowing limits* required in infinite horizon economies -* a counterpart of the **law of iterated expectations** known as a **law of iterated values** +* a counterpart of the *law of iterated expectations* known as a *law of iterated values* -* a **state-variable degeneracy** that prevails within a competitive equilibrium and that opens the way to various appearances of resolvent operators +* a *state-variable degeneracy* that prevails within a competitive equilibrium and that opens the way to various appearances of resolvent operators ## The setting @@ -76,7 +76,7 @@ probability of observing a particular sequence of events $s^t$ is given by a probability measure $\pi_t(s^t)$. For $t > \tau$, we write the probability -of observing $s^t$ conditional on the realization of $s^\tau$as $\pi_t(s^t\vert s^\tau)$. +of observing $s^t$ conditional on the realization of $s^\tau$ as $\pi_t(s^t\vert s^\tau)$. We assume that trading occurs after observing $s_0$, @@ -345,7 +345,7 @@ $$ ```{note} The matrix geometric sum $(I - Q)^{-1} = I + Q + Q^2 + \cdots $ -is an example of a **resolvent operator**. +is an example of a *resolvent operator*. ``` Below, we describe an equilibrium model with trading of one-period Arrow securities in which the pricing kernel is endogenous. @@ -372,12 +372,12 @@ We'll use these objects to state a useful property in asset pricing theory. ### Laws of Iterated Expectations and Iterated Values -A **law of iterated values** has a mathematical structure that parallels a -**law of iterated expectations** +A *law of iterated values* has a mathematical structure that parallels a +*law of iterated expectations* We can describe its structure readily in the Markov setting of this lecture -Recall the following recursion satisfied by $j$ step ahead transition probabilites +Recall the following recursion satisfied by $j$ step ahead transition probabilities for our finite state Markov chain: $$ @@ -405,7 +405,7 @@ Q_j(s_{t+j}| s_t) = \sum_{s_{t+1}} Q_{j-1}(s_{t+j}| s_{t+1}) Q(s_{t+1} | s_t) $$ -The time $t$ **value** in Markov state $s_t$ of a time $t+j$ payout $d(s_{t+j})$ +The time $t$ *value* in Markov state $s_t$ of a time $t+j$ payout $d(s_{t+j})$ is @@ -413,7 +413,7 @@ $$ V(d(s_{t+j})|s_t) = \sum_{s_{t+j}} d(s_{t+j}) Q_j(s_{t+j}| s_t) $$ -The **law of iterated values** states +The *law of iterated values* states $$ V \left[ V (d(s_{t+j}) | s_{t+1}) | s_t \right] = V(d(s_{t+j})| s_t ) @@ -436,7 +436,7 @@ $$ Now we are ready to do some fun calculations. -We find it interesting to think in terms of analytical **inputs** into and **outputs** from our general equilibrium theorizing. +We find it interesting to think in terms of analytical *inputs* into and *outputs* from our general equilibrium theorizing. ### Inputs @@ -501,7 +501,7 @@ This follows from agent $k$'s first-order necessary conditions. But with the CRRA preferences that we have assumed, individual consumptions vary proportionately with aggregate consumption and therefore with the aggregate endowment. - * This is a consequence of our preference specification implying that **Engle curves** are affine in wealth and therefore satisfy conditions for **Gorman aggregation** + * This is a consequence of our preference specification implying that *Engel curves* are affine in wealth and therefore satisfy conditions for *Gorman aggregation* Thus, @@ -509,7 +509,7 @@ $$ c^k \left(s\right) = \alpha_k c\left(s\right) = \alpha_k y\left(s\right) $$ -for an arbitrary **distribution of wealth** in the form of an $K \times 1$ vector $\alpha$ +for an arbitrary *distribution of wealth* in the form of an $K \times 1$ vector $\alpha$ that satisfies $$ \alpha_k \in \left(0, 1\right), \quad \sum_{k=1}^K \alpha_k = 1 $$ @@ -525,12 +525,12 @@ Note that $Q_{ij}$ is independent of vector $\alpha$. -**Key finding:** We can compute competitive equilibrium **prices** prior to computing a **distribution of wealth**. +*Key finding:* We can compute competitive equilibrium *prices* prior to computing a *distribution of wealth*. ### Values -Having computed an equilibrium pricing kernel $Q$, we can compute several **values** that are required +Having computed an equilibrium pricing kernel $Q$, we can compute several *values* that are required to pose or represent the solution of an individual household's optimum problem. @@ -582,7 +582,7 @@ even if he consumes zero goods forevermore. ```{prf:remark} If we have an Inada condition at zero consumption or just impose that consumption -be nonnegative, then in a *finite horizon* economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. See the section on a [finite orizon economy](#finite-horizon) below. +be nonnegative, then in a *finite horizon* economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. See the section on a [finite horizon economy](#finite-horizon) below. ``` @@ -846,6 +846,7 @@ import jax import jax.numpy as jnp import numpy as np from typing import NamedTuple +from functools import partial import matplotlib.pyplot as plt ``` @@ -869,13 +870,14 @@ class RecurCompetitive(NamedTuple): s: jax.Array # state vector P: jax.Array # transition matrix ys: jax.Array # endowments ys = [y1, y2, .., yT] - y: jax.Array # total endownment under each state - n: float # number of states - K: float # number of agents + y: jax.Array # total endowment under each state + n: int # number of states + K: int # number of agents γ: float # risk aversion β: float # discount rate T: float # time horizon, 0 if infinite - Q: jax.Array # pricing kernek + Q: jax.Array # pricing kernel + V: jax.Array # resolvent / partial-sum matrices PRF: jax.Array # price of risk-free bond R: jax.Array # risk-free rate A: jax.Array # natural debt limit @@ -883,51 +885,94 @@ class RecurCompetitive(NamedTuple): ψ: jax.Array # continuation value J: jax.Array # optimal value - -def compute_rc_model(s, P, ys, 0_idx=0, γ=0.5, β=0.98, T=0): +@partial(jax.jit, static_argnums=(6,)) +def compute_rc_model(s, P, ys, s0_idx=0, γ=0.5, β=0.98, T=0): + """Complete equilibrium objects under the endogenous pricing kernel. + + Args + ---- + s : array-like + Markov states. + P : array-like + Transition matrix. + ys : array-like + Endowment matrix; rows index states, columns index agents. + s0_idx : int, optional + Index of the initial zero-asset-holding state. + γ : float, optional + Risk aversion parameter. + β : float, optional + Discount factor. + T : int, optional + Number of periods; 0 means an infinite-horizon economy. + + Returns + ------- + RecurCompetitive + Instance containing all parameters and computed equilibrium results. + """ n, K = ys.shape - y = jnp.sum(ys, 1) + y = jnp.sum(ys, axis=1) def u(c): - "The CRRA utility" + "CRRA utility evaluated elementwise." return c ** (1 - γ) / (1 - γ) def u_prime(c): - "The first derivative of CRRA utility" + "Marginal utility for the CRRA specification." return c ** (-γ) def pricing_kernel(c): - "Compute the pricing kernel matrix Q" + "Build the Arrow-security pricing kernel matrix." Q = jnp.empty((n, n)) - for i in range(n): - for j in range(n): + # fori_loop iterates over each state i while carrying the partially + # filled matrix Q as the loop carry. + def body_i(i, Q): + # fills row i entry-by-entry. + def body_j(j, q): ratio = u_prime(c[j]) / u_prime(c[i]) - Q[i, j] = β * ratio * P[i, j] + # Return a (n,) array + return q.at[j].set(β * ratio * P[i, j]) + + q = jax.lax.fori_loop( + 0, n, body_j, jnp.zeros((n,)) + ) + return Q.at[i, :].set(q) + Q = jax.lax.fori_loop( + 0, n, body_i, jnp.zeros((n, n)) + ) return Q - def V(Q): + def resolvent_operator(Q): + "Compute the resolvent or finite partial sums of Q depending on T." if T==0: V = jnp.empty((1, n, n)) - V.at[0].set(jnp.linalg.inv(jnp.eye(n) - Q)) + V = V.at[0].set(jnp.linalg.inv(jnp.eye(n) - Q)) # V = [I + Q + Q^2 + ... + Q^T] (finite case) else: V = jnp.empty((T+1, n, n)) - V.at[0].set(jnp.eye(n)) + V = V.at[0].set(jnp.eye(n)) Qt = jnp.eye(n) - for t in range(1, T+1): + + # Loop body advances the Q power and accumulates the geometric sum. + def body(t, carry): + Qt, V = carry Qt = Qt @ Q - V.at[t].set(V[t-1] + Qt) + V = V.at[t].set(V[t-1] + Qt) + return Qt, V + + _, V = jax.lax.fori_loop(1, T+1, body, (Qt, V)) return V - def A(ys, V): - "Give endownments and pricing kernel, compute natural debt limit" + def natural_debt_limit(ys, V): + "Compute natural debt limits from the terminal resolvent block." return V[-1] @ ys - def wealth_distribution(V, ys, y, 0_idx): - "Solve for wealth distribution α" + def wealth_distribution(V, ys, y, s0_idx): + "Recover equilibrium wealth shares α from the initial state row." # row of V corresponding to s0 Vs0 = V[-1, s0_idx, :] @@ -936,61 +981,71 @@ def compute_rc_model(s, P, ys, 0_idx=0, γ=0.5, β=0.98, T=0): return α def continuation_wealths(V, α): - "Given α, compute the continuation wealths ψ" + "Back out continuation wealths for each agent." diff = jnp.empty((n, K)) - for k in range(K): - diff.at[:, k].set(α[k] * y - ys[:, k]) + + # Loop scatters each agent's state-dependent surplus into the column k. + def body(k, carry): + diff = carry + return diff.at[:, k].set(α[k] * y - ys[:, k]) + + # Applies body sequentially while threading diff. + diff = jax.lax.scan(0, K, body, diff) ψ = V @ diff return ψ def price_risk_free_bond(Q): - "Give Q, compute price of one-period risk free bond" - - PRF = jnp.sum(Q, axis=1) - - return PRF + "Given Q, compute price of one-period risk-free bond" + return jnp.sum(Q, axis=1) def risk_free_rate(Q): "Given Q, compute one-period gross risk-free interest rate R" + return jnp.reciprocal(price_risk_free_bond(Q)) - R = jnp.sum(Q, axis=1) - R = jnp.reciprocal(R) - - return R def value_functions(α, y): - "Given α, compute the optimal value functions J in equilibrium" + "Assemble lifetime value functions for each agent." # compute (I - βP)^(-1) in infinite case if T==0: P_seq = jnp.empty((1, n, n)) - P_seq.at[0].set( + P_seq = P_seq.at[0].set( jnp.linalg.inv(jnp.eye(n) - β * P) ) # and (I + βP + ... + β^T P^T) in finite case else: P_seq = jnp.empty((T+1, n, n)) - P_seq.at[0].set(np.eye(n)) + P_seq = P_seq.at[0].set(np.eye(n)) Pt = jnp.eye(n) - for t in range(1, T+1): + + def body(t, carry): + Pt, P_seq = carry Pt = Pt @ P - P_seq.at[t].set(P_seq[t-1] + Pt * β ** t) + P_seq = P_seq.at[t].set(P_seq[t-1] + Pt * β ** t) + return Pt, P_seq + + _, P_seq = jax.lax.fori_loop( + 1, T+1, body, (Pt, P_seq) + ) # compute the matrix [u(α_1 y), ..., u(α_K, y)] flow = jnp.empty((n, K)) - for k in range(K): - flow.at[:, k].set(u(α[k] * y)) - + def body(k, carry): + flow = carry + return flow.at[:, k].set(u(α[k] * y)) + + flow = jax.lax.fori_loop(0, K, body, flow) + J = P_seq @ flow return J - Q = pricing_kernek(y) - V = V(Q) - A = A(ys, Q) + Q = pricing_kernel(y) + V = resolvent_operator(Q) + A = natural_debt_limit(ys, V) α = wealth_distribution(V, ys, y, s0_idx) ψ = continuation_wealths(V, α) PRF = price_risk_free_bond(Q) @@ -998,7 +1053,7 @@ def compute_rc_model(s, P, ys, 0_idx=0, γ=0.5, β=0.98, T=0): J = value_functions(α, y) return RecurCompetitive( - s=s, P=P, ys=ys, y=y n=n, K=K, γ=γ, β=β, T=T, + s=s, P=P, ys=ys, y=y, n=n, K=K, γ=γ, β=β, T=T, Q=Q, V=V, A=A, α=α, ψ=ψ, PRF=PRF, R=R, J=J ) ``` @@ -1029,8 +1084,8 @@ P = jnp.array([[.5, .5], [.5, .5]]) # endowments ys = jnp.empty((n, K)) -ys.at[:, 0].set(1 - s) # y1 -ys.at[:, 1].set(s) # y2 +ys = ys.at[:, 0].set(1 - s) # y1 +ys = ys.at[:, 1].set(s) # y2 ``` ```{code-cell} ipython3 @@ -1043,7 +1098,7 @@ ex1.ys ``` ```{code-cell} ipython3 -# pricing kernal +# pricing kernel ex1.Q ``` @@ -1068,7 +1123,7 @@ print(f'J = \n{ex1.J}') # when the initial state is state 2 ex1 = compute_rc_model(s, P, ys, s0_idx=1) print(f'α = {ex1.α}') -print(f'ψ = \n{wx1.ψ}') +print(f'ψ = \n{ex1.ψ}') print(f'J = \n{ex1.J}') ``` @@ -1086,8 +1141,8 @@ P = jnp.array([[.5, .5], [.5, .5]]) # endowments ys = jnp.empty((n, K)) -ys.at[:, 0].set(1.5) # y1 -ys.at[:, 1].set(s) # y2 +ys = ys.at[:, 0].set(1.5) # y1 +ys = ys.at[:, 1].set(s) # y2 ``` ```{code-cell} ipython3 @@ -1099,7 +1154,7 @@ ex2 = compute_rc_model(s, P, ys) print("ys = \n", ex2.ys) -# pricing kernal +# pricing kernel print ("Q = \n", ex2.Q) # Risk free rate R @@ -1107,11 +1162,11 @@ print("R = ", ex2.R) ``` ```{code-cell} ipython3 -# pricing kernal +# pricing kernel ex2.Q ``` -Note that the pricing kernal in example economies 1 and 2 differ. +Note that the pricing kernel in example economies 1 and 2 differ. This comes from differences in the aggregate endowments in state 1 and 2 in example 1. @@ -1169,8 +1224,8 @@ P = jnp.array([[1-λ, λ], [0, 1]]) # endowments ys = jnp.empty((n, K)) -ys.at[:, 0].set([1, 0]) # y1 -ys.at[:, 1].set([0, 1]) # y2 +ys = ys.at[:, 0].set([1, 0]) # y1 +ys = ys.at[:, 1].set([0, 1]) # y2 ``` ```{code-cell} ipython3 @@ -1225,16 +1280,22 @@ For the specification of the Markov chain in example 3, let's take a look at how αs0_seq = jnp.empty((len(λ_seq), 2)) αs1_seq = jnp.empty((len(λ_seq), 2)) -for i, λ in enumerate(λ_seq): +def body(i, carry): + αs0_seq, αs1_seq = carry + λ = λ_seq[i] P = jnp.array([[1-λ, λ], [0, 1]]) - ex3 = compute_rc_model(s, P, ys) - # initial state s0 = 1 - αs0_seq.at[i, :].set(ex3.α) + ex3_s0 = compute_rc_model(s, P, ys) # initial state s0 = 2 - ex3 = compute_rc_model(s, P, ys, s0_idx=1) - αs1_seq.at[i, :].set(ex3.α) + ex3_s1 = compute_rc_model(s, P, ys, s0_idx=1) + + return (αs0_seq.at[i, :].set(ex3_s0.α), + αs1_seq.at[i, :].set(ex3_s1.α)) + +αs0_seq, αs1_seq = jax.lax.fori_loop( + 0, 100, body, (αs0_seq, αs1_seq) + ) ``` ```{code-cell} ipython3 @@ -1269,8 +1330,8 @@ P = jnp.array([[1-λ, λ, 0], [μ/2, μ, μ/2], [(1-δ)/2, (1-δ)/2, δ]]) # endowments ys = jnp.empty((n, K)) -ys.at[:, 0].set([.25, .75, .2]) # y1 -ys.at[:, 1].set([1.25, .25, .2]) # y2 +ys = ys.at[:, 0].set([.25, .75, .2]) # y1 +ys = ys.at[:, 1].set([1.25, .25, .2]) # y2 ``` ```{code-cell} ipython3 @@ -1281,7 +1342,7 @@ ex4 = compute_rc_model(s, P, ys) # endowments print("ys = \n", ex4.ys) -# pricing kernal +# pricing kernel print ("Q = \n", ex4.Q) # Risk free rate R @@ -1318,8 +1379,8 @@ P = jnp.array([[.5, .5], [.5, .5]]) # endowments ys = jnp.empty((n, K)) -ys.at[:, 0].set(1 - s) # y1 -ys.at[:, 1].set(s) # y2 +ys = ys.at[:, 0].set(1 - s) # y1 +ys = ys.at[:, 1].set(s) # y2 ``` ```{code-cell} ipython3 @@ -1337,7 +1398,7 @@ ex1_finite.ys ``` ```{code-cell} ipython3 -# pricing kernal +# pricing kernel ex1_finite.Q ``` From 283a5c476590724dfa7b350a04b5bd317eaa449a Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Mon, 17 Nov 2025 16:55:24 +0800 Subject: [PATCH 09/16] modified: lectures/ge_arrow.md --- lectures/ge_arrow.md | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index f7543991a..1ed31ec6f 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -145,7 +145,7 @@ $$ for all $t$ and for all $s^t$. -## Recursive Formulation +## Recursive formulation Following descriptions in section 9.3.3 of Ljungqvist and Sargent {cite}`Ljungqvist2012` chapter 9, we set up a competitive equilibrium of a pure exchange economy with complete markets in one-period Arrow securities. @@ -239,7 +239,7 @@ are zero net aggregate claims. -## State Variable Degeneracy +## State variable degeneracy Please see Ljungqvist and Sargent {cite}`Ljungqvist2012` for a description of timing protocol for trades consistent with an Arrow-Debreu vision in which @@ -284,7 +284,7 @@ This outcome depends critically on there being complete markets in Arrow securit For example, it does not prevail in the incomplete markets setting of this lecture {doc}`The Aiyagari Model ` -## Markov Asset Prices +## Markov asset prices Let's start with a brief summary of formulas for computing asset prices in @@ -315,7 +315,7 @@ $$ * The gross rate of return on a one-period risk-free bond Markov state $\bar s_i$ is $R_i = (\sum_j Q_{i,j})^{-1}$ -### Exogenous Pricing Kernel +### Exogenous pricing kernel At this point, we'll take the pricing kernel $Q$ as exogenous, i.e., determined outside the model @@ -352,7 +352,7 @@ Below, we describe an equilibrium model with trading of one-period Arrow securit In constructing our model, we'll repeatedly encounter formulas that remind us of our asset pricing formulas. -### Multi-Step-Forward Transition Probabilities and Pricing Kernels +### Multi-step-forward transition probabilities and pricing kernels The $(i,j)$ component of the $\ell$-step ahead transition probability $P^\ell$ is @@ -370,7 +370,7 @@ $$ We'll use these objects to state a useful property in asset pricing theory. -### Laws of Iterated Expectations and Iterated Values +### Laws of iterated expectations and iterated values A *law of iterated values* has a mathematical structure that parallels a *law of iterated expectations* @@ -432,7 +432,7 @@ V \left[ V ( d(s_{t+j}) | s_{t+1} ) | s_t \right] \end{aligned} $$ -## General Equilibrium +## General equilibrium Now we are ready to do some fun calculations. @@ -483,7 +483,7 @@ $$ * A collection of $n \times 1$ vectors of individual $k$ consumptions: $c^k\left(s\right), k=1,\ldots, K$ -### $Q$ is the Pricing Kernel +### $Q$ is the pricing kernel For any agent $k \in \left[1, \ldots, K\right]$, at the equilibrium allocation, @@ -587,7 +587,7 @@ be nonnegative, then in a *finite horizon* economy with sequential trading of on -### Continuation Wealth +### Continuation wealth Continuation wealth plays an important role in Bellmanizing a competitive equilibrium with sequential trading of a complete set of one-period Arrow securities. @@ -645,7 +645,7 @@ the economy begins with all agents being debt-free and financial-asset-free at Note that all agents' continuation wealths recurrently return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. ``` -### Optimal Portfolios +### Optimal portfolios A nifty feature of the model is that an optimal portfolio of a type $k$ agent equals the continuation wealth that we just computed. @@ -656,7 +656,7 @@ $$ a_k(s) = \psi^k(s), \quad s \in \left[\bar s_1, \ldots, \bar s_n \right] $$ (eqn:optport) -### Equilibrium Wealth Distribution $\alpha$ +### Equilibrium wealth distribution $\alpha$ With the initial state being a particular state $s_0 \in \left[\bar{s}_1, \ldots, \bar{s}_n\right]$, @@ -703,7 +703,7 @@ $$ J^k = (I - \beta P)^{-1} u(\alpha_k y) , \quad u(c) = \frac{c^{1-\gamma}}{1- where it is understood that $ u(\alpha_k y)$ is a vector. -## Finite Horizon +## Finite horizon We now describe a finite-horizon version of the economy that operates for $T+1$ periods $t \in {\bf T} = \{ 0, 1, \ldots, T\}$. @@ -717,7 +717,7 @@ one-period utility function $u(c)$ satisfies an Inada condition that sets the ma limits borrowing. -### Continuation Wealths +### Continuation wealths We denote a $K \times 1$ vector of state-dependent continuation wealths in Markov state $s$ at time $t$ as @@ -834,7 +834,7 @@ where it is understood that $ u(\alpha_k y)$ is a vector. -## Python Code +## Python code We are ready to dive into some Python code. @@ -885,6 +885,7 @@ class RecurCompetitive(NamedTuple): ψ: jax.Array # continuation value J: jax.Array # optimal value + @partial(jax.jit, static_argnums=(6,)) def compute_rc_model(s, P, ys, s0_idx=0, γ=0.5, β=0.98, T=0): """Complete equilibrium objects under the endogenous pricing kernel. From 7af6f1fe30ef7d1d284372e1ddf255f6388ef879 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Mon, 17 Nov 2025 17:03:17 +0800 Subject: [PATCH 10/16] modified: lectures/ge_arrow.md --- lectures/ge_arrow.md | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index 1ed31ec6f..ed74f2e5d 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -328,7 +328,9 @@ Two examples would be We'll write down implications of Markov asset pricing in a nutshell for two types of assets - * the price in Markov state $s$ at time $t$ of a **cum dividend** stock that entitles the owner at the beginning of time $t$ to the time $t$ dividend and the option to sell the asset at time $t+1$. The price evidently satisfies $p^h(\bar s_i) = d^h(\bar s_i) + \sum_j Q_{ij} p^h(\bar s_j) $, which implies that the vector $p^h$ satisfies $p^h = d^h + Q p^h$ which implies the formula + * the price in Markov state $s$ at time $t$ of a **cum dividend** stock that entitles the owner at the beginning of time $t$ to the time $t$ dividend and the option to sell the asset at time $t+1$. + + * The price evidently satisfies $p^h(\bar s_i) = d^h(\bar s_i) + \sum_j Q_{ij} p^h(\bar s_j) $, which implies that the vector $p^h$ satisfies $p^h = d^h + Q p^h$ which implies the formula $$ p^h = (I - Q)^{-1} d^h @@ -337,7 +339,9 @@ $$ -* the price in Markov state $s$ at time $t$ of an **ex dividend** stock that entitles the owner at the end of time $t$ to the time $t+1$ dividend and the option to sell the stock at time $t+1$. The price is +* the price in Markov state $s$ at time $t$ of an **ex dividend** stock that entitles the owner at the end of time $t$ to the time $t+1$ dividend and the option to sell the stock at time $t+1$. + +The price is $$ p^h = (I - Q)^{-1} Q d^h @@ -582,7 +586,9 @@ even if he consumes zero goods forevermore. ```{prf:remark} If we have an Inada condition at zero consumption or just impose that consumption -be nonnegative, then in a *finite horizon* economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. See the section on a [finite horizon economy](#finite-horizon) below. +be nonnegative, then in a *finite horizon* economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. + +See the section on a [finite horizon economy](#finite-horizon) below. ``` @@ -637,7 +643,9 @@ Note that $\sum_{k=1}^K \psi^k = {0}_{n \times 1}$. ```{prf:remark} At the initial state $s_0 \in \begin{bmatrix} \bar s_1, \ldots, \bar s_n \end{bmatrix}$, -the continuation wealth $\psi^k(s_0) = 0$ for all agents $k = 1, \ldots, K$. This indicates that +the continuation wealth $\psi^k(s_0) = 0$ for all agents $k = 1, \ldots, K$. + +This indicates that the economy begins with all agents being debt-free and financial-asset-free at time $0$, state $s_0$. ``` @@ -771,13 +779,17 @@ Note that $\sum_{k=1}^K \psi_t^k = {0}_{n \times 1}$ for all $t \in {\bf T}$. ```{prf:remark} At the initial state $s_0 \in \begin{bmatrix} \bar s_1, \ldots, \bar s_n \end{bmatrix}$, - for all agents $k = 1, \ldots, K$, continuation wealth $\psi_0^k(s_0) = 0$. This indicates that + for all agents $k = 1, \ldots, K$, continuation wealth $\psi_0^k(s_0) = 0$. + + This indicates that the economy begins with all agents being debt-free and financial-asset-free at time $0$, state $s_0$. ``` ```{prf:remark} -Note that all agents' continuation wealths return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. This will recur if the Markov chain makes the initial state $s_0$ recurrent. +Note that all agents' continuation wealths return to zero when the Markov state returns to whatever value $s_0$ it had at time $0$. + +This will recur if the Markov chain makes the initial state $s_0$ recurrent. ``` From 868fb4172aa3685aa5914955fe7acc75add796e4 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Mon, 17 Nov 2025 17:14:10 +0800 Subject: [PATCH 11/16] recover markov_asset.md file --- lectures/markov_asset.md | 254 +++++++++++++++------------------------ 1 file changed, 98 insertions(+), 156 deletions(-) diff --git a/lectures/markov_asset.md b/lectures/markov_asset.md index ef6c8d3c2..1739dde0f 100644 --- a/lectures/markov_asset.md +++ b/lectures/markov_asset.md @@ -35,16 +35,6 @@ kernelspec: "Asset pricing is all about covariances" -- Lars Peter Hansen ``` -```{admonition} GPU -:class: warning - -This lecture is accelerated via [hardware](status:machine-details) that has access to a GPU and JAX for GPU programming. - -Free GPUs are available on Google Colab. To use this option, please click on the play icon top right, select Colab, and set the runtime environment to include a GPU. - -Alternatively, if you have your own GPU, you can follow the [instructions](https://github.com/google/jax) for installing JAX with GPU support. If you would like to install JAX running on the `cpu` only you can use `pip install jax[cpu]` -``` - In addition to what's in Anaconda, this lecture will need the following libraries: ```{code-cell} ipython @@ -85,10 +75,7 @@ Let's start with some imports: import matplotlib.pyplot as plt import numpy as np import quantecon as qe -import jax -import jax.numpy as jnp -from jax.numpy.linalg import eigvals, solve -from typing import NamedTuple +from numpy.linalg import eigvals, solve ``` ## {index}`Pricing Models ` @@ -104,7 +91,7 @@ Let $\{d_t\}_{t \geq 0}$ be a stream of dividends Let's look at some equations that we expect to hold for prices of assets under ex-dividend contracts (we will consider cum-dividend pricing in the exercises). -### Risk-neutral pricing +### Risk-Neutral Pricing ```{index} single: Pricing Models; Risk-Neutral ``` @@ -129,7 +116,7 @@ Here ${\mathbb E}_t [y]$ denotes the best forecast of $y$, conditioned on inform More precisely, ${\mathbb E}_t [y]$ is the mathematical expectation of $y$ conditional on information available at time $t$. -### Pricing with random discount factor +### Pricing with Random Discount Factor ```{index} single: Pricing Models; Risk Aversion ``` @@ -158,7 +145,7 @@ This is because such assets pay well when funds are more urgently wanted. We give examples of how the stochastic discount factor has been modeled below. -### Asset pricing and covariances +### Asset Pricing and Covariances Recall that, from the definition of a conditional covariance ${\rm cov}_t (x_{t+1}, y_{t+1})$, we have @@ -185,9 +172,9 @@ It is useful to regard equation {eq}`lteeqs102` as a generalization of equatio Equation {eq}`lteeqs102` asserts that the covariance of the stochastic discount factor with the one period payout $d_{t+1} + p_{t+1}$ is an important determinant of the price $p_t$. -We give examples of some models of stochastic discount factors that have been proposed later in this lecture and also in a {doc}`later lecture`. +We give examples of some models of stochastic discount factors that have been proposed later in this lecture and also in a [later lecture](https://python-advanced.quantecon.org/lucas_model.html). -### The price-dividend ratio +### The Price-Dividend Ratio Aside from prices, another quantity of interest is the **price-dividend ratio** $v_t := p_t / d_t$. @@ -203,7 +190,7 @@ v_t = {\mathbb E}_t \left[ m_{t+1} \frac{d_{t+1}}{d_t} (1 + v_{t+1}) \right] Below we'll discuss the implication of this equation. -## Prices in the risk-neutral case +## Prices in the Risk-Neutral Case What can we say about price dynamics on the basis of the models described above? @@ -216,7 +203,7 @@ For now we'll study the risk-neutral case in which the stochastic discount fac We'll focus on how an asset price depends on a dividend process. -### Example 1: constant dividends +### Example 1: Constant Dividends The simplest case is risk-neutral price of a constant, non-random dividend stream $d_t = d > 0$. @@ -247,7 +234,7 @@ This is the equilibrium price in the constant dividend case. Indeed, simple algebra shows that setting $p_t = \bar p$ for all $t$ satisfies the difference equation $p_t = \beta (d + p_{t+1})$. -### Example 2: dividends with deterministic growth paths +### Example 2: Dividends with Deterministic Growth Paths Consider a growing, non-random dividend process $d_{t+1} = g d_t$ where $0 < g \beta < 1$. @@ -280,7 +267,7 @@ $$ This is called the *Gordon formula*. (mass_mg)= -### Example 3: Markov growth, risk-neutral pricing +### Example 3: Markov Growth, Risk-Neutral Pricing Next, we consider a dividend process @@ -323,21 +310,14 @@ The next figure shows a simulation, where * $\{X_t\}$ evolves as a discretized AR1 process produced using {ref}`Tauchen's method `. * $g_t = \exp(X_t)$, so that $\ln g_t = X_t$ is the growth rate. -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: | - State, growth, and dividend simulation - name: fig_markov_sim ---- +```{code-cell} ipython n = 7 mc = qe.tauchen(n, 0.96, 0.25) sim_length = 80 -x_series = mc.simulate(sim_length, init=jnp.median(mc.state_values)) -g_series = jnp.exp(x_series) -d_series = jnp.cumprod(g_series) # Assumes d_0 = 1 +x_series = mc.simulate(sim_length, init=np.median(mc.state_values)) +g_series = np.exp(x_series) +d_series = np.cumprod(g_series) # Assumes d_0 = 1 series = [x_series, g_series, d_series, np.log(d_series)] labels = ['$X_t$', '$g_t$', '$d_t$', r'$\log \, d_t$'] @@ -350,7 +330,7 @@ plt.tight_layout() plt.show() ``` -#### Pricing formula +#### Pricing Formula To obtain asset prices in this setting, let's adapt our analysis from the case of deterministic growth. @@ -420,25 +400,18 @@ As before, we'll generate $\{X_t\}$ as a {ref}`discretized AR1 process Here's the code, including a test of the spectral radius condition -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: | - Price-dividend ratio risk-neutral case - name: fig_pdv_neutral ---- +```{code-cell} python3 n = 25 # Size of state space β = 0.9 mc = qe.tauchen(n, 0.96, 0.02) -K = mc.P * jnp.exp(mc.state_values) +K = mc.P * np.exp(mc.state_values) warning_message = "Spectral radius condition fails" -assert jnp.max(jnp.abs(eigvals(K))) < 1 / β, warning_message +assert np.max(np.abs(eigvals(K))) < 1 / β, warning_message -I = jnp.identity(n) -v = solve(I - β * K, β * K @ jnp.ones(n)) +I = np.identity(n) +v = solve(I - β * K, β * K @ np.ones(n)) fig, ax = plt.subplots() ax.plot(mc.state_values, v, 'g-o', lw=2, alpha=0.7, label='$v$') @@ -467,7 +440,7 @@ We'll price several distinct assets, including * A consol (a type of bond issued by the UK government in the 19th century) * Call options on a consol -### Pricing a Lucas tree +### Pricing a Lucas Tree ```{index} single: Finite Markov Asset Pricing; Lucas Tree ``` @@ -496,7 +469,7 @@ m_{t+1} = \beta \frac{u'(c_{t+1})}{u'(c_t)} where $u$ is a concave utility function and $c_t$ is time $t$ consumption of a representative consumer. -(A derivation of this expression is given in a {doc}`later lecture`) +(A derivation of this expression is given in a [later lecture](https://python-advanced.quantecon.org/lucas_model.html)) Assume the existence of an endowment that follows growth process {eq}`mass_fmce`. @@ -566,51 +539,46 @@ v = (I - \beta J)^{-1} \beta J {\mathbb 1} We will define a function tree_price to compute $v$ given parameters stored in the class AssetPriceModel -```{code-cell} ipython3 -class AssetPriceModel(NamedTuple): +```{code-cell} python3 +class AssetPriceModel: """ A class that stores the primitives of the asset pricing model. Parameters ---------- + β : scalar, float + Discount factor mc : MarkovChain - Contains the transition matrix and set of state values + Contains the transition matrix and set of state values for the state + process + γ : scalar(float) + Coefficient of risk aversion g : callable The function mapping states to growth rates - β : float - Discount factor - γ : float - Coefficient of risk aversion - n: int - The number of states - """ - mc: qe.MarkovChain - g: callable - β: float - γ: float - n: int - - -def create_ap_model(mc=None, g=jnp.exp, β=0.96, γ=2.0): - """Create an AssetPriceModel class""" - if mc is None: - n, ρ, σ = 25, 0.9, 0.02 - mc = qe.tauchen(n, ρ, σ) - else: - mc = mc - n = mc.P.shape[0] - - return AssetPriceModel(mc=mc, g=g, β=β, γ=γ, n=n) - - -def test_stability(Q, β): - """ - Stability test for a given matrix Q. + """ - sr = np.max(np.abs(eigvals(Q))) - if not sr < 1 / β: - msg = f"Spectral radius condition failed with radius = {sr}" - raise ValueError(msg) + def __init__(self, β=0.96, mc=None, γ=2.0, g=np.exp): + self.β, self.γ = β, γ + self.g = g + + # A default process for the Markov chain + if mc is None: + self.ρ = 0.9 + self.σ = 0.02 + self.mc = qe.tauchen(n, self.ρ, self.σ) + else: + self.mc = mc + + self.n = self.mc.P.shape[0] + + def test_stability(self, Q): + """ + Stability test for a given matrix Q. + """ + sr = np.max(np.abs(eigvals(Q))) + if not sr < 1 / self.β: + msg = f"Spectral radius condition failed with radius = {sr}" + raise ValueError(msg) def tree_price(ap): @@ -633,11 +601,11 @@ def tree_price(ap): J = P * ap.g(y)**(1 - γ) # Make sure that a unique solution exists - test_stability(J, β) + ap.test_stability(J) # Compute v - I = jnp.identity(ap.n) - Ones = jnp.ones(ap.n) + I = np.identity(ap.n) + Ones = np.ones(ap.n) v = solve(I - β * J, β * J @ Ones) return v @@ -646,25 +614,19 @@ def tree_price(ap): Here's a plot of $v$ as a function of the state for several values of $\gamma$, with a positively correlated Markov process and $g(x) = \exp(x)$ -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: | - Lucas tree prices for varying risk aversion - name: fig_lucas_gamma ---- +```{code-cell} python3 γs = [1.2, 1.4, 1.6, 1.8, 2.0] -ap = create_ap_model() +ap = AssetPriceModel() states = ap.mc.state_values fig, ax = plt.subplots() for γ in γs: - tem_ap = create_ap_model(mc=ap.mc, g=ap.g, β=ap.β, γ=γ) - v = tree_price(tem_ap) + ap.γ = γ + v = tree_price(ap) ax.plot(states, v, lw=2, alpha=0.6, label=rf"$\gamma = {γ}$") +ax.set_title('Price-dividend ratio as a function of the state') ax.set_ylabel("price-dividend ratio") ax.set_xlabel("state") ax.legend(loc='upper right') @@ -744,7 +706,7 @@ p = (I - \beta M)^{-1} \beta M \zeta {\mathbb 1} The above is implemented in the function consol_price. -```{code-cell} ipython3 +```{code-cell} python3 def consol_price(ap, ζ): """ Computes price of a consol bond with payoff ζ @@ -753,7 +715,7 @@ def consol_price(ap, ζ): ---------- ap: AssetPriceModel An instance of AssetPriceModel containing primitives - + ζ : scalar(float) Coupon of the console @@ -761,17 +723,18 @@ def consol_price(ap, ζ): ------- p : array_like(float) Console bond prices + """ # Simplify names, set up matrices β, γ, P, y = ap.β, ap.γ, ap.mc.P, ap.mc.state_values M = P * ap.g(y)**(- γ) # Make sure that a unique solution exists - test_stability(M, β) + ap.test_stability(M) # Compute price - I = jnp.identity(ap.n) - Ones = jnp.ones(ap.n) + I = np.identity(ap.n) + Ones = np.ones(ap.n) p = solve(I - β * M, β * ζ * M @ Ones) return p @@ -849,7 +812,7 @@ Start at some initial $w$ and iterate with $T$ to convergence . We can find the solution with the following function call_option -```{code-cell} ipython3 +```{code-cell} python3 def call_option(ap, ζ, p_s, ϵ=1e-7): """ Computes price of a call option on a consol bond. @@ -865,7 +828,7 @@ def call_option(ap, ζ, p_s, ϵ=1e-7): p_s : scalar(float) Strike price - ϵ : scalar(float), optional(default=1e-7) + ϵ : scalar(float), optional(default=1e-8) Tolerance for infinite horizon problem Returns @@ -879,42 +842,26 @@ def call_option(ap, ζ, p_s, ϵ=1e-7): M = P * ap.g(y)**(- γ) # Make sure that a unique consol price exists - test_stability(M, β) + ap.test_stability(M) # Compute option price p = consol_price(ap, ζ) - w = jnp.zeros(ap.n) + w = np.zeros(ap.n) error = ϵ + 1 - - def step(state): - w, error = state + while error > ϵ: # Maximize across columns - w_new = jnp.maximum(β * M @ w, p - p_s) + w_new = np.maximum(β * M @ w, p - p_s) # Find maximal difference of each component and update - error_new = jnp.amax(jnp.abs(w - w_new)) - return (w_new, error_new) - - # Check whether converged - def cond(state): - _, error = state - return error > ϵ + error = np.amax(np.abs(w - w_new)) + w = w_new - final_w, _ = jax.lax.while_loop(cond, step, (w, error)) - - return final_w + return w ``` Here's a plot of $w$ compared to the consol price when $P_S = 40$ -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: | - Consol price and call option value - name: fig_consol_call ---- -ap = create_ap_model(β=0.9) +```{code-cell} python3 +ap = AssetPriceModel(β=0.9) ζ = 1.0 strike_price = 40 @@ -1013,12 +960,10 @@ $$ Consider the following primitives -```{code-cell} ipython3 +```{code-cell} python3 n = 5 # Size of State Space -P = jnp.full((n, n), 0.0125) -P = P.at[jnp.arange(n), jnp.arange(n)].set( - P[jnp.arange(n), jnp.arange(n)] + 1 - P.sum(1) - ) +P = np.full((n, n), 0.0125) +P[range(n), range(n)] += 1 - P.sum(1) # State values of the Markov chain s = np.array([0.95, 0.975, 1.0, 1.025, 1.05]) γ = 2.0 @@ -1043,13 +988,11 @@ Do the same for First, let's enter the parameters: -```{code-cell} ipython3 +```{code-cell} python3 n = 5 -P = jnp.full((n, n), 0.0125) -P = P.at[jnp.arange(n), jnp.arange(n)].set( - P[jnp.arange(n), jnp.arange(n)] + 1 - P.sum(1) - ) -s = jnp.array([0.95, 0.975, 1.0, 1.025, 1.05]) # State values +P = np.full((n, n), 0.0125) +P[range(n), range(n)] += 1 - P.sum(1) +s = np.array([0.95, 0.975, 1.0, 1.025, 1.05]) # State values mc = qe.MarkovChain(P, state_values=s) γ = 2.0 @@ -1061,27 +1004,27 @@ p_s = 150.0 Next, we'll create an instance of `AssetPriceModel` to feed into the functions -```{code-cell} ipython3 -apm = create_ap_model(mc=mc, g=lambda x: x, β=β, γ=γ) +```{code-cell} python3 +apm = AssetPriceModel(β=β, mc=mc, γ=γ, g=lambda x: x) ``` Now we just need to call the relevant functions on the data: -```{code-cell} ipython3 +```{code-cell} python3 tree_price(apm) ``` -```{code-cell} ipython3 +```{code-cell} python3 consol_price(apm, ζ) ``` -```{code-cell} ipython3 +```{code-cell} python3 call_option(apm, ζ, p_s) ``` Let's show the last two functions as a plot -```{code-cell} ipython3 +```{code-cell} python3 fig, ax = plt.subplots() ax.plot(s, consol_price(apm, ζ), label='consol') ax.plot(s, call_option(apm, ζ, p_s), label='call option') @@ -1142,7 +1085,7 @@ Is one higher than the other? Can you give intuition? Here's a suitable function: -```{code-cell} ipython3 +```{code-cell} python3 def finite_horizon_call_option(ap, ζ, p_s, k): """ Computes k period option value. @@ -1152,23 +1095,22 @@ def finite_horizon_call_option(ap, ζ, p_s, k): M = P * ap.g(y)**(- γ) # Make sure that a unique solution exists - test_stability(M, β) + ap.test_stability(M) + # Compute option price p = consol_price(ap, ζ) - def step(i, w): + w = np.zeros(ap.n) + for i in range(k): # Maximize across columns - w = jnp.maximum(β * M @ w, p - p_s) - return w - - w = jax.lax.fori_loop(0, k, step, jnp.zeros(ap.n)) + w = np.maximum(β * M @ w, p - p_s) return w ``` Now let's compute the option values at `k=5` and `k=25` -```{code-cell} ipython3 +```{code-cell} python3 fig, ax = plt.subplots() for k in [5, 25]: w = finite_horizon_call_option(apm, ζ, p_s, k) From ba2fd6475d1e878fa5025d70f0dfbe74fb6d8c19 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 2 Dec 2025 09:10:02 +0900 Subject: [PATCH 12/16] fix bug --- lectures/ge_arrow.md | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index ed74f2e5d..a861f1fc0 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -998,12 +998,11 @@ def compute_rc_model(s, P, ys, s0_idx=0, γ=0.5, β=0.98, T=0): diff = jnp.empty((n, K)) # Loop scatters each agent's state-dependent surplus into the column k. - def body(k, carry): - diff = carry + def body(k, diff): return diff.at[:, k].set(α[k] * y - ys[:, k]) # Applies body sequentially while threading diff. - diff = jax.lax.scan(0, K, body, diff) + diff = jax.lax.fori_loop(0, K, body, diff) ψ = V @ diff From b58cafedc25db2a30bc925a909067f5ed30f21c7 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 2 Dec 2025 09:52:37 +0900 Subject: [PATCH 13/16] fix reference --- lectures/ge_arrow.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index a861f1fc0..f2b1034b8 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -505,7 +505,7 @@ This follows from agent $k$'s first-order necessary conditions. But with the CRRA preferences that we have assumed, individual consumptions vary proportionately with aggregate consumption and therefore with the aggregate endowment. - * This is a consequence of our preference specification implying that *Engel curves* are affine in wealth and therefore satisfy conditions for *Gorman aggregation* +* This is a consequence of our preference specification implying that *Engel curves* are affine in wealth and therefore satisfy conditions for *Gorman aggregation* Thus, @@ -588,7 +588,7 @@ even if he consumes zero goods forevermore. If we have an Inada condition at zero consumption or just impose that consumption be nonnegative, then in a *finite horizon* economy with sequential trading of one-period Arrow securities there is no need to impose natural debt limits. -See the section on a [finite horizon economy](#finite-horizon) below. +See the section on a [finite horizon economy](finite_horizon) below. ``` @@ -710,7 +710,7 @@ $$ J^k = (I - \beta P)^{-1} u(\alpha_k y) , \quad u(c) = \frac{c^{1-\gamma}}{1- where it is understood that $ u(\alpha_k y)$ is a vector. - +(finite_horizon)= ## Finite horizon We now describe a finite-horizon version of the economy that operates for $T+1$ periods From bf55c41770625cc7c35ea250b28b541a6d87beca Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 2 Dec 2025 16:49:57 +0900 Subject: [PATCH 14/16] add definition block --- lectures/ge_arrow.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index f2b1034b8..3698e5a3d 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -206,7 +206,7 @@ where it is understood that $c_t^k = c^k(s_t)$ and $c_{t+1}^k = c^k(s_{t+1})$. - +```{prf:definition} A **recursive competitive equilibrium** is an initial distribution of wealth $\vec a_0$, a set of borrowing limits $\{\bar A^k(s)\}_{k=1}^K$, a pricing kernel $Q(s' | s)$, sets of value functions $\{v^k(a,s)\}_{k=1}^K$, and @@ -230,7 +230,7 @@ $\sum_k \hat a_{t+1}^k(s') = 0$ for all $t$ and $s'$. * The initial financial wealth vector $\vec a_0$ satisfies $\sum_{k=1}^K a_0^k = 0 $. - +``` The third condition asserts that there are zero net aggregate claims in all Markov states. From 253dea657e3714bdcfef5373b242ad4dcda38961 Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 2 Dec 2025 16:55:27 +0900 Subject: [PATCH 15/16] add figure metadata cell --- lectures/ge_arrow.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index 3698e5a3d..18aafad70 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -1311,6 +1311,13 @@ def body(i, carry): ``` ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: | + Wealth distribution under different transition + name: fig_wealth +--- fig, axs = plt.subplots(1, 2, figsize=(12, 4)) for i, αs_seq in enumerate([αs0_seq, αs1_seq]): From be072763aec841b2a1c6b9bd2e4389d7ad8f1efe Mon Sep 17 00:00:00 2001 From: Kenko LI Date: Tue, 2 Dec 2025 17:01:09 +0900 Subject: [PATCH 16/16] refine styling --- lectures/ge_arrow.md | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/lectures/ge_arrow.md b/lectures/ge_arrow.md index 18aafad70..6d42b6b52 100644 --- a/lectures/ge_arrow.md +++ b/lectures/ge_arrow.md @@ -541,21 +541,21 @@ to pose or represent the solution of an individual household's optimum problem. We denote an $K \times 1$ vector of state-dependent values of agents' endowments in Markov state $s$ as $$ -A\left(s\right)=\left[\begin{array}{c} +A\left(s\right)=\begin{bmatrix} A^{1}\left(s\right)\\ \vdots\\ A^{K}\left(s\right) -\end{array}\right], \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] +\end{bmatrix}, \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] $$ and an $n \times 1$ vector-form function of continuation endowment values for each individual $k$ as $$ -A^{k}=\left[\begin{array}{c} +A^{k}=\begin{bmatrix} A^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ A^{k}\left(\bar{s}_{n}\right) -\end{array}\right], \quad k \in \left[1, \ldots, K\right] +\end{bmatrix}, \quad k \in \left[1, \ldots, K\right] $$ $A^k$ of consumer $k$ satisfies @@ -567,11 +567,11 @@ $$ where $$ -y^{k}=\left[\begin{array}{c} +y^{k}=\begin{bmatrix} y^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ y^{k}\left(\bar{s}_{n}\right) -\end{array}\right] \equiv \begin{bmatrix} y^k_1 \cr \vdots \cr y^k_n \end{bmatrix} +\end{bmatrix} \equiv \begin{bmatrix} y^k_1 \cr \vdots \cr y^k_n \end{bmatrix} $$ @@ -602,21 +602,21 @@ trading of a complete set of one-period Arrow securities. We denote an $K \times 1$ vector of state-dependent continuation wealths in Markov state $s$ as $$ -\psi\left(s\right)=\left[\begin{array}{c} +\psi\left(s\right)=\begin{bmatrix} \psi^{1}\left(s\right)\\ \vdots\\ \psi^{K}\left(s\right) -\end{array}\right], \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] +\end{bmatrix}, \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] $$ and an $n \times 1$ vector-form function of continuation wealths for each individual $k\in \left[1, \ldots, K\right]$ as $$ -\psi^{k}=\left[\begin{array}{c} +\psi^{k}=\begin{bmatrix} \psi^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ \psi^{k}\left(\bar{s}_{n}\right) -\end{array}\right] +\end{bmatrix} $$ Continuation wealth $\psi^k$ of consumer $k$ satisfies @@ -628,15 +628,15 @@ $$ (eq:continwealth) where $$ -y^{k}=\left[\begin{array}{c} +y^{k}=\begin{bmatrix} y^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ y^{k}\left(\bar{s}_{n}\right) -\end{array}\right],\quad y=\left[\begin{array}{c} +\end{bmatrix},\quad y=\begin{bmatrix} y\left(\bar{s}_{1}\right)\\ \vdots\\ y\left(\bar{s}_{n}\right) -\end{array}\right] +\end{bmatrix} $$ Note that $\sum_{k=1}^K \psi^k = {0}_{n \times 1}$. @@ -731,21 +731,21 @@ limits borrowing. We denote a $K \times 1$ vector of state-dependent continuation wealths in Markov state $s$ at time $t$ as $$ -\psi_t\left(s\right)=\left[\begin{array}{c} +\psi_t\left(s\right)=\begin{bmatrix} \psi^{1}\left(s\right)\\ \vdots\\ \psi^{K}\left(s\right) -\end{array}\right], \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] +\end{bmatrix}, \quad s \in \left[\bar{s}_1, \ldots, \bar{s}_n\right] $$ and an $n \times 1$ vector of continuation wealths for each individual $k$ as $$ -\psi_t^{k}=\left[\begin{array}{c} +\psi_t^{k}=\begin{bmatrix} \psi_t^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ \psi_t^{k}\left(\bar{s}_{n}\right) -\end{array}\right], \quad k \in \left[1, \ldots, K\right] +\end{bmatrix}, \quad k \in \left[1, \ldots, K\right] $$ @@ -764,15 +764,15 @@ $$ (eq:vv) where $$ -y^{k}=\left[\begin{array}{c} +y^{k}=\begin{bmatrix} y^{k}\left(\bar{s}_{1}\right)\\ \vdots\\ y^{k}\left(\bar{s}_{n}\right) -\end{array}\right],\quad y=\left[\begin{array}{c} +\end{bmatrix},\quad y=\begin{bmatrix} y\left(\bar{s}_{1}\right)\\ \vdots\\ y\left(\bar{s}_{n}\right) -\end{array}\right] +\end{bmatrix} $$ Note that $\sum_{k=1}^K \psi_t^k = {0}_{n \times 1}$ for all $t \in {\bf T}$.