diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 309ba66c..6a50b834 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,4 +1,4 @@ style = "sciml" -margin = 84 +margin = 80 yas_style_nesting = true annotate_untyped_fields_with_any = false \ No newline at end of file diff --git a/README.md b/README.md index bf018b82..13c53b4b 100644 --- a/README.md +++ b/README.md @@ -58,10 +58,14 @@ julia format_myst.jl lectures/getting_started_julia/getting_started.md ``` As a helper, you can call a shell script to do it for an entire folder + ```bash bash format_all_directory.sh lectures/dynamic_programming ``` + or to also do the unicode substitutions + +```bash bash format_all_directory.sh lectures/dynamic_programming true ``` diff --git a/lectures/Manifest.toml b/lectures/Manifest.toml index 5202eeb7..40cf09b2 100644 --- a/lectures/Manifest.toml +++ b/lectures/Manifest.toml @@ -417,9 +417,9 @@ version = "1.9.1" [[deps.DiffEqBase]] deps = ["ArrayInterface", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "SparseArrays", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"] -git-tree-sha1 = "53159da9fef2dd92815499e18f93b699ed54f397" +git-tree-sha1 = "044648af911974c3928058c1f8c83f159dece274" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.145.5" +version = "6.145.6" [deps.DiffEqBase.extensions] DiffEqBaseChainRulesCoreExt = "ChainRulesCore" @@ -624,9 +624,9 @@ version = "2.0.0" [[deps.FileIO]] deps = ["Pkg", "Requires", "UUIDs"] -git-tree-sha1 = "299dc33549f68299137e51e6d49a13b5b1da9673" +git-tree-sha1 = "c5c28c245101bd59154f649e19b038d15901b5dc" uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" -version = "1.16.1" +version = "1.16.2" [[deps.FilePathsBase]] deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"] @@ -1265,9 +1265,9 @@ version = "1.0.3" [[deps.LoopVectorization]] deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "PrecompileTools", "SIMDTypes", "SLEEFPirates", "Static", "StaticArrayInterface", "ThreadingUtilities", "UnPack", "VectorizationBase"] -git-tree-sha1 = "0f5648fbae0d015e3abe5867bca2b362f67a5894" +git-tree-sha1 = "fc712dd664097440f19a91a704299cca02134ca0" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.12.166" +version = "0.12.168" weakdeps = ["ChainRulesCore", "ForwardDiff", "SpecialFunctions"] [deps.LoopVectorization.extensions] @@ -1313,12 +1313,6 @@ git-tree-sha1 = "d2a140e551c9ec9028483e3c7d1244f417567ac0" uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" version = "1.24.0" -[[deps.MathProgBase]] -deps = ["LinearAlgebra", "SparseArrays"] -git-tree-sha1 = "9abbe463a1e9fc507f12a69e7f29346c2cdc472c" -uuid = "fdba3010-5040-5b88-9595-932c9decdf73" -version = "0.7.8" - [[deps.MatrixFactorizations]] deps = ["ArrayLayouts", "LinearAlgebra", "Printf", "Random"] git-tree-sha1 = "78f6e33434939b0ac9ba1df81e6d005ee85a7396" @@ -1396,10 +1390,14 @@ uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" version = "7.8.3" [[deps.NLopt]] -deps = ["MathOptInterface", "MathProgBase", "NLopt_jll"] -git-tree-sha1 = "5a7e32c569200a8a03c3d55d286254b0321cd262" +deps = ["NLopt_jll"] +git-tree-sha1 = "19d2a1c8a3c5b5a459f54a10e54de630c4a05701" uuid = "76087f3c-5699-56af-9a33-bf431cd00edd" -version = "0.6.5" +version = "1.0.0" +weakdeps = ["MathOptInterface"] + + [deps.NLopt.extensions] + NLoptMathOptInterfaceExt = ["MathOptInterface"] [[deps.NLopt_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] @@ -1663,9 +1661,9 @@ version = "0.2.4" [[deps.PreallocationTools]] deps = ["Adapt", "ArrayInterface", "ForwardDiff"] -git-tree-sha1 = "13078a8afc9f88ede295c17c3c673eb52c05b0b4" +git-tree-sha1 = "64bb68f76f789f5fe5930a80af310f19cdafeaed" uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46" -version = "0.4.16" +version = "0.4.17" [deps.PreallocationTools.extensions] PreallocationToolsReverseDiffExt = "ReverseDiff" @@ -1731,9 +1729,9 @@ version = "2.9.1" [[deps.QuantEcon]] deps = ["DSP", "DataStructures", "Distributions", "FFTW", "Graphs", "LinearAlgebra", "Markdown", "NLopt", "Optim", "Pkg", "Primes", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "Test"] -git-tree-sha1 = "4f989614763519b089f4127063c1795d7a412753" +git-tree-sha1 = "034293b29fdbcae73aeb7ca0b2755e693f04701b" uuid = "fcd29c91-0bd7-5a09-975d-7ac3f643a60c" -version = "0.16.5" +version = "0.16.6" [[deps.Query]] deps = ["DataValues", "IterableTables", "MacroTools", "QueryOperators", "Statistics"] @@ -2095,9 +2093,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"] [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "42d5373c10272d14ef49cc68ffc22df3b93c549a" +git-tree-sha1 = "4e17a790909b17f7bf1496e3aec138cf01b60b3b" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.8.2" +version = "1.9.0" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] diff --git a/lectures/continuous_time/covid_sde.md b/lectures/continuous_time/covid_sde.md index fae1e19b..9ba314db 100644 --- a/lectures/continuous_time/covid_sde.md +++ b/lectures/continuous_time/covid_sde.md @@ -605,4 +605,3 @@ In this case, there are significant differences between the early and late death This bleak simulation has assumed that no individuals has long-term immunity and that there will be no medical advancements on that time horizon - both of which are unlikely to be true. Nevertheless, it suggests that the timing of lifting lockdown has a more profound impact after 18 months if we allow stochastic shocks imperfect immunity. - diff --git a/lectures/continuous_time/seir_model.md b/lectures/continuous_time/seir_model.md index 5cc01c4c..c0a36aa9 100644 --- a/lectures/continuous_time/seir_model.md +++ b/lectures/continuous_time/seir_model.md @@ -517,4 +517,3 @@ Despite its richness, the model above is fully deterministic. The policy $\bar{ One way that randomness can lead to aggregate fluctuations is the granularity that comes through the discreteness of individuals. This topic, the connection between SDEs and the Langevin equations typically used in the approximation of chemical reactions in well-mixed media is explored in further lectures on continuous time Markov chains. Instead, in the {doc}`next lecture `, we will concentrate on randomness that comes from aggregate changes in behavior or policy. - diff --git a/lectures/multi_agent_models/aiyagari.md b/lectures/multi_agent_models/aiyagari.md index 7cfd8054..0253cfaa 100644 --- a/lectures/multi_agent_models/aiyagari.md +++ b/lectures/multi_agent_models/aiyagari.md @@ -195,6 +195,7 @@ The object also includes a default set of parameters that we'll adopt unless oth ```{code-cell} julia using LinearAlgebra, Statistics +using LaTeXStrings, Plots, QuantEcon ``` ```{code-cell} julia @@ -205,30 +206,27 @@ using Test, Random ``` ```{code-cell} julia -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)) - -function setup_Q!(Q, s_i_vals, z_chain) +function 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), + # -Inf is the utility of dying (0 consumption) + u = sigma == 1 ? x -> log(x) : + x -> (x^(1 - sigma) - 1) / (1 - sigma)) + + # Create grids + 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) + + # Fill in the Q and R + Q = zeros(n, a_size, n) for next_s_i in 1:size(Q, 3) for a_i in 1:size(Q, 2) for s_i in 1:size(Q, 1) @@ -241,10 +239,8 @@ function setup_Q!(Q, s_i_vals, z_chain) end end end - return Q -end -function setup_R!(R, a_vals, s_vals, r, w, u) + R = fill(-Inf, n, a_size) for new_a_i in 1:size(R, 2) a_new = a_vals[new_a_i] for s_i in 1:size(R, 1) @@ -256,7 +252,8 @@ function setup_R!(R, a_vals, s_vals, r, w, u) end end end - return R + return (; 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) end ``` @@ -271,16 +268,16 @@ Random.seed!(42); ```{code-cell} julia # Create an instance of Household -am = Household(a_max = 20.0, r = 0.03, w = 0.956) +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) # Simplify names -(;z_size, a_size, n, a_vals) = am +(; z_size, a_size, n, a_vals) = am z_vals = am.z_chain.state_values # Get all optimal actions across the set of @@ -298,10 +295,10 @@ plot!(xlabel = "current assets", ylabel = "next period assets", grid = false) tags: [remove-cell] --- @testset begin - #test a_vals[4] ≈ 0.3015075377869347 - #test a_star[4] ≈ 0.2010050252246231 - #test results.v[4] ≈ -27.48291672016239 - #test z_vals ≈ [0.1, 1.0] + @test a_vals[4] ≈ 0.3015075377869347 + @test a_star[4] ≈ 0.2010050252246231 + @test results.v[4] ≈ -27.48291672016239 + @test z_vals ≈ [0.1, 1.0] end ``` @@ -323,29 +320,16 @@ Random.seed!(42); ``` ```{code-cell} julia -# Firms' parameters -const A = 1 -const N = 1 -const α = 0.33 -const β = 0.96 -const δ = 0.05 - -function r_to_w(r) - return A * (1 - α) * (A * α / (r + δ)) ^ (α / (1 - α)) -end -function rd(K) - return A * α * (N / K) ^ (1 - α) - δ -end - -function prices_to_capital_stock(am, r) +# Calculate supply of capital for a given r +function prices_to_capital_stock(r; beta, A, N, alpha, delta, a_max) + # Create an instance of Household given the parameters - # Set up problem - w = r_to_w(r) - (;a_vals, s_vals, u) = am - setup_R!(am.R, a_vals, s_vals, r, w, u) + # Calculate the equilibrium wages + w = A * (1 - alpha) * (A * alpha / (r + delta))^(alpha / (1 - alpha)) + am = Household(; beta, a_max, w, r) - 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) @@ -354,33 +338,39 @@ function prices_to_capital_stock(am, r) stationary_probs = stationary_distributions(results.mc)[:, 1][1] # Return K - return dot(am.s_vals[:, 1], stationary_probs) + K = dot(am.s_vals[:, 1], stationary_probs) + + # Return capital + return K end -# Create an instance of Household -am = Household(β = β, a_max = 20.0) +# Inverse Demand for capital +function r_inverse_demand(K; A, N, alpha, delta) + return A * alpha * (N / K)^(1 - alpha) - delta +end # Create a grid of r values at which to compute demand and supply of capital r_vals = range(0.005, 0.04, length = 20) -# Compute supply of capital -k_vals = prices_to_capital_stock.(Ref(am), r_vals) +# Firms' parameters +A = 1 +N = 1 +alpha = 0.33 +beta = 0.96 +delta = 0.05 +a_max = 20.0 -# Plot against demand for capital by firms -demand = rd.(k_vals) -labels = ["demand for capital" "supply of capital"] -plot(k_vals, [demand r_vals], label = labels, lw = 2, alpha = 0.6) -plot!(xlabel = "capital", ylabel = "interest rate", xlim = (2, 14), ylim = (0.0, 0.1)) -``` +prices_to_capital_stock(r_vals[1]; A, N, alpha, beta, delta, a_max) -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test k_vals[4] ≈ 3.920775511050653 - #test demand[4] ≈ 0.08211578674946372 - #test r_vals[4] ≈ 0.010526315789473684 -end -``` +# Compute supply of capital +k_vals = prices_to_capital_stock.(r_vals; A, N, alpha, beta, delta, a_max) + +r_inverse_demand_vals = r_inverse_demand.(k_vals; A, N, alpha, delta) +# Plot against demand for capital by firms +labels = ["demand for capital" "supply of capital"] +plot(k_vals, [r_inverse_demand_vals r_vals], label = labels, lw = 2, + alpha = 0.6) +plot!(xlabel = "capital", ylabel = "interest rate", xlim = (2, 14), + ylim = (0.0, 0.1)) +``` \ No newline at end of file diff --git a/lectures/multi_agent_models/arellano.md b/lectures/multi_agent_models/arellano.md index 28dd86bc..41805bcc 100644 --- a/lectures/multi_agent_models/arellano.md +++ b/lectures/multi_agent_models/arellano.md @@ -301,7 +301,7 @@ The code can be found below: tags: [hide-output] --- using LinearAlgebra, Statistics -using LaTeXStrings, Parameters, QuantEcon, DataFrames, Plots, Random +using LaTeXStrings, QuantEcon, DataFrames, Plots, Random ``` ```{code-cell} julia @@ -312,25 +312,25 @@ using Test ``` ```{code-cell} julia -function ArellanoEconomy(;β = .953, - γ = 2., - r = 0.017, - ρ = 0.945, - η = 0.025, - θ = 0.282, - ny = 21, - nB = 251) +function ArellanoEconomy(; beta = 0.953, + gamma = 2.0, + r = 0.017, + 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 + Bgrid = collect(range(-0.4, 0.4, length = nB)) + mc = tauchen(ny, rho, eta) + Pi = mc.p ygrid = exp.(mc.state_values) - ydefgrid = min.(.969 * mean(ygrid), ygrid) + ydefgrid = min.(0.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,16 +351,17 @@ function one_step_update!(ae, EVc) # unpack stuff - @unpack β, γ, r, ρ, η, θ, ny, nB = ae - @unpack ygrid, ydefgrid, Bgrid, Π, vf, vd, vc, policy, q, defprob = ae - zero_ind = searchsortedfirst(Bgrid, 0.) + (; 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.0) for iy in 1:ny y = ae.ygrid[iy] 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 @@ -369,15 +369,14 @@ function one_step_update!(ae, current_max = -1e14 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] + for ib_next in 1:nB + c = max(y - ae.q[ib_next, iy] * Bgrid[ib_next] + B, 1e-14) + m = u(ae, c) + beta * EV[ib_next, iy] if m > current_max current_max = m pol_ind = ib_next end - end # update value and policy functions @@ -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,13 +404,13 @@ 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 - dist = 10. + dist = 10.0 # allocate memory for update V_upd = similar(ae.vf) @@ -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) @@ -443,20 +442,19 @@ end function QuantEcon.simulate(ae, capT = 5000; y_init = mean(ae.ygrid), - B_init = mean(ae.Bgrid), - ) + B_init = mean(ae.Bgrid),) # get initial indices - zero_index = searchsortedfirst(ae.Bgrid, 0.) + zero_index = searchsortedfirst(ae.Bgrid, 0.0) y_init_ind = searchsortedfirst(ae.ygrid, y_init) 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 - y_sim_val = zeros(capT+1) + y_sim_val = zeros(capT + 1) B_sim_val, q_sim_val = similar(y_sim_val), similar(y_sim_val) B_sim_indices = fill(0, capT + 1) default_status = fill(false, capT + 1) @@ -478,7 +476,7 @@ function QuantEcon.simulate(ae, default_status[t + 1] = true y_sim_val[t] = ae.ydefgrid[y_sim_indices[t]] B_sim_indices[t + 1] = zero_index - B_sim_val[t+1] = 0. + B_sim_val[t + 1] = 0.0 q_sim_val[t] = ae.q[zero_index, y_sim_indices[t]] else default_status[t] = false @@ -488,15 +486,15 @@ function QuantEcon.simulate(ae, q_sim_val[t] = ae.q[B_sim_indices[t + 1], y_sim_indices[t]] end - # if you are in default + # if you are in default else B_sim_indices[t + 1] = zero_index - B_sim_val[t+1] = 0. + B_sim_val[t + 1] = 0.0 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,14 +589,14 @@ using DataFrames, Plots Compute the value function, policy and equilibrium prices ```{code-cell} julia -ae = ArellanoEconomy(β = .953, # time discount rate - γ = 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 - ny = 21, # number of points in y grid - nB = 251) # number of points in B grid +ae = ArellanoEconomy(beta = 0.953, # time discount rate + gamma = 2.0, # risk aversion + r = 0.017, # international interest rate + rho = 0.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 # now solve the model on the grid. vfi!(ae) @@ -625,9 +623,9 @@ iy_high, iy_low = map(x -> searchsortedfirst(ae.ygrid, x), (high, low)) x = zeros(0) q_low = zeros(0) q_high = zeros(0) -for i in 1:ae.nB +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]) @@ -638,7 +636,8 @@ end plot(x, q_low, label = "Low") plot!(x, q_high, label = "High") plot!(title = L"Bond price schedule $q(y, B^\prime)$", - xlabel = L"B^\prime", ylabel = L"q", legend_title = L"y", legend = :topleft) + xlabel = L"B^\prime", ylabel = L"q", legend_title = L"y", + legend = :topleft) ``` ```{code-cell} julia @@ -658,17 +657,18 @@ Draw a plot of the value functions plot(ae.Bgrid, ae.vf[:, iy_low], label = "Low") plot!(ae.Bgrid, ae.vf[:, iy_high], label = "High") plot!(xlabel = L"B", ylabel = L"V(y,B)", title = "Value functions", - legend_title=L"y", legend = :topleft) + legend_title = L"y", legend = :topleft) ``` Draw a heat map for default probability ```{code-cell} julia -heatmap(ae.Bgrid[1:end-1], - ae.ygrid[2:end], - reshape(clamp.(vec(ae.defprob[1:end - 1, 1:end - 1]), 0, 1), 250, 20)') +heatmap(ae.Bgrid[1:(end - 1)], + ae.ygrid[2:end], + reshape(clamp.(vec(ae.defprob[1:(end - 1), 1:(end - 1)]), 0, 1), 250, + 20)') plot!(xlabel = L"B^\prime", ylabel = L"y", title = "Probability of default", - legend = :topleft) + legend = :topleft) ``` Plot a time series of major variables simulated from the model @@ -696,10 +696,12 @@ plots = plot(layout = (3, 1), size = (700, 800)) # Plot the three variables, and for each each variable shading the period(s) of default # in grey for i in 1:3 - plot!(plots[i], 1:T, y_vals[i], title = titles[i], xlabel = "time", label = "", lw = 2) + plot!(plots[i], 1:T, y_vals[i], title = titles[i], xlabel = "time", + label = "", lw = 2) for j in 1:length(def_start) plot!(plots[i], [def_start[j], def_end[j]], fill(maximum(y_vals[i]), 2), - fillrange = [extrema(y_vals[i])...], fcolor = :grey, falpha = 0.3, label = "") + fillrange = [extrema(y_vals[i])...], fcolor = :grey, falpha = 0.3, + label = "") end end diff --git a/lectures/multi_agent_models/harrison_kreps.md b/lectures/multi_agent_models/harrison_kreps.md index 70e35213..2517dfdf 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 = 0.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,28 +416,28 @@ 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, - tol=1e-16) + beta = 0.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] + 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 = [maximum((q * p_old) + (q * dividend_payoff)) - for q in transitions] - p_new = β * temp + for q in transitions] + p_new = beta * temp # If we succed in converging, break out of for loop - if maximum(sqrt, ((p_new - p_old).^2)) < 1e-12 + if maximum(sqrt, ((p_new - p_old) .^ 2)) < 1e-12 break end end - temp=[minimum((q * p_old) + (q * dividend_payoff)) for q in transitions] - ptwiddle = β * temp + temp = [minimum((q * p_old) + (q * dividend_payoff)) for q in transitions] + ptwiddle = beta * temp phat_a = [p_new[1], ptwiddle[2]] phat_b = [ptwiddle[1], p_new[2]] @@ -484,22 +484,21 @@ Constraints on short sales prevent that. Here's code to solve for $\check p$ using iteration ```{code-cell} julia -function price_pessimistic_beliefs(transitions, - dividend_payoff; - β=.75, max_iter=50000, - tol=1e-16) +function price_pessimistic_beliefs(transitions, dividend_payoff; beta = 0.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 + temp = [minimum((q * p_old) + (q * dividend_payoff)) + for q in transitions] + p_new = beta * temp # If we succed in converging, break out of for loop - if maximum(sqrt, ((p_new - p_old).^2)) < 1e-12 + if maximum(sqrt, ((p_new - p_old) .^ 2)) < 1e-12 break end end @@ -567,7 +566,8 @@ labels = ["p_a", "p_b", "p_optimistic", "p_pessimistic"] for (transition, label) in zip(transitions, labels) println(label) println(repeat("=", 20)) - s0, s1 = round.(price_single_beliefs(transition, dividendreturn), digits=2) + s0, s1 = round.(price_single_beliefs(transition, dividendreturn), + digits = 2) println("State 0: $s0") println("State 1: $s1") println(repeat("-", 20)) @@ -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 b0e43fdd..c8aee14f 100644 --- a/lectures/multi_agent_models/lake_model.md +++ b/lectures/multi_agent_models/lake_model.md @@ -196,9 +196,8 @@ Here's the code: --- tags: [hide-output] --- -using LinearAlgebra, Statistics -using Distributions, Expectations, NLsolve, Parameters, Plots -using QuantEcon, Roots, Random +using LinearAlgebra, Statistics, Distributions +using NLsolve, Plots, QuantEcon, Roots, Random ``` ```{code-cell} julia @@ -208,49 +207,36 @@ tags: [remove-cell] using Test ``` -```{code-cell} julia - -``` +Reusable functions for simulating linear maps and finding their fixed points ```{code-cell} julia -LakeModel = @with_kw (λ = 0.283, α = 0.013, b = 0.0124, d = 0.00822) - -function transition_matrices(lm) - (;λ, α, 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, Â = Â) +function simulate_linear(A, x_0, T) + X = zeros(length(x_0), T + 1) + X[:, 1] = x_0 + for t in 2:(T + 1) + X[:, t] = A * X[:, t - 1] + end + return X end -function rate_steady_state(lm) - (;Â) = transition_matrices(lm) - sol = fixedpoint(x -> Â * x, fill(0.5, 2)) - converged(sol) || error("Failed to converge in $(result.iterations) iterations") +function linear_steady_state(A, x_0 = ones(size(A, 1)) / size(A, 1)) + sol = fixedpoint(x -> A * x, x_0) + converged(sol) || error("Failed to converge in $(sol.iterations) iter") return sol.zero end +``` -function simulate_stock_path(lm, X0, T) - (;A) = transition_matrices(lm) - X_path = zeros(eltype(X0), 2, T) - X = copy(X0) - for t in 1:T - X_path[:, t] = X - X = A * X - end - return X_path -end +```{code-cell} julia +function LakeModel(; lambda = 0.283, alpha = 0.013, b = 0.0124, d = 0.00822) + # calculate transition matrices + g = b - d + A = [(1 - lambda) * (1 - d)+b (1 - d) * alpha+b + (1 - d)*lambda (1 - d)*(1 - alpha)] + A_hat = A ./ (1 + g) -function simulate_rate_path(lm, x0, T) - (;Â) = transition_matrices(lm) - x_path = zeros(eltype(x0), 2, T) - x = copy(x0) - for t in 1:T - x_path[:, t] = x - x = Â * x - end - return x_path + # Solve for fixed point to find the steady-state u and e + x_bar = linear_steady_state(A_hat) + return (; lambda, alpha, b, d, A, A_hat, x_bar) end ``` @@ -258,24 +244,22 @@ Let's observe these matrices for the baseline model ```{code-cell} julia lm = LakeModel() -A, Â = transition_matrices(lm) -A +lm.A ``` ```{code-cell} julia -Â +lm.A_hat ``` And a revised model ```{code-cell} julia -lm = LakeModel(α = 2.0) -A, Â = transition_matrices(lm) -A +lm = LakeModel(; alpha = 0.2) +lm.A ``` ```{code-cell} julia -Â +lm.A_hat ``` ```{code-cell} julia @@ -283,8 +267,8 @@ Â tags: [remove-cell] --- @testset begin - @test lm.α ≈ 2.0 - @test A[1][1] ≈ 0.7235062600000001 + @test lm.alpha ≈ 0.2 + @test lm.A[1][1] ≈ 0.7235062600000001 end ``` @@ -303,15 +287,18 @@ U_0 = u_0 * N_0 E_0 = e_0 * N_0 X_0 = [U_0; E_0] -X_path = simulate_stock_path(lm, X_0, T) +X_path = simulate_linear(lm.A, X_0, T - 1) x1 = X_path[1, :] x2 = X_path[2, :] -x3 = dropdims(sum(X_path, dims = 1), dims = 1) +x3 = sum(X_path, dims = 1)' -plt_unemp = plot(title = "Unemployment", 1:T, x1, color = :blue, lw = 2, grid = true, label = "") -plt_emp = plot(title = "Employment", 1:T, x2, color = :blue, lw = 2, grid = true, label = "") -plt_labor = plot(title = "Labor force", 1:T, x3, color = :blue, lw = 2, grid = true, label = "") +plt_unemp = plot(title = "Unemployment", 1:T, x1, color = :blue, lw = 2, + grid = true, label = "") +plt_emp = plot(title = "Employment", 1:T, x2, color = :blue, lw = 2, + grid = true, label = "") +plt_labor = plot(title = "Labor force", 1:T, x3, color = :blue, lw = 2, + grid = true, label = "") plot(plt_unemp, plt_emp, plt_labor, layout = (3, 1), size = (800, 600)) ``` @@ -343,8 +330,7 @@ This is the case for our default parameters: ```{code-cell} julia lm = LakeModel() -A, Â = transition_matrices(lm) -e, f = eigvals(Â) +e, f = eigvals(lm.A_hat) abs(e), abs(f) ``` @@ -362,21 +348,23 @@ Let's look at the convergence of the unemployment and employment rate to steady ```{code-cell} julia lm = LakeModel() -e_0 = 0.92 # initial employment rate -u_0 = 1 - e_0 # initial unemployment rate -T = 50 # simulation length +e_0 = 0.92 # initial employment rate +u_0 = 1 - e_0 # initial unemployment rate +T = 50 # simulation length -xbar = rate_steady_state(lm) +u_bar, e_bar = lm.x_bar x_0 = [u_0; e_0] -x_path = simulate_rate_path(lm, x_0, T) - -plt_unemp = plot(title ="Unemployment rate", 1:T, x_path[1, :],color = :blue, lw = 2, - alpha = 0.5, grid = true, label = "") -plot!(plt_unemp, [xbar[1]], color=:red, linetype = :hline, linestyle = :dash, lw = 2, label = "") -plt_emp = plot(title = "Employment rate", 1:T, x_path[2, :],color = :blue, lw = 2, alpha = 0.5, - grid = true, label = "") -plot!(plt_emp, [xbar[2]], color=:red, linetype = :hline, linestyle = :dash, lw = 2, label = "") -plot(plt_unemp, plt_emp, layout = (2, 1), size=(700,500)) +x_path = simulate_linear(lm.A_hat, x_0, T - 1) + +plt_unemp = plot(title = "Unemployment rate", 1:T, x_path[1, :], color = :blue, + lw = 2, alpha = 0.5, grid = true, label = "") +plot!(plt_unemp, [u_bar], color = :red, linetype = :hline, linestyle = :dash, + lw = 2, label = "") +plt_emp = plot(title = "Employment rate", 1:T, x_path[2, :], color = :blue, + lw = 2, alpha = 0.5, grid = true, label = "") +plot!(plt_emp, [e_bar], color = :red, linetype = :hline, linestyle = :dash, + lw = 2, label = "") +plot(plt_unemp, plt_emp, layout = (2, 1), size = (700, 500)) ``` ```{code-cell} julia @@ -480,35 +468,33 @@ MarkovChain type to investigate this. Let's plot the path of the sample averages over 5,000 periods ```{code-cell} julia -using QuantEcon, Roots, Random -``` - -```{code-cell} julia -lm = LakeModel(d = 0, b = 0) +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 Random.seed!(42) +u_bar, e_bar = lm.x_bar 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] -plt_unemp = plot(title = "Percent of time unemployed", 1:T, s_bars[:,1],color = :blue, lw = 2, - alpha = 0.5, label = "", grid = true) -plot!(plt_unemp, [xbar[1]], linetype = :hline, linestyle = :dash, color=:red, lw = 2, label = "") -plt_emp = plot(title = "Percent of time employed", 1:T, s_bars[:,2],color = :blue, lw = 2, - alpha = 0.5, label = "", grid = true) -plot!(plt_emp, [xbar[2]], linetype = :hline, linestyle = :dash, color=:red, lw = 2, label = "") -plot(plt_unemp, plt_emp, layout = (2, 1), size=(700,500)) +s_path = simulate(mc, T; init = 2) +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) +plot!(plt_unemp, [u_bar], linetype = :hline, linestyle = :dash, color = :red, + lw = 2, label = "") +plt_emp = plot(title = "Percent of time employed", 1:T, s_bars[:, 2], + color = :blue, lw = 2, alpha = 0.5, label = "", grid = true) +plot!(plt_emp, [e_bar], linetype = :hline, linestyle = :dash, color = :red, + lw = 2, label = "") +plot(plt_unemp, plt_emp, layout = (2, 1), size = (700, 500)) ``` ```{code-cell} julia @@ -516,8 +502,8 @@ plot(plt_unemp, plt_emp, layout = (2, 1), size=(700,500)) tags: [remove-cell] --- @testset begin - #test xbar[1] ≈ 0.04391891891891919 - #test s_bars[end,end] ≈ 0.957 + @test u_bar ≈ 0.04391891891891919 + @test s_bars[end,end] ≈ 0.9488 end ``` @@ -630,160 +616,190 @@ 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, +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 + (; alpha, beta, sigma, c, gamma, w, w_probs, u) = mcm - # necessary objects - u_w = u.(w, σ) - u_c = u(c, σ) + # pre-calculate utilities + 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] + V = x[1:(end - 1)] U = x[end] - [u_w + β * ((1 - α) * V .+ α * U); u_c + β * (1 - γ) * U + β * γ * E * max.(U, V)] + return [u_w + beta * ((1 - alpha) * V .+ alpha * U); + u_c + beta * (1 - gamma) * U + + beta * gamma * dot(w_probs, max.(U, V))] end # value function iteration x_iv = [V_iv; U_iv] # initial x val xstar = fixedpoint(T, x_iv, iterations = iter, xtol = tol, m = 0).zero - V = xstar[1:end-1] + V = xstar[1:(end - 1)] U = xstar[end] # 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 ``` And the McCall object ```{code-cell} julia -# a default utility function -u(c, σ) = c > 0 ? (c^(1 - σ) - 1) / (1 - σ) : -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 +function McCallModel(; alpha, beta, gamma, c, sigma, w, w_probs, + u = (c, sigma) -> c > 0 ? + (c^(1 - sigma) - 1) / (1 - sigma) : + -10e-6) + return (; alpha, beta, gamma, c, sigma, u, w, w_probs) +end ``` Now let's compute and plot welfare, employment, unemployment, and tax revenue as a function of the unemployment compensation rate ```{code-cell} julia -# some global variables that will stay constant -α = 0.013 -α_q = (1 - (1 - α)^3) -b_param = 0.0124 -d_param = 0.00822 -β = 0.98 -γ = 1.0 -σ = 2.0 - -# the default wage distribution: a discretized log normal -log_wage_mean, wage_grid_size, max_wage = 20, 200, 170 -w_vec = range(1e-3, max_wage, length = wage_grid_size + 1) - -logw_dist = Normal(log(log_wage_mean), 1) -cdf_logw = cdf.(logw_dist, log.(w_vec)) -pdf_logw = cdf_logw[2:end] - cdf_logw[1:end-1] - -p_vec = pdf_logw ./ sum(pdf_logw) -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 - E = E) # expectation operator - - (;V, U, w̄) = solve_mccall_model(mcm) - indicator = wage -> wage > w̄ - λ = γ * E * indicator.(w_vec .- τ) - - return w̄, λ, V, U +function compute_optimal_quantities(c_pretax, tau; w_probs, sigma, gamma, beta, + alpha, w_pretax) + mcm = McCallModel(; alpha, beta, gamma, sigma, w_probs, + c = c_pretax - tau, # post-tax compensation + w = w_pretax .- tau) + + (; V, U, w_bar) = solve_mccall_model(mcm) + accept_wage = w_pretax .- tau .> w_bar + + # sum up proportion accepting the wages + lambda = gamma * dot(w_probs, accept_wage) + 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_pretax, tau; w_probs, sigma, gamma, + beta, alpha, w_pretax, b, d) + w_bar, lambda, V, U = compute_optimal_quantities(c_pretax, tau; w_probs, + sigma, + gamma, beta, alpha, + w_pretax) # compute steady state employment and unemployment rates - lm = LakeModel(λ = λ_param, α = α_q, b = b_param, d = d_param) - x = rate_steady_state(lm) - u_rate, e_rate = x + lm = LakeModel(; lambda, alpha, b, d) + u_rate, e_rate = lm.x_bar # compute steady state welfare - indicator(wage) = wage > w̄ - decisions = indicator.(w_vec .- τ) - w = (E * (V .* decisions)) / (E * decisions) + accept_wage = w_pretax .- tau .> w_bar + w = (dot(w_probs, V .* accept_wage)) / dot(w_probs, accept_wage) welfare = e_rate .* w + u_rate .* U return u_rate, e_rate, welfare end -function find_balanced_budget_tax(c) +function find_balanced_budget_tax(c_pretax; w_probs, sigma, gamma, beta, alpha, + w_pretax, b, d) function steady_state_budget(t) - u_rate, e_rate, w = compute_steady_state_quantities(c, t) - return t - u_rate * c + u_rate, e_rate, w = compute_steady_state_quantities(c_pretax, t; + w_probs, sigma, + gamma, beta, alpha, + w_pretax, b, d) + return t - u_rate * c_pretax end - τ = find_zero(steady_state_budget, (0.0, 0.9c)) - return τ + tau = find_zero(steady_state_budget, (0.0, 0.9 * c_pretax)) + return tau end +``` -# levels of unemployment insurance we wish to study -Nc = 60 -c_vec = range(5, 140, length = Nc) - -tax_vec = zeros(Nc) -unempl_vec = similar(tax_vec) -empl_vec = similar(tax_vec) -welfare_vec = similar(tax_vec) - -for i in 1:Nc - t = find_balanced_budget_tax(c_vec[i]) - u_rate, e_rate, welfare = compute_steady_state_quantities(c_vec[i], t) - tax_vec[i] = t - unempl_vec[i] = u_rate - empl_vec[i] = e_rate - welfare_vec[i] = welfare +Helper function to calculate for various unemployment insurance levels + +```{code-cell} julia +function calculate_equilibriums(c_pretax; w_probs, sigma, gamma, beta, alpha, + w_pretax, b, d) + tau_vec = similar(c_pretax) + u_vec = similar(c_pretax) + e_vec = similar(c_pretax) + welfare_vec = similar(c_pretax) + + for (i, c_pre) in enumerate(c_pretax) + tau = find_balanced_budget_tax(c_pre; w_probs, sigma, gamma, beta, + alpha, w_pretax, b, d) + u_rate, e_rate, welfare = compute_steady_state_quantities(c_pre, tau; + w_probs, + sigma, + gamma, beta, + alpha, + w_pretax, + b, d) + tau_vec[i] = tau + u_vec[i] = u_rate + e_vec[i] = e_rate + welfare_vec[i] = welfare + end + return tau_vec, u_vec, e_vec, welfare_vec end +``` -plt_unemp = plot(title = "Unemployment", c_vec, unempl_vec, color = :blue, lw = 2, alpha=0.7, - label = "",grid = true) -plt_tax = plot(title = "Tax", c_vec, tax_vec, color = :blue, lw = 2, alpha=0.7, label = "", - grid = true) -plt_emp = plot(title = "Employment", c_vec, empl_vec, color = :blue, lw = 2, alpha=0.7, label = "", - grid = true) -plt_welf = plot(title = "Welfare", c_vec, welfare_vec, color = :blue, lw = 2, alpha=0.7, label = "", - grid = true) +Parameters for our experiment + +```{code-cell} julia +alpha_base = 0.013 +alpha = (1 - (1 - alpha_base)^3) +b = 0.0124 +d = 0.00822 +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 +w_pretax = range(1e-3, max_wage, length = wage_grid_size + 1) +logw_dist = Normal(log(log_wage_mean), 1) +cdf_logw = cdf.(logw_dist, log.(w_pretax)) +pdf_logw = cdf_logw[2:end] - cdf_logw[1:(end - 1)] +w_probs = pdf_logw ./ sum(pdf_logw) # probabilities +w_pretax = (w_pretax[1:(end - 1)] + w_pretax[2:end]) / 2 -plot(plt_unemp, plt_emp, plt_tax, plt_welf, layout = (2,2), size = (800, 700)) +# levels of unemployment insurance we wish to study +c_pretax = range(5, 140, length = 60) +tau_vec, u_vec, e_vec, welfare_vec = calculate_equilibriums(c_pretax; w_probs, + sigma, gamma, beta, + alpha, w_pretax, b, + d) + +# plots +plt_unemp = plot(title = "Unemployment", c_pretax, u_vec, color = :blue, + lw = 2, alpha = 0.7, label = "", grid = true) +plt_tax = plot(title = "Tax", c_pretax, tau_vec, color = :blue, lw = 2, + alpha = 0.7, label = "", grid = true) +plt_emp = plot(title = "Employment", c_pretax, e_vec, color = :blue, lw = 2, + alpha = 0.7, label = "", grid = true) +plt_welf = plot(title = "Welfare", c_pretax, welfare_vec, color = :blue, lw = 2, + alpha = 0.7, label = "", grid = true) + +plot(plt_unemp, plt_emp, plt_tax, plt_welf, layout = (2, 2), size = (800, 700)) ``` Welfare first increases and then decreases as unemployment benefits rise. The level that maximizes steady state welfare is approximately 62. +```{code-cell} julia +--- +tags: [remove-cell] +--- +@testset begin + @test tau_vec[2] ≈ 0.859664806251665 + @test u_vec[10] ≈ 0.22587607040020583 + @test e_vec[4] ≈ 0.8527180941230578 +end +``` + ## Exercises ### Exercise 1 @@ -830,7 +846,7 @@ steady state values to x0 ```{code-cell} julia lm = LakeModel() -x0 = rate_steady_state(lm) +x0 = lm.x_bar println("Initial Steady State: $x0") ``` @@ -844,14 +860,10 @@ T = 50 New legislation changes $\lambda$ to $0.2$ ```{code-cell} julia -lm = LakeModel(λ = 0.2) -``` - -```{code-cell} julia -xbar = rate_steady_state(lm) # new steady state -X_path = simulate_stock_path(lm, x0 * N0, T) -x_path = simulate_rate_path(lm, x0, T) -println("New Steady State: $xbar") +lm = LakeModel(; lambda = 0.2) +X_path = simulate_linear(lm.A, x0 * N0, T - 1) +x_path = simulate_linear(lm.A_hat, x0, T - 1) +println("New Steady State: $(lm.x_bar)") ``` Now plot stocks @@ -859,14 +871,14 @@ Now plot stocks ```{code-cell} julia x1 = X_path[1, :] x2 = X_path[2, :] -x3 = dropdims(sum(X_path, dims = 1), dims = 1) +x3 = sum(X_path, dims = 1)' -plt_unemp = plot(title = "Unemployment", 1:T, x1, color = :blue, grid = true, label = "", - bg_inside = :lightgrey) -plt_emp = plot(title = "Employment", 1:T, x2, color = :blue, grid = true, label = "", - bg_inside = :lightgrey) -plt_labor = plot(title = "Labor force", 1:T, x3, color = :blue, grid = true, label = "", - bg_inside = :lightgrey) +plt_unemp = plot(title = "Unemployment", 1:T, x1, color = :blue, grid = true, + label = "", bg_inside = :lightgrey) +plt_emp = plot(title = "Employment", 1:T, x2, color = :blue, grid = true, + label = "", bg_inside = :lightgrey) +plt_labor = plot(title = "Labor force", 1:T, x3, color = :blue, grid = true, + label = "", bg_inside = :lightgrey) plot(plt_unemp, plt_emp, plt_labor, layout = (3, 1), size = (800, 600)) ``` @@ -876,23 +888,23 @@ plot(plt_unemp, plt_emp, plt_labor, layout = (3, 1), size = (800, 600)) tags: [remove-cell] --- @testset begin - #test x1[1] ≈ 8.266626766923284 - #test x2[2] ≈ 91.43632870031433 - #test x3[3] ≈ 100.83774723999996 + @test x3[21] ≈ 108.70045145781441 + @test x2[12] ≈ 93.05432951806618 end ``` And how the rates evolve ```{code-cell} julia -plt_unemp = plot(title = "Unemployment rate", 1:T, x_path[1,:], color = :blue, grid = true, - label = "", bg_inside = :lightgrey) -plot!(plt_unemp, [xbar[1]], linetype = :hline, linestyle = :dash, color =:red, label = "") - -plt_emp = plot(title = "Employment rate", 1:T, x_path[2,:], color = :blue, grid = true, - label = "", bg_inside = :lightgrey) -plot!(plt_emp, [xbar[2]], linetype = :hline, linestyle = :dash, color =:red, label = "") +plt_unemp = plot(title = "Unemployment rate", 1:T, x_path[1, :], color = :blue, + grid = true, label = "", bg_inside = :lightgrey) +plot!(plt_unemp, [u_bar], linetype = :hline, linestyle = :dash, color = :red, + label = "") +plt_emp = plot(title = "Employment rate", 1:T, x_path[2, :], color = :blue, + grid = true, label = "", bg_inside = :lightgrey) +plot!(plt_emp, [e_bar], linetype = :hline, linestyle = :dash, color = :red, + label = "") plot(plt_unemp, plt_emp, layout = (2, 1), size = (800, 600)) ``` @@ -901,8 +913,8 @@ plot(plt_unemp, plt_emp, layout = (2, 1), size = (800, 600)) tags: [remove-cell] --- @testset begin - #test x_path[1,3] ≈ 0.09471014989625384 - #test x_path[2,7] ≈ 0.8936171021324064 + @test x_path[1,3] ≈ 0.09471014989625384 + @test x_path[2,7] ≈ 0.8936171021324064 end ``` @@ -921,7 +933,7 @@ state ```{code-cell} julia lm = LakeModel() -x0 = rate_steady_state(lm) +x0 = lm.x_bar ``` ```{code-cell} julia @@ -929,23 +941,23 @@ x0 = rate_steady_state(lm) tags: [remove-cell] --- @testset begin - #test x0[1] ≈ 0.08266626766923285 + @test x0[1] ≈ 0.08266626766923285 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_linear(lm.A, x0 * N0, T_hat - 1) +x_path1 = simulate_linear(lm.A_hat, x0, T_hat - 1) ``` Now we reset $b$ to the original value and then, using the state @@ -953,9 +965,9 @@ after 20 periods for the new initial conditions, we simulate for the 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 +lm = LakeModel(; b = 0.0124) +X_path2 = simulate_linear(lm.A, X_path1[:, end - 1], T - T_hat) +x_path2 = simulate_linear(lm.A_hat, x_path1[:, end - 1], T - T_hat) ``` Finally we combine these two paths and plot @@ -966,20 +978,20 @@ X_path = hcat(X_path1, X_path2[:, 2:end]) ``` ```{code-cell} julia -x1 = X_path[1,:] -x2 = X_path[2,:] -x3 = dropdims(sum(X_path, dims = 1), dims = 1) +x1 = X_path[1, :] +x2 = X_path[2, :] +x3 = sum(X_path, dims = 1)' -plt_unemp = plot(title = "Unemployment", 1:T, x1, color = :blue, lw = 2, alpha = 0.7, - grid = true, label = "", bg_inside = :lightgrey) +plt_unemp = plot(title = "Unemployment", 1:T, x1, color = :blue, lw = 2, + alpha = 0.7, grid = true, label = "", bg_inside = :lightgrey) plot!(plt_unemp, ylims = extrema(x1) .+ (-1, 1)) -plt_emp = plot(title = "Employment", 1:T, x2, color = :blue, lw = 2, alpha = 0.7, grid = true, - label = "", bg_inside = :lightgrey) +plt_emp = plot(title = "Employment", 1:T, x2, color = :blue, lw = 2, + alpha = 0.7, grid = true, label = "", bg_inside = :lightgrey) plot!(plt_emp, ylims = extrema(x2) .+ (-1, 1)) -plt_labor = plot(title = "Labor force", 1:T, x3, color = :blue, alpha = 0.7, grid = true, - label = "", bg_inside = :lightgrey) +plt_labor = plot(title = "Labor force", 1:T, x3, color = :blue, alpha = 0.7, + grid = true, label = "", bg_inside = :lightgrey) plot!(plt_labor, ylims = extrema(x3) .+ (-1, 1)) plot(plt_unemp, plt_emp, plt_labor, layout = (3, 1), size = (800, 600)) ``` @@ -989,23 +1001,24 @@ plot(plt_unemp, plt_emp, plt_labor, layout = (3, 1), size = (800, 600)) tags: [remove-cell] --- @testset begin - #test x1[1] ≈ 8.266626766923284 - #test x2[2] ≈ 92.11681873319097 - #test x3[3] ≈ 98.95872483999996 + @test x1[1] ≈ 8.266626766923284 + @test x2[2] ≈ 92.11681873319097 + @test x3[3] ≈ 98.95872483999996 end ``` And the rates ```{code-cell} julia -plt_unemp = plot(title = "Unemployment Rate", 1:T, x_path[1,:], color = :blue, grid = true, - label = "", bg_inside = :lightgrey, lw = 2) -plot!(plt_unemp, [x0[1]], linetype = :hline, linestyle = :dash, color =:red, label = "", lw = 2) - -plt_emp = plot(title = "Employment Rate", 1:T, x_path[2,:], color = :blue, grid = true, - label = "", bg_inside = :lightgrey, lw = 2) -plot!(plt_emp, [x0[2]], linetype = :hline, linestyle = :dash, color =:red, label = "", lw = 2) - +plt_unemp = plot(title = "Unemployment Rate", 1:T, x_path[1, :], color = :blue, + grid = true, label = "", bg_inside = :lightgrey, lw = 2) +plot!(plt_unemp, [x0[1]], linetype = :hline, linestyle = :dash, color = :red, + label = "", lw = 2) + +plt_emp = plot(title = "Employment Rate", 1:T, x_path[2, :], color = :blue, + grid = true, label = "", bg_inside = :lightgrey, lw = 2) +plot!(plt_emp, [x0[2]], linetype = :hline, linestyle = :dash, color = :red, + label = "", lw = 2) plot(plt_unemp, plt_emp, layout = (2, 1), size = (800, 600)) ``` @@ -1014,8 +1027,8 @@ plot(plt_unemp, plt_emp, layout = (2, 1), size = (800, 600)) tags: [remove-cell] --- @testset begin - #test x_path[1,3] ≈ 0.06791408368459205 - #test x_path[2,7] ≈ 0.9429334437639298 + @test x_path[1,3] ≈ 0.06791408368459205 + @test x_path[2,7] ≈ 0.9429334437639298 end ``` diff --git a/lectures/multi_agent_models/lucas_model.md b/lectures/multi_agent_models/lucas_model.md index ac282a27..50082321 100644 --- a/lectures/multi_agent_models/lucas_model.md +++ b/lectures/multi_agent_models/lucas_model.md @@ -386,69 +386,54 @@ using Test ``` ```{code-cell} julia -using LinearAlgebra, Statistics -using Distributions, Interpolations, LaTeXStrings, Parameters, Plots, Random +using LinearAlgebra, Statistics, Random +using Distributions, Interpolations, LaTeXStrings, Plots, NLsolve ``` ```{code-cell} julia -# model -function LucasTree(;γ = 2.0, - β = 0.95, - α = 0.9, - σ = 0.1, - grid_size = 100) - - ϕ = LogNormal(0.0, σ) - shocks = rand(ϕ, 500) +function LucasTree(; gamma = 2.0, + beta = 0.95, + alpha = 0.9, + sigma = 0.1, + grid_size = 100, + num_z = 500) + phi = LogNormal(0.0, sigma) + z = rand(phi, num_z) # build a grid with mass around stationary distribution - ssd = σ / sqrt(1 - α^2) - grid_min, grid_max = exp(-4ssd), exp(4ssd) + ssd = sigma / sqrt(1 - alpha^2) + grid_min = exp(-4 * ssd) + grid_max = exp(4 * ssd) 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 .* z) .^ (1 - gamma)) end - return (γ = γ, β = β, α = α, σ = σ, ϕ = ϕ, 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 - z = lt.shocks - - Af = LinearInterpolation(grid, f, extrapolation_bc=Line()) - - Tf = [ h[i] + β * mean(Af.(grid[i]^α .* z)) for i in 1:length(grid) ] - return Tf + return (; gamma, beta, alpha, sigma, phi, grid, z, h) end # get equilibrium price for Lucas tree -function solve_lucas_model(lt; - tol = 1e-6, - max_iter = 500) - - (;grid, γ) = lt - - i = 0 - f = zero(grid) # Initial guess of f - error = tol + 1 - - while (error > tol) && (i < max_iter) - f_new = lucas_operator(lt, f) - error = maximum(abs, f_new - f) - f = f_new - i += 1 +function solve_lucas_model(lt; ftol = 1e-8, iterations = 500) + (; grid, gamma, alpha, beta, h, z) = lt + + # approximate Lucas operator, which returns the updated function Tf on the grid + function T(f) + Af = linear_interpolation(grid, f, extrapolation_bc = Line()) + # Using z for monte-carlo integration + Tf = [h[i] + beta * mean(Af.(grid[i]^alpha .* z)) + for i in 1:length(grid)] + return Tf end - # p(y) = f(y) * y ^ γ - price = f .* grid.^γ + sol = fixedpoint(T, zero(grid); ftol, iterations) + converged(sol) || error("Failed to converge in $(sol.iterations) iter") + f = sol.zero + + price = f .* grid .^ gamma # f(y)*y^gamma return price end @@ -459,7 +444,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); ``` @@ -468,9 +453,9 @@ price_vals = solve_lucas_model(tree); tags: [remove-cell] --- @testset begin - #test price_vals[57] ≈ 44.5077566916004 - #test price_vals[78] ≈ 68.42956586308563 - #test price_vals[13] ≈ 9.880376662058682 + @test price_vals[57] ≈ 41.35601991726328 + @test price_vals[78] ≈ 63.474988006734925 + @test price_vals[13] ≈ 9.24450665849126 end ``` @@ -519,11 +504,11 @@ Random.seed!(42); ```{code-cell} julia plot() -for β in (.95, 0.98) - tree = LucasTree(;β = β) +for beta in (0.95, 0.98) + tree = LucasTree(; 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,9 +520,9 @@ 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)) - #test price_vals[20] ≈ 35.00073581199659 - #test price_vals[57] ≈ 124.32987344509688 + price_vals = solve_lucas_model(LucasTree(beta = 0.98)) + @test price_vals[20] ≈ 32.17514173843709 + @test price_vals[57] ≈ 113.89100588660085 end ``` diff --git a/lectures/multi_agent_models/markov_asset.md b/lectures/multi_agent_models/markov_asset.md index 68723095..ccf296ba 100644 --- a/lectures/multi_agent_models/markov_asset.md +++ b/lectures/multi_agent_models/markov_asset.md @@ -290,13 +290,12 @@ The next figure shows a simulation, where --- tags: [remove-cell] --- -using Test, Random +using Test ``` ```{code-cell} julia -using LinearAlgebra, Statistics -using LaTeXStrings, Parameters, Plots, QuantEcon - +using LinearAlgebra, Statistics, Random +using LaTeXStrings, Plots, QuantEcon, NLsolve ``` ```{code-cell} julia @@ -317,7 +316,7 @@ d_series = cumprod(g_series) # assumes d_0 = 1 series = [x_series g_series d_series log.(d_series)] labels = [L"X_t" L"g_t" L"d_t" L"ln(d_t)"] -plot(series, layout = 4, labels = labels) +plot(series; layout = 4, labels) ``` ```{code-cell} julia @@ -325,10 +324,9 @@ plot(series, layout = 4, labels = labels) tags: [remove-cell] --- @testset begin - #test x_series[4] ≈ -0.669642857142857 - #test g_series[5] ≈ 0.4094841251523643 - #test d_series[9] ≈ 0.03514706070454392 # Near the inflection point. - #test log.(d_series)[72] ≈ -29.24107142857142 # Something near the end. + @test x_series[4] ≈ -0.22321428571428567 + @test g_series[6] ≈ 1.2500884211272658 + @test d_series[9] ≈ 0.5118913634883987 end ``` @@ -404,20 +402,15 @@ 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, - lw = 2, - ylabel = "price-dividend ratio", - xlabel = "state", - alpha = 0.7, - label = L"v") +plot(mc.state_values, v; lw = 2, ylabel = "price-dividend ratio", + xlabel = L"X_t", alpha = 0.7, label = L"v") ``` ```{code-cell} julia @@ -425,10 +418,8 @@ plot(mc.state_values, tags: [remove-cell] --- @testset begin - #test v[2] ≈ 3.4594684257743284 atol = 1e-7 - #test v[1] ≈ 3.2560393349907755 - #test v[5] ≈ 4.526909446326235 - #test K[8] ≈ 8.887213530262768e-10 + @test v[2] ≈ 3.4594684257743284 + @test K[8] ≈ 8.887213530262768e-10 end ``` @@ -536,45 +527,27 @@ Assuming that the spectral radius of $J$ is strictly less than $\beta^{-1}$, thi v = (I - \beta J)^{-1} \beta J {\mathbb 1} ``` -We will define a function tree_price to solve for $v$ given parameters stored in +We will define a function `tree_price` to solve for $v$ given parameters stored in the AssetPriceModel objects +The default Markov Chain for will be a discretized AR(1) with $\rho = 0.9, \sigma = 0.02$ and discretized into 25 states using Tauchen's method. + ```{code-cell} julia -# A default Markov chain for the state process -ρ = 0.9 -σ = 0.02 -n = 25 -default_mc = tauchen(n, ρ, σ) - -AssetPriceModel = @with_kw (β = 0.96, - γ = 2.0, - mc = default_mc, - n = size(mc.p)[1], - g = exp) - -# test stability of matrix Q -function test_stability(ap, Q) - sr = maximum(abs, eigvals(Q)) - if sr ≥ 1 / ap.β - msg = "Spectral radius condition failed with radius = $sr" - throw(ArgumentError(msg)) - end +function AssetPriceModel(; beta = 0.96, gamma = 2.0, g = exp, + mc = tauchen(25, 0.9, 0.02)) + return (; beta, gamma, mc, g) end # price/dividend ratio of the Lucas tree -function tree_price(ap; γ = ap.γ) - # Simplify names, set up matrices - (;β, mc) = ap - P, y = mc.p, mc.state_values - y = reshape(y, 1, ap.n) - J = P .* ap.g.(y).^(1 - γ) - - # Make sure that a unique solution exists - test_stability(ap, J) +function tree_price(ap) + (; beta, mc, gamma, g) = ap + P = mc.p + y = mc.state_values' + J = P .* g.(y) .^ (1 - gamma) + @assert maximum(abs, eigvals(J)) < 1 / beta # check stability # Compute v - v = (I - β * J) \ sum(β * J, dims = 2) - + v = (I - beta * J) \ sum(beta * J, dims = 2) return v end ``` @@ -583,25 +556,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] -ap = AssetPriceModel() -states = ap.mc.state_values - -lines = [] -labels = [] - -for γ in γs - v = tree_price(ap, γ = γ) - label = L"\gamma = %$γ" - push!(labels, label) - push!(lines, v) +gammas = [1.2, 1.4, 1.6, 1.8, 2.0] +p = plot(title = "Price-dividend ratio as a function of the state", + xlabel = L"X_t", ylabel = "price-dividend ratio") + +for gamma in gammas + ap = AssetPriceModel(; gamma) + states = ap.mc.state_values + plot!(states, tree_price(ap); label = L"\gamma = %$gamma") end - -plot(lines, - labels = reshape(labels, 1, length(labels)), - title = "Price-dividend ratio as a function of the state", - ylabel = "price-dividend ratio", - xlabel = "state") +p ``` ```{code-cell} julia @@ -609,10 +573,10 @@ plot(lines, tags: [remove-cell] --- @testset begin - #test lines[2][4] ≈ 33.36574362637905 - #test lines[3][12] ≈ 28.52560591264372 - #test lines[4][18] ≈ 22.38597470787489 - #test lines[5][24] ≈ 15.81947255704859 + ap = AssetPriceModel() + v = tree_price(ap) + @test v[5] ≈ 74.54830456854316 + @test v[16] ≈ 29.426651157109678 end ``` @@ -687,24 +651,23 @@ yields the solution p = (I - \beta M)^{-1} \beta M \zeta {\mathbb 1} ``` -The above is implemented in the function consol_price +The above is implemented in the function `consol_price` ```{code-cell} julia -function consol_price(ap, ζ) - # Simplify names, set up matrices - (;β, γ, mc, g, n) = ap - P, y = mc.p, mc.state_values - y = reshape(y, 1, n) - M = P .* g.(y).^(-γ) - - # Make sure that a unique solution exists - test_stability(ap, M) +function consol_price(ap, zeta) + (; beta, gamma, mc, g) = ap + P = mc.p + y = mc.state_values' + M = P .* g.(y) .^ (-gamma) + @assert maximum(abs, eigvals(M)) < 1 / beta # Compute price - return (I - β * M) \ sum(β * ζ * M, dims = 2) + return (I - beta * M) \ sum(beta * zeta * M, dims = 2) end ``` +Note that the `sum(Q, dims=2)` is equivalent to `Q * ones(size(Q)[1], 1)`. + ### Pricing an Option to Purchase the Consol Let's now price options of varying maturity that give the right to purchase a consol at a price $p_S$. @@ -771,52 +734,46 @@ T w = \max \{ \beta M w,\; p - p_S {\mathbb 1} \} $$ -Start at some initial $w$ and iterate to convergence with $T$. +Start at some initial $w$ and iterate to convergence with $T$, or use a fixed point algorithm. 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) - - # Simplify names, set up matrices - (;β, γ, mc, g, n) = ap - P, y = mc.p, mc.state_values - y = reshape(y, 1, n) - M = P .* g.(y).^(-γ) - - # Make sure that a unique console price exists - test_stability(ap, M) - - # Compute option price - p = consol_price(ap, ζ) - w = zeros(ap.n, 1) - error = ϵ + 1 - while (error > ϵ) - # Maximize across columns - w_new = max.(β * M * w, p .- p_s) - # Find maximal difference of each component and update - error = maximum(abs, w - w_new) - w = w_new - end - - return w +function call_option(ap, zeta, p_s) + (; beta, gamma, mc, g) = ap + P = mc.p + y = mc.state_values' + M = P .* g.(y) .^ (-gamma) + @assert maximum(abs, eigvals(M)) < 1 / beta + + # Find consol prices + p = consol_price(ap, zeta) + + # Operator for fixed point, using consol prices + T(w) = max.(beta * M * w, p .- p_s) + + # Compute option price as fixed point + sol = fixedpoint(T, zeros(length(y), 1)) + converged(sol) || error("Failed to converge in $(sol.iterations) iter") + return sol.zero 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") +plot(x, p, color = "blue", lw = 2, xlabel = L"X_t", label = "consol price") +plot!(x, w, color = "green", lw = 2, + label = "value of call option with strike at $strike_price") ``` ```{code-cell} julia @@ -824,8 +781,8 @@ plot!(x, w, color = "green", lw = 2, label = "value of call option") tags: [remove-cell] --- @testset begin - #test p[17] ≈ 9.302197030956606 - #test w[20] ≈ 0.46101660813737866 + @test p[17] ≈ 9.302197030956606 + @test w[20] ≈ 0.4610168409491948 end ``` @@ -836,15 +793,6 @@ where the consol prices is high --- will eventually be visited. The reason is that $\beta=0.9$, so the future is discounted relatively rapidly -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test x[2] ≈ -0.126178653628809 - #test p[5] ≈ 52.85568616593254 -end -``` ### Risk Free Rates @@ -899,9 +847,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). @@ -961,25 +909,25 @@ Is one higher than the other? Can you give intuition? ```{code-cell} julia n = 5 -P = fill(0.0125, n, n) + (0.95 - 0.0125)I +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 +g = x -> x # identity +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, mc, gamma, g) ``` +Lucas tree prices are ```{code-cell} julia -v = tree_price(ap) -println("Lucas Tree Prices: $v\n") +tree_price(ap) ``` ```{code-cell} julia @@ -987,17 +935,18 @@ println("Lucas Tree Prices: $v\n") tags: [remove-cell] --- @testset begin - #test v[2] ≈ 21.935706611219704 + v = tree_price(ap) + @test v[2] ≈ 21.935706611219704 end ``` +Consol Bond Prices ```{code-cell} julia -v_consol = consol_price(ap, 1.0) -println("Consol Bond Prices: $(v_consol)\n") +consol_price(ap, 1.0) ``` ```{code-cell} julia -w = call_option(ap, ζ, p_s) +w = call_option(ap, zeta, p_s) ``` ```{code-cell} julia @@ -1005,8 +954,8 @@ w = call_option(ap, ζ, p_s) tags: [remove-cell] --- @testset begin - #test v_consol[1] ≈ 753.8710047641985 - #test w[2][1] ≈ 176.83933430191294 + @test consol_price(ap, 1.0)[1] ≈ 753.8710047641985 + @test w[2] ≈ 176.83933430191294 end ``` @@ -1015,23 +964,20 @@ end Here's a suitable function: ```{code-cell} julia -function finite_horizon_call_option(ap, ζ, p_s, k) - - # Simplify names, set up matrices - (;β, γ, mc) = ap - P, y = mc.p, mc.state_values - y = y' - M = P .* ap.g.(y).^(- γ) - - # Make sure that a unique console price exists - test_stability(ap, M) +function finite_horizon_call_option(ap, zeta, p_s, k) + (; beta, gamma, mc) = ap + P = mc.p + y = mc.state_values' + M = P .* ap.g.(y) .^ (-gamma) + @assert maximum(abs, eigvals(M)) < 1 / beta # Compute option price - p = consol_price(ap, ζ) - w = zeros(ap.n, 1) + p = consol_price(ap, zeta) + + w = zeros(length(y), 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 @@ -1043,30 +989,18 @@ end tags: [remove-cell] --- @testset begin - #test p[3] ≈ 70.00064625026326 - #test w[2] ≈ 176.83933430191294 + @test finite_horizon_call_option(ap, zeta, p_s, 5)[3] ≈ 31.798938780647198 end ``` ```{code-cell} julia -lines = [] -labels = [] +p = plot(title = "Value Finite Horizon Call Option", xlabel = L"t", + ylabel = "value") for k in [5, 25] - w = finite_horizon_call_option(ap, ζ, p_s, k) - push!(lines, w) - push!(labels, L"k = %$k") -end -plot(lines, labels = reshape(labels, 1, length(labels))) -``` - -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test lines[1][4] ≈ 29.859285398252347 - #test lines[2][2] ≈ 147.00074548801277 + w = finite_horizon_call_option(ap, zeta, p_s, k) + plot!(w; label = L"k = %$k") end +p ``` Not surprisingly, the option has greater value with larger $k$. diff --git a/lectures/multi_agent_models/markov_perf.md b/lectures/multi_agent_models/markov_perf.md index e70da347..3f062158 100644 --- a/lectures/multi_agent_models/markov_perf.md +++ b/lectures/multi_agent_models/markov_perf.md @@ -432,28 +432,28 @@ 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) +A = I + zeros(3, 3) B1 = [0.0, 1.0, 0.0] B2 = [0.0, 0.0, 1.0] -R1 = [ 0.0 -a0 / 2.0 0.0; - -a0 / 2.0 a1 a1 / 2.0; - 0.0 a1 / 2.0 0.0] +R1 = [0.0 -a0/2.0 0.0; + -a0/2.0 a1 a1/2.0; + 0.0 a1/2.0 0.0] -R2 = [ 0.0 0.0 -a0 / 2.0; - 0.0 0.0 a1 / 2.0; - -a0 / 2.0 a1 / 2.0 a1] +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=β) +Lambda1 = A - (B2 * F2) +lq1 = QuantEcon.LQ(Q1, R1, Lambda1, 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 Lambda1[1, 1] ≈ 1.0 && Lambda1[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 @@ -504,7 +504,7 @@ This is close enough for rock and roll, as they say in the trade. Indeed, isapprox agrees with our assessment ```{code-cell} julia -isapprox(F1, F1_ih, atol=1e-7) +isapprox(F1, F1_ih, atol = 1e-7) ``` ### Dynamics @@ -522,22 +522,21 @@ The following program ```{code-cell} julia using LaTeXStrings, Plots - AF = A - B1 * F1 - B2 * F2 n = 20 x = zeros(3, n) x[:, 1] = [1 1 1] -for t in 1:n-1 - x[:, t+1] = AF * x[:, t] +for t in 1:(n - 1) + x[:, t + 1] = AF * x[:, t] end q1 = x[2, :] q2 = x[3, :] q = q1 + q2 # total output, MPE p = a0 .- a1 * q # price, MPE -plt = plot(q, color=:blue, lw=2, alpha=0.75, label="total output") -plot!(plt, p, color=:green, lw=2, alpha=0.75, label="price") -plot!(plt, title="Output and prices, duopoly MPE") +plt = plot(q, color = :blue, lw = 2, alpha = 0.75, label = "total output") +plot!(plt, p, color = :green, lw = 2, alpha = 0.75, label = "price") +plot!(plt, title = "Output and prices, duopoly MPE") ``` ```{code-cell} julia @@ -651,9 +650,9 @@ 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 -D = [ -1 0.5; - 0.5 -1] +delta = 0.02 +D = [-1 0.5; + 0.5 -1] b = [25, 25] c1 = c2 = [1, -2, 1] e1 = e2 = [10, 10, 3] @@ -683,28 +682,28 @@ 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) B1 = [0.0, 1.0, 0.0] B2 = [0.0, 0.0, 1.0] -R1 = [ 0.0 -a0 / 2.0 0.0; - -a0 / 2.0 a1 a1 / 2.0; - 0.0 a1 / 2.0 0.0] +R1 = [0.0 -a0/2.0 0.0; + -a0/2.0 a1 a1/2.0; + 0.0 a1/2.0 0.0] -R2 = [ 0.0 0.0 -a0 / 2.0; - 0.0 0.0 a1 / 2.0; - -a0 / 2.0 a1 / 2.0 a1] +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 +782,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) @@ -814,15 +813,19 @@ end Let's have a look at the different time paths ```{code-cell} julia -plt_q = plot(qm, color=:blue, lw=2, alpha=0.75, label="monopolist output") -plot!(plt_q, q, color=:green, lw=2, alpha=0.75, label="MPE total output") -plot!(plt_q, xlabel="time", ylabel="output", ylim=(2,4),legend=:topright) +plt_q = plot(qm, color = :blue, lw = 2, alpha = 0.75, + label = "monopolist output") +plot!(plt_q, q, color = :green, lw = 2, alpha = 0.75, + label = "MPE total output") +plot!(plt_q, xlabel = "time", ylabel = "output", ylim = (2, 4), + legend = :topright) -plt_p = plot(pm, color=:blue, lw=2, alpha=0.75, label="monopolist price") -plot!(plt_p, p, color=:green, lw=2, alpha=0.75, label="MPE price") -plot!(plt_p, xlabel="time", ylabel="price",legend=:topright) +plt_p = plot(pm, color = :blue, lw = 2, alpha = 0.75, + label = "monopolist price") +plot!(plt_p, p, color = :green, lw = 2, alpha = 0.75, label = "MPE price") +plot!(plt_p, xlabel = "time", ylabel = "price", legend = :topright) -plot(plt_q, plt_p, layout=(2,1), size=(700,600)) +plot(plt_q, plt_p, layout = (2, 1), size = (700, 600)) ``` ### Exercise 2 @@ -830,13 +833,13 @@ plot(plt_q, plt_p, layout=(2,1), size=(700,600)) We treat the case $\delta = 0.02$ ```{code-cell} julia -δ = 0.02 -D = [-1 0.5; +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,42 +863,42 @@ 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]; - 0 0 1] - -B1 = δ_1 * [1 -D[1, 1]; - 0 -D[2, 1]; - 0 0] -B2 = δ_1 * [0 -D[1, 2]; - 1 -D[2, 2]; - 0 0] - -R1 = -[0.5 * c1[3] 0 0.5 * c1[2]; - 0 0 0; - 0.5 * c1[2] 0 c1[1]] - -R2 = -[0 0 0; - 0 0.5 * c2[3] 0.5*c2[2]; - 0 0.5 * c2[2] c2[1]] - -Q1 = [-0.5*e1[3] 0; - 0 D[1, 1]] -Q2 = [-0.5*e2[3] 0; - 0 D[2, 2]] +A = [delta_1 0 -delta_1*b[1]; + 0 delta_1 -delta_1*b[2]; + 0 0 1] + +B1 = delta_1 * [1 -D[1, 1]; + 0 -D[2, 1]; + 0 0] +B2 = delta_1 * [0 -D[1, 2]; + 1 -D[2, 2]; + 0 0] + +R1 = -[0.5*c1[3] 0 0.5*c1[2]; + 0 0 0; + 0.5*c1[2] 0 c1[1]] + +R2 = -[0 0 0; + 0 0.5*c2[3] 0.5*c2[2]; + 0 0.5*c2[2] c2[1]] + +Q1 = [-0.5*e1[3] 0; + 0 D[1, 1]] +Q2 = [-0.5*e2[3] 0; + 0 D[2, 2]] S1 = zeros(2, 2) S2 = copy(S1) -W1 = [ 0.0 0.0; - 0.0 0.0; - -0.5 * e1[2] b[1] / 2.0] -W2 = [ 0.0 0.0; - 0.0 0.0; - -0.5 * e2[2] b[2] / 2.0] +W1 = [0.0 0.0; + 0.0 0.0; + -0.5*e1[2] b[1]/2.0] +W2 = [0.0 0.0; + 0.0 0.0; + -0.5*e2[2] b[2]/2.0] -M1 = [0.0 0.0; - 0.0 D[1, 2] / 2.0] +M1 = [0.0 0.0; + 0.0 D[1, 2]/2.0] M2 = copy(M1) ``` @@ -943,16 +946,16 @@ corresponding to $\delta = 0.02$ AF = A - B1 * F1 - B2 * F2 n = 25 x = zeros(3, n) -x[:, 1] = [2 0 1] -for t in 1:(n-1) - x[:, t+1] = AF * x[:, t] +x[:, 1] = [2 0 1] +for t in 1:(n - 1) + x[:, t + 1] = AF * x[:, t] end I1 = x[1, :] I2 = x[2, :] -plot(I1, color=:blue, lw=2, alpha=0.75, label="inventories, firm 1") -plot!(I2, color=:green, lw=2, alpha=0.75, label="inventories, firm 2") -plot!(title=L"\delta = 0.02") +plot(I1, color = :blue, lw = 2, alpha = 0.75, label = "inventories, firm 1") +plot!(I2, color = :green, lw = 2, alpha = 0.75, label = "inventories, firm 2") +plot!(title = L"\delta = 0.02") ``` ```{code-cell} julia diff --git a/lectures/multi_agent_models/matsuyama.md b/lectures/multi_agent_models/matsuyama.md index b0ac9108..ab0a1611 100644 --- a/lectures/multi_agent_models/matsuyama.md +++ b/lectures/multi_agent_models/matsuyama.md @@ -336,13 +336,12 @@ using Test ``` ```{code-cell} julia -using LinearAlgebra, Statistics -using Plots, Parameters +using LinearAlgebra, Statistics, Plots ``` ```{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 +353,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 +362,49 @@ function h_j(j, nk, s1, s2, θ, δ, ρ) return root end -DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≤ s1_ρ) && (n2 ≤ s2_ρ) +function DLL(n1, n2, s1_rho, s2_rho, s1, s2, theta, delta, rho) + (n1 <= s1_rho) && (n2 <= s2_rho) +end -DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≥ h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 ≥ h_j(2, n1, s1, s2, θ, δ, ρ)) +function 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)) +end -DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≥ s1_ρ) && (n2 ≤ h_j(2, n1, s1, s2, θ, δ, ρ)) +function 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)) +end -DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ) = - (n1 ≤ h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 ≥ s2_ρ) +function 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) +end -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, θ, δ, ρ) +function 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) +end -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,16 +414,17 @@ 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 nsync += 1 - # If not, then reset the nsync counter + # If not, then reset the nsync counter else nsync = 0 end @@ -432,20 +440,21 @@ 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 - unit_range = range(0.0, 1.0, length = npts) + unit_range = range(0.0, 1.0, length = npts) # Allocate space to store time to sync time_2_sync = zeros(npts, npts) # 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 +462,17 @@ 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,30 +482,28 @@ 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 end -function pers_till_sync(model, n1_0, n2_0, - maxiter = 500, npers = 3) - # Unpack parameters - @unpack s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = model +function pers_till_sync(model, n1_0, n2_0, maxiter = 500, npers = 3) + (; 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; maxiter = 250, 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 +519,9 @@ 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,20 +576,21 @@ 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, θ, δ, ρ) - ab = create_attraction_basis(model,npts=npts) + unitrange = range(0, 1, length = npts) + model = MSGSync(s1, theta, delta, rho) + ab = create_attraction_basis(model, npts = npts) plt = Plots.heatmap(ab, legend = false) end ``` ```{code-cell} julia params = [[0.5, 2.5, 0.7, 0.2], - [0.5, 2.5, 0.7, 0.4], - [0.5, 2.5, 0.7, 0.6], - [0.5, 2.5, 0.7, 0.8]] + [0.5, 2.5, 0.7, 0.4], + [0.5, 2.5, 0.7, 0.6], + [0.5, 2.5, 0.7, 0.8]] plots = (plot_attraction_basis(p...) for p in params) plot(plots..., size = (1000, 1000)) @@ -593,13 +601,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 9a752c52..19654114 100644 --- a/lectures/multi_agent_models/rational_expectations.md +++ b/lectures/multi_agent_models/rational_expectations.md @@ -595,7 +595,7 @@ using Test ``` ```{code-cell} julia -using QuantEcon, Printf, LinearAlgebra +using QuantEcon, LinearAlgebra ``` To map a problem into a {doc}`discounted optimal linear control problem <../dynamic_programming/lqcontrol>`, we need to define. @@ -662,34 +662,36 @@ 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 0 1] +A = [1 0 0 + 0 kappa1 kappa0 + 0 0 1] B = [1.0, 0.0, 0.0] -R = [ 0 a1/2 -a0/2 - a1/2 0 0 - -a0/2 0 0] +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] - -@printf("F = [%.3f, %.3f, %.3f]\n", F[1], F[2], F[3]) -@printf("(h0, h1, h2) = [%.3f, %.3f, %.3f]\n", h0, h1, h2) +# unpack results +h0 = -F[3] +h1 = 1 - F[1] +h2 = -F[2] +@show F +@show h0, h1, h2; ``` ```{code-cell} julia @@ -748,16 +750,18 @@ candidates = ([94.0886298678, 0.923409232937], [95.0818452486, 0.952459076301]) for (k0, k1) in candidates - A = [1 0 0 + A = [1 0 0 0 k1 k0 - 0 0 1] - lq = QuantEcon.LQ(Q, R, A, B; bet=β) + 0 0 1] + 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] + h0 = -F[3] + h1 = 1 - F[1] + h2 = -F[2] if isapprox(k0, h0; atol = 1e-4) && isapprox(k1, h1 + h2; atol = 1e-4) - @printf("Equilibrium pair = (%.6f, %.6f)\n", k0, k1) - @printf("(h0, h1, h2) = [%.6f, %.6f, %.6f]\n", h0, h1, h2) + println("Equilibrium pair = ($k0, $k1)") + @show h0, h1, h2 end end ``` @@ -830,18 +834,19 @@ $\kappa_0 = -F_1$ and $\kappa_1 = 1-F_0$ A = I + zeros(2, 2) B = [1.0, 0.0] -R = [ a1 / 2.0 -a0 / 2.0 - -a0 / 2.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 = -F[2] +kappa1 = 1 - F[1] +@show kappa0, kappa1; ``` ```{code-cell} julia @@ -849,8 +854,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 ``` @@ -876,18 +881,18 @@ The problem can be solved as follows A = I + zeros(2, 2) B = [1.0, 0.0] -R = [ a1 -a0 / 2.0 - -a0 / 2.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 m0, m1 = -F[2], 1 - F[1] -println("m0=$m0\tm1=$m1") +@show m0, m1; ``` ```{code-cell} julia diff --git a/lectures/multi_agent_models/schelling.md b/lectures/multi_agent_models/schelling.md index 3ef8929a..3bfc53b9 100644 --- a/lectures/multi_agent_models/schelling.md +++ b/lectures/multi_agent_models/schelling.md @@ -8,7 +8,6 @@ kernelspec: language: julia name: julia-1.10 --- - (schelling)= ```{raw} html
@@ -148,7 +147,7 @@ using Test, Random ``` ```{code-cell} julia -using Parameters, Plots, LinearAlgebra, Statistics +using Plots, LinearAlgebra, Statistics ``` @@ -160,7 +159,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 +178,7 @@ function is_happy(a) # partialsortperm(get_distance.(Ref(a), agents), # 1:neighborhood_size)) - return share ≥ preference + return share >= preference end function update!(a) @@ -205,7 +204,8 @@ function plot_distribution(agents) end end - p = scatter(x_vals_0, y_vals_0, color = :orange, markersize = 8, alpha = 0.6) + p = scatter(x_vals_0, y_vals_0, color = :orange, markersize = 8, + alpha = 0.6) scatter!(x_vals_1, y_vals_1, color = :green, markersize = 8, alpha = 0.6) return plot!(legend = :none) end diff --git a/lectures/multi_agent_models/uncertainty_traps.md b/lectures/multi_agent_models/uncertainty_traps.md index dbfc2343..b2a8d799 100644 --- a/lectures/multi_agent_models/uncertainty_traps.md +++ b/lectures/multi_agent_models/uncertainty_traps.md @@ -231,22 +231,26 @@ using Test, Random ```{code-cell} julia using LinearAlgebra, Statistics -using DataFrames, LaTeXStrings, Parameters, Plots +using DataFrames, LaTeXStrings, 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 +function # 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) +end ``` In the results below we use this code to simulate time series for the major variables. @@ -360,18 +364,19 @@ 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!(xlabel = L"\gamma", ylabel = L"\gamma^\prime", legend_title = L"M", legend = :bottomright) +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) ``` ```{code-cell} julia @@ -379,9 +384,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 +401,53 @@ 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,19 +468,20 @@ 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!(xlabel = L"x", ylabel = L"y", legend_title = "Variable", legend = :bottomright) +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)) +plt = plot(layout = (4, 1), size = (600, 600)) for i in 1:4 plot!(plt[i], len, yvals[i], xlabel = L"t", ylabel = vars[i], label = "") @@ -484,7 +494,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 diff --git a/lectures/tools_and_techniques/numerical_linear_algebra.md b/lectures/tools_and_techniques/numerical_linear_algebra.md index f4cda6f9..b4f81c25 100644 --- a/lectures/tools_and_techniques/numerical_linear_algebra.md +++ b/lectures/tools_and_techniques/numerical_linear_algebra.md @@ -54,7 +54,7 @@ The theme of this lecture, and numerical linear algebra in general, comes down t --- tags: [hide-output] --- -using LinearAlgebra, Statistics, BenchmarkTools, SparseArrays, Random, Parameters +using LinearAlgebra, Statistics, BenchmarkTools, SparseArrays, Random Random.seed!(42); # seed random numbers for reproducibility ``` @@ -491,14 +491,11 @@ To manually use the QR decomposition in solving linear least squares: ```{code-cell} julia Af = qr(A) -Q = Af.Q -R = [Af.R; zeros(N - M, M)] # Stack with zeros -@show Q * R ≈ A -x = R \ Matrix(Q)' * b # simplified QR solution for least squares +@show Af.Q * Af.R ≈ A +x = Af.R \ (Matrix(Af.Q)' * b) # simplified QR solution for least squares ``` -This stacks the `R` with zeros, but the more specialized algorithm would not multiply directly -in that way. +See [here](https://discourse.julialang.org/t/qr-decomposition-with-julia-how-to/92508/8) for more, thought keep in mind that more specialized algorithm would be more efficient. In some cases, if an LU is not available for a particular matrix structure, the QR factorization can also be used to solve systems of equations (i.e., not just LLS). This tends to be about 2 times slower than the LU