Multi agent models update (#276)

* Updated formatting and code changes for the multi_agent models
* Julia 1.10 support

---------

Co-authored-by: Jesse Perla <[email protected]>

Author: annabellasd

.JuliaFormatter.toml

@@ -1,4 +1,4 @@
 style = "sciml"
-margin = 84
+margin = 80
 yas_style_nesting = true
 annotate_untyped_fields_with_any = false

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
 ```
 

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]

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.
-

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 <covid_sde>`, we will concentrate on randomness that comes from aggregate changes in behavior or policy.
-

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))
+```