diff --git a/lectures/mccall_fitted_vfi.md b/lectures/mccall_fitted_vfi.md index 42ed0a7f3..9395cab87 100644 --- a/lectures/mccall_fitted_vfi.md +++ b/lectures/mccall_fitted_vfi.md @@ -74,6 +74,7 @@ import jax import jax.numpy as jnp from jax import lax from typing import NamedTuple +from functools import partial import quantecon as qe ``` @@ -81,12 +82,12 @@ import quantecon as qe The model is the same as in the {doc}`discrete case `, with the following features: -- Each period, an unemployed agent receives a wage offer $w$ +- Each period, an unemployed agent receives a wage offer $W_t$ - Wage offers follow a continuous Markov process: $W_t = \exp(X_t)$ where $X_{t+1} = \rho X_t + \nu Z_{t+1}$ - $\{Z_t\}$ is IID and standard normal - Jobs terminate with probability $\alpha$ each period (separation rate) - Unemployed workers receive compensation $c$ per period -- Workers have CRRA utility $u(c) = \frac{c^{1-\gamma} - 1}{1-\gamma}$ +- Workers have CRRA utility $u(x) = \frac{x^{1-\gamma} - 1}{1-\gamma}$ - Future payoffs are discounted by factor $\beta \in (0,1)$ ## The algorithm @@ -132,6 +133,12 @@ $$ where $\psi$ is the standard normal density. +To understand this expression, recall that $W_t = \exp(X_t)$ where $X_{t+1} = \rho X_t + \nu Z_{t+1}$. + +If the current wage is $w = \exp(x)$, then $x = \log(w)$ and the next period's log-wage is $X_{t+1} = \rho \log(w) + \nu Z_{t+1}$. + +Hence the next period's wage is $W_{t+1} = \exp(X_{t+1}) = \exp(\rho \log(w) + \nu Z_{t+1}) = w^\rho \exp(\nu Z_{t+1})$. + Here we are thinking of $v_u$ as a function on all of $\mathbb{R}_+$. @@ -216,7 +223,7 @@ The first step is to build a JAX-compatible structure for the McCall model with The key computational challenge is evaluating the conditional expectation $(Pv_u)(w) = \int v_u(w') p(w, w') dw'$ at each wage grid point. -From the equation above, we have: +Recall that we have: $$ (Pv_u)(w) = \int v_u(w^\rho \exp(\nu z)) \psi(z) dz @@ -227,14 +234,14 @@ where $\psi$ is the standard normal density. We approximate this integral using Monte Carlo integration with draws from the standard normal distribution: $$ -(Pv_u)(w) \approx \frac{1}{N} \sum_{i=1}^N v_u(w^\rho \exp(\nu z_i)) +(Pv_u)(w) \approx \frac{1}{N} \sum_{i=1}^N v_u(w^\rho \exp(\nu Z_i)) $$ We use the same CRRA utility function as in the discrete case: ```{code-cell} ipython3 -def u(c, γ): - return (c**(1 - γ) - 1) / (1 - γ) +def u(x, γ): + return (x**(1 - γ) - 1) / (1 - γ) ``` Here's our model structure using a NamedTuple. @@ -270,9 +277,11 @@ def create_mccall_model( mc = qe.markov.tauchen(grid_size, ρ, ν) w_grid = jnp.exp(jnp.array(mc.state_values)) - return Model(c=c, α=α, β=β, ρ=ρ, ν=ν, γ=γ, w_grid=w_grid, z_draws=z_draws) + return Model(c, α, β, ρ, ν, γ, w_grid, z_draws) ``` +Here is the Bellman operator, where we use Monte Carlo integration to evaluate the expectation. + ```{code-cell} ipython3 def T(model, v): """Update the value function.""" @@ -293,9 +302,9 @@ def T(model, v): Pv = compute_exp_all(w_grid) d = 1 / (1 - β * (1 - α)) - accept = d * (u(w_grid, γ) + α * β * Pv) - reject = u(c, γ) + β * Pv - return jnp.maximum(accept, reject) + v_e = d * (u(w_grid, γ) + α * β * Pv) + continuation_values = u(c, γ) + β * Pv + return jnp.maximum(v_e, continuation_values) ``` Here's the solver: @@ -349,9 +358,9 @@ def get_greedy(v: jnp.ndarray, model: Model) -> jnp.ndarray: Pv = compute_exp_all(w_grid) d = 1 / (1 - β * (1 - α)) - accept = d * (u(w_grid, γ) + α * β * Pv) - reject = u(c, γ) + β * Pv - σ = accept >= reject + v_e = d * (u(w_grid, γ) + α * β * Pv) + continuation_values = u(c, γ) + β * Pv + σ = v_e >= continuation_values return σ ``` @@ -409,26 +418,330 @@ compute_exp_all = jax.vmap(compute_expectation) Pv = compute_exp_all(w_grid) d = 1 / (1 - β * (1 - α)) -accept = d * (u(w_grid, γ) + α * β * Pv) -h_star = u(c, γ) + β * Pv -w_star = get_reservation_wage(σ_star, model) +v_e = d * (u(w_grid, γ) + α * β * Pv) +h = u(c, γ) + β * Pv +w_bar = get_reservation_wage(σ_star, model) ``` Let's plot our results. ```{code-cell} ipython3 fig, ax = plt.subplots(figsize=(9, 5.2)) -ax.plot(w_grid, h_star, linewidth=4, ls="--", alpha=0.4, - label="continuation value") -ax.plot(w_grid, accept, linewidth=4, ls="--", alpha=0.4, - label="stopping value") -ax.plot(w_grid, v_star, "k-", alpha=0.7, label=r"$v_u^*(w)$") +ax.plot(w_grid, h, 'g-', linewidth=2, + label="continuation value function $h$") +ax.plot(w_grid, v_e, 'b-', linewidth=2, + label="employment value function $v_e$") ax.legend(frameon=False) ax.set_xlabel(r"$w$") plt.show() ``` -The exercises ask you to explore the solution and how it changes with parameters. +The reservation wage is at the intersection of the employment value function $v_e$ and the continuation value function $h$. + +## Simulation + +Let's simulate the employment path of a single agent under the optimal policy. + +We need a function to update the agent's state by one period. + +```{code-cell} ipython3 +def update_agent(key, status, wage, model, w_bar): + """ + Updates an agent's employment status and current wage. + + Parameters: + - key: JAX random key + - status: Current employment status (0 or 1) + - wage: Current wage if employed, current offer if unemployed + - model: Model instance + - w_bar: Reservation wage + + """ + c, α, β, ρ, ν, γ, w_grid, z_draws = model + + # Draw new wage offer based on current wage + key1, key2 = jax.random.split(key) + z = jax.random.normal(key1) + new_wage = wage**ρ * jnp.exp(ν * z) + + # Check if separation occurs (for employed workers) + separation_occurs = jax.random.uniform(key2) < α + + # Accept if current wage meets or exceeds reservation wage + accepts = wage >= w_bar + + # If employed: status = 1 if no separation, 0 if separation + # If unemployed: status = 1 if accepts, 0 if rejects + next_status = jnp.where( + status, + 1 - separation_occurs.astype(jnp.int32), # employed path + accepts.astype(jnp.int32) # unemployed path + ) + + # If employed: wage = current if no separation, new if separation + # If unemployed: wage = current if accepts, new if rejects + next_wage = jnp.where( + status, + jnp.where(separation_occurs, new_wage, wage), # employed path + jnp.where(accepts, wage, new_wage) # unemployed path + ) + + return next_status, next_wage +``` + +Here's a function to simulate the employment path of a single agent. + +```{code-cell} ipython3 +def simulate_employment_path( + model: Model, # Model details + w_bar: float, # Reservation wage + T: int = 2_000, # Simulation length + seed: int = 42 # Set seed for simulation + ): + """ + Simulate employment path for T periods starting from unemployment. + + """ + key = jax.random.PRNGKey(seed) + c, α, β, ρ, ν, γ, w_grid, z_draws = model + + # Initial conditions: start unemployed with initial wage draw + status = 0 + key, subkey = jax.random.split(key) + wage = jnp.exp(jax.random.normal(subkey) * ν) + + wage_path = [] + status_path = [] + + for t in range(T): + wage_path.append(wage) + status_path.append(status) + + key, subkey = jax.random.split(key) + status, wage = update_agent( + subkey, status, wage, model, w_bar + ) + + return jnp.array(wage_path), jnp.array(status_path) +``` + +Let's create a comprehensive plot of the employment simulation: + +```{code-cell} ipython3 +model = create_mccall_model() + +# Calculate reservation wage for plotting +v_star = vfi(model) +σ_star = get_greedy(v_star, model) +w_bar = get_reservation_wage(σ_star, model) + +wage_path, employment_status = simulate_employment_path(model, w_bar) + +fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 6)) + +# Plot employment status +ax1.plot(employment_status, 'b-', alpha=0.7, linewidth=1) +ax1.fill_between( + range(len(employment_status)), employment_status, alpha=0.3, color='blue' +) +ax1.set_ylabel('employment status') +ax1.set_title('Employment path (0=unemployed, 1=employed)') +ax1.set_yticks((0, 1)) +ax1.set_ylim(-0.1, 1.1) + +# Plot wage path with reservation wage +ax2.plot(wage_path, 'b-', alpha=0.7, linewidth=1) +ax2.axhline(y=w_bar, color='black', linestyle='--', alpha=0.8, + label=f'Reservation wage: {w_bar:.2f}') +ax2.set_xlabel('time') +ax2.set_ylabel('wage') +ax2.set_title('Wage path (actual and offers)') +ax2.legend() + +# Plot cumulative fraction of time unemployed +unemployed_indicator = (employment_status == 0).astype(int) +cumulative_unemployment = ( + jnp.cumsum(unemployed_indicator) / + jnp.arange(1, len(employment_status) + 1) +) + +ax3.plot(cumulative_unemployment, 'r-', alpha=0.8, linewidth=2) +ax3.axhline(y=jnp.mean(unemployed_indicator), color='black', + linestyle='--', alpha=0.7, + label=f'Final rate: {jnp.mean(unemployed_indicator):.3f}') +ax3.set_xlabel('time') +ax3.set_ylabel('cumulative unemployment rate') +ax3.set_title('Cumulative fraction of time spent unemployed') +ax3.legend() +ax3.set_ylim(0, 1) + +plt.tight_layout() +plt.show() +``` + +The simulation shows the agent cycling between employment and unemployment. + +The agent starts unemployed and receives wage offers according to the Markov process. + +When unemployed, the agent accepts offers that exceed the reservation wage. + +When employed, the agent faces job separation with probability $\alpha$ each period. + +### Cross-Sectional Analysis + +Now let's simulate many agents simultaneously to examine the cross-sectional unemployment rate. + +We first create a vectorized version of `update_agent` to efficiently update all agents in parallel: + +```{code-cell} ipython3 +# Create vectorized version of update_agent +update_agents_vmap = jax.vmap( + update_agent, in_axes=(0, 0, 0, None, None) +) +``` + +Next we define the core simulation function, which uses `lax.fori_loop` to efficiently iterate many agents forward in time: + +```{code-cell} ipython3 +@partial(jax.jit, static_argnums=(3, 4)) +def _simulate_cross_section_compiled( + key: jnp.ndarray, + model: Model, + w_bar: float, + n_agents: int, + T: int + ): + """JIT-compiled core simulation loop using lax.fori_loop. + Returns only the final employment state to save memory.""" + c, α, β, ρ, ν, γ, w_grid, z_draws = model + + # Initialize arrays + key, subkey = jax.random.split(key) + wages = jnp.exp(jax.random.normal(subkey, (n_agents,)) * ν) + status = jnp.zeros(n_agents, dtype=jnp.int32) + + def update(t, loop_state): + key, status, wages = loop_state + + # Shift loop state forwards + key, subkey = jax.random.split(key) + agent_keys = jax.random.split(subkey, n_agents) + + status, wages = update_agents_vmap( + agent_keys, status, wages, model, w_bar + ) + + return key, status, wages + + # Run simulation using fori_loop + initial_loop_state = (key, status, wages) + final_loop_state = lax.fori_loop(0, T, update, initial_loop_state) + + # Return only final employment state + _, final_is_employed, _ = final_loop_state + return final_is_employed + + +def simulate_cross_section( + model: Model, + n_agents: int = 100_000, + T: int = 200, + seed: int = 42 + ) -> float: + """ + Simulate employment paths for many agents and return final unemployment rate. + + Parameters: + - model: Model instance with parameters + - n_agents: Number of agents to simulate + - T: Number of periods to simulate + - seed: Random seed for reproducibility + + Returns: + - unemployment_rate: Fraction of agents unemployed at time T + """ + key = jax.random.PRNGKey(seed) + + # Solve for optimal reservation wage + v_star = vfi(model) + σ_star = get_greedy(v_star, model) + w_bar = get_reservation_wage(σ_star, model) + + # Run JIT-compiled simulation + final_status = _simulate_cross_section_compiled( + key, model, w_bar, n_agents, T + ) + + # Calculate unemployment rate at final period + unemployment_rate = 1 - jnp.mean(final_status) + + return unemployment_rate +``` + +This function generates a histogram showing the distribution of employment status across many agents: + +```{code-cell} ipython3 +def plot_cross_sectional_unemployment(model: Model, t_snapshot: int = 200, + n_agents: int = 20_000): + """ + Generate histogram of cross-sectional unemployment at a specific time. + + Parameters: + - model: Model instance with parameters + - t_snapshot: Time period at which to take the cross-sectional snapshot + - n_agents: Number of agents to simulate + """ + # Get final employment state directly + key = jax.random.PRNGKey(42) + v_star = vfi(model) + σ_star = get_greedy(v_star, model) + w_bar = get_reservation_wage(σ_star, model) + final_status = _simulate_cross_section_compiled( + key, model, w_bar, n_agents, t_snapshot + ) + + # Calculate unemployment rate + unemployment_rate = 1 - jnp.mean(final_status) + + fig, ax = plt.subplots(figsize=(8, 5)) + + # Plot histogram as density (bars sum to 1) + weights = jnp.ones_like(final_status) / len(final_status) + ax.hist(final_status, bins=[-0.5, 0.5, 1.5], + alpha=0.7, color='blue', edgecolor='black', + density=True, weights=weights) + + ax.set_xlabel('employment status (0=unemployed, 1=employed)') + ax.set_ylabel('density') + ax.set_title(f'Cross-sectional distribution at t={t_snapshot}, ' + + f'unemployment rate = {unemployment_rate:.3f}') + ax.set_xticks([0, 1]) + + plt.tight_layout() + plt.show() +``` + +Now let's compare the time-average unemployment rate (from a single agent's long simulation) with the cross-sectional unemployment rate (from many agents at a single point in time). + +```{code-cell} ipython3 +model = create_mccall_model() +cross_sectional_unemp = simulate_cross_section( + model, n_agents=20_000, T=200 +) + +time_avg_unemp = jnp.mean(unemployed_indicator) +print(f"Time-average unemployment rate (single agent): " + f"{time_avg_unemp:.4f}") +print(f"Cross-sectional unemployment rate (at t=200): " + f"{cross_sectional_unemp:.4f}") +print(f"Difference: {abs(time_avg_unemp - cross_sectional_unemp):.4f}") +``` + +Now let's visualize the cross-sectional distribution: + +```{code-cell} ipython3 +plot_cross_sectional_unemployment(model) +``` ## Exercises @@ -489,21 +802,21 @@ We compute the reservation wage for different values of the risk aversion parame ```{code-cell} ipython3 γ_vals = jnp.linspace(1.2, 2.5, 15) -w_star_vec = jnp.empty_like(γ_vals) +w_bar_vec = jnp.empty_like(γ_vals) for i, γ in enumerate(γ_vals): model = create_mccall_model(γ=γ) v_star = vfi(model) σ_star = get_greedy(v_star, model) - w_star = get_reservation_wage(σ_star, model) - w_star_vec = w_star_vec.at[i].set(w_star) + w_bar = get_reservation_wage(σ_star, model) + w_bar_vec = w_bar_vec.at[i].set(w_bar) fig, ax = plt.subplots(figsize=(9, 5.2)) -ax.plot(γ_vals, w_star_vec, linewidth=2, alpha=0.6, +ax.plot(γ_vals, w_bar_vec, linewidth=2, alpha=0.6, label='reservation wage') ax.legend(frameon=False) ax.set_xlabel(r'$\gamma$') -ax.set_ylabel(r'$w^*$') +ax.set_ylabel(r'$\bar{w}$') ax.set_title('Reservation wage as a function of risk aversion') plt.show() ``` diff --git a/lectures/mccall_model.md b/lectures/mccall_model.md index 43469f062..febec370e 100644 --- a/lectures/mccall_model.md +++ b/lectures/mccall_model.md @@ -77,19 +77,19 @@ from quantecon.distributions import BetaBinomial ```{index} single: Models; McCall ``` -An unemployed agent receives in each period a job offer at wage $w_t$. +An unemployed agent receives in each period a job offer at wage $W_t$. In this lecture, we adopt the following simple environment: -* The offer sequence $\{w_t\}_{t \geq 0}$ is IID, with $q(w)$ being the probability of observing wage $w$ in finite set $\mathbb{W}$. -* The agent observes $w_t$ at the start of $t$. -* The agent knows that $\{w_t\}$ is IID with common distribution $q$ and can use this when computing expectations. +* The offer sequence $\{W_t\}_{t \geq 0}$ is IID, with $q(w)$ being the probability of observing wage $w$ in finite set $\mathbb{W}$. +* The agent observes $W_t$ at the start of $t$. +* The agent knows that $\{W_t\}$ is IID with common distribution $q$ and can use this when computing expectations. (In later lectures, we will relax these assumptions.) At time $t$, our agent has two choices: -1. Accept the offer and work permanently at constant wage $w_t$. +1. Accept the offer and work permanently at constant wage $W_t$. 1. Reject the offer, receive unemployment compensation $c$, and reconsider next period. The agent is infinitely lived and aims to maximize the expected discounted @@ -107,7 +107,7 @@ The smaller is $\beta$, the more the agent discounts future earnings relative to The variable $y_t$ is income, equal to -* his/her wage $w_t$ when employed +* his/her wage $W_t$ when employed * unemployment compensation $c$ when unemployed diff --git a/lectures/mccall_model_with_sep_markov.md b/lectures/mccall_model_with_sep_markov.md index fc9f1cb0b..177fa4a28 100644 --- a/lectures/mccall_model_with_sep_markov.md +++ b/lectures/mccall_model_with_sep_markov.md @@ -208,8 +208,8 @@ The optimal policy turns out to be a reservation wage strategy: accept all wages The default utility function is a CRRA utility function ```{code-cell} ipython3 -def u(c, γ): - return (c**(1 - γ) - 1) / (1 - γ) +def u(x, γ): + return (x**(1 - γ) - 1) / (1 - γ) ``` Let's set up a `Model` class to store information needed to solve the model. @@ -406,8 +406,8 @@ This is implemented via `jnp.searchsorted` on the precomputed cumulative sum The function `update_agent` advances the agent's state by one period. -The agent's state is a pair $(s_t, w_t)$, where $s_t$ is employment status (0 if -unemployed, 1 if employed) and $w_t$ is +The agent's state is a pair $(S_t, W_t)$, where $S_t$ is employment status (0 if +unemployed, 1 if employed) and $W_t$ is * their current wage offer, if unemployed, or * their current wage, if employed. @@ -569,10 +569,10 @@ fraction of time an agent spends unemployed over a long time series. We will see that these two values are approximately equal -- if fact they are exactly equal in the limit. -The reason is that the process $(s_t, w_t)$, where +The reason is that the process $(S_t, W_t)$, where -- $s_t$ is the employment status and -- $w_t$ is the wage +- $S_t$ is the employment status and +- $W_t$ is the wage is Markovian, since the next pair depends only on the current pair and iid randomness, and ergodic. @@ -595,7 +595,7 @@ In particular, the fraction of time a single agent spends unemployed (across all wage states) converges to the cross-sectional unemployment rate: $$ - \lim_{T \to \infty} \frac{1}{T} \sum_{t=1}^{T} \mathbb{1}\{s_t = \text{unemployed}\} = \sum_{w=1}^{n} \pi(\text{unemployed}, w) + \lim_{T \to \infty} \frac{1}{T} \sum_{t=1}^{T} \mathbb{1}\{S_t = \text{unemployed}\} = \sum_{w=1}^{n} \pi(\text{unemployed}, w) $$ This holds regardless of initial conditions -- provided that we burn in the diff --git a/lectures/mccall_model_with_separation.md b/lectures/mccall_model_with_separation.md index adb711ab3..245aa8202 100644 --- a/lectures/mccall_model_with_separation.md +++ b/lectures/mccall_model_with_separation.md @@ -90,7 +90,7 @@ introducing a utility function $u$. It satisfies $u'> 0$ and $u'' < 0$. -Wage offers $\{ w_t \}$ are IID with common distribution $q$. +Wage offers $\{ W_t \}$ are IID with common distribution $q$. The set of possible wage values is denoted by $\mathbb W$. @@ -106,9 +106,9 @@ If currently employed at wage $w$, the worker 1. receives utility $u(w)$ from their current wage and 1. is fired with some (small) probability $\alpha$, becoming unemployed next period. -If currently unemployed, the worker receives random wage offer $w_t$ and either accepts or rejects. +If currently unemployed, the worker receives random wage offer $W_t$ and either accepts or rejects. -If he accepts, then he begins work immediately at wage $w_t$. +If he accepts, then he begins work immediately at wage $W_t$. If he rejects, then he receives unemployment compensation $c$. @@ -239,8 +239,8 @@ Let's now implement a solution method based on the two Bellman equations The default utility function is a CRRA utility function ```{code-cell} ipython3 -def u(c, γ): - return (c**(1 - γ) - 1) / (1 - γ) +def u(x, γ): + return (x**(1 - γ) - 1) / (1 - γ) ``` Also, here's a default wage distribution, based around the BetaBinomial diff --git a/lectures/mccall_persist_trans.md b/lectures/mccall_persist_trans.md index d1606ee80..53bd462ef 100644 --- a/lectures/mccall_persist_trans.md +++ b/lectures/mccall_persist_trans.md @@ -68,20 +68,20 @@ from typing import NamedTuple Wages at each point in time are given by $$ -w_t = \exp(z_t) + y_t +W_t = \exp(Z_t) + Y_t $$ where $$ -y_t \sim \exp(\mu + s \zeta_t) +Y_t \sim \exp(\mu + s \zeta_t) \quad \text{and} \quad -z_{t+1} = d + \rho z_t + \sigma \epsilon_{t+1} +Z_{t+1} = d + \rho Z_t + \sigma \epsilon_{t+1} $$ Here $\{ \zeta_t \}$ and $\{ \epsilon_t \}$ are both IID and standard normal. -Here $\{y_t\}$ is a transitory component and $\{z_t\}$ is persistent. +Here $\{Y_t\}$ is a transitory component and $\{Z_t\}$ is persistent. As before, the worker can either @@ -146,7 +146,7 @@ $$ \frac{u(w)}{1-\beta} \geq f^*(z) $$ -For utility, we take $u(c) = \ln(c)$. +For utility, we take $u(x) = \ln(x)$. The reservation wage is the wage where equality holds in the last expression. @@ -327,7 +327,7 @@ at all state values. Next, we study how mean unemployment duration varies with unemployment compensation. -For simplicity, we'll fix the initial state at $z_t = 0$. +For simplicity, we'll fix the initial state at $Z_0 = 0$. ```{code-cell} ipython @jax.jit