From befd388cc847d1d11caa0b329a9375f03e07c378 Mon Sep 17 00:00:00 2001 From: annabellasd Date: Tue, 10 Oct 2023 15:34:07 -0800 Subject: [PATCH] updates --- lectures/multi_agent_models/aiyagari.md | 50 +++--- lectures/multi_agent_models/arellano.md | 69 ++++----- lectures/multi_agent_models/harrison_kreps.md | 22 +-- lectures/multi_agent_models/lake_model.md | 146 +++++++++--------- lectures/multi_agent_models/lucas_model.md | 42 ++--- lectures/multi_agent_models/markov_asset.md | 92 ++++++----- lectures/multi_agent_models/markov_perf.md | 40 ++--- lectures/multi_agent_models/matsuyama.md | 104 ++++++------- .../rational_expectations.md | 32 ++-- lectures/multi_agent_models/schelling.md | 6 +- .../multi_agent_models/uncertainty_traps.md | 89 +++++------ 11 files changed, 343 insertions(+), 349 deletions(-) diff --git a/lectures/multi_agent_models/aiyagari.md b/lectures/multi_agent_models/aiyagari.md index 11ae13bb..de0f8a5e 100644 --- a/lectures/multi_agent_models/aiyagari.md +++ b/lectures/multi_agent_models/aiyagari.md @@ -210,23 +210,23 @@ using LaTeXStrings, Parameters, Plots, QuantEcon ``` ```{code-cell} julia -Household = @with_kw (r = 0.01, - w = 1.0, - σ = 1.0, - β = 0.96, - z_chain = MarkovChain([0.9 0.1; 0.1 0.9], [0.1; 1.0]), - a_min = 1e-10, - a_max = 18.0, - a_size = 200, - a_vals = range(a_min, a_max, length = a_size), - z_size = length(z_chain.state_values), - n = a_size * z_size, - s_vals = gridmake(a_vals, z_chain.state_values), - s_i_vals = gridmake(1:a_size, 1:z_size), - u = σ == 1 ? x -> log(x) : x -> (x^(1 - σ) - 1) / (1 - σ), - R = setup_R!(fill(-Inf, n, a_size), a_vals, s_vals, r, w, u), - # -Inf is the utility of dying (0 consumption) - Q = setup_Q!(zeros(n, a_size, n), s_i_vals, z_chain)) +Household(;r = 0.01, + w = 1.0, + sigma = 1.0, + beta = 0.96, + z_chain = MarkovChain([0.9 0.1; 0.1 0.9], [0.1; 1.0]), + a_min = 1e-10, + a_max = 18.0, + a_size = 200, + a_vals = range(a_min, a_max, length = a_size), + z_size = length(z_chain.state_values), + n = a_size * z_size, + s_vals = gridmake(a_vals, z_chain.state_values), + s_i_vals = gridmake(1:a_size, 1:z_size), + u = sigma == 1 ? x -> log(x) : x -> (x^(1 - sigma) - 1) / (1 - sigma), + R = setup_R!(fill(-Inf, n, a_size), a_vals, s_vals, r, w, u), + # -Inf is the utility of dying (0 consumption) + Q = setup_Q!(zeros(n, a_size, n), s_i_vals, z_chain)) = (; r, w, sigma, beta, z_chain, a_min, a_max, a_size, a_vals, z_size, n, s_vals, s_i_vals, u, R, Q) function setup_Q!(Q, s_i_vals, z_chain) for next_s_i in 1:size(Q, 3) @@ -274,7 +274,7 @@ Random.seed!(42); am = Household(a_max = 20.0, r = 0.03, w = 0.956) # Use the instance to build a discrete dynamic program -am_ddp = DiscreteDP(am.R, am.Q, am.β) +am_ddp = DiscreteDP(am.R, am.Q, am.beta) # Solve using policy function iteration results = solve(am_ddp, PFI) @@ -326,16 +326,16 @@ Random.seed!(42); # Firms' parameters const A = 1 const N = 1 -const α = 0.33 -const β = 0.96 -const δ = 0.05 +const alpha = 0.33 +const beta = 0.96 +const delta = 0.05 function r_to_w(r) - return A * (1 - α) * (A * α / (r + δ)) ^ (α / (1 - α)) + return A * (1 - alpha) * (A * alpha / (r + delta)) ^ (alpha / (1 - alpha)) end function rd(K) - return A * α * (N / K) ^ (1 - α) - δ + return A * alpha * (N / K) ^ (1 - alpha) - delta end function prices_to_capital_stock(am, r) @@ -345,7 +345,7 @@ function prices_to_capital_stock(am, r) (;a_vals, s_vals, u) = am setup_R!(am.R, a_vals, s_vals, r, w, u) - aiyagari_ddp = DiscreteDP(am.R, am.Q, am.β) + aiyagari_ddp = DiscreteDP(am.R, am.Q, am.beta) # Compute the optimal policy results = solve(aiyagari_ddp, PFI) @@ -358,7 +358,7 @@ function prices_to_capital_stock(am, r) end # Create an instance of Household -am = Household(β = β, a_max = 20.0) +am = Household(beta = beta, a_max = 20.0) # Create a grid of r values at which to compute demand and supply of capital r_vals = range(0.005, 0.04, length = 20) diff --git a/lectures/multi_agent_models/arellano.md b/lectures/multi_agent_models/arellano.md index 8b0fb7cb..e9e92a40 100644 --- a/lectures/multi_agent_models/arellano.md +++ b/lectures/multi_agent_models/arellano.md @@ -312,25 +312,25 @@ using Test ``` ```{code-cell} julia -function ArellanoEconomy(;β = .953, - γ = 2., +function ArellanoEconomy(;beta = .953, + gamma = 2., r = 0.017, - ρ = 0.945, - η = 0.025, - θ = 0.282, + rho = 0.945, + eta = 0.025, + theta = 0.282, ny = 21, nB = 251) # create grids Bgrid = collect(range(-.4, .4, length = nB)) - mc = tauchen(ny, ρ, η) - Π = mc.p + mc = tauchen(ny, rho, eta) + Pi = mc.p ygrid = exp.(mc.state_values) ydefgrid = min.(.969 * mean(ygrid), ygrid) # define value functions # notice ordered different than Python to take - # advantage of column major layout of Julia) + # advantage of column major layout of Julia vf = zeros(nB, ny) vd = zeros(1, ny) vc = zeros(nB, ny) @@ -338,13 +338,12 @@ function ArellanoEconomy(;β = .953, q = ones(nB, ny) .* (1 / (1 + r)) defprob = zeros(nB, ny) - return (β = β, γ = γ, r = r, ρ = ρ, η = η, θ = θ, ny = ny, - nB = nB, ygrid = ygrid, ydefgrid = ydefgrid, - Bgrid = Bgrid, Π = Π, vf = vf, vd = vd, vc = vc, - policy = policy, q = q, defprob = defprob) + return (;beta, gamma, r, rho, eta, theta, ny, + nB, ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, + policy, q, defprob) end -u(ae, c) = c^(1 - ae.γ) / (1 - ae.γ) +u(ae, c) = c^(1 - ae.gamma) / (1 - ae.gamma) function one_step_update!(ae, EV, @@ -352,8 +351,8 @@ function one_step_update!(ae, EVc) # unpack stuff - @unpack β, γ, r, ρ, η, θ, ny, nB = ae - @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae + (;beta, gamma, r, rho, eta, theta, ny, nB) = ae + (;ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob) = ae zero_ind = searchsortedfirst(Bgrid, 0.) for iy in 1:ny @@ -361,7 +360,7 @@ function one_step_update!(ae, ydef = ae.ydefgrid[iy] # value of being in default with income y - defval = u(ae, ydef) + β * (θ * EVc[zero_ind, iy] + (1-θ) * EVd[1, iy]) + defval = u(ae, ydef) + beta * (theta * EVc[zero_ind, iy] + (1-theta) * EVd[1, iy]) ae.vd[1, iy] = defval for ib in 1:nB @@ -371,7 +370,7 @@ function one_step_update!(ae, pol_ind = 0 for ib_next=1:nB c = max(y - ae.q[ib_next, iy]*Bgrid[ib_next] + B, 1e-14) - m = u(ae, c) + β * EV[ib_next, iy] + m = u(ae, c) + beta * EV[ib_next, iy] if m > current_max current_max = m @@ -390,14 +389,14 @@ end function compute_prices!(ae) # unpack parameters - @unpack β, γ, r, ρ, η, θ, ny, nB = ae + (;beta, gamma, r, rho, eta, theta, ny, nB) = ae # create default values with a matching size vd_compat = repeat(ae.vd, nB) default_states = vd_compat .> ae.vc # update default probabilities and prices - copyto!(ae.defprob, default_states * ae.Π') + copyto!(ae.defprob, default_states * ae.Pi') copyto!(ae.q, (1 .- ae.defprob) / (1 + r)) return end @@ -405,9 +404,9 @@ end function vfi!(ae; tol = 1e-8, maxit = 10000) # unpack stuff - @unpack β, γ, r, ρ, η, θ, ny, nB = ae - @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae - Πt = Π' + (;beta, gamma, r, rho, eta, theta, ny, nB) = ae + (;ygrid, ydefgrid, Bgrid, Pi, vf, vd, vc, policy, q, defprob) = ae + Pit = Pi' # Iteration stuff it = 0 @@ -420,11 +419,11 @@ function vfi!(ae; tol = 1e-8, maxit = 10000) it += 1 # compute expectations for this iterations - # (we need Π' because of order value function dimensions) + # (we need Pi' because of order value function dimensions) copyto!(V_upd, ae.vf) - EV = ae.vf * Πt - EVd = ae.vd * Πt - EVc = ae.vc * Πt + EV = ae.vf * Pit + EVd = ae.vd * Pit + EVc = ae.vc * Pit # update value function one_step_update!(ae, EV, EVd, EVc) @@ -452,7 +451,7 @@ function QuantEcon.simulate(ae, B_init_ind = searchsortedfirst(ae.Bgrid, B_init) # create a QE MarkovChain - mc = MarkovChain(ae.Π) + mc = MarkovChain(ae.Pi) y_sim_indices = simulate(mc, capT + 1; init = y_init_ind) # allocate and fill output @@ -495,8 +494,8 @@ function QuantEcon.simulate(ae, y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]] q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]] - # with probability θ exit default status - default_status[t + 1] = rand() ≥ ae.θ + # with probability theta exit default status + default_status[t + 1] = rand() >= ae.theta end end @@ -591,12 +590,12 @@ using DataFrames, Plots Compute the value function, policy and equilibrium prices ```{code-cell} julia -ae = ArellanoEconomy(β = .953, # time discount rate - γ = 2., # risk aversion +ae = ArellanoEconomy(beta = .953, # time discount rate + gamma = 2., # risk aversion r = 0.017, # international interest rate - ρ = .945, # persistence in output - η = 0.025, # st dev of output shock - θ = 0.282, # prob of regaining access + rho = .945, # persistence in output + eta = 0.025, # st dev of output shock + theta = 0.282, # prob of regaining access ny = 21, # number of points in y grid nB = 251) # number of points in B grid @@ -627,7 +626,7 @@ q_low = zeros(0) q_high = zeros(0) for i in 1:ae.nB b = ae.Bgrid[i] - if -0.35 ≤ b ≤ 0 # to match fig 3 of Arellano + if -0.35 <= b <= 0 # to match fig 3 of Arellano push!(x, b) push!(q_low, ae.q[i, iy_low]) push!(q_high, ae.q[i, iy_high]) diff --git a/lectures/multi_agent_models/harrison_kreps.md b/lectures/multi_agent_models/harrison_kreps.md index 156ac748..d23d9997 100644 --- a/lectures/multi_agent_models/harrison_kreps.md +++ b/lectures/multi_agent_models/harrison_kreps.md @@ -290,12 +290,12 @@ Here's a function that can be used to compute these values using LinearAlgebra function price_single_beliefs(transition, dividend_payoff; - β=.75) + beta=.75) # First compute inverse piece - imbq_inv = inv(I - β * transition) + imbq_inv = inv(I - beta * transition) # Next compute prices - prices = β * ((imbq_inv * transition) * dividend_payoff) + prices = beta * ((imbq_inv * transition) * dividend_payoff) return prices end @@ -416,7 +416,7 @@ Here's code to solve for $\bar p$, $\hat p_a$ and $\hat p_b$ using the iterative ```{code-cell} julia function price_optimistic_beliefs(transitions, dividend_payoff; - β=.75, max_iter=50000, + beta=.75, max_iter=50000, tol=1e-16) # We will guess an initial price vector of [0, 0] @@ -424,11 +424,11 @@ function price_optimistic_beliefs(transitions, p_old = [10.0,10.0] # We know this is a contraction mapping, so we can iterate to conv - for i ∈ 1:max_iter + for i in 1:max_iter p_old = p_new temp = [maximum((q * p_old) + (q * dividend_payoff)) for q in transitions] - p_new = β * temp + p_new = beta * temp # If we succed in converging, break out of for loop if maximum(sqrt, ((p_new - p_old).^2)) < 1e-12 @@ -437,7 +437,7 @@ function price_optimistic_beliefs(transitions, end temp=[minimum((q * p_old) + (q * dividend_payoff)) for q in transitions] - ptwiddle = β * temp + ptwiddle = beta * temp phat_a = [p_new[1], ptwiddle[2]] phat_b = [ptwiddle[1], p_new[2]] @@ -486,17 +486,17 @@ Here's code to solve for $\check p$ using iteration ```{code-cell} julia function price_pessimistic_beliefs(transitions, dividend_payoff; - β=.75, max_iter=50000, + beta=.75, max_iter=50000, tol=1e-16) # We will guess an initial price vector of [0, 0] p_new = [0, 0] p_old = [10.0, 10.0] # We know this is a contraction mapping, so we can iterate to conv - for i ∈ 1:max_iter + for i in 1:max_iter p_old = p_new temp=[minimum((q * p_old) + (q* dividend_payoff)) for q in transitions] - p_new = β * temp + p_new = beta * temp # If we succed in converging, break out of for loop if maximum(sqrt, ((p_new - p_old).^2)) < 1e-12 @@ -594,7 +594,7 @@ heterogeneous beliefs. opt_beliefs = price_optimistic_beliefs([qa, qb], dividendreturn) labels = ["p_optimistic", "p_hat_a", "p_hat_b"] -for (p, label) ∈ zip(opt_beliefs, labels) +for (p, label) in zip(opt_beliefs, labels) println(label) println(repeat("=", 20)) s0, s1 = round.(p, digits = 2) diff --git a/lectures/multi_agent_models/lake_model.md b/lectures/multi_agent_models/lake_model.md index 1e87f825..1866c63e 100644 --- a/lectures/multi_agent_models/lake_model.md +++ b/lectures/multi_agent_models/lake_model.md @@ -213,20 +213,20 @@ using Test ``` ```{code-cell} julia -LakeModel = @with_kw (λ = 0.283, α = 0.013, b = 0.0124, d = 0.00822) +LakeModel(;lambda = 0.283, alpha = 0.013, b = 0.0124, d = 0.00822) = (;lambda, alpha, b, d) function transition_matrices(lm) - (;λ, α, b, d) = lm + (;lambda, alpha, b, d) = lm g = b - d - A = [(1 - λ) * (1 - d) + b (1 - d) * α + b - (1 - d) * λ (1 - d) * (1 - α)] - Â = A ./ (1 + g) - return (A = A, Â = Â) + A = [(1 - lambda) * (1 - d) + b (1 - d) * alpha + b + (1 - d) * lambda (1 - d) * (1 - alpha)] + A_hat = A ./ (1 + g) + return (A, A_hat) end function rate_steady_state(lm) - (;Â) = transition_matrices(lm) - sol = fixedpoint(x -> Â * x, fill(0.5, 2)) + (;A_hat) = transition_matrices(lm) + sol = fixedpoint(x -> A_hat * x, fill(0.5, 2)) converged(sol) || error("Failed to converge in $(result.iterations) iterations") return sol.zero end @@ -243,12 +243,12 @@ function simulate_stock_path(lm, X0, T) end function simulate_rate_path(lm, x0, T) - (;Â) = transition_matrices(lm) + (;A_hat) = transition_matrices(lm) x_path = zeros(eltype(x0), 2, T) x = copy(x0) for t in 1:T x_path[:, t] = x - x = Â * x + x = A_hat * x end return x_path end @@ -258,24 +258,24 @@ Let's observe these matrices for the baseline model ```{code-cell} julia lm = LakeModel() -A, Â = transition_matrices(lm) +A, A_hat = transition_matrices(lm) A ``` ```{code-cell} julia -Â +A_hat ``` And a revised model ```{code-cell} julia -lm = LakeModel(α = 2.0) -A, Â = transition_matrices(lm) +lm = LakeModel(alpha = 2.0) +A, A_hat = transition_matrices(lm) A ``` ```{code-cell} julia -Â +A_hat ``` ```{code-cell} julia @@ -283,7 +283,7 @@ Â tags: [remove-cell] --- @testset begin - @test lm.α ≈ 2.0 + @test lm.alpha ≈ 2.0 @test A[1][1] ≈ 0.7235062600000001 end ``` @@ -343,8 +343,8 @@ This is the case for our default parameters: ```{code-cell} julia lm = LakeModel() -A, Â = transition_matrices(lm) -e, f = eigvals(Â) +A, A_hat = transition_matrices(lm) +e, f = eigvals(A_hat) abs(e), abs(f) ``` @@ -487,9 +487,9 @@ using QuantEcon, Roots, Random lm = LakeModel(d = 0, b = 0) T = 5000 # Simulation length -(;α, λ) = lm -P = [(1 - λ) λ - α (1 - α)] +(;alpha, lambda) = lm +P = [(1 - lambda) lambda + alpha (1 - alpha)] ``` ```{code-cell} julia @@ -498,9 +498,9 @@ mc = MarkovChain(P, [0; 1]) # 0=unemployed, 1=employed xbar = rate_steady_state(lm) s_path = simulate(mc, T; init=2) -s̄_e = cumsum(s_path) ./ (1:T) -s̄_u = 1 .- s̄_e -s_bars = [s̄_u s̄_e] +s_bar_e = cumsum(s_path) ./ (1:T) +s_bar_u = 1 .- s_bar_e +s_bars = [s_bar_u s_bar_e] plt_unemp = plot(title = "Percent of time unemployed", 1:T, s_bars[:,1],color = :blue, lw = 2, alpha = 0.5, label = "", grid = true) @@ -630,19 +630,18 @@ Following {cite}`davis2006flow`, we set $\alpha$, the hazard rate of leaving emp We will make use of (with some tweaks) the code we wrote in the {doc}`McCall model lecture <../dynamic_programming/mccall_model>`, embedded below for convenience. ```{code-cell} julia -function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1e-5, - iter = 2_000) - (;α, β, σ, c, γ, w, E, u) = mcm +function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1e-5,iter = 2_000) + (;alpha, beta, sigma, c, gamma, w, E, u) = mcm # necessary objects - u_w = u.(w, σ) - u_c = u(c, σ) + u_w = u.(w, sigma) + u_c = u(c, sigma) # Bellman operator T. Fixed point is x* s.t. T(x*) = x* function T(x) V = x[1:end-1] U = x[end] - [u_w + β * ((1 - α) * V .+ α * U); u_c + β * (1 - γ) * U + β * γ * E * max.(U, V)] + [u_w + beta * ((1 - alpha) * V .+ alpha * U); u_c + beta * (1 - gamma) * U + beta * gamma * E * max.(U, V)] end # value function iteration @@ -654,13 +653,13 @@ function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)), tol = 1 # compute the reservation wage w_barindex = searchsortedfirst(V .- U, 0.0) if w_barindex >= length(w) # if this is true, you never want to accept - w̄ = Inf + w_bar = Inf else - w̄ = w[w_barindex] # otherwise, return the number + w_bar = w[w_barindex] # otherwise, return the number end # return a NamedTuple, so we can select values by name - return (V = V, U = U, w̄ = w̄) + return (;V, U, w_bar) end ``` @@ -668,17 +667,18 @@ And the McCall object ```{code-cell} julia # a default utility function -u(c, σ) = c > 0 ? (c^(1 - σ) - 1) / (1 - σ) : -10e-6 +u(c, sigma) = c > 0 ? (c^(1 - sigma) - 1) / (1 - sigma) : -10e-6 # model constructor -McCallModel = @with_kw (α = 0.2, - β = 0.98, # discount rate - γ = 0.7, - c = 6.0, # unemployment compensation - σ = 2.0, - u = u, # utility function - w = range(10, 20, length = 60), # wage values - E = Expectation(BetaBinomial(59, 600, 400))) # distribution over wage values +McCallModel(;alpha = 0.2, + beta = 0.98, # discount rate + gamma = 0.7, + c = 6.0, # unemployment compensation + sigma = 2.0, + u = u, # utility function + w = range(10, 20, length = 60), # w= wage values; E= distribution over wage values + E = Expectation(BetaBinomial(59, 600, 400))) = (;alpha, beta, gamma, c, sigma, u, w, E) + ``` Now let's compute and plot welfare, employment, unemployment, and tax revenue as a @@ -686,13 +686,13 @@ function of the unemployment compensation rate ```{code-cell} julia # some global variables that will stay constant -α = 0.013 -α_q = (1 - (1 - α)^3) +alpha = 0.013 +alpha_q = (1 - (1 - alpha)^3) b_param = 0.0124 d_param = 0.00822 -β = 0.98 -γ = 1.0 -σ = 2.0 +beta = 0.98 +gamma = 1.0 +sigma = 2.0 # the default wage distribution: a discretized log normal log_wage_mean, wage_grid_size, max_wage = 20, 200, 170 @@ -707,33 +707,33 @@ w_vec = (w_vec[1:end-1] + w_vec[2:end]) / 2 E = expectation(Categorical(p_vec)) # expectation object -function compute_optimal_quantities(c, τ) - mcm = McCallModel(α = α_q, - β = β, - γ = γ, - c = c - τ, # post-tax compensation - σ = σ, - w = w_vec .- τ, # post-tax wages +function compute_optimal_quantities(c, tau) + mcm = McCallModel(alpha = alpha_q, + beta = beta, + gamma = gamma, + c = c - tau, # post-tax compensation + sigma = sigma, + w = w_vec .- tau, # post-tax wages E = E) # expectation operator - (;V, U, w̄) = solve_mccall_model(mcm) - indicator = wage -> wage > w̄ - λ = γ * E * indicator.(w_vec .- τ) + (;V, U, w_bar) = solve_mccall_model(mcm) + indicator = wage -> wage > w_bar + lambda = gamma * E * indicator.(w_vec .- tau) - return w̄, λ, V, U + return w_bar, lambda, V, U end -function compute_steady_state_quantities(c, τ) - w̄, λ_param, V, U = compute_optimal_quantities(c, τ) +function compute_steady_state_quantities(c, tau) + w_bar, lambda_param, V, U = compute_optimal_quantities(c, tau) # compute steady state employment and unemployment rates - lm = LakeModel(λ = λ_param, α = α_q, b = b_param, d = d_param) + lm = LakeModel(lambda = lambda_param, alpha = alpha_q, b = b_param, d = d_param) x = rate_steady_state(lm) u_rate, e_rate = x # compute steady state welfare - indicator(wage) = wage > w̄ - decisions = indicator.(w_vec .- τ) + indicator(wage) = wage > w_bar + decisions = indicator.(w_vec .- tau) w = (E * (V .* decisions)) / (E * decisions) welfare = e_rate .* w + u_rate .* U @@ -746,8 +746,8 @@ function find_balanced_budget_tax(c) return t - u_rate * c end - τ = find_zero(steady_state_budget, (0.0, 0.9c)) - return τ + tau = find_zero(steady_state_budget, (0.0, 0.9c)) + return tau end # levels of unemployment insurance we wish to study @@ -844,7 +844,7 @@ T = 50 New legislation changes $\lambda$ to $0.2$ ```{code-cell} julia -lm = LakeModel(λ = 0.2) +lm = LakeModel(lambda = 0.2) ``` ```{code-cell} julia @@ -936,16 +936,16 @@ end Here are the other parameters: ```{code-cell} julia -b̂ = 0.003 -T̂ = 20 +b_hat = 0.003 +T_hat = 20 ``` Let's increase $b$ to the new value and simulate for 20 periods ```{code-cell} julia -lm = LakeModel(b=b̂) -X_path1 = simulate_stock_path(lm, x0 * N0, T̂) # simulate stocks -x_path1 = simulate_rate_path(lm, x0, T̂) # simulate rates +lm = LakeModel(b=b_hat) +X_path1 = simulate_stock_path(lm, x0 * N0, T_hat) # simulate stocks +x_path1 = simulate_rate_path(lm, x0, T_hat) # simulate rates ``` Now we reset $b$ to the original value and then, using the state @@ -954,8 +954,8 @@ additional 30 periods ```{code-cell} julia lm = LakeModel(b = 0.0124) -X_path2 = simulate_stock_path(lm, X_path1[:, end-1], T-T̂+1) # simulate stocks -x_path2 = simulate_rate_path(lm, x_path1[:, end-1], T-T̂+1) # simulate rates +X_path2 = simulate_stock_path(lm, X_path1[:, end-1], T-T_hat+1) # simulate stocks +x_path2 = simulate_rate_path(lm, x_path1[:, end-1], T-T_hat+1) # simulate rates ``` Finally we combine these two paths and plot diff --git a/lectures/multi_agent_models/lucas_model.md b/lectures/multi_agent_models/lucas_model.md index b5239a57..2b26c385 100644 --- a/lectures/multi_agent_models/lucas_model.md +++ b/lectures/multi_agent_models/lucas_model.md @@ -393,39 +393,39 @@ using Distributions, Interpolations, LaTeXStrings, Parameters, Plots, Random ```{code-cell} julia # model -function LucasTree(;γ = 2.0, - β = 0.95, - α = 0.9, - σ = 0.1, +function LucasTree(;gamma = 2.0, + beta = 0.95, + alpha = 0.9, + sigma = 0.1, grid_size = 100) - ϕ = LogNormal(0.0, σ) - shocks = rand(ϕ, 500) + phi = LogNormal(0.0, sigma) + shocks = rand(phi, 500) # build a grid with mass around stationary distribution - ssd = σ / sqrt(1 - α^2) + ssd = sigma / sqrt(1 - alpha^2) grid_min, grid_max = exp(-4ssd), exp(4ssd) grid = range(grid_min, grid_max, length = grid_size) - # set h(y) = β * int u'(G(y,z)) G(y,z) ϕ(dz) + # set h(y) = beta * int u'(G(y,z)) G(y,z) phi(dz) h = similar(grid) for (i, y) in enumerate(grid) - h[i] = β * mean((y^α .* shocks).^(1 - γ)) + h[i] = beta * mean((y^alpha .* shocks).^(1 - gamma)) end - return (γ = γ, β = β, α = α, σ = σ, ϕ = ϕ, grid = grid, shocks = shocks, h = h) + return (gamma = gamma, beta = beta, alpha = alpha, sigma = sigma, phi = phi, grid = grid, shocks = shocks, h = h) end # approximate Lucas operator, which returns the updated function Tf on the grid function lucas_operator(lt, f) # unpack input - (;grid, α, β, h) = lt + (;grid, alpha, beta, h) = lt z = lt.shocks - Af = LinearInterpolation(grid, f, extrapolation_bc=Line()) + Af = linear_interpolation(grid, f, extrapolation_bc=Line()) - Tf = [ h[i] + β * mean(Af.(grid[i]^α .* z)) for i in 1:length(grid) ] + Tf = [ h[i] + beta * mean(Af.(grid[i]^alpha .* z)) for i in 1:length(grid) ] return Tf end @@ -434,7 +434,7 @@ function solve_lucas_model(lt; tol = 1e-6, max_iter = 500) - (;grid, γ) = lt + (;grid, gamma) = lt i = 0 f = zero(grid) # Initial guess of f @@ -447,8 +447,8 @@ function solve_lucas_model(lt; i += 1 end - # p(y) = f(y) * y ^ γ - price = f .* grid.^γ + # p(y) = f(y) * y ^ gamma + price = f .* grid.^gamma return price end @@ -459,7 +459,7 @@ An example of usage is given in the docstring and repeated here ```{code-cell} julia Random.seed!(42) # For reproducible results. -tree = LucasTree(γ = 2.0, β = 0.95, α = 0.90, σ = 0.1) +tree = LucasTree(gamma = 2.0, beta = 0.95, alpha = 0.90, sigma = 0.1) price_vals = solve_lucas_model(tree); ``` @@ -519,11 +519,11 @@ Random.seed!(42); ```{code-cell} julia plot() -for β in (.95, 0.98) - tree = LucasTree(;β = β) +for beta in (.95, 0.98) + tree = LucasTree(;beta = beta) grid = tree.grid price_vals = solve_lucas_model(tree) - plot!(grid, price_vals, lw = 2, label = L"\beta = %$β") + plot!(grid, price_vals, lw = 2, label = L"\beta = %$beta") end plot!(xlabel = L"y", ylabel = "price", legend = :topleft) @@ -535,7 +535,7 @@ tags: [remove-cell] --- @testset begin # For the 0.98, since the other one is overwritten. Random.seed!(42) - price_vals = solve_lucas_model(LucasTree(β = 0.98)) + price_vals = solve_lucas_model(LucasTree(beta = 0.98)) #test price_vals[20] ≈ 35.00073581199659 #test price_vals[57] ≈ 124.32987344509688 end diff --git a/lectures/multi_agent_models/markov_asset.md b/lectures/multi_agent_models/markov_asset.md index 27bc9878..821c702c 100644 --- a/lectures/multi_agent_models/markov_asset.md +++ b/lectures/multi_agent_models/markov_asset.md @@ -404,12 +404,12 @@ Here's the code, including a test of the spectral radius condition ```{code-cell} julia n = 25 # size of state space -β = 0.9 +beta = 0.9 mc = tauchen(n, 0.96, 0.02) K = mc.p .* exp.(mc.state_values)' -v = (I - β * K) \ (β * K * ones(n, 1)) +v = (I - beta * K) \ (beta * K * ones(n, 1)) plot(mc.state_values, v, @@ -541,39 +541,35 @@ the AssetPriceModel objects ```{code-cell} julia # A default Markov chain for the state process -ρ = 0.9 -σ = 0.02 +rho = 0.9 +sigma = 0.02 n = 25 -default_mc = tauchen(n, ρ, σ) +default_mc = tauchen(n, rho, sigma) -AssetPriceModel = @with_kw (β = 0.96, - γ = 2.0, - mc = default_mc, - n = size(mc.p)[1], - g = exp) +AssetPriceModel(;beta = 0.96, gamma = 2.0, mc = default_mc, n = size(mc.p)[1], g = exp) = (;beta, gamma, mc, n, g) # test stability of matrix Q function test_stability(ap, Q) sr = maximum(abs, eigvals(Q)) - if sr ≥ 1 / ap.β + if sr >= 1 / ap.beta msg = "Spectral radius condition failed with radius = $sr" throw(ArgumentError(msg)) end end # price/dividend ratio of the Lucas tree -function tree_price(ap; γ = ap.γ) +function tree_price(ap; gamma = ap.gamma) # Simplify names, set up matrices - (;β, mc) = ap + (;beta, mc) = ap P, y = mc.p, mc.state_values y = reshape(y, 1, ap.n) - J = P .* ap.g.(y).^(1 - γ) + J = P .* ap.g.(y).^(1 - gamma) # Make sure that a unique solution exists test_stability(ap, J) # Compute v - v = (I - β * J) \ sum(β * J, dims = 2) + v = (I - beta * J) \ sum(beta * J, dims = 2) return v end @@ -583,16 +579,16 @@ Here's a plot of $v$ as a function of the state for several values of $\gamma$, with a positively correlated Markov process and $g(x) = \exp(x)$ ```{code-cell} julia -γs = [1.2, 1.4, 1.6, 1.8, 2.0] +gammas = [1.2, 1.4, 1.6, 1.8, 2.0] ap = AssetPriceModel() states = ap.mc.state_values lines = [] labels = [] -for γ in γs - v = tree_price(ap, γ = γ) - label = L"\gamma = %$γ" +for gamma in gammas + v = tree_price(ap, gamma = gamma) + label = L"\gamma = %$gamma" push!(labels, label) push!(lines, v) end @@ -690,18 +686,18 @@ p = (I - \beta M)^{-1} \beta M \zeta {\mathbb 1} The above is implemented in the function consol_price ```{code-cell} julia -function consol_price(ap, ζ) +function consol_price(ap, zeta) # Simplify names, set up matrices - (;β, γ, mc, g, n) = ap + (;beta, gamma, mc, g, n) = ap P, y = mc.p, mc.state_values y = reshape(y, 1, n) - M = P .* g.(y).^(-γ) + M = P .* g.(y).^(-gamma) # Make sure that a unique solution exists test_stability(ap, M) # Compute price - return (I - β * M) \ sum(β * ζ * M, dims = 2) + return (I - beta * M) \ sum(beta * zeta * M, dims = 2) end ``` @@ -777,24 +773,24 @@ We can find the solution with the following function call_option ```{code-cell} julia # price of perpetual call on consol bond -function call_option(ap, ζ, p_s, ϵ = 1e-7) +function call_option(ap, zeta, p_s, epsilon = 1e-7) # Simplify names, set up matrices - (;β, γ, mc, g, n) = ap + (;beta, gamma, mc, g, n) = ap P, y = mc.p, mc.state_values y = reshape(y, 1, n) - M = P .* g.(y).^(-γ) + M = P .* g.(y).^(-gamma) # Make sure that a unique console price exists test_stability(ap, M) # Compute option price - p = consol_price(ap, ζ) + p = consol_price(ap, zeta) w = zeros(ap.n, 1) - error = ϵ + 1 - while (error > ϵ) + error = epsilon + 1 + while (error > epsilon) # Maximize across columns - w_new = max.(β * M * w, p .- p_s) + w_new = max.(beta * M * w, p .- p_s) # Find maximal difference of each component and update error = maximum(abs, w - w_new) w = w_new @@ -807,13 +803,13 @@ end Here's a plot of $w$ compared to the consol price when $P_S = 40$ ```{code-cell} julia -ap = AssetPriceModel(β=0.9) -ζ = 1.0 +ap = AssetPriceModel(beta=0.9) +zeta = 1.0 strike_price = 40.0 x = ap.mc.state_values -p = consol_price(ap, ζ) -w = call_option(ap, ζ, strike_price) +p = consol_price(ap, zeta) +w = call_option(ap, zeta, strike_price) plot(x, p, color = "blue", lw = 2, xlabel = "state", label = "consol price") plot!(x, w, color = "green", lw = 2, label = "value of call option") @@ -899,9 +895,9 @@ Consider the following primitives n = 5 P = fill(0.0125, n, n) + (0.95 - 0.0125)I s = [1.05, 1.025, 1.0, 0.975, 0.95] -γ = 2.0 -β = 0.94 -ζ = 1.0 +gamma = 2.0 +beta = 0.94 +zeta = 1.0 ``` Let $g$ be defined by $g(x) = x$ (that is, $g$ is the identity map). @@ -965,16 +961,16 @@ P = fill(0.0125, n, n) + (0.95 - 0.0125)I s = [0.95, 0.975, 1.0, 1.025, 1.05] # state values mc = MarkovChain(P, s) -γ = 2.0 -β = 0.94 -ζ = 1.0 +gamma = 2.0 +beta = 0.94 +zeta = 1.0 p_s = 150.0 ``` Next we'll create an instance of AssetPriceModel to feed into the functions. ```{code-cell} julia -ap = AssetPriceModel(β = β, mc = mc, γ = γ, g = x -> x) +ap = AssetPriceModel(beta = beta, mc = mc, gamma = gamma, g = x -> x) ``` ```{code-cell} julia @@ -997,7 +993,7 @@ println("Consol Bond Prices: $(v_consol)\n") ``` ```{code-cell} julia -w = call_option(ap, ζ, p_s) +w = call_option(ap, zeta, p_s) ``` ```{code-cell} julia @@ -1015,23 +1011,23 @@ end Here's a suitable function: ```{code-cell} julia -function finite_horizon_call_option(ap, ζ, p_s, k) +function finite_horizon_call_option(ap, zeta, p_s, k) # Simplify names, set up matrices - (;β, γ, mc) = ap + (;beta, gamma, mc) = ap P, y = mc.p, mc.state_values y = y' - M = P .* ap.g.(y).^(- γ) + M = P .* ap.g.(y).^(- gamma) # Make sure that a unique console price exists test_stability(ap, M) # Compute option price - p = consol_price(ap, ζ) + p = consol_price(ap, zeta) w = zeros(ap.n, 1) for i in 1:k # Maximize across columns - w = max.(β * M * w, p .- p_s) + w = max.(beta * M * w, p .- p_s) end return w @@ -1052,7 +1048,7 @@ end lines = [] labels = [] for k in [5, 25] - w = finite_horizon_call_option(ap, ζ, p_s, k) + w = finite_horizon_call_option(ap, zeta, p_s, k) push!(lines, w) push!(labels, L"k = %$k") end diff --git a/lectures/multi_agent_models/markov_perf.md b/lectures/multi_agent_models/markov_perf.md index a7b61d10..06fcfcf9 100644 --- a/lectures/multi_agent_models/markov_perf.md +++ b/lectures/multi_agent_models/markov_perf.md @@ -432,8 +432,8 @@ using QuantEcon, LinearAlgebra # parameters a0 = 10.0 a1 = 2.0 -β = 0.96 -γ = 12.0 +beta = 0.96 +gamma = 12.0 # in LQ form A = I + zeros(3, 3) @@ -448,12 +448,12 @@ R2 = [ 0.0 0.0 -a0 / 2.0; 0.0 0.0 a1 / 2.0; -a0 / 2.0 a1 / 2.0 a1] -Q1 = Q2 = γ +Q1 = Q2 = gamma S1 = S2 = W1 = W2 = M1 = M2 = 0.0 # solve using QE's nnash function F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2, - beta=β) + beta=beta) # display policies println("Computed policies for firm 1 and firm 2:") @@ -480,8 +480,8 @@ In particular, let's take F2 as computed above, plug it into {eq}`eq_mpe_p1p` an We hope that the resulting policy will agree with F1 as computed above ```{code-cell} julia -Λ1 = A - (B2 * F2) -lq1 = QuantEcon.LQ(Q1, R1, Λ1, B1, bet=β) +Lamda1 = A - (B2 * F2) +lq1 = QuantEcon.LQ(Q1, R1, Lamda1, B1, bet=beta) P1_ih, F1_ih, d = stationary_values(lq1) F1_ih ``` @@ -493,7 +493,7 @@ tags: [remove-cell] @testset begin @test P1_ih[2, 2] ≈ 5.441368459897164 @test d ≈ 0.0 - @test Λ1[1, 1] ≈ 1.0 && Λ1[3, 2] ≈ -0.07584666305807419 + @test Lamda1[1, 1] ≈ 1.0 && Lamda1[3, 2] ≈ -0.07584666305807419 @test F1_ih ≈ [-0.6684661291052371 0.29512481789806305 0.07584666292394007] @test isapprox(F1, F1_ih, atol=1e-7) # Make sure the test below comes up true. end @@ -651,7 +651,7 @@ The exercise is to calculate these matrices and compute the following figures. The first figure shows the dynamics of inventories for each firm when the parameters are ```{code-cell} julia -δ = 0.02 +delta = 0.02 D = [ -1 0.5; 0.5 -1] b = [25, 25] @@ -683,8 +683,8 @@ First let's compute the duopoly MPE under the stated parameters # parameters a0 = 10.0 a1 = 2.0 -β = 0.96 -γ = 12.0 +beta = 0.96 +gamma = 12.0 # in LQ form A = I + zeros(3, 3) @@ -699,12 +699,12 @@ R2 = [ 0.0 0.0 -a0 / 2.0; 0.0 0.0 a1 / 2.0; -a0 / 2.0 a1 / 2.0 a1] -Q1 = Q2 = γ +Q1 = Q2 = gamma S1 = S2 = W1 = W2 = M1 = M2 = 0.0 # solve using QE's nnash function F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2, - beta=β) + beta=beta) ``` ```{code-cell} julia @@ -783,9 +783,9 @@ resulting dynamics of $\{q_t\}$, starting at $q_0 = 2.0$. tags: [hide-output] --- R = a1 -Q = γ +Q = gamma A = B = 1 -lq_alt = QuantEcon.LQ(Q, R, A, B, bet=β) +lq_alt = QuantEcon.LQ(Q, R, A, B, bet=beta) P, F, d = stationary_values(lq_alt) q̄ = a0 / (2.0 * a1) qm = zeros(n) @@ -830,13 +830,13 @@ plot(plt_q, plt_p, layout=(2,1), size=(700,600)) We treat the case $\delta = 0.02$ ```{code-cell} julia -δ = 0.02 +delta = 0.02 D = [-1 0.5; 0.5 -1] b = [25, 25] c1 = c2 = [1, -2, 1] e1 = e2 = [10, 10, 3] -δ_1 = 1-δ +delta_1 = 1-delta ``` Recalling that the control and state are @@ -860,14 +860,14 @@ we set up the matrices as follows: ```{code-cell} julia # create matrices needed to compute the Nash feedback equilibrium -A = [δ_1 0 -δ_1 * b[1]; - 0 δ_1 -δ_1 * b[2]; +A = [delta_1 0 -delta_1 * b[1]; + 0 delta_1 -delta_1 * b[2]; 0 0 1] -B1 = δ_1 * [1 -D[1, 1]; +B1 = delta_1 * [1 -D[1, 1]; 0 -D[2, 1]; 0 0] -B2 = δ_1 * [0 -D[1, 2]; +B2 = delta_1 * [0 -D[1, 2]; 1 -D[2, 2]; 0 0] diff --git a/lectures/multi_agent_models/matsuyama.md b/lectures/multi_agent_models/matsuyama.md index 3153335a..88ee3b50 100644 --- a/lectures/multi_agent_models/matsuyama.md +++ b/lectures/multi_agent_models/matsuyama.md @@ -342,7 +342,7 @@ using Plots, Parameters ``` ```{code-cell} julia -function h_j(j, nk, s1, s2, θ, δ, ρ) +function h_j(j, nk, s1, s2, theta, delta, rho) # Find out who's h we are evaluating if j == 1 sj = s1 @@ -354,8 +354,8 @@ function h_j(j, nk, s1, s2, θ, δ, ρ) # Coefficients on the quadratic a x^2 + b x + c = 0 a = 1.0 - b = ((ρ + 1 / ρ) * nk - sj - sk) - c = (nk * nk - (sj * nk) / ρ - sk * ρ * nk) + b = ((rho + 1 / rho) * nk - sj - sk) + c = (nk * nk - (sj * nk) / rho - sk * rho * nk) # Positive solution of quadratic form root = (-b + sqrt(b * b - 4 * a * c)) / (2 * a) @@ -363,41 +363,41 @@ function h_j(j, nk, s1, s2, θ, δ, ρ) return root end -DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≤ s1_ρ) && (n2 ≤ s2_ρ) +DLL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 <= s1_rho) && (n2 <= s2_rho) -DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≥ h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 ≥ h_j(2, n1, s1, s2, θ, δ, ρ)) +DHH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 >= h_j(1, n2, s1, s2, theta, delta, rho)) && (n2 >= h_j(2, n1, s1, s2, theta, delta, rho)) -DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≥ s1_ρ) && (n2 ≤ h_j(2, n1, s1, s2, θ, δ, ρ)) +DHL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 >= s1_rho) && (n2 <= h_j(2, n1, s1, s2, theta, delta, rho)) -DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≤ h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 ≥ s2_ρ) +DLH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + (n1 <= h_j(1, n2, s1, s2, theta, delta, rho)) && (n2 >= s2_rho) -function one_step(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) +function one_step(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) # Depending on where we are, evaluate the right branch - if DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * (θ * s1_ρ + (1 - θ) * n1) - n2_tp1 = δ * (θ * s2_ρ + (1 - θ) * n2) - elseif DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * n1 - n2_tp1 = δ * n2 - elseif DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * n1 - n2_tp1 = δ * (θ * h_j(2, n1, s1, s2, θ, δ, ρ) + (1 - θ) * n2) - elseif DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) - n1_tp1 = δ * (θ * h_j(1, n2, s1, s2, θ, δ, ρ) + (1 - θ) * n1) - n2_tp1 = δ * n2 + if DLL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * (theta * s1_rho + (1 - theta) * n1) + n2_tp1 = delta * (theta * s2_rho + (1 - theta) * n2) + elseif DHH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * n1 + n2_tp1 = delta * n2 + elseif DHL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * n1 + n2_tp1 = delta * (theta * h_j(2, n1, s1, s2, theta, delta, rho) + (1 - theta) * n2) + elseif DLH(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + n1_tp1 = delta * (theta * h_j(1, n2, s1, s2, theta, delta, rho) + (1 - theta) * n1) + n2_tp1 = delta * n2 end return n1_tp1, n2_tp1 end -new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - one_step(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) +new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) = + one_step(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) -function pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, +function pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers) # Initialize the status of synchronization @@ -407,11 +407,11 @@ function pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, nsync = 0 - while (~synchronized) && (iters < maxiter) + while (!synchronized) && (iters < maxiter) # Increment the number of iterations and get next values iters += 1 - n1_t, n2_t = new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) + n1_t, n2_t = new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) # Check whether same in this period if abs(n1_t - n2_t) < 1e-8 @@ -432,7 +432,7 @@ function pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, return synchronized, pers_2_sync end -function create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, +function create_attraction_basis(s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers, npts) # Create unit range with npts synchronized, pers_2_sync = false, 0 @@ -443,9 +443,7 @@ function create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, # Iterate over initial conditions for (i, n1_0) in enumerate(unit_range) for (j, n2_0) in enumerate(unit_range) - synchronized, pers_2_sync = pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, - s1, s2, θ, δ, ρ, - maxiter, npers) + synchronized, pers_2_sync = pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho, maxiter, npers) time_2_sync[i, j] = pers_2_sync end end @@ -453,18 +451,18 @@ function create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, ρ, end # model -function MSGSync(s1 = 0.5, θ = 2.5, δ = 0.7, ρ = 0.2) +function MSGSync(s1 = 0.5, theta = 2.5, delta = 0.7, rho = 0.2) # Store other cutoffs and parameters we use s2 = 1 - s1 - s1_ρ = min((s1 - ρ * s2) / (1 - ρ), 1) - s2_ρ = 1 - s1_ρ + s1_rho = min((s1 - rho * s2) / (1 - rho), 1) + s2_rho = 1 - s1_rho - return (s1 = s1, s2 = s2, s1_ρ = s1_ρ, s2_ρ = s2_ρ, θ = θ, δ = δ, ρ = ρ) + return (;s1, s2, s1_rho, s2_rho, theta, delta, rho) end function simulate_n(model, n1_0, n2_0, T) # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model + (;s1, s2, theta, delta, rho, s1_rho, s2_rho) = model # Allocate space n1 = zeros(T) @@ -474,7 +472,7 @@ function simulate_n(model, n1_0, n2_0, T) for t in 1:T # Get next values n1[t], n2[t] = n1_0, n2_0 - n1_0, n2_0 = new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) + n1_0, n2_0 = new_n1n2(n1_0, n2_0, s1_rho, s2_rho, s1, s2, theta, delta, rho) end return n1, n2 @@ -483,10 +481,10 @@ end function pers_till_sync(model, n1_0, n2_0, maxiter = 500, npers = 3) # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model + (;s1, s2, theta, delta, rho, s1_rho, s2_rho) = model - return pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, - θ, δ, ρ, maxiter, npers) + return pers_till_sync(n1_0, n2_0, s1_rho, s2_rho, s1, s2, + theta, delta, rho, maxiter, npers) end function create_attraction_basis(model; @@ -494,10 +492,10 @@ function create_attraction_basis(model; npers = 3, npts = 50) # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model + (;s1, s2, theta, delta, rho, s1_rho, s2_rho) = model - ab = create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ, - ρ, maxiter, npers, npts) + ab = create_attraction_basis(s1_rho, s2_rho, s1, s2, theta, delta, + rho, maxiter, npers, npts) return ab end ``` @@ -513,8 +511,8 @@ The time series share parameters but differ in their initial condition. Here's the function ```{code-cell} julia -function plot_timeseries(n1_0, n2_0, s1 = 0.5, θ = 2.5, δ = 0.7, ρ = 0.2) - model = MSGSync(s1, θ, δ, ρ) +function plot_timeseries(n1_0, n2_0, s1 = 0.5, theta = 2.5, delta = 0.7, rho = 0.2) + model = MSGSync(s1, theta, delta, rho) n1, n2 = simulate_n(model, n1_0, n2_0, 25) return [n1 n2] end @@ -569,10 +567,10 @@ Replicate the figure {ref}`shown above ` by coloring initial conditions ### Exercise 1 ```{code-cell} julia -function plot_attraction_basis(s1 = 0.5, θ = 2.5, δ = 0.7, ρ = 0.2; npts = 250) +function plot_attraction_basis(s1 = 0.5, theta = 2.5, delta = 0.7, rho = 0.2; npts = 250) # Create attraction basis unitrange = range(0, 1, length = npts) - model = MSGSync(s1, θ, δ, ρ) + model = MSGSync(s1, theta, delta, rho) ab = create_attraction_basis(model,npts=npts) plt = Plots.heatmap(ab, legend = false) end @@ -593,13 +591,13 @@ plot(plots..., size = (1000, 1000)) tags: [remove-cell] --- function plot_attraction_basis(s1 = 0.5, - θ = 2.5, - δ = 0.7, - ρ = 0.2; + theta = 2.5, + delta = 0.7, + rho = 0.2; npts = 250) # Create attraction basis unitrange = range(0, 1, length = npts) - model = MSGSync(s1, θ, δ, ρ) + model = MSGSync(s1, theta, delta, rho) ab = create_attraction_basis(model, npts = npts) end diff --git a/lectures/multi_agent_models/rational_expectations.md b/lectures/multi_agent_models/rational_expectations.md index 4c3d9fe1..2215b90e 100644 --- a/lectures/multi_agent_models/rational_expectations.md +++ b/lectures/multi_agent_models/rational_expectations.md @@ -662,16 +662,16 @@ Here's our solution # model parameters a0 = 100 a1 = 0.05 -β = 0.95 -γ = 10.0 +beta = 0.95 +gamma = 10.0 # beliefs -κ0 = 95.5 -κ1 = 0.95 +kappa0 = 95.5 +kappa1 = 0.95 # formulate the LQ problem A = [1 0 0 - 0 κ1 κ0 + 0 kappa1 kappa0 0 0 1] B = [1.0, 0.0, 0.0] @@ -680,10 +680,10 @@ R = [ 0 a1/2 -a0/2 a1/2 0 0 -a0/2 0 0] -Q = 0.5 * γ +Q = 0.5 * gamma # solve for the optimal policy -lq = QuantEcon.LQ(Q, R, A, B; bet = β) +lq = QuantEcon.LQ(Q, R, A, B; bet = beta) P, F, d = stationary_values(lq) hh = h0, h1, h2 = -F[3], 1 - F[1], -F[2] @@ -751,7 +751,7 @@ for (k0, k1) in candidates A = [1 0 0 0 k1 k0 0 0 1] - lq = QuantEcon.LQ(Q, R, A, B; bet=β) + lq = QuantEcon.LQ(Q, R, A, B; bet=beta) P, F, d = stationary_values(lq) hh = h0, h1, h2 = -F[3], 1 - F[1], -F[2] @@ -833,15 +833,15 @@ B = [1.0, 0.0] R = [ a1 / 2.0 -a0 / 2.0 -a0 / 2.0 0.0] -Q = γ / 2.0 +Q = gamma / 2.0 # solve for the optimal policy -lq = QuantEcon.LQ(Q, R, A, B; bet=β) +lq = QuantEcon.LQ(Q, R, A, B; bet=beta) P, F, d = stationary_values(lq) # print the results -κ0, κ1 = -F[2], 1 - F[1] -println("κ0=$κ0\tκ1=$κ1") +kappa0, kappa1 = -F[2], 1 - F[1] +println("kappa0=$kappa0\tkappa1=$kappa1") ``` ```{code-cell} julia @@ -849,8 +849,8 @@ println("κ0=$κ0\tκ1=$κ1") tags: [remove-cell] --- @testset begin - @test κ0 ≈ 95.081845248600 rtol = 4 - @test κ1 ≈ 0.952459076301 rtol = 4 + @test kappa0 ≈ 95.081845248600 rtol = 4 + @test kappa1 ≈ 0.952459076301 rtol = 4 end ``` @@ -879,10 +879,10 @@ B = [1.0, 0.0] R = [ a1 -a0 / 2.0 -a0 / 2.0 0.0] -Q = γ / 2.0 +Q = gamma / 2.0 # solve for the optimal policy -lq = QuantEcon.LQ(Q, R, A, B; bet=β) +lq = QuantEcon.LQ(Q, R, A, B; bet=beta) P, F, d = stationary_values(lq) # print the results diff --git a/lectures/multi_agent_models/schelling.md b/lectures/multi_agent_models/schelling.md index 868aed51..4bb0f874 100644 --- a/lectures/multi_agent_models/schelling.md +++ b/lectures/multi_agent_models/schelling.md @@ -8,7 +8,7 @@ kernelspec: language: julia name: julia-1.9 --- - +##AT THIS ONE BUT LUCAS (schelling)= ```{raw} html
@@ -160,7 +160,7 @@ Random.seed!(42); ``` ```{code-cell} julia -Agent = @with_kw (kind, location = rand(2)) +Agent(;kind, location = rand(2)) = (;kind, location) draw_location!(a) = a.location .= rand(2) @@ -179,7 +179,7 @@ function is_happy(a) # partialsortperm(get_distance.(Ref(a), agents), # 1:neighborhood_size)) - return share ≥ preference + return share >= preference end function update!(a) diff --git a/lectures/multi_agent_models/uncertainty_traps.md b/lectures/multi_agent_models/uncertainty_traps.md index a637b118..477ef665 100644 --- a/lectures/multi_agent_models/uncertainty_traps.md +++ b/lectures/multi_agent_models/uncertainty_traps.md @@ -236,17 +236,18 @@ using DataFrames, LaTeXStrings, Parameters, Plots ``` ```{code-cell} julia -UncertaintyTrapEcon = @with_kw (a = 1.5, # risk aversion - γ_x = 0.5, # production shock precision - ρ = 0.99, # correlation coefficient for θ - σ_θ = 0.5, # standard dev. of θ shock - num_firms = 100, # number of firms - σ_F = 1.5, # standard dev. of fixed costs - c = -420.0, # external opportunity cost - μ_init = 0.0, # initial value for μ - γ_init = 4.0, # initial value for γ - θ_init = 0.0, # initial value for θ - σ_x = sqrt(a / γ_x)) # standard dev. of shock +UncertaintyTrapEcon(;a = 1.5, # risk aversion + gamma_x = 0.5, # production shock precision + rho = 0.99, # correlation coefficient for theta + sigma_theta = 0.5, # standard dev. of theta shock + num_firms = 100, # number of firms + sigma_F = 1.5, # standard dev. of fixed costs + c = -420.0, # external opportunity cost + mu_init = 0.0, # initial value for mu + gamma_init = 4.0, # initial value for gamma + theta_init = 0.0, # initial value for theta + sigma_x = sqrt(a / gamma_x)) = # standard dev. of shock + (;a, gamma_x, rho, sigma_theta, num_firms, sigma_F, c, mu_init, gamma_init, theta_init, sigma_x) ``` In the results below we use this code to simulate time series for the major variables. @@ -360,17 +361,17 @@ different values of $M$ ```{code-cell} julia econ = UncertaintyTrapEcon() -@unpack ρ, σ_θ, γ_x = econ # simplify names +(;rho, sigma_theta, gamma_x) = econ # simplify names -# grid for γ and γ_{t+1} -γ = range(1e-10, 3, length = 200) +# grid for gamma and gamma_{t+1} +gamma = range(1e-10, 3, length = 200) M_range = 0:6 -γp = 1 ./ (ρ^2 ./ (γ .+ γ_x .* M_range') .+ σ_θ^2) +gammap = 1 ./ (rho^2 ./ (gamma .+ gamma_x .* M_range') .+ sigma_theta^2) labels = ["0" "1" "2" "3" "4" "5" "6"] -plot(γ, γ, lw = 2, label = "45 Degree") -plot!(γ, γp, lw = 2, label = labels) +plot(gamma, gamma, lw = 2, label = "45 Degree") +plot!(gamma, gammap, lw = 2, label = labels) plot!(xlabel = L"\gamma", ylabel = L"\gamma^\prime", legend_title = L"M", legend = :bottomright) ``` @@ -379,9 +380,9 @@ plot!(xlabel = L"\gamma", ylabel = L"\gamma^\prime", legend_title = L"M", legend tags: [remove-cell] --- @testset begin - @test γp[2,2] ≈ 0.46450522950184053 - @test γp[3,3] ≈ 0.8323524432613787 - @test γp[5,5] ≈ 1.3779664509290432 + @test gammap[2,2] ≈ 0.46450522950184053 + @test gammap[3,3] ≈ 0.8323524432613787 + @test gammap[5,5] ≈ 1.3779664509290432 end ``` @@ -396,49 +397,49 @@ is, the number of active firms and average output ```{code-cell} julia function simulate(uc, capT = 2_000) # unpack parameters - @unpack a, γ_x, ρ, σ_θ, num_firms, σ_F, c, μ_init, γ_init, θ_init, σ_x = uc + (;a, gamma_x, rho, sigma_theta, num_firms, sigma_F, c, mu_init, gamma_init, theta_init, sigma_x) = uc # draw standard normal shocks w_shocks = randn(capT) # aggregate functions - # auxiliary function ψ - function ψ(γ, μ, F) - temp1 = -a * (μ - F) - temp2 = 0.5 * a^2 / (γ + γ_x) + # auxiliary function psi + function psi(gamma, mu, F) + temp1 = -a * (mu - F) + temp2 = 0.5 * a^2 / (gamma + gamma_x) return (1 - exp(temp1 + temp2)) / a - c end # compute X, M - function gen_aggregates(γ, μ, θ) - F_vals = σ_F * randn(num_firms) - M = sum(ψ.(Ref(γ), Ref(μ), F_vals) .> 0) # counts number of active firms - if any(ψ(γ, μ, f) > 0 for f in F_vals) # ∃ an active firm - x_vals = θ .+ σ_x * randn(M) + function gen_aggregates(gamma, mu, theta) + F_vals = sigma_F * randn(num_firms) + M = sum(psi.(Ref(gamma), Ref(mu), F_vals) .> 0) # counts number of active firms + if any(psi(gamma, mu, f) > 0 for f in F_vals) # there is an active firm + x_vals = theta .+ sigma_x * randn(M) X = mean(x_vals) else X = 0.0 end - return (X = X, M = M) + return (;X, M) end # initialize dataframe - X_init, M_init = gen_aggregates(γ_init, μ_init, θ_init) - df = DataFrame(γ = γ_init, μ = μ_init, θ = θ_init, X = X_init, M = M_init) + X_init, M_init = gen_aggregates(gamma_init, mu_init, theta_init) + df = DataFrame(gamma = gamma_init, mu = mu_init, theta = theta_init, X = X_init, M = M_init) # update dataframe for t in 2:capT # unpack old variables - θ_old, γ_old, μ_old, X_old, M_old = (df.θ[end], df.γ[end], df.μ[end], df.X[end], df.M[end]) + theta_old, gamma_old, mu_old, X_old, M_old = (df.theta[end], df.gamma[end], df.mu[end], df.X[end], df.M[end]) # define new beliefs - θ = ρ * θ_old + σ_θ * w_shocks[t-1] - μ = (ρ * (γ_old * μ_old + M_old * γ_x * X_old))/(γ_old + M_old * γ_x) - γ = 1 / (ρ^2 / (γ_old + M_old * γ_x) + σ_θ^2) + theta = rho * theta_old + sigma_theta * w_shocks[t-1] + mu = (rho * (gamma_old * mu_old + M_old * gamma_x * X_old))/(gamma_old + M_old * gamma_x) + gamma = 1 / (rho^2 / (gamma_old + M_old * gamma_x) + sigma_theta^2) # compute new aggregates - X, M = gen_aggregates(γ, μ, θ) - push!(df, (γ = γ, μ = μ, θ = θ, X = X, M = M)) + X, M = gen_aggregates(gamma, mu, theta) + push!(df, (;gamma, mu, theta, X, M)) end # return @@ -459,16 +460,16 @@ Random.seed!(42); # set random seed for reproducible results ```{code-cell} julia df = simulate(econ) -plot(eachindex(df.μ), df.μ, lw = 2, label = L"\mu") -plot!(eachindex(df.θ), df.θ, lw = 2, label = L"\theta") +plot(eachindex(df.mu), df.mu, lw = 2, label = L"\mu") +plot!(eachindex(df.theta), df.theta, lw = 2, label = L"\theta") plot!(xlabel = L"x", ylabel = L"y", legend_title = "Variable", legend = :bottomright) ``` Now let's plot the whole thing together ```{code-cell} julia -len = eachindex(df.θ) -yvals = [df.θ, df.μ, df.γ, df.M] +len = eachindex(df.theta) +yvals = [df.theta, df.mu, df.gamma, df.M] vars = [L"\theta", L"\mu", L"\gamma", L"M"] plt = plot(layout = (4,1), size = (600, 600)) @@ -484,7 +485,7 @@ plot(plt) --- tags: [remove-cell] --- -mdf = DataFrame(t = eachindex(df.θ), θ = df.θ, μ = df.μ, γ = df.γ, M = df.M) +mdf = DataFrame(t = eachindex(df.theta), theta = df.theta, mu = df.mu, gamma = df.gamma, M = df.M) @testset begin @test stack(mdf, collect(2:5))[:value][3] ≈ -0.49742498224730913 atol = 1e-3