Skip to content

Commit

Permalink
Multi agent models update (#276)
Browse files Browse the repository at this point in the history
* Updated formatting and code changes for the multi_agent models
* Julia 1.10 support

---------

Co-authored-by: Jesse Perla <jesseperla@gmail.com>
  • Loading branch information
annabellasd and jlperla authored Jan 4, 2024
1 parent bf0ae39 commit 6d3b093
Show file tree
Hide file tree
Showing 17 changed files with 860 additions and 913 deletions.
2 changes: 1 addition & 1 deletion .JuliaFormatter.toml
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
style = "sciml"
margin = 84
margin = 80
yas_style_nesting = true
annotate_untyped_fields_with_any = false
4 changes: 4 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
```
Expand Down
40 changes: 19 additions & 21 deletions lectures/Manifest.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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"]
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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"]
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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"]
Expand Down Expand Up @@ -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]
Expand Down
1 change: 0 additions & 1 deletion lectures/continuous_time/covid_sde.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.

1 change: 0 additions & 1 deletion lectures/continuous_time/seir_model.md
Original file line number Diff line number Diff line change
Expand Up @@ -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 <covid_sde>`, we will concentrate on randomness that comes from aggregate changes in behavior or policy.

142 changes: 66 additions & 76 deletions lectures/multi_agent_models/aiyagari.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
```

Expand All @@ -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
Expand All @@ -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
```

Expand All @@ -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)
Expand All @@ -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))
```
Loading

0 comments on commit 6d3b093

Please sign in to comment.