diff --git a/code/jl/expo_curve.py b/code/jl/expo_curve.py index 4f56d29..dae2b58 100644 --- a/code/jl/expo_curve.py +++ b/code/jl/expo_curve.py @@ -1,41 +1,18 @@ -# --- -# jupyter: -# jupytext: -# text_representation: -# extension: .py -# format_name: light -# format_version: '1.5' -# jupytext_version: 1.14.6 -# kernelspec: -# display_name: Python 3 -# language: python -# name: python3 -# --- - import numpy as np import matplotlib.pyplot as plt import matplotlib +from scipy.linalg import expm, eigvals matplotlib.rc("text", usetex=True) # allow tex rendering -fontsize=16 -from scipy.linalg import expm, eigvals -# %matplotlib inline +fontsize = 16 -# + A = ((-2, -0.4, 0), (-1.4, -1, 2.2), (0, -2, -0.6)) A = np.array(A) -# - - ev = eigvals(A) - -np.imag(ev) - -np.max(np.real(ev)) - h = 0.01 s0 = 0.01 * np.array((1, 1, 1)) @@ -48,7 +25,6 @@ y.append(b) z.append(c) -# + ax = plt.figure().add_subplot(projection='3d') ax.plot(x, y, z, label='$t \mapsto \mathrm{e}^{t A} u_0$') @@ -67,8 +43,3 @@ plt.savefig("../figures/expo_curve_1.pdf") plt.show() -# - - - - - diff --git a/code/jl/markov_js_with_sep.jl b/code/jl/markov_js_with_sep.jl index 7997201..45cd133 100644 --- a/code/jl/markov_js_with_sep.jl +++ b/code/jl/markov_js_with_sep.jl @@ -117,7 +117,7 @@ function plot_w_stars(; α_vals=LinRange(0.0, 1.0, 10), ax.plot(α_vals, w_star_vec, lw=2, alpha=0.6, label="reservation wage") ax.legend(frameon=false, fontsize=fontsize) ax.set_xlabel(L"\alpha", fontsize=fontsize) - ax.set_xlabel(L"w", fontsize=fontsize) + ax.set_ylabel(L"w", fontsize=fontsize) plt.show() if savefig fig.savefig(figname) diff --git a/code/jl/stoch_dom_finite.jl b/code/jl/stoch_dom_finite.jl index 3e5060d..9a34a6a 100644 --- a/code/jl/stoch_dom_finite.jl +++ b/code/jl/stoch_dom_finite.jl @@ -6,7 +6,7 @@ PyPlot.matplotlib[:rc]("text", usetex=true) # allow tex rendering p, q = 0.75, 0.25 fig, axes = plt.subplots(1, 2, figsize=(10, 5.2)) ax = axes[1] -ax.bar(1:2, (p, 1-p), label=L"\phi") +ax.bar(1:2, (p, 1-p), label=L"\varphi") ax = axes[2] ax.bar(1:2, (q, 1-q), label=L"\psi") diff --git a/code/py/american_option.py b/code/py/american_option.py index 5859cac..091e677 100644 --- a/code/py/american_option.py +++ b/code/py/american_option.py @@ -7,6 +7,7 @@ from quantecon import compute_fixed_point import numpy as np +import matplotlib.pyplot as plt from collections import namedtuple from numba import njit, prange @@ -36,7 +37,7 @@ def create_american_option_model( z_vals, Q = mc.state_values + μ, mc.P w_vals, φ, β = np.array([-s, s]), np.array([0.5, 0.5]), 1 / (1 + r) return Model(t_vals=t_vals, z_vals=z_vals, w_vals=w_vals, Q=Q, - φ=φ, T=T, β=β, K=K) + φ=φ, T=T, β=β, K=K) @njit(parallel=True) @@ -72,12 +73,8 @@ def compute_cvf(model): # Plots - -import matplotlib.pyplot as plt - - def plot_contours(savefig=False, - figname="./figures/american_option_1.pdf"): + figname="../figures_py/american_option_1.png"): model = create_american_option_model() t_vals, z_vals, w_vals, Q, φ, T, β, K = model @@ -110,7 +107,7 @@ def plot_contours(savefig=False, def plot_strike(savefig=False, fontsize=9, - figname="./figures/american_option_2.pdf"): + figname="../figures_py/american_option_2.png"): model = create_american_option_model() t_vals, z_vals, w_vals, Q, φ, T, β, K = model h_star = compute_cvf(model) diff --git a/code/py/ar1_spec_rad.py b/code/py/ar1_spec_rad.py index 37fb9c3..3020a15 100644 --- a/code/py/ar1_spec_rad.py +++ b/code/py/ar1_spec_rad.py @@ -9,6 +9,7 @@ """ import numpy as np +import matplotlib.pyplot as plt from quantecon.markov import tauchen @@ -69,12 +70,9 @@ def compute_mc_spec_rad(n, ρ, σ, μ, m, b): # Plots -import matplotlib.pyplot as plt - - def plot_beta_sim(T=80, savefig=True, - figname="./figures/ar1_spec_rad.png"): + figname="../figures_py/ar1_spec_rad.png"): β_vals = np.zeros(T) Z = 1 for t in range(T): @@ -92,5 +90,3 @@ def plot_beta_sim(T=80, if savefig: fig.savefig(figname) plt.show() - -plot_beta_sim() diff --git a/code/py/bellman_envelope.py b/code/py/bellman_envelope.py new file mode 100644 index 0000000..5f01e8c --- /dev/null +++ b/code/py/bellman_envelope.py @@ -0,0 +1,56 @@ +import numpy as np +import matplotlib.pyplot as plt + + +xmin = -0.5 +xmax = 2.0 + +xgrid = np.linspace(xmin, xmax, 1000) + +a1, b1 = 0.15, 0.5 # first T_σ +a2, b2 = 0.5, 0.4 # second T_σ +a3, b3 = 0.75, 0.2 # third T_σ + +v1 = b1/(1-a1) +v2 = b2/(1-a2) +v3 = b3/(1-a3) + +T1 = a1 * xgrid + b1 +T2 = a2 * xgrid + b2 +T3 = a3 * xgrid + b3 +T = np.maximum.reduce([T1, T2, T3]) + +fig, ax = plt.subplots(figsize=(6, 5)) +for spine in ["left", "bottom"]: + ax.spines[spine].set_position("zero") + +for spine in ["right", "top"]: + ax.spines[spine].set_color("none") + +ax.plot(xgrid, T1, "k-", lw=1) +ax.plot(xgrid, T2, "k-", lw=1) +ax.plot(xgrid, T3, "k-", lw=1) + +ax.plot(xgrid, T, lw=6, alpha=0.3, color="blue", + label=r"$T = \bigvee_{\sigma \in \Sigma} T_\sigma$") + + +ax.text(2.1, 0.6, r"$T_{\sigma'}$") +ax.text(2.1, 1.4, r"$T_{\sigma''}$") +ax.text(2.1, 1.9, r"$T_{\sigma'''}$") + +ax.legend(frameon=False, loc="upper center") + + +ax.set_xlim(xmin, xmax+0.5) +ax.set_ylim(-0.2, 2.5) +ax.text(2.4, -0.15, r"$v$") + +ax.set_xticks([]) +ax.set_yticks([]) + +plt.show() + +file_name = "../figures_py/bellman_envelope.png" +fig.savefig(file_name) + diff --git a/code/py/binom_stoch_dom.py b/code/py/binom_stoch_dom.py new file mode 100644 index 0000000..f9579b4 --- /dev/null +++ b/code/py/binom_stoch_dom.py @@ -0,0 +1,25 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import binom + + +n, m, p = 10, 18, 0.5 + +phi = binom(n, p) +psi = binom(m, p) + +x = np.arange(0, m+1) + +fig, ax = plt.subplots(figsize=(9, 5.2)) +lb = r"$\phi = B({}, {})$".format(n, p) +ax.plot(x, np.concatenate((phi.pmf(np.arange(0, n+1)), + np.zeros(m-n))), "-o", alpha=0.6, label=lb) +lb = r"$\psi = B({}, {})$".format(m, p) +ax.plot(x, psi.pmf(x), "-o", alpha=0.6, label=lb) + +ax.legend(fontsize=16, frameon=False) + +plt.show() + +file_name = "../figures_py/binom_stoch_dom.png" +fig.savefig(file_name) diff --git a/code/py/concave_map_fp.py b/code/py/concave_map_fp.py new file mode 100644 index 0000000..7db96ac --- /dev/null +++ b/code/py/concave_map_fp.py @@ -0,0 +1,42 @@ +import numpy as np +import matplotlib.pyplot as plt + +x0 = 0.25 +xmin, xmax = 0, 3 +fs = 18 + +x_grid = np.linspace(xmin, xmax, 1200) + +def g(x): + return 1 + 0.5 * x**0.5 + +xstar = 1.64 + +fig, ax = plt.subplots(figsize=(10, 5.5)) + +# Plot the functions +lb = r"$g$" +ax.plot(x_grid, g(x_grid), lw=2, alpha=0.6, label=lb) +ax.plot(x_grid, x_grid, "k--", lw=1, alpha=0.7, label=r"$45$") + +# Show and annotate the fixed point +fps = (xstar,) +ax.plot(fps, fps, "go", ms=10, alpha=0.6) +ax.set_xlabel(r"$x$", fontsize=fs) +ax.annotate(r"$x^*$", + xy=(xstar, xstar), + xycoords="data", + xytext=(-20, 20), + textcoords="offset points", + fontsize=fs) + +ax.legend(loc="upper left", frameon=False, fontsize=fs) +ax.set_xticks((0, 1, 2, 3)) +ax.set_yticks((0, 1, 2, 3)) +ax.set_ylim(0, 3) +ax.set_xlim(0, 3) + +plt.show() + +file_name = "../figures_py/concave_map_fp.png" +fig.savefig(file_name) diff --git a/code/py/cont_time_js.py b/code/py/cont_time_js.py new file mode 100644 index 0000000..86c3927 --- /dev/null +++ b/code/py/cont_time_js.py @@ -0,0 +1,190 @@ +""" + +Continuous time job search. + +The job status is + + s = 0 for unemployed and s = 1 for employed + +The policy function has the form + + σ[j] = optimal choice in when s = 0 + +We use σ[j] = 0 for reject and 1 for accept + +""" + +import numpy as np +import matplotlib.pyplot as plt +from quantecon.markov import tauchen +from collections import namedtuple +from numba import njit, prange +import time + +# Model +Model = namedtuple("Model", ("n", "w_vals", "P", "Q", "δ", "κ", "c", "α")) + +def create_js_model(α=0.1, κ=1.0, δ=0.1, n=100, ρ=0.9, ν=0.2, c=1.0): + mc = tauchen(n, ρ, ν) + w_vals, P = np.exp(mc.state_values), mc.P + + def Π(s, j, a, s_prime, j_prime): + # a -= 1 + if s == 0 and s_prime == 0: + return P[j, j_prime] * (1 - a) + elif s == 0 and s_prime == 1: + return P[j, j_prime] * a + elif s == 1 and s_prime == 0: + return P[j, j_prime] + else: + return 0.0 + + Q = np.zeros((2, n, 2, 2, n)) + for s, j, a, s_prime, j_prime in np.ndindex(Q.shape): + λ = κ if s == 0 else α + Q[s, j, a, s_prime, j_prime] = λ * (Π(s, j, a, s_prime, j_prime) + - (s == s_prime and j == j_prime)) + + return Model(n=n, w_vals=w_vals, P=P, Q=Q, δ=δ, + κ=κ, c=c, α=α) + +@njit +def B(s, j, a, v, model): + n, w_vals, P, Q, δ, κ, c, α = model + r = c if s == 0 else w_vals[j] + continuation_value = 0 + for s_prime, j_prime in np.ndindex(2, n): + continuation_value += v[s_prime, j_prime] * Q[s, j, a, s_prime, j_prime] + return r + continuation_value + + +@njit(parallel=True) +def get_greedy(v, model): + n = model.n + σ = np.zeros(n, dtype=np.int32) + for j in prange(n): + σ[j] = np.argmax(np.array([B(0, j, a, v, model) + for a in range(2)])) + return σ + + +@njit(parallel=True) +def get_value(σ, model): + n, w_vals, P, Q, δ, κ, c, α = model + A = np.zeros((2 * n, 2 * n)) + r_σ = np.zeros(2 * n) + for s, j in np.ndindex(2, n): + r_σ[s * n + j] = c if s == 0 else w_vals[j] + + for s, j, s_prime, j_prime in np.ndindex(2, n, 2, n): + A[s * n + j, s_prime * n + j_prime] = δ * (s == s_prime and j == j_prime) - Q[s, j, σ[j], s_prime, j_prime] + + v_σ = np.linalg.solve(A, r_σ) + return v_σ.reshape(2, n) + + +@njit +def policy_iteration(v_init, model, tolerance=1e-9, max_iter=1000, verbose=False): + v = v_init + error = tolerance + 1 + k = 1 + while error > tolerance and k < max_iter: + last_v = v + σ = get_greedy(v, model) + v = get_value(σ, model) + error = np.max(np.abs(v - last_v)) + if verbose: + print(f"Completed iteration {k}.") + k += 1 + return v, get_greedy(v, model) + + +def plot_policy(savefig=False, figname="../figures_py/cont_time_js_pol.png"): + model = create_js_model() + n, w_vals = model.n, model.w_vals + v_init = np.ones((2, n)) + start = time.time() + v_star, σ_star = policy_iteration(v_init, model, verbose=True) + end = time.time() + print(f"Elapsed time: {end - start:.2f} seconds.") + + fig, ax = plt.subplots(figsize=(9, 5)) + ax.plot(w_vals, σ_star) + ax.set_xlabel("wage offer", fontsize=14) + ax.set_yticks((0, 1)) + ax.set_ylabel("action (reject/accept)", fontsize=14) + if savefig: + plt.savefig(figname) + plt.show() + + +def plot_reswage(savefig=False, figname="../figures_py/cont_time_js_res.png"): + α_vals = np.linspace(0.05, 1.0, 100) + res_wages_alpha = [] + for α in α_vals: + model = create_js_model(α=α) + n, w_vals = model.n, model.w_vals + v_init = np.ones((2, n)) + v_star, σ_star = policy_iteration(v_init, model) + w_idx = np.searchsorted(σ_star, 1) + w_bar = w_vals[w_idx] + res_wages_alpha.append(w_bar) + + κ_vals = np.linspace(0.5, 1.5, 100) + res_wages_kappa = [] + for κ in κ_vals: + model = create_js_model(κ=κ) + n, w_vals = model.n, model.w_vals + v_init = np.ones((2, n)) + v_star, σ_star = policy_iteration(v_init, model) + w_idx = np.searchsorted(σ_star, 1) + w_bar = w_vals[w_idx] + res_wages_kappa.append(w_bar) + + δ_vals = np.linspace(0.05, 1.0, 100) + res_wages_delta = [] + for δ in δ_vals: + model = create_js_model(δ=δ) + v_init = np.ones((2, n)) + v_star, σ_star = policy_iteration(v_init, model) + w_idx = np.searchsorted(σ_star, 1) + w_bar = w_vals[w_idx] + res_wages_delta.append(w_bar) + + c_vals = np.linspace(0.5, 1.5, 100) + res_wages_c = [] + for c in c_vals: + model = create_js_model(c=c) + v_init = np.ones((2, n)) + v_star, σ_star = policy_iteration(v_init, model) + w_idx = np.searchsorted(σ_star, 1) + w_bar = w_vals[w_idx] + res_wages_c.append(w_bar) + + fig, axes = plt.subplots(2, 2, figsize=(9, 5)) + + ax = axes[0, 0] + ax.plot(α_vals, res_wages_alpha) + ax.set_xlabel("separation rate", fontsize=14) + ax.set_ylabel("res. wage", fontsize=14) + + ax = axes[0, 1] + ax.plot(κ_vals, res_wages_kappa) + ax.set_xlabel("offer rate", fontsize=14) + ax.set_ylabel("res. wage", fontsize=14) + + ax = axes[1, 0] + ax.plot(δ_vals, res_wages_delta) + ax.set_xlabel("discount rate", fontsize=14) + ax.set_ylabel("res. wage", fontsize=14) + + ax = axes[1, 1] + ax.plot(c_vals, res_wages_c) + ax.set_xlabel("unempl. compensation", fontsize=14) + ax.set_ylabel("res. wage", fontsize=14) + + fig.tight_layout() + + if savefig: + plt.savefig(figname) + plt.show() diff --git a/code/py/discount_spec_rad.py b/code/py/discount_spec_rad.py new file mode 100644 index 0000000..55a86d3 --- /dev/null +++ b/code/py/discount_spec_rad.py @@ -0,0 +1,45 @@ +import numpy as np +from numba import njit +import matplotlib.pyplot as plt +from quantecon import tauchen + +fontsize = 18 + +# Spectral radius +@njit +def ρ(A): + return np.max(np.abs(np.linalg.eigvals(A))) + +def plot_contours(savefig=False, figname="../figures_py/discount_spec_rad.png"): + fig, ax = plt.subplots(figsize=(8, 4.2)) + grid_size = 50 + n = 6 + a_vals = np.linspace(0.1, 0.95, grid_size) + s_vals = np.linspace(0.1, 0.32, grid_size) + μ = 0.96 + + R = np.zeros((grid_size, grid_size)) + L = np.zeros((n, n)) + + for i_a, a in enumerate(a_vals): + for i_s, s in enumerate(s_vals): + mc = tauchen(n, a, np.sqrt(1 - a**2) * s, (1 - a) * μ) + z_vals, Q = mc.state_values, mc.P + β_vals = z_vals # np.maximum(z_vals, 0) + L = β_vals * Q + R[i_a, i_s] = ρ(L) + + cs1 = ax.contourf(a_vals, s_vals, R.T, alpha=0.5) + ctr1 = ax.contour(a_vals, s_vals, R.T, levels=[1.0]) + plt.clabel(ctr1, inline=1, fontsize=13) + plt.colorbar(cs1, ax=ax) + + ax.set_xlabel(r"$a$", fontsize=fontsize) + ax.set_ylabel(r"$s$", fontsize=fontsize) + + if savefig: + fig.savefig(figname) + plt.show() + + +plot_contours(True) diff --git a/code/py/ez_dp_code.py b/code/py/ez_dp_code.py new file mode 100644 index 0000000..2f7044e --- /dev/null +++ b/code/py/ez_dp_code.py @@ -0,0 +1,85 @@ +import numpy as np +from numba import njit, prange +from ez_model import B, B2 + + +@njit(parallel=True) +def T_σ(v, σ, model): + w_n, e_n = v.shape + v_new = np.zeros_like(v) + for i in prange(w_n): + for j in range(e_n): + v_new[i, j] = B(i, j, σ[i, j], v, model) + return v_new + + +@njit(parallel=True) +def T2_σ(h, σ, model): + w_n = len(h) + h_new = np.zeros_like(h) + for i in prange(w_n): + h_new[i] = B2(i, σ[i], h, model) + return h_new + + +@njit(parallel=True) +def get_greedy(v, model): + w_n, e_n = v.shape + σ = np.zeros_like(v, dtype=np.int32) + for i in prange(w_n): + for j in range(e_n): + B_values = np.array([B(i, j, k, v, model) for k in range(w_n)]) + σ[i, j] = np.argmax(B_values) + return σ + + +@njit(parallel=True) +def get_greedy2(h, model): + w_n = len(h) + σ = np.zeros(w_n, dtype=np.int32) + for i in prange(w_n): + B_values = np.array([B2(i, k, h, model) for k in range(w_n)]) + σ[i] = np.argmax(B_values) + return σ + + +@njit +def get_value(v_init, σ, m, model): + v = v_init + for _ in range(m): + v = T_σ(v, σ, model) + return v + +@njit +def get_value2(h_init, σ, m, model): + h = h_init + for _ in range(m): + h = T2_σ(h, σ, model) + return h + + +def optimistic_policy_iteration(v_init, model, tolerance=1e-9, max_iter=1000, m=100): + v = v_init + error = tolerance + 1 + k = 1 + while error > tolerance and k < max_iter: + last_v = v + σ = get_greedy(v, model) + v = get_value(v, σ, m, model) + error = np.max(np.abs(v - last_v)) + print(f"Completed iteration {k} with error {error}.") + k += 1 + return v, get_greedy(v, model) + +def optimistic_policy_iteration2(h_init, model, tolerance=1e-9, max_iter=1000, m=100): + h = h_init + error = tolerance + 1 + k = 1 + while error > tolerance and k < max_iter: + last_h = h + σ = get_greedy2(h, model) + h = get_value2(h, σ, m, model) + error = np.max(np.abs(h - last_h)) + print(f"Completed iteration {k} with error {error}.") + k += 1 + return h, get_greedy2(h, model) diff --git a/code/py/ez_f_shapes.py b/code/py/ez_f_shapes.py new file mode 100644 index 0000000..17d08d3 --- /dev/null +++ b/code/py/ez_f_shapes.py @@ -0,0 +1,26 @@ +import numpy as np +import matplotlib.pyplot as plt +from numba import njit + + +@njit +def F(w, r=1, β=0.5, θ=5): + return (r + β * w**(1/θ))**θ + +w_grid = np.linspace(0.1, 2.0, 200) + +fig, axes = plt.subplots(2, 2) + +θ_vals = [-2, -0.5, 0.5, 2] + +for i in range(2): + for j in range(2): + θ = θ_vals[i*2 + j] + f = lambda w: F(w, θ=θ) + axes[i, j].plot(w_grid, w_grid, "k--", alpha=0.6, label=r"$45$") + axes[i, j].plot(w_grid, f(w_grid), label=r"$U$") + axes[i, j].legend() + axes[i, j].set_title(f"$\\theta = {θ}$") + +fig.tight_layout() +plt.show() diff --git a/code/py/ez_model.py b/code/py/ez_model.py new file mode 100644 index 0000000..9302762 --- /dev/null +++ b/code/py/ez_model.py @@ -0,0 +1,65 @@ +import numpy as np +from scipy.stats import binom +from numba import njit, prange +from collections import namedtuple + +Model = namedtuple('Model', ['α', 'β', 'γ', 'θ', 'φ', 'e_grid', 'w_grid']) + +def create_ez_model(ψ=1.97, + β=0.96, + γ=-7.89, + n=80, + p=0.5, + e_max=0.5, + w_size=50, + w_max=2): + α = 1 - 1/ψ + θ = γ / α + b = binom(n - 1, p) + φ = b.pmf(np.arange(n)) + e_grid = np.linspace(1e-5, e_max, n) + w_grid = np.linspace(0, w_max, w_size) + return Model(α=α, β=β, γ=γ, θ=θ, φ=φ, + e_grid=e_grid, w_grid=w_grid) + +@njit +def B(i, j, k, v, model): + """ + Action-value aggregator for the original model. + """ + α, β, γ, θ, φ, e_grid, w_grid = model + w, e, s = w_grid[i], e_grid[j], w_grid[k] + if s<= w: + Rv = np.dot(v[k, :]**γ, φ)**(1/γ) + return ((w - s + e)**α + β * Rv**α)**(1/α) + return -np.inf + +@njit +def B2(i, k, h, model): + α, β, γ, θ, φ, e_grid, w_grid = model + w, s = w_grid[i], w_grid[k] + if s <= w: + Ge = ((w - s + e_grid)**α + β * h[k]**α)**(1/α) + return np.dot(Ge**γ, φ)**(1/γ) + return -np.inf + + +@njit +def G_obj(i, j, k, h, model): + α, β, γ, θ, φ, e_grid, w_grid = model + w, e, s = w_grid[i], e_grid[j], w_grid[k] + if s <= w: + return ((w - s + e)**α + β * h[k]**α)**(1/α) + return -np.inf + + +@njit +def G_max(h, model): + w_n, e_n = len(model.w_grid), len(model.e_grid) + σ_star_mod = np.zeros((w_n, e_n), dtype=np.int32) + for i in prange(w_n): + for j in range(e_n): + G_values = np.array([G_obj(i, j, k, h, model) + for k in range(w_n)]) + σ_star_mod[i, j] = np.argmax(G_values) + return σ_star_mod diff --git a/code/py/ez_noncontraction.py b/code/py/ez_noncontraction.py new file mode 100644 index 0000000..3db9a20 --- /dev/null +++ b/code/py/ez_noncontraction.py @@ -0,0 +1,23 @@ +import numpy as np +import matplotlib.pyplot as plt + +def F(w, r=0.5, β=0.5, θ=5): + return (r + β * w**(1/θ))**θ + +w_grid = np.linspace(0.001, 2.0, 200) + +def plot_F(savefig=False, figname="../figures_py/ez_noncontraction.png", fs=16): + fig, ax = plt.subplots(figsize=(9, 5.2)) + f = lambda w: F(w, θ=-10) + ax.plot(w_grid, w_grid, "k--", alpha=0.6, label=r"$45$") + ax.plot(w_grid, f(w_grid), label=r"$\hat{K} = F$") + ax.set_xticks((0, 1, 2)) + ax.set_yticks((0, 1, 2)) + ax.legend(fontsize=fs, frameon=False) + + plt.show() + + if savefig: + fig.savefig(figname) + +plot_F() diff --git a/code/py/ez_plot_functions.py b/code/py/ez_plot_functions.py new file mode 100644 index 0000000..a75d75d --- /dev/null +++ b/code/py/ez_plot_functions.py @@ -0,0 +1,35 @@ +import numpy as np +import matplotlib.pyplot as plt + +def plot_policy(σ, model, title="", savefig=False, figname="../figures_py/ez_policies.png", fontsize=16): + w_grid = model.w_grid + fig, ax = plt.subplots(figsize=(9, 5.2)) + ax.plot(w_grid, w_grid, "k--", label=r"$45$") + ax.plot(w_grid, w_grid[σ[:, 0]], label=r"$\sigma^*(\cdot, e_1)$") + ax.plot(w_grid, w_grid[σ[:, -1]], label=r"$\sigma^*(\cdot, e_N)$") + ax.legend(fontsize=fontsize) + ax.set_xlabel("$w$", fontsize=fontsize) + ax.set_ylabel("$\sigma^*$", fontsize=fontsize) + ax.set_title(title, fontsize=fontsize) + if savefig: + plt.savefig(figname) + plt.show() + +def plot_value_orig(v, model, fontsize=16): + w_grid = model.w_grid + fig, ax = plt.subplots(figsize=(9, 5.2)) + ax.plot(w_grid, v[:, 0], label=r"$v^*(\cdot, e_1)$") + ax.plot(w_grid, v[:, -1], label=r"$v^*(\cdot, e_N)$") + ax.legend(fontsize=fontsize) + ax.set_xlabel("$w$", fontsize=fontsize) + ax.set_ylabel("$v^*$", fontsize=fontsize) + plt.show() + +def plot_value_mod(h, model, fontsize=16): + w_grid = model.w_grid + fig, ax = plt.subplots(figsize=(9, 5.2)) + ax.plot(w_grid, h, label=r"$h^*$") + ax.legend(fontsize=fontsize) + ax.set_xlabel("$w$", fontsize=fontsize) + ax.set_ylabel("$h^*$", fontsize=fontsize) + plt.show() diff --git a/code/py/ez_policy_plot.py b/code/py/ez_policy_plot.py new file mode 100644 index 0000000..aeedd56 --- /dev/null +++ b/code/py/ez_policy_plot.py @@ -0,0 +1,13 @@ +import numpy as np + +from ez_model import create_ez_model, G_max +from ez_dp_code import optimistic_policy_iteration2 +from ez_plot_functions import plot_policy + +model = create_ez_model() + +h_init = np.ones(len(model.w_grid)) +h_star, _ = optimistic_policy_iteration2(h_init, model) + +σ_star_mod = G_max(h_star, model) +plot_policy(σ_star_mod, model, title="optimal savings", savefig=True) diff --git a/code/py/ez_sub_test.py b/code/py/ez_sub_test.py new file mode 100644 index 0000000..212c26a --- /dev/null +++ b/code/py/ez_sub_test.py @@ -0,0 +1,34 @@ +""" +Quick test and plots + +""" + +import numpy as np +import time + +from ez_model import create_ez_model, G_max +from ez_dp_code import optimistic_policy_iteration2, optimistic_policy_iteration +from ez_plot_functions import plot_policy + + + +model = create_ez_model() + +print("Solving unmodified model.") +v_init = np.ones((len(model.w_grid), len(model.e_grid))) +in_time = time.time() +v_star, σ_star = optimistic_policy_iteration(v_init, model) +out_time = time.time() +print("Time to solve unmodified model: ", out_time - in_time) +print("Solving modified model.") + +h_init = np.ones(len(model.w_grid)) +in_time = time.time() +h_star, _ = optimistic_policy_iteration2(h_init, model) +out_time = time.time() +print("Time to solve modified model: ", out_time - in_time) + +σ_star_mod = G_max(h_star, model) + +plot_policy(σ_star, model, title="original") +plot_policy(σ_star_mod, model, title="transformed") diff --git a/code/py/ez_timings.py b/code/py/ez_timings.py new file mode 100644 index 0000000..b8e3fa3 --- /dev/null +++ b/code/py/ez_timings.py @@ -0,0 +1,43 @@ +import numpy as np +import matplotlib.pyplot as plt +import time + +from ez_model import create_ez_model +from ez_dp_code import optimistic_policy_iteration2, optimistic_policy_iteration + +n_vals = np.arange(20, 101, 10) +β_vals = [0.96, 0.98] +gains = np.zeros((len(β_vals), len(n_vals))) + +for β_i, β in enumerate(β_vals): + for n_i, n in enumerate(n_vals): + model = create_ez_model(n=n, β=β) + v_init = np.ones((len(model.w_grid), len(model.e_grid))) + + in_time = time.time() + optimistic_policy_iteration(v_init, model) + unmod_time = time.time() - in_time + + h_init = np.ones(len(model.w_grid)) + in_time = time.time() + optimistic_policy_iteration2(h_init, model) + mod_time = time.time() - in_time + + gains[β_i, n_i] = unmod_time / mod_time + +def plot_function(savefig=False, figname="../figures_py/ez_rel_timing.png"): + fig, ax = plt.subplots(figsize=(9, 5)) + global gains, n_vals, β_vals + for β_i, β in enumerate(β_vals): + label = f"speed gain with $\\beta$ = {β}" + ax.plot(n_vals, gains[β_i, :], "-o", label=label) + + ax.legend(loc="lower right", fontsize=16) + ax.set_xticks(n_vals) + ax.set_xlabel("size of $\\mathsf{E}$", fontsize=16) + ax.set_ylabel("Speed Gain", fontsize=16) + if savefig: + plt.savefig(figname) + plt.show() + +plot_function(True) diff --git a/code/py/ez_utility.py b/code/py/ez_utility.py index 9a6a28b..8b14849 100644 --- a/code/py/ez_utility.py +++ b/code/py/ez_utility.py @@ -52,7 +52,7 @@ def compute_ez_utility(model): def plot_convergence(savefig=False, num_iter=100, - figname="./figures/ez_utility_c.pdf"): + figname="../figures_py/ez_utility_c.png"): fig, ax = plt.subplots(figsize=(10, 5.2)) model = create_ez_utility_model() @@ -82,7 +82,7 @@ def plot_convergence(savefig=False, def plot_v(savefig=False, - figname="./figures/ez_utility_1.pdf"): + figname="../figures_py/ez_utility_1.png"): fig, ax = plt.subplots(figsize=(10, 5.2)) model = create_ez_utility_model() @@ -99,7 +99,7 @@ def plot_v(savefig=False, def vary_gamma(gamma_vals=[1.0, -8.0], savefig=False, - figname="./figures/ez_utility_2.pdf"): + figname="../figures_py/ez_utility_2.png"): fig, ax = plt.subplots(figsize=(10, 5.2)) @@ -119,7 +119,7 @@ def vary_gamma(gamma_vals=[1.0, -8.0], def vary_alpha(alpha_vals=[0.5, 0.6], savefig=False, - figname="./figures/ez_utility_3.pdf"): + figname="../figures_py/ez_utility_3.png"): fig, ax = plt.subplots(figsize=(10, 5.2)) @@ -135,3 +135,8 @@ def vary_alpha(alpha_vals=[0.5, 0.6], plt.show() if savefig: fig.savefig(figname) + +plot_convergence(savefig=True) +plot_v(savefig=True) +vary_gamma(savefig=True) +vary_alpha(savefig=True) diff --git a/code/py/finite_lq.py b/code/py/finite_lq.py index 4039d9c..acc7361 100644 --- a/code/py/finite_lq.py +++ b/code/py/finite_lq.py @@ -157,7 +157,7 @@ def optimistic_policy_iteration(model, tol=1e-5, m=100): import matplotlib.pyplot as plt -def plot_policy(savefig=False, figname="./figures/finite_lq_0.pdf"): +def plot_policy(savefig=False, figname="../figures_py/finite_lq_0.png"): model = create_investment_model() β, a_0, a_1, γ, c, y_grid, z_grid, Q = model σ_star = optimistic_policy_iteration(model) @@ -171,7 +171,7 @@ def plot_policy(savefig=False, figname="./figures/finite_lq_0.pdf"): fig.savefig(figname) -def plot_sim(savefig=False, figname="./figures/finite_lq_1.pdf"): +def plot_sim(savefig=False, figname="../figures_py/finite_lq_1.png"): ts_length = 200 fig, axes = plt.subplots(4, 1, figsize=(9, 11.2)) @@ -209,16 +209,14 @@ def plot_sim(savefig=False, figname="./figures/finite_lq_1.pdf"): def plot_timing(m_vals=np.arange(1, 601, 10), savefig=False, - figname="./figures/finite_lq_time.pdf" + figname="../figures_py/finite_lq_time.png" ): - # NOTE: Uncomment the following lines in this function to - # include Policy iteration plot model = create_investment_model() - # print("Running Howard policy iteration.") - # t1 = time.time() - # σ_pi = policy_iteration(model) - # pi_time = time.time() - t1 - # print(f"PI completed in {pi_time} seconds.") + print("Running Howard policy iteration.") + t1 = time.time() + σ_pi = policy_iteration(model) + pi_time = time.time() - t1 + print(f"PI completed in {pi_time} seconds.") print("Running value function iteration.") t1 = time.time() σ_vfi = value_iteration(model) @@ -236,8 +234,8 @@ def plot_timing(m_vals=np.arange(1, 601, 10), fig, ax = plt.subplots(figsize=(9, 5.2)) ax.plot(m_vals, [vfi_time]*len(m_vals), linewidth=2, label="value function iteration") - # ax.plot(m_vals, [pi_time]*len(m_vals), - # linewidth=2, label="Howard policy iteration") + ax.plot(m_vals, [pi_time]*len(m_vals), + linewidth=2, label="Howard policy iteration") ax.plot(m_vals, opi_times, linewidth=2, label="optimistic policy iteration") ax.legend(frameon=False) ax.set_xlabel(r"$m$") @@ -246,3 +244,7 @@ def plot_timing(m_vals=np.arange(1, 601, 10), if savefig: fig.savefig(figname) return (vfi_time, opi_times) + +plot_policy(savefig=True) +plot_sim(savefig=True) +plot_timing(savefig=True) diff --git a/code/py/finite_opt_saving_2.py b/code/py/finite_opt_saving_2.py index 72ae431..aa3a05a 100644 --- a/code/py/finite_opt_saving_2.py +++ b/code/py/finite_opt_saving_2.py @@ -118,7 +118,7 @@ def plot_timing(m_vals=np.arange(1, 601, 10), ax.set_ylabel("time") plt.show() if savefig: - fig.savefig("./figures/finite_opt_saving_2_1.png") + fig.savefig("../figures_py/finite_opt_saving_2_1.png") return (pi_time, vfi_time, opi_times) @@ -140,7 +140,7 @@ def plot_policy(method="pi", savefig=False): plt.title(f"Method: {method}") plt.show() if savefig: - fig.savefig(f"./figures/finite_opt_saving_2_2_{method}.png") + fig.savefig(f"../figures_py/finite_opt_saving_2_2_{method}.png") def plot_time_series(m=2_000, savefig=False): @@ -151,7 +151,7 @@ def plot_time_series(m=2_000, savefig=False): ax.legend() plt.show() if savefig: - fig.savefig("./figures/finite_opt_saving_ts.pdf") + fig.savefig("../figures_py/finite_opt_saving_ts.png") def plot_histogram(m=1_000_000, savefig=False): @@ -165,7 +165,7 @@ def plot_histogram(m=1_000_000, savefig=False): plt.show() if savefig: - fig.savefig("./figures/finite_opt_saving_hist.pdf") + fig.savefig("../figures_py/finite_opt_saving_hist.png") def plot_lorenz(m=1_000_000, savefig=False): @@ -180,4 +180,10 @@ def plot_lorenz(m=1_000_000, savefig=False): plt.show() if savefig: - fig.savefig("./figures/finite_opt_saving_lorenz.pdf") + fig.savefig("../figures_py/finite_opt_saving_lorenz.png") + +plot_timing(savefig=True) +plot_histogram(savefig=True) +plot_lorenz(savefig=True) +plot_policy(savefig=True) +plot_time_series(savefig=True) diff --git a/code/py/firm_exit.py b/code/py/firm_exit.py index 38730c1..1d2c577 100644 --- a/code/py/firm_exit.py +++ b/code/py/firm_exit.py @@ -68,7 +68,7 @@ def vfi(model): def plot_val(savefig=False, - figname="./figures/firm_exit_1.pdf"): + figname="../figures_py/firm_exit_1.png"): fig, ax = plt.subplots(figsize=(9, 5.2)) @@ -91,7 +91,7 @@ def plot_val(savefig=False, def plot_comparison(savefig=False, - figname="./figures/firm_exit_2.pdf"): + figname="../figures_py/firm_exit_2.png"): fig, ax = plt.subplots(figsize=(9, 5.2)) @@ -110,3 +110,6 @@ def plot_comparison(savefig=False, plt.show() if savefig: fig.savefig(figname) + +plot_val(savefig=True) +plot_comparison(savefig=True) diff --git a/code/py/firm_hiring.py b/code/py/firm_hiring.py index 43b246c..aa430aa 100644 --- a/code/py/firm_hiring.py +++ b/code/py/firm_hiring.py @@ -88,7 +88,7 @@ def optimistic_policy_iteration(model, tolerance=1e-5, m=100): def plot_policy(savefig=False, - figname="./figures/firm_hiring_pol.pdf"): + figname="../figures_py/firm_hiring_pol.png"): model = create_hiring_model() β, κ, α, p, w, l_grid, z_grid, Q = model σ_star = optimistic_policy_iteration(model) @@ -130,7 +130,7 @@ def sim_dynamics(model, ts_length): def plot_sim(savefig=False, - figname="./figures/firm_hiring_ts.pdf", + figname="../figures_py/firm_hiring_ts.png", ts_length = 250): model = create_hiring_model() β, κ, α, p, w, l_grid, z_grid, Q = model @@ -149,7 +149,7 @@ def plot_sim(savefig=False, def plot_growth(savefig=False, - figname="./figures/firm_hiring_g.pdf", + figname="../figures_py/firm_hiring_g.png", ts_length = 10_000_000): model = create_hiring_model() @@ -170,3 +170,7 @@ def plot_growth(savefig=False, plt.show() if savefig: fig.savefig(figname) + +plot_sim(savefig=True) +plot_policy(savefig=True) +plot_growth(savefig=True) diff --git a/code/py/fosd_tauchen.py b/code/py/fosd_tauchen.py new file mode 100644 index 0000000..146bfb2 --- /dev/null +++ b/code/py/fosd_tauchen.py @@ -0,0 +1,32 @@ +import numpy as np +import matplotlib.pyplot as plt +from quantecon.markov import tauchen + +n = 25 +nu = 1.0 +a = 0.5 + +# Generating Markov Chain using Tauchen method +mc = tauchen(n, a, nu) +state_values, P = mc.state_values, mc.P +i, j = 8, 12 + +fig, axes = plt.subplots(2, 1, figsize=(10, 5.2)) +fontsize = 16 + +# Plotting the first subplot +ax = axes[0] +ax.plot(state_values, P[i, :], "b-o", alpha=0.4, lw=2, label=r"$\varphi$") +ax.plot(state_values, P[j, :], "g-o", alpha=0.4, lw=2, label=r"$\psi$") +ax.legend(frameon=False, fontsize=fontsize) + +# Plotting the second subplot +ax = axes[1] +F = [np.sum(P[i, k:]) for k in range(n)] +G = [np.sum(P[j, k:]) for k in range(n)] +ax.plot(state_values, F, "b-o", alpha=0.4, lw=2, label=r"$G^\varphi$") +ax.plot(state_values, G, "g-o", alpha=0.4, lw=2, label=r"$G^\psi$") +ax.legend(frameon=False, fontsize=fontsize) + +plt.show() +fig.savefig("../figures_py/fosd_tauchen_1.png") diff --git a/code/py/howard_newton.py b/code/py/howard_newton.py new file mode 100644 index 0000000..260ce7e --- /dev/null +++ b/code/py/howard_newton.py @@ -0,0 +1,50 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import newton +from autograd import grad + +v0 = 0.5 +T = lambda x: 1 + 0.2 * x**2 +DT = grad(T) +Tsp = lambda v, v0=v0: T(v0) + DT(v0) * (v - v0) # T_\sigma' + +vs = newton(lambda v: T(v) - v, 0.5) # find fixed point of T +v1 = (T(v0) - DT(v0) * v0) / (1 - DT(v0)) # v_\sigma' + +fs = 16 + +xmin, xmax = -0.2, 2.5 +xgrid = np.linspace(xmin, xmax, 1000) + +fig, ax = plt.subplots() + +lb_T = r"$T$" +ax.plot(xgrid, T(xgrid), lw=2, alpha=0.6, label=lb_T) + +lb_Tsp = r"$T_{\sigma'}$" +ax.plot(xgrid, Tsp(xgrid), lw=2, alpha=0.6, label=lb_Tsp) + +ax.plot(xgrid, xgrid, "k--", lw=1, alpha=0.7, label=r"$45^{\circ}$") + +fp1 = [v1] +ax.plot(fp1, fp1, "go", ms=5, alpha=0.6) +ax.plot([v0], [Tsp(v0)], "go", ms=5, alpha=0.6) +ax.plot([vs], [vs], "go", ms=5, alpha=0.6) + +ax.vlines([v0, vs, v1], [0, 0, 0], [Tsp(v0), vs, v1], + color="k", linestyle="-.", lw=0.4) + +ax.legend(frameon=False, fontsize=fs) + +ax.set_xticks([v0, vs, v1]) +ax.set_xticklabels([r"$v_\sigma$", r"$v^*$", r"$v_{\sigma'}$"], fontsize=fs) +ax.set_yticks([0]) + +ax.set_xlim(0, 2.6) +ax.set_ylim(0, 2.6) + +plt.show() + +filename = "../figures_py/howard_newton_1.png" + +fig.savefig(filename) diff --git a/code/py/iid_job_search_cv.py b/code/py/iid_job_search_cv.py new file mode 100644 index 0000000..e8cacd4 --- /dev/null +++ b/code/py/iid_job_search_cv.py @@ -0,0 +1,117 @@ +import numpy as np +import matplotlib.pyplot as plt +from numba import njit +from s_approx import successive_approx +from two_period_job_search import create_job_search_model + + +@njit +def g(h, model): + n, w_vals, ϕ, β, c = model + return c + β * np.dot(np.maximum(w_vals / (1 - β), h), ϕ) + +def compute_hstar_wstar(model, h_init=0.0): + n, w_vals, ϕ, β, c = model + h_star = successive_approx(lambda h: g(h, model), h_init) + return h_star, (1 - β) * h_star + +# Plot functions +def fig_g(model, savefig=False, fs=16, + figname="../figures_py/iid_job_search_g.png"): + n, w_vals, ϕ, β, c = model + h_grid = np.linspace(600, max(c, n) / (1 - β), 100) + g_vals = np.array([g(h, model) for h in h_grid]) + + fig, ax = plt.subplots(figsize=(9, 5.5)) + ax.plot(h_grid, g_vals, lw=2.0, label=r"$g$") + ax.plot(h_grid, h_grid, "k--", lw=1.0, label="45") + + ax.legend(frameon=False, fontsize=fs, loc="lower right") + + h_star, w_star = compute_hstar_wstar(model) + ax.plot(h_star, h_star, "go", ms=10, alpha=0.6) + + ax.annotate(r"$h^*$", + xy=(h_star, h_star), + xycoords="data", + xytext=(40, -40), + textcoords="offset points", + fontsize=fs) + + if savefig: + fig.savefig(figname) + + plt.show() + +def fig_tg(betas=[0.95, 0.96], savefig=False, fs=16, + figname="../figures_py/iid_job_search_tg.png"): + h_grid = np.linspace(600, 1200, 100) + fig, ax = plt.subplots(figsize=(9, 5.5)) + ax.plot(h_grid, h_grid, "k--", lw=1.0, label="45") + + for i, β in enumerate(betas): + model = create_job_search_model(β=β) + n, w_vals, ϕ, β, c = model + g_vals = np.array([g(h, model) for h in h_grid]) + + lb = f"$g_{i+1} \\; (\\beta_{i+1} = {β})$" + ax.plot(h_grid, g_vals, lw=2.0, label=lb) + + ax.legend(frameon=False, fontsize=fs, loc="lower right") + + h_star, w_star = compute_hstar_wstar(model) + ax.plot(h_star, h_star, "go", ms=10, alpha=0.6) + + lb = f"$h^*_{i+1}$" + ax.annotate(lb, + xy=(h_star, h_star), + xycoords="data", + xytext=(40, -40), + textcoords="offset points", + fontsize=fs) + + if savefig: + fig.savefig(figname) + + plt.show() + +def fig_cv(model, fs=16, savefig=False, + figname="../figures_py/iid_job_search_4.png"): + n, w_vals, ϕ, β, c = model + h_star, w_star = compute_hstar_wstar(model) + vhat = np.maximum(w_vals / (1 - β), h_star) + + fig, ax = plt.subplots() + ax.plot(w_vals, vhat, "k-", lw=2.0, label="value function") + ax.legend(fontsize=fs) + ax.set_ylim(0, np.max(vhat)) + + plt.show() + if savefig: + fig.savefig(figname) + +def fig_bf(betas=np.linspace(0.9, 0.99, 20), savefig=False, fs=16, + figname="../figures_py/iid_job_search_bf.png"): + h_vals = np.zeros_like(betas) + for i, β in enumerate(betas): + model = create_job_search_model(β=β) + h, w = compute_hstar_wstar(model) + h_vals[i] = h + + fig, ax = plt.subplots() + ax.plot(betas, h_vals, lw=2.0, alpha=0.7, label=r"$h^*(\beta)$") + ax.legend(frameon=False, fontsize=fs) + ax.set_xlabel(r"$\beta$", fontsize=fs) + ax.set_ylabel("continuation value", fontsize=fs) + + if savefig: + fig.savefig(figname) + + plt.show() + +default_model = create_job_search_model() + +fig_g(default_model, savefig=True) +fig_tg(savefig=True) +fig_cv(default_model, savefig=True) +fig_bf(savefig=True) diff --git a/code/py/inventory_cont_time.py b/code/py/inventory_cont_time.py new file mode 100644 index 0000000..38ce213 --- /dev/null +++ b/code/py/inventory_cont_time.py @@ -0,0 +1,50 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import expon, geom + +def sim_path(T=10, seed=123, λ=0.5, α=0.7, b=10): + J, Y = 0.0, b + J_vals, Y_vals = [J], [Y] + np.random.seed(seed) + + # Wait times are exponential + φ = expon(scale=1/λ) + # Orders are geometric + G = geom(p=α) + + while True: + W = φ.rvs() + J += W + J_vals.append(J) + if Y == 0: + Y = b + else: + U = G.rvs() + 1 # Geometric on 1, 2,... + Y = Y - min(Y, U) + Y_vals.append(Y) + if J > T: + break + + def X(t): + k = np.searchsorted(J_vals, t, side='right') - 1 + return Y_vals[k+1] + + return X + +T = 50 +X = sim_path(T=T) + +# Plotting +grid = np.linspace(0, T - 0.001, 100) + +fig, ax = plt.subplots(figsize=(9, 5.2)) +ax.step(grid, [X(t) for t in grid], label=r"$X_t$", alpha=0.7) + +ax.set_xticks([0, 10, 20, 30, 40, 50]) + +ax.set_xlabel("time", fontsize=12) +ax.set_ylabel("inventory", fontsize=12) +ax.legend(fontsize=12) + +plt.savefig("../figures_py/inventory_cont_time_1.png") +plt.show() diff --git a/code/py/inventory_dp.py b/code/py/inventory_dp.py index f14069a..931871b 100644 --- a/code/py/inventory_dp.py +++ b/code/py/inventory_dp.py @@ -3,6 +3,7 @@ import numpy as np from collections import namedtuple from numba import njit +import matplotlib.pyplot as plt # NamedTuple Model @@ -70,9 +71,6 @@ def solve_inventory_model(v_init, model): # == Plots == # -import matplotlib.pyplot as plt - - # Create an instance of the model and solve it model = create_inventory_model() β, K, c, κ, p = model @@ -93,7 +91,7 @@ def sim_inventories(ts_length=400, X_init=0): def plot_vstar_and_opt_policy(fontsize=10, - figname="./figures/inventory_dp_vs.pdf", + figname="../figures_py/inventory_dp_vs.png", savefig=False): fig, axes = plt.subplots(2, 1, figsize=(8, 6.5)) @@ -113,7 +111,7 @@ def plot_vstar_and_opt_policy(fontsize=10, def plot_ts(fontsize=10, - figname="./figures/inventory_dp_ts.pdf", + figname="../figures_py/inventory_dp_ts.png", savefig=False): X = sim_inventories() fig, ax = plt.subplots(figsize=(9, 5.5)) @@ -125,3 +123,6 @@ def plot_ts(fontsize=10, plt.show() if savefig: fig.savefig(figname) + +plot_vstar_and_opt_policy(savefig=True) +plot_ts(savefig=True) diff --git a/code/py/inventory_sdd.py b/code/py/inventory_sdd.py index 20f719a..2f698e2 100644 --- a/code/py/inventory_sdd.py +++ b/code/py/inventory_sdd.py @@ -15,6 +15,7 @@ from time import time from numba import njit, prange from collections import namedtuple +import matplotlib.pyplot as plt # NamedTuple Model Model = namedtuple("Model", ("K", "c", "κ", "p", "r", @@ -147,8 +148,6 @@ def optimistic_policy_iteration(v_init, # == Plots == # -import matplotlib.pyplot as plt - # Create an instance of the model and solve it model = create_sdd_inventory_model() K, c, κ, p, r, R, y_vals, z_vals, Q = model @@ -171,7 +170,7 @@ def sim_inventories(ts_length, X_init=0): def plot_ts(ts_length=400, fontsize=10, - figname="./figures/inventory_sdd_ts.pdf", + figname="../figures_py/inventory_sdd_ts.png", savefig=False): X, Z = sim_inventories(ts_length) @@ -221,5 +220,8 @@ def plot_timing(m_vals=np.arange(1, 400, 10), ax.set_ylabel("time", fontsize=fontsize) plt.show() if savefig: - fig.savefig("./figures/inventory_sdd_timing.pdf") - return (opi_time, vfi_time, opi_times) \ No newline at end of file + fig.savefig("../figures_py/inventory_sdd_timing.png") + return (opi_time, vfi_time, opi_times) + +plot_ts(savefig=True) +plot_timing(savefig=True) diff --git a/code/py/inventory_sim.py b/code/py/inventory_sim.py index ec8f3a4..78eb877 100644 --- a/code/py/inventory_sim.py +++ b/code/py/inventory_sim.py @@ -3,6 +3,7 @@ from itertools import product from quantecon import MarkovChain from collections import namedtuple +import matplotlib.pyplot as plt # NamedTuple Model Model = namedtuple("Model", ("S", "s", "p", "φ", "h")) @@ -45,15 +46,12 @@ def compute_stationary_dist(model): # Plots -import matplotlib.pyplot as plt - - default_model = create_inventory_model() def plot_ts(model, fontsize=16, - figname="./figures/inventory_sim_1.pdf", - savefig=False): + figname="../figures_py/inventory_sim_1.png", + savefig=False): S, s, p, φ, h = model X = sim_inventories(model) fig, ax = plt.subplots(figsize=(9, 5.2)) @@ -70,8 +68,8 @@ def plot_ts(model, fontsize=16, def plot_hist(model, fontsize=16, - figname="./figures/inventory_sim_2.pdf", - savefig=False): + figname="../figures_py/inventory_sim_2.png", + savefig=False): S, s, p, φ, h = model state_values, ψ_star = compute_stationary_dist(model) X = sim_inventories(model, 1_000_000) @@ -90,3 +88,6 @@ def plot_hist(model, fontsize=16, plt.show() if savefig: fig.savefig(figname) + +plot_ts(default_model, savefig=True) +plot_hist(default_model, savefig=True) diff --git a/code/py/js_with_sep_sim.py b/code/py/js_with_sep_sim.py new file mode 100644 index 0000000..b391b3c --- /dev/null +++ b/code/py/js_with_sep_sim.py @@ -0,0 +1,63 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import rv_discrete +from markov_js_with_sep import create_js_with_sep_model, vfi +from numba import njit + +# Create and solve model +model = create_js_with_sep_model() +n, w_vals, P, β, c, α = model +v_star, σ_star = vfi(model) + +# Create Markov distributions to draw from +P_dists = [rv_discrete(values=(np.arange(n), P[i, :])) for i in range(n)] + +def update_wages_idx(w_idx): + return P_dists[w_idx].rvs() + +@njit +def sim_wages(ts_length=100): + w_idx = np.random.randint(0, n-1) + W = np.zeros(ts_length) + for t in range(ts_length): + W[t] = w_vals[w_idx] + w_idx = update_wages_idx(w_idx) + return W + +def sim_outcomes(ts_length=100): + status = 0 + E, W = [], [] + w_idx = np.random.randint(0, n-1) + for t in range(ts_length): + if status == 0: + status = σ_star[w_idx] if σ_star[w_idx] else 0 + else: + status = 0 if np.random.rand() < α else 1 + W.append(w_vals[w_idx]) + E.append(status) + w_idx = update_wages_idx(w_idx) + return W, E + +# == Plots == # + +def plot_status(ts_length=100, savefig=False, + figname="../figures_py/js_with_sep_sim_1.png"): + W, E = sim_outcomes(ts_length=ts_length) + fs = 16 + fig, axes = plt.subplots(2, 1, figsize=(10, 8)) + + ax = axes[0] + ax.plot(W, label="wage offers") + ax.legend(fontsize=fs, frameon=False) + + ax = axes[1] + ax.set_yticks([0, 1]) + ax.set_yticklabels(["unempl.", "employed"]) + ax.plot(E, label="status") + ax.legend(fontsize=fs, frameon=False) + + if savefig: + fig.savefig(figname) + plt.show() + +plot_status(savefig=True) diff --git a/code/py/lake.py b/code/py/lake.py index 0855de7..de611a1 100644 --- a/code/py/lake.py +++ b/code/py/lake.py @@ -1,4 +1,5 @@ import numpy as np +import matplotlib.pyplot as plt α, λ, d, b = 0.01, 0.1, 0.02, 0.025 g = b - d @@ -14,10 +15,7 @@ # == Plots == # -import matplotlib.pyplot as plt - - -def plot_paths(figname="./figures/lake_1.pdf", savefig=False): +def plot_paths(figname="../figures_py/lake_1.png", savefig=False): path_length = 100 x_path_1 = np.zeros((2, path_length)) @@ -88,7 +86,7 @@ def plot_paths(figname="./figures/lake_1.pdf", savefig=False): -def plot_growth(savefig=False, figname="./figures/lake_2.pdf"): +def plot_growth(savefig=False, figname="../figures_py/lake_2.png"): path_length = 100 x_0 = 2.1, 1.2 @@ -114,3 +112,6 @@ def plot_growth(savefig=False, figname="./figures/lake_2.pdf"): plt.show() if savefig: fig.savefig(figname) + +plot_paths(savefig=True) +plot_growth(savefig=True) diff --git a/code/py/linear_iter_fig.py b/code/py/linear_iter_fig.py index 675f21f..bed580e 100644 --- a/code/py/linear_iter_fig.py +++ b/code/py/linear_iter_fig.py @@ -3,7 +3,7 @@ from linear_iter import x_star, T -def plot_main(savefig=False, figname="./figures/linear_iter_fig_1.pdf"): +def plot_main(savefig=False, figname="../figures_py/linear_iter_fig_1.png"): fig, ax = plt.subplots() @@ -42,3 +42,5 @@ def plot_main(savefig=False, figname="./figures/linear_iter_fig_1.pdf"): plt.show() if savefig: fig.savefig(figname) + +plot_main(savefig=True) diff --git a/code/py/loc_stable.py b/code/py/loc_stable.py new file mode 100644 index 0000000..9c69e1d --- /dev/null +++ b/code/py/loc_stable.py @@ -0,0 +1,48 @@ +import numpy as np +import matplotlib.pyplot as plt + +fontsize = 16 + +def subplots(): + """Custom subplots with axes through the origin""" + fig, ax = plt.subplots() + + # Set the axes through the origin + for spine in ["left", "bottom"]: + ax.spines[spine].set_position("zero") + for spine in ["right", "top"]: + ax.spines[spine].set_color("none") + + return fig, ax + +x_grid = np.linspace(0, 2, 200) + +# Function +def g(x): + return np.sqrt(x) + +# First order approximation at 1 +def hat_g(x): + return 1 + (1/2) * (x - 1) + +fig, ax = subplots() +ax.plot(x_grid, x_grid, 'k--', label=r"$45^\circ$") +ax.plot(x_grid, g(x_grid), label=r"$g$") +ax.plot(x_grid, hat_g(x_grid), label=r"$\hat g$") +ax.set_xticks([]) +ax.set_yticks([]) +ax.legend(loc="upper center", fontsize=16) + +fp = [1] +ax.plot(fp, fp, 'ko', ms=8, alpha=0.6) + +ax.annotate(r"$x^*$", + xy=(1, 1), + xycoords="data", + xytext=(40, -80), + textcoords="offset points", + fontsize=16, + arrowprops=dict(arrowstyle="->")) + +plt.savefig("../figures_py/loc_stable_1.png") +plt.show() diff --git a/code/py/markov_js.py b/code/py/markov_js.py index c41eeeb..2f60a4e 100644 --- a/code/py/markov_js.py +++ b/code/py/markov_js.py @@ -7,6 +7,7 @@ import numpy as np from collections import namedtuple from s_approx import successive_approx +import matplotlib.pyplot as plt # NamedTuple Model @@ -85,16 +86,13 @@ def policy_iteration(model): # == Plots == # -import matplotlib.pyplot as plt - - default_model = create_markov_js_model() def plot_main(model=default_model, method="vfi", savefig=False, - figname="./figures/markov_js_vfix.png"): + figname="../figures_py/markov_js_vfix.png"): n, w_vals, P, β, c = model if method == "vfi": @@ -115,3 +113,5 @@ def plot_main(model=default_model, plt.show() if savefig: fig.savefig(figname) + +plot_main(savefig=True) diff --git a/code/py/markov_js_with_sep.py b/code/py/markov_js_with_sep.py index 02b5caa..191a4f9 100644 --- a/code/py/markov_js_with_sep.py +++ b/code/py/markov_js_with_sep.py @@ -7,6 +7,7 @@ import numpy as np from collections import namedtuple from s_approx import successive_approx +import matplotlib.pyplot as plt # NamedTuple Model @@ -54,15 +55,12 @@ def vfi(model): # == Plots == # -import matplotlib.pyplot as plt - - default_model = create_js_with_sep_model() def plot_main(model=default_model, savefig=False, - figname="./figures/markov_js_with_sep_1.pdf"): + figname="../figures_py/markov_js_with_sep_1.png"): n, w_vals, P, β, c, α = model v_star, σ_star = vfi(model) @@ -93,7 +91,7 @@ def plot_main(model=default_model, def plot_w_stars(α_vals=np.linspace(0.0, 1.0, 10), savefig=False, - figname="./figures/markov_js_with_sep_2.pdf"): + figname="../figures_py/markov_js_with_sep_2.png"): w_star_vec = np.empty_like(α_vals) for (i_α, α) in enumerate(α_vals): @@ -124,3 +122,6 @@ def plot_w_stars(α_vals=np.linspace(0.0, 1.0, 10), plt.show() if savefig: fig.savefig(figname) + +plot_main(savefig=True) +plot_w_stars(savefig=True) diff --git a/code/py/modified_opt_savings.py b/code/py/modified_opt_savings.py index c025fe5..5ccfbb9 100644 --- a/code/py/modified_opt_savings.py +++ b/code/py/modified_opt_savings.py @@ -4,6 +4,7 @@ from collections import namedtuple from numba import njit, prange from math import floor +import matplotlib.pyplot as plt # NamedTuple Model @@ -184,11 +185,8 @@ def lorenz(v): # assumed sorted vector # Plots -import matplotlib.pyplot as plt - - def plot_contours(savefig=False, - figname="./figures/modified_opt_savings_1.pdf"): + figname="../figures_py/modified_opt_savings_1.png"): model = create_savings_model() β, γ, η_grid, φ, w_grid, y_grid, Q = model @@ -221,6 +219,7 @@ def plot_contours(savefig=False, fig.savefig(figname) plt.show() + def plot_policies(savefig=False): model = create_savings_model() β, γ, η_grid, φ, w_grid, y_grid, Q = model @@ -237,7 +236,8 @@ def plot_policies(savefig=False): ax.legend() plt.show() if savefig: - fig.savefig(f"./figures/modified_opt_saving_2.pdf") + fig.savefig("../figures_py/modified_opt_saving_2.png") + def plot_time_series(m=2_000, savefig=False): @@ -248,7 +248,8 @@ def plot_time_series(m=2_000, savefig=False): ax.legend() plt.show() if savefig: - fig.savefig("./figures/modified_opt_saving_ts.pdf") + fig.savefig("../figures_py/modified_opt_saving_ts.png") + def plot_histogram(m=1_000_000, savefig=False): @@ -262,7 +263,8 @@ def plot_histogram(m=1_000_000, savefig=False): plt.show() if savefig: - fig.savefig("./figures/modified_opt_saving_hist.pdf") + fig.savefig("../figures_py/modified_opt_saving_hist.png") + def plot_lorenz(m=1_000_000, savefig=False): @@ -277,4 +279,10 @@ def plot_lorenz(m=1_000_000, savefig=False): plt.show() if savefig: - fig.savefig("./figures/modified_opt_saving_lorenz.pdf") \ No newline at end of file + fig.savefig("../figures_py/modified_opt_saving_lorenz.png") + +plot_contours(savefig=True) +plot_histogram(savefig=True) +plot_lorenz(savefig=True) +plot_policies(savefig=True) +plot_time_series(savefig=True) diff --git a/code/py/newton.py b/code/py/newton.py new file mode 100644 index 0000000..b4ab03a --- /dev/null +++ b/code/py/newton.py @@ -0,0 +1,63 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.optimize import root_scalar +from autograd import grad + +fs = 16 + +# Define the function T and its derivative +def T(x): + return 1 + (x / (x + 1)) + +# Autograd to compute the derivative of T +DT = grad(T) + +x0 = 0.5 + +def T_hat(x, x0=x0): + return T(x0) + DT(x0) * (x - x0) + +# Find the fixed point of T +res = root_scalar(lambda x: T(x) - x, x0=0.5) +xs = res.root + +# Compute the Newton approximation x1 +x1 = (T(x0) - DT(x0) * x0) / (1 - DT(x0)) + +# Plot +def plot_45(file_name="../figures_py/newton_1.png", xmin=0.0, xmax=2.6, + savefig=False): + xgrid = np.linspace(xmin, xmax, 1000) + + fig, ax = plt.subplots() + + lb_T = r"$T$" + ax.plot(xgrid, T(xgrid), lw=2, alpha=0.6, label=lb_T) + + lb_T_hat = r"$\hat{T}$" + ax.plot(xgrid, T_hat(xgrid), lw=2, alpha=0.6, label=lb_T_hat) + + ax.plot(xgrid, xgrid, "k--", lw=1, alpha=0.7, label=r"$45^\circ$") + + fp1 = [x1] + ax.plot(fp1, fp1, "go", ms=5, alpha=0.6) + ax.plot([x0], [T_hat(x0)], "go", ms=5, alpha=0.6) + ax.plot([xs], [xs], "go", ms=5, alpha=0.6) + + ax.vlines([x0, xs, x1], [0, 0, 0], [T_hat(x0), xs, x1], + color="k", linestyle="-.", lw=0.4) + + ax.legend(frameon=False, fontsize=fs) + + ax.set_xticks([x0, xs, x1]) + ax.set_xticklabels([r"$u_0$", r"$u^*$", r"$u_1$"], fontsize=fs) + ax.set_yticks([0]) + + ax.set_xlim(0, 2.6) + ax.set_ylim(0, 2.6) + + plt.show() + if savefig: + fig.savefig(file_name) + +plot_45(savefig=True) diff --git a/code/py/newton_solow.py b/code/py/newton_solow.py new file mode 100644 index 0000000..5246557 --- /dev/null +++ b/code/py/newton_solow.py @@ -0,0 +1,94 @@ +import numpy as np +import matplotlib.pyplot as plt +from autograd import elementwise_grad + +A, s, alpha, delta = 2, 0.3, 0.3, 0.4 +x0 = 0.25 +n = 14 + +def g(k): + return A * s * k**alpha + (1 - delta) * k + +Dg = elementwise_grad(g) + +def q(x): + return (g(x) - Dg(x) * x) / (1 - Dg(x)) + +fs = 14 +kstar = ((s * A) / delta)**(1/(1 - alpha)) + +def plot_45(file_name="../figures_py/newton_solow_45.png", xmin=0.0, + xmax=4, savefig=False): + xgrid = np.linspace(xmin, xmax, 1200) + + fig, ax = plt.subplots() + + lb_g = r"$g$" + ax.plot(xgrid, g(xgrid), lw=2, alpha=0.6, label=lb_g) + + lb_q = r"$Q$" + ax.plot(xgrid, q(xgrid), lw=2, alpha=0.6, label=lb_q) + + ax.plot(xgrid, xgrid, "k--", lw=1, alpha=0.7, label=r"$45^\circ$") + + fps = [kstar] + ax.plot(fps, fps, "go", ms=10, alpha=0.6) + + ax.legend(frameon=False, fontsize=fs) + + ax.set_xlabel(r"$k_t$", fontsize=fs) + ax.set_ylabel(r"$k_{t+1}$", fontsize=fs) + + ax.set_ylim(-3, 4) + ax.set_xlim(0, 4) + + plt.show() + if savefig: + fig.savefig(file_name) + + +def compute_iterates(k0, f): + k = k0 + k_iterates = [] + for _ in range(n): + k_iterates.append(k) + k = f(k) + return k_iterates + + +def plot_trajectories(file_name="../figures_py/newton_solow_traj.png", savefig=False): + x_grid = np.arange(1, n + 1) + + fig, axes = plt.subplots(2, 1) + ax1, ax2 = axes + + k0_a, k0_b = 0.8, 3.1 + + ks1 = compute_iterates(k0_a, g) + ax1.plot(x_grid, ks1, "-o", label="successive approximation") + + ks2 = compute_iterates(k0_b, g) + ax2.plot(x_grid, ks2, "-o", label="successive approximation") + + ks3 = compute_iterates(k0_a, q) + ax1.plot(x_grid, ks3, "-o", label="newton steps") + + ks4 = compute_iterates(k0_b, q) + ax2.plot(x_grid, ks4, "-o", label="newton steps") + + for ax in axes: + ax.plot(x_grid, kstar * np.ones(n), "k--") + ax.legend(fontsize=fs, frameon=False) + ax.set_ylim(0.6, 3.2) + xticks = [2, 4, 6, 8, 10, 12] + ax.set_xticks(xticks) + ax.set_xticklabels([str(s) for s in xticks], fontsize=fs) + ax.set_yticks([kstar]) + ax.set_yticklabels([r"$k^*$"], fontsize=fs) + + plt.show() + if savefig: + fig.savefig(file_name) + +plot_45(savefig=True) +plot_trajectories(savefig=True) diff --git a/code/py/optimality_illustration.py b/code/py/optimality_illustration.py new file mode 100644 index 0000000..cccafb4 --- /dev/null +++ b/code/py/optimality_illustration.py @@ -0,0 +1,56 @@ +import numpy as np +import matplotlib.pyplot as plt + +fs = 18 + +xmin = -0.5 +xmax = 2.0 + +xgrid = np.linspace(xmin, xmax, 1000) + +a1, b1 = 0.1, 0.5 # first T_σ +a2, b2 = 0.65, 0.4 # second T_σ + +v1 = b1 / (1 - a1) +v2 = b2 / (1 - a2) + +T1 = a1 * xgrid + b1 +T2 = a2 * xgrid + b2 +T = np.maximum(T1, T2) + +fig, ax = plt.subplots() + +# Set the axes through the origin +for spine in ["left", "bottom"]: + ax.spines[spine].set_position("zero") +for spine in ["right", "top"]: + ax.spines[spine].set_color("none") + +ax.plot(xgrid, xgrid, "k--", lw=1, alpha=0.7, label=r"$45^{\circ}$") +ax.plot(xgrid, T1, "k-", lw=1) +ax.plot(xgrid, T2, "k-", lw=1) + +ax.plot(xgrid, T, lw=6, alpha=0.3, color="blue", label=r"$T$") + +ax.plot((v1,), (v1,), "go", ms=5, alpha=0.8) +ax.plot((v2,), (v2,), "go", ms=5, alpha=0.8) + +ax.vlines((v1, v2), (0, 0), (v1, v2), color="k", linestyle="-.", lw=0.4) + +ax.text(2.1, 0.6, r"$T_{\sigma'}$", fontsize=fs) +ax.text(2.1, 1.4, r"$T_{\sigma''}$", fontsize=fs) + +ax.legend(frameon=False, loc="upper center", fontsize=fs) + +ax.set_xticks([v1, v2]) +ax.set_xticklabels([r"$v_{\sigma'}$", r"$v_{\sigma''} = v^*$"], fontsize=fs) +ax.set_yticks([]) + +ax.set_xlim(xmin, xmax + 0.5) +ax.set_ylim(-0.2, 2) +ax.text(2.4, -0.15, r"$v$", fontsize=22) + +plt.show() + +file_name = "../figures_py/optimality_illustration_1.png" +fig.savefig(file_name) diff --git a/code/py/pd_ratio.py b/code/py/pd_ratio.py index 16047e6..8ca0694 100644 --- a/code/py/pd_ratio.py +++ b/code/py/pd_ratio.py @@ -6,6 +6,7 @@ from quantecon.markov import tauchen import numpy as np from collections import namedtuple +import matplotlib.pyplot as plt # NamedTuple Model @@ -50,16 +51,12 @@ def pd_ratio(model): # == Plots == # - -import matplotlib.pyplot as plt - - default_model = create_asset_pricing_model() def plot_main(μ_d_vals=(0.02, 0.08), savefig=False, - figname="./figures/pd_ratio_1.pdf"): + figname="../figures_py/pd_ratio_1.png"): fig, ax = plt.subplots(figsize=(9, 5.2)) for μ_d in μ_d_vals: @@ -74,3 +71,5 @@ def plot_main(μ_d_vals=(0.02, 0.08), plt.show() if savefig: fig.savefig(figname) + +plot_main(savefig=True) diff --git a/code/py/plot_interest_rates.py b/code/py/plot_interest_rates.py index 9e02d17..cdab4a5 100644 --- a/code/py/plot_interest_rates.py +++ b/code/py/plot_interest_rates.py @@ -8,8 +8,8 @@ import numpy as np import matplotlib.pyplot as plt -df_nominal = pd.read_csv("./data/GS1.csv") -df_real = pd.read_csv("./data/WFII10.csv") +df_nominal = pd.read_csv("data/GS1.csv") +df_real = pd.read_csv("data/WFII10.csv") def plot_rates(df, fontsize=16, savefig=False): r_type = 'nominal' if df.equals(df_nominal) else 'real' @@ -20,4 +20,7 @@ def plot_rates(df, fontsize=16, savefig=False): ax.legend(fontsize=fontsize, frameon=False) plt.show() if savefig: - fig.savefig(f'./figures/plot_interest_rates_{r_type}.pdf') \ No newline at end of file + fig.savefig(f'../figures_py/plot_interest_rates_{r_type}.png') + +plot_rates(df_nominal, savefig=True) +plot_rates(df_real, savefig=True) \ No newline at end of file diff --git a/code/py/quantile_function.py b/code/py/quantile_function.py index 45d59ff..331111c 100644 --- a/code/py/quantile_function.py +++ b/code/py/quantile_function.py @@ -1,20 +1,22 @@ from scipy.stats import rv_discrete import numpy as np - from numba import njit -"Compute the τ-th quantile of v(X) when X ∼ ϕ and v = sort(v)." + @njit def quantile(τ, v, ϕ): + """Compute the τ-th quantile of v(X) when X ∼ ϕ and v = sort(v).""" for (i, v_value) in enumerate(v): p = sum(ϕ[:i+1]) # sum all ϕ[j] s.t. v[j] ≤ v_value if p >= τ: # exit and return v_value if prob ≥ τ return v_value -"For each i, compute the τ-th quantile of v(Y) when Y ∼ P(i, ⋅)" + def R(τ, v, P): + """For each i, compute the τ-th quantile of v(Y) when Y ∼ P(i, ⋅)""" return np.array([quantile(τ, v, P[i, :]) for i in range(len(v))]) + def quantile_test(τ): ϕ = [0.1, 0.2, 0.7] v = [10, 20, 30] diff --git a/code/py/quantile_js.py b/code/py/quantile_js.py index 545ac07..048facd 100644 --- a/code/py/quantile_js.py +++ b/code/py/quantile_js.py @@ -4,7 +4,8 @@ """ from quantecon import tauchen, MarkovChain import numpy as np -from quantile_function import * +from quantile_function import R +import matplotlib.pyplot as plt "Creates an instance of the job search model." def create_markov_js_model( @@ -31,15 +32,16 @@ def T_σ(v, σ, model): e = w_vals / (1 - β) return σ * e + (1 - σ) * h -" Get a v-greedy policy." + def get_greedy(v, model): + """Get a v-greedy policy.""" n, w_vals, P, β, c, τ = model σ = w_vals / (1 - β) >= c + β * R(τ, v, P) return σ -"Optimistic policy iteration routine." def optimistic_policy_iteration(model, tolerance=1e-5, m=100): + """Optimistic policy iteration routine.""" n, w_vals, P, β, c, τ = model v = np.ones(n) error = tolerance + 1 @@ -57,12 +59,9 @@ def optimistic_policy_iteration(model, tolerance=1e-5, m=100): # == Plots == # -import matplotlib.pyplot as plt - - def plot_main(tau_vals=(0.1, 0.25, 0.5, 0.6, 0.7, 0.8), savefig=False, - figname="./figures/quantile_js.pdf"): + figname="../figures_py/quantile_js.png"): w_star_vals = np.zeros(len(tau_vals)) @@ -91,3 +90,4 @@ def plot_main(tau_vals=(0.1, 0.25, 0.5, 0.6, 0.7, 0.8), if savefig: fig.savefig(figname) +plot_main(savefig=True) diff --git a/code/py/random_walk.py b/code/py/random_walk.py new file mode 100644 index 0000000..012e64f --- /dev/null +++ b/code/py/random_walk.py @@ -0,0 +1,19 @@ +import numpy as np +import matplotlib.pyplot as plt + +fontsize = 16 + +fig, ax = plt.subplots(figsize=(9, 5.2)) + +n, m = 100, 12 +cols = ["k-", "b-", "g-"] + +for _ in range(m): + s = np.random.choice(cols) + ax.plot(np.cumsum(np.random.randn(n)), s, alpha=0.5) + +ax.set_xlabel("time") +plt.show() + +file_name = "../figures_py/random_walk_1.png" +fig.savefig(file_name) diff --git a/code/py/rates.py b/code/py/rates.py new file mode 100644 index 0000000..f1f7042 --- /dev/null +++ b/code/py/rates.py @@ -0,0 +1,25 @@ +import matplotlib.pyplot as plt + +fontsize = 16 + +qs = [1, 2] +labels = ["linear", "quadratic"] +x0 = 0.9 +N = 26 +beta = 0.9 + +fig, ax = plt.subplots() + +for q, label in zip(qs, labels): + current_x = x0 + x = [] + for t in range(1, N + 1): + x.append(current_x) + current_x = beta * current_x**q + ax.plot(x, "o-", label=label, alpha=0.6) + +ax.legend(fontsize=16) + +file_name = "../figures_py/rates_1.png" +plt.savefig(file_name) +plt.show() diff --git a/code/py/risk_sensitive_js.py b/code/py/risk_sensitive_js.py index b85c50a..6f99f9d 100644 --- a/code/py/risk_sensitive_js.py +++ b/code/py/risk_sensitive_js.py @@ -9,6 +9,7 @@ import numpy as np from numba import njit from collections import namedtuple +import matplotlib.pyplot as plt # NamedTuple Model Model = namedtuple("Model", ("n", "w_vals", "P", "β", "c", "θ")) @@ -63,13 +64,9 @@ def vfi(model): # == Plots == # - -import matplotlib.pyplot as plt - - def plot_main(theta_vals=(-10, 0.0001, 0.1), savefig=False, - figname="./figures/risk_sensitive_js.pdf"): + figname="../figures_py/risk_sensitive_js.png"): fig, axes = plt.subplots(len(theta_vals), 1, figsize=(9, 22)) @@ -93,3 +90,5 @@ def plot_main(theta_vals=(-10, 0.0001, 0.1), if savefig: fig.savefig(figname) + +plot_main(savefig=True) diff --git a/code/py/rs_utility.py b/code/py/rs_utility.py index f091d28..9d6612a 100644 --- a/code/py/rs_utility.py +++ b/code/py/rs_utility.py @@ -4,6 +4,7 @@ import numpy as np from numba import njit from collections import namedtuple +import matplotlib.pyplot as plt # NamedTuple Model @@ -38,12 +39,8 @@ def compute_rs_utility(model): # Plots - -import matplotlib.pyplot as plt - - def plot_v(savefig=False, - figname="./figures/rs_utility_1.pdf"): + figname="../figures_py/rs_utility_1.png"): fig, ax = plt.subplots(figsize=(10, 5.2)) model = create_rs_utility_model() @@ -67,7 +64,7 @@ def plot_v(savefig=False, def plot_multiple_v(savefig=False, - figname="./figures/rs_utility_2.pdf"): + figname="../figures_py/rs_utility_2.png"): fig, ax = plt.subplots(figsize=(10, 5.2)) σ_vals = 0.05, 0.1 @@ -85,3 +82,6 @@ def plot_multiple_v(savefig=False, plt.show() if savefig: fig.savefig(figname) + +plot_v(savefig=True) +plot_multiple_v(savefig=True) diff --git a/code/py/solow_fp.py b/code/py/solow_fp.py index c780aa3..c11735e 100644 --- a/code/py/solow_fp.py +++ b/code/py/solow_fp.py @@ -61,4 +61,4 @@ def plot_45(ax, k0=0.5, plot_45(ax, A=2.0, s=0.3, alpha=0.4, delta=0.4) fig.tight_layout() plt.show() -fig.savefig("./figures/solow_fp.pdf") +fig.savefig("../figures_py/solow_fp.png") diff --git a/code/py/solow_fp_adjust.py b/code/py/solow_fp_adjust.py new file mode 100644 index 0000000..80901d8 --- /dev/null +++ b/code/py/solow_fp_adjust.py @@ -0,0 +1,62 @@ +import matplotlib.pyplot as plt +import numpy as np + + +def subplots(fs): + """Custom subplots with axes through the origin""" + fig, ax = plt.subplots(2, 2, figsize=fs) + + for i in range(2): + for j in range(2): + # Set the axes through the origin + for spine in ["left", "bottom"]: + ax[i, j].spines[spine].set_position("zero") + ax[i, j].spines[spine].set_color("black") + for spine in ["right", "top"]: + ax[i, j].spines[spine].set_color("none") + + return fig, ax + +A = [2, 2.5, 2, 2] +s = [0.3, 0.3, 0.2, 0.3] +alpha = [0.3, 0.3, 0.3, 0.3] +delta = [0.4, 0.4, 0.4, 0.6] +x0 = 0.25 +num_arrows = 8 +ts_length = 12 +xmin, xmax = 0, 3 + +def g(k, s, A, delta, alpha): + return A * s * k**alpha + (1 - delta) * k + +def kstar(s, A, delta, alpha): + return ((s * A) / delta)**(1/(1 - alpha)) + +xgrid = np.linspace(xmin, xmax, 120) + +fig, ax = subplots((10, 7)) + +# (0,0) is the default parameters +# (0,1) increases A +# (1,0) decreases s +# (1,1) increases delta + +lb = ["default", r"$A=2.5$", r"$s=.2$", r"$\delta=.6$"] +count = 0 + +for i in range(2): + for j in range(2): + ax[i, j].set_xlim(xmin, xmax) + ax[i, j].set_ylim(xmin, xmax) + ax[i, j].plot(xgrid, g(xgrid, s[count], A[count], delta[count], alpha[count]), + "b-", lw=2, alpha=0.6, label=lb[count]) + ks = kstar(s[count], A[count], delta[count], alpha[count]) + ax[i, j].plot(ks, ks, "go") + ax[i, j].plot(xgrid, xgrid, "k-", lw=1, alpha=0.7) + count += 1 + ax[i, j].legend(loc="lower right", frameon=False, fontsize=14) + +plt.show() + +file_name = "../figures_py/solow_fp_adjust.png" +fig.savefig(file_name) diff --git a/code/py/stoch_dom_finite.py b/code/py/stoch_dom_finite.py new file mode 100644 index 0000000..a843e0e --- /dev/null +++ b/code/py/stoch_dom_finite.py @@ -0,0 +1,23 @@ +import matplotlib.pyplot as plt +from matplotlib.ticker import MaxNLocator + +p, q = 0.75, 0.25 +fig, axes = plt.subplots(1, 2, figsize=(10, 5.2)) + +# First subplot +ax = axes[0] +ax.bar([1, 2], [p, 1-p], label=r"$\phi$") +ax.set_xticks([1, 2]) +ax.xaxis.set_major_locator(MaxNLocator(integer=True)) + +# Second subplot +ax = axes[1] +ax.bar([1, 2], [q, 1-q], label=r"$\psi$") +ax.set_xticks([1, 2]) +ax.xaxis.set_major_locator(MaxNLocator(integer=True)) + +# Add legends +for ax in axes: + ax.legend() + +plt.show() diff --git a/code/py/tauchen.py b/code/py/tauchen.py new file mode 100644 index 0000000..c6250ca --- /dev/null +++ b/code/py/tauchen.py @@ -0,0 +1,33 @@ +import numpy as np +import matplotlib.pyplot as plt +from quantecon.markov import tauchen + +ρ, b, ν = 0.9, 0.0, 1.0 +μ_x = b / (1 - ρ) +σ_x = np.sqrt(ν**2 / (1.0 - ρ**2)) + +# Number of states +n = 15 + +mc = tauchen(n, ρ, ν) +approx_sd = mc.stationary_distributions[0] + +def psi_star(y): + c = 1 / (np.sqrt(2 * np.pi) * σ_x) + return c * np.exp(-(y - μ_x)**2 / (2 * σ_x**2)) + +# == Plots == # +fontsize = 10 + +fig, ax = plt.subplots() + +ax.bar(mc.state_values, approx_sd, width=0.6, alpha=0.6, label="approximation") + +x_grid = np.linspace(min(mc.state_values) - 2, max(mc.state_values) + 2, 100) +ax.plot(x_grid, psi_star(x_grid), '-k', lw=2, alpha=0.6, label=r"$\psi^*$") + +ax.legend(fontsize=fontsize) + +plt.show() + +plt.savefig("../figures_py/tauchen_1.png") diff --git a/code/py/three_fixed_points.py b/code/py/three_fixed_points.py new file mode 100644 index 0000000..5b67ba8 --- /dev/null +++ b/code/py/three_fixed_points.py @@ -0,0 +1,38 @@ +import numpy as np +import matplotlib.pyplot as plt + +fs = 14 +xmin, xmax = 0.0000001, 2.0 + +def g(u): + return 2.125 / (1 + u**-4) + +xgrid = np.linspace(xmin, xmax, 200) + +fig, ax = plt.subplots(figsize=(6.5, 6)) + +ax.set_xlim(xmin, xmax) +ax.set_ylim(xmin, xmax) + +ax.plot(xgrid, g(xgrid), "b-", lw=2, alpha=0.6, label=r"$T$") +ax.plot(xgrid, xgrid, "k-", lw=1, alpha=0.7, label=r"$45^\circ$") + +ax.legend(fontsize=fs) + +fps = [0.01, 0.94, 1.98] +fps_labels = [r"$u_\ell$", r"$u_m$", r"$u_h$"] +coords = [(40, 80), (40, -40), (-40, -80)] + +ax.plot(fps, fps, "ro", ms=8, alpha=0.6) + +for fp, lb, coord in zip(fps, fps_labels, coords): + ax.annotate(lb, + xy=(fp, fp), + xycoords="data", + xytext=coord, + textcoords="offset points", + fontsize=fs, + arrowprops=dict(arrowstyle="->")) + +plt.savefig("../figures_py/three_fixed_points.png") +plt.show() diff --git a/code/py/two_period_job_search.py b/code/py/two_period_job_search.py index f8ebc17..4dd2d63 100644 --- a/code/py/two_period_job_search.py +++ b/code/py/two_period_job_search.py @@ -3,6 +3,7 @@ import numpy as np from numba import njit from collections import namedtuple +import matplotlib.pyplot as plt # NamedTuple Model @@ -52,9 +53,6 @@ def res_wage(model): ##### Plots ##### -import matplotlib.pyplot as plt - - default_model = create_job_search_model() @@ -69,7 +67,7 @@ def fig_dist(model=default_model, fs=10): def fig_v1(model=default_model, savefig=False, - figname="./figures/iid_job_search_0_py.pdf", fs=18): + figname="../figures_py/iid_job_search_0_py.png", fs=12): """ Plot two-period value function and res wage """ @@ -100,3 +98,5 @@ def fig_v1(model=default_model, savefig=False, if savefig: fig.savefig(figname) plt.show() + +fig_v1(savefig=True) diff --git a/code/py/up_down_stable.py b/code/py/up_down_stable.py new file mode 100644 index 0000000..72fe908 --- /dev/null +++ b/code/py/up_down_stable.py @@ -0,0 +1,48 @@ +import numpy as np +import matplotlib.pyplot as plt + +fs = 12 +xmin, xmax = 0., 1.0 + +def g(x): + return 0.2 + 0.6 * x**1.2 + +xgrid = np.linspace(xmin, xmax, 200) + +fig, ax = plt.subplots(figsize=(8.0, 6)) + +for spine in ["left", "bottom"]: + ax.spines[spine].set_position("zero") +for spine in ["right", "top"]: + ax.spines[spine].set_color("none") + +ax.set_xlim(xmin, xmax) +ax.set_ylim(xmin, xmax) + +ax.plot(xgrid, g(xgrid), "b-", lw=2, alpha=0.6, label=r"$T$") +ax.plot(xgrid, xgrid, "k-", lw=1, alpha=0.7, label=r"$45^\circ$") + +ax.legend(frameon=False, fontsize=fs) + +fp = (0.4,) +fps_label = r"$\bar{v}$" +coords = (40, -20) + +ax.plot(fp, fp, "ro", ms=8, alpha=0.6) + +ax.annotate(fps_label, + xy=(fp[0], fp[0]), + xycoords="data", + xytext=coords, + textcoords="offset points", + fontsize=fs, + arrowprops=dict(arrowstyle="->")) + +ax.set_xticks([0, 1]) +ax.set_xticklabels([r"$0$", r"$1$"], fontsize=fs) +ax.set_yticks([]) + +ax.set_xlabel(r"$V$", fontsize=fs) + +plt.savefig("../figures_py/up_down_stable.png") +plt.show() diff --git a/code/py/v_star_illus.py b/code/py/v_star_illus.py new file mode 100644 index 0000000..bbfae7f --- /dev/null +++ b/code/py/v_star_illus.py @@ -0,0 +1,41 @@ +import numpy as np +import matplotlib.pyplot as plt + +fs = 12 + +xmin, xmax = 0.01, 2.0 + +xgrid = np.linspace(xmin, xmax, 1000) + +v1 = xgrid ** 0.7 +v2 = xgrid ** 0.1 + 0.05 +v = np.maximum(v1, v2) + +fig, ax = plt.subplots() + +for spine in ["left", "bottom"]: + ax.spines[spine].set_position("zero") +for spine in ["right", "top"]: + ax.spines[spine].set_color("none") + +ax.plot(xgrid, v1, "k-", lw=1) +ax.plot(xgrid, v2, "k-", lw=1) +ax.plot(xgrid, v, lw=6, alpha=0.3, color="blue", + label=r"$v^* = \bigvee_{\sigma \in \Sigma} v_\sigma$") + +ax.text(2.1, 1.1, r"$v_{\sigma'}$", fontsize=fs) +ax.text(2.1, 1.6, r"$v_{\sigma''}$", fontsize=fs) +ax.text(1.2, 0.3, r"$\Sigma = \{\sigma', \sigma''\}$", fontsize=fs) + +ax.legend(frameon=False, loc="upper left", fontsize=fs) + +ax.set_xlim(xmin, xmax + 0.5) +ax.set_ylim(0.0, 2) +ax.text(2.4, -0.15, r"$x$", fontsize=20) + +ax.set_xticks([]) +ax.set_yticks([]) + +plt.show() +file_name = "../figures_py/v_star_illus.png" +fig.savefig(file_name) diff --git a/code/py/val_consumption.py b/code/py/val_consumption.py new file mode 100644 index 0000000..921c032 --- /dev/null +++ b/code/py/val_consumption.py @@ -0,0 +1,30 @@ +import numpy as np +import matplotlib.pyplot as plt +from quantecon.markov import tauchen +from numpy.linalg import solve + + +def compute_v(n=25, β=0.98, ρ=0.96, ν=0.05, γ=2.0): + mc = tauchen(n, ρ, ν) + x_vals = mc.state_values + P = mc.P + r = np.exp((1 - γ) * x_vals) / (1 - γ) + v = solve(np.eye(n) - β * P, r) + return x_vals, v + + +def plot_v(savefig=False, figname="../figures_py/val_consumption_1.png"): + fontsize = 12 + + fig, ax = plt.subplots(figsize=(10, 5.2)) + x_vals, v = compute_v() + ax.plot(x_vals, v, lw=2, alpha=0.7, label=r"$v$") + ax.set_xlabel(r"$x$", fontsize=fontsize) + ax.legend(frameon=False, fontsize=fontsize, loc="upper left") + ax.set_ylim(-65, -40) + plt.show() + if savefig: + fig.savefig(figname) + + +plot_v(savefig=True) diff --git a/pdf/dp.pdf b/pdf/dp.pdf index 875e41a..ad82764 100644 Binary files a/pdf/dp.pdf and b/pdf/dp.pdf differ