From 8a91e665224033d7f0cd8d9d9bdc20a157f1840c Mon Sep 17 00:00:00 2001 From: Annabella Stoll-Dansereau <38112312+annabellasd@users.noreply.github.com> Date: Thu, 4 Jan 2024 11:07:10 -0800 Subject: [PATCH] Dynamic programming squared update (#277) * updated lectures --------- Co-authored-by: mmcky Co-authored-by: mmcky Co-authored-by: Trang Truong <91804044+htrangtr@users.noreply.github.com> Co-authored-by: Jesse Perla --- format_myst.jl | 4 +- lectures/dynamic_programming_squared/amss.md | 518 +++++++++--------- .../dynamic_programming_squared/dyn_stack.md | 219 ++++---- .../dynamic_programming_squared/lqramsey.md | 248 ++++----- .../opt_tax_recur.md | 351 ++++++------ 5 files changed, 630 insertions(+), 710 deletions(-) diff --git a/format_myst.jl b/format_myst.jl index 0e24e1dd..fcf44d03 100644 --- a/format_myst.jl +++ b/format_myst.jl @@ -58,7 +58,9 @@ function format_myst(input_file_path, output_file_path, extra_replacements = fal replacements = Dict("α" => "alpha", "β" => "beta", "γ" => "gamma", "≤" => "<=", "≥" => ">=", "Σ" => "Sigma", "σ" => "sigma","μ"=>"mu","ϕ"=>"phi","ψ"=>"psi","ϵ"=>"epsilon", "δ"=>"delta","θ" => "theta","ζ"=>"zeta","X̄" => "X_bar","p̄" => "p_bar","x̂" => "x_hat","λ"=>"lambda", - "ρ"=>"rho","u′" => "u_prime" , "f′"=>"f_prime"," ∂u∂c"=>"dudc","Π"=>"Pi","π"=>"pi"," ξ"=>"Xi","c̄"=>"c_bar","w̄"=>"w_bar") + "ρ"=>"rho","u′" => "u_prime" , "f′"=>"f_prime"," ∂u∂c"=>"dudc","Π"=>"Pi","π"=>"pi","ξ"=>"xi","c̄"=>"c_bar", + "w̄"=>"w_bar","Θ" => "Theta", "Ξ" =>"Xi", "Q̃" => "Q_tilde","R̃" => "R_tilde","Ã" => "A_tilde", "B̃" => "B_tilde", + "P̃" => "P_tilde","F̃" => "F_tilde","d̃" => "d_tilde") # Replace the code blocks in the content and handle exceptions try diff --git a/lectures/dynamic_programming_squared/amss.md b/lectures/dynamic_programming_squared/amss.md index 1bc4fee1..5ae14997 100644 --- a/lectures/dynamic_programming_squared/amss.md +++ b/lectures/dynamic_programming_squared/amss.md @@ -44,11 +44,6 @@ In this lecture, we We begin with an introduction to the model. - -```{code-cell} julia -using LinearAlgebra, Statistics -``` - ## Competitive Equilibrium with Distorting Taxes Many but not all features of the economy are identical to those of {doc}`the Lucas-Stokey economy <../dynamic_programming_squared/opt_tax_recur>`. @@ -390,13 +385,15 @@ on optimal taxation with state-contingent debt sequential allocation implementa --- tags: [remove-cell] --- -using Test, Random +using Test ``` ```{code-cell} julia --- tags: [output_scroll] --- + +using LinearAlgebra, Statistics, Random using QuantEcon, NLsolve, NLopt import QuantEcon.simulate @@ -404,10 +401,10 @@ import QuantEcon.simulate mutable struct Model{TF <: AbstractFloat, TM <: AbstractMatrix{TF}, TV <: AbstractVector{TF}} - β::TF - Π::TM + beta::TF + Pi::TM G::TV - Θ::TV + Theta::TV transfers::Bool U::Function Uc::Function @@ -417,7 +414,6 @@ mutable struct Model{TF <: AbstractFloat, n_less_than_one::Bool end - struct SequentialAllocation{TP <: Model, TI <: Integer, TV <: AbstractVector} @@ -426,19 +422,18 @@ struct SequentialAllocation{TP <: Model, S::TI cFB::TV nFB::TV - ΞFB::TV + XiFB::TV zFB::TV end - function SequentialAllocation(model::Model) - β, Π, G, Θ = model.β, model.Π, model.G, model.Θ - mc = MarkovChain(Π) - S = size(Π, 1) # Number of states + beta, Pi, G, Theta = model.beta, model.Pi, model.G, model.Theta + mc = MarkovChain(Pi) + S = size(Pi, 1) # Number of states # Now find the first best allocation - cFB, nFB, ΞFB, zFB = find_first_best(model, S, 1) + cFB, nFB, XiFB, zFB = find_first_best(model, S, 1) - return SequentialAllocation(model, mc, S, cFB, nFB, ΞFB, zFB) + return SequentialAllocation(model, mc, S, cFB, nFB, XiFB, zFB) end @@ -446,13 +441,13 @@ function find_first_best(model::Model, S::Integer, version::Integer) if version != 1 && version != 2 throw(ArgumentError("version must be 1 or 2")) end - β, Θ, Uc, Un, G, Π = - model.β, model.Θ, model.Uc, model.Un, model.G, model.Π + beta, Theta, Uc, Un, G, Pi = + model.beta, model.Theta, model.Uc, model.Un, model.G, model.Pi function res!(out, z) c = z[1:S] n = z[S+1:end] - out[1:S] = Θ .* Uc.(c, n) + Un.(c, n) - out[S+1:end] = Θ .* n .- c .- G + out[1:S] = Theta .* Uc.(c, n) + Un.(c, n) + out[S+1:end] = Theta .* n .- c .- G end res = nlsolve(res!, 0.5 * ones(2 * S)) @@ -463,32 +458,32 @@ function find_first_best(model::Model, S::Integer, version::Integer) if version == 1 cFB = res.zero[1:S] nFB = res.zero[S+1:end] - ΞFB = Uc(cFB, nFB) # Multiplier on the resource constraint - zFB = vcat(cFB, nFB, ΞFB) - return cFB, nFB, ΞFB, zFB + XiFB = Uc(cFB, nFB) # Multiplier on the resource constraint + zFB = vcat(cFB, nFB, XiFB) + return cFB, nFB, XiFB, zFB elseif version == 2 cFB = res.zero[1:S] nFB = res.zero[S+1:end] IFB = Uc(cFB, nFB) .* cFB + Un(cFB, nFB) .* nFB - xFB = \(LinearAlgebra.I - β * Π, IFB) + xFB = \(LinearAlgebra.I - beta * Pi, IFB) zFB = [vcat(cFB[s], xFB[s], xFB) for s in 1:S] return cFB, nFB, IFB, xFB, zFB end end -function time1_allocation(pas::SequentialAllocation, μ::Real) +function time1_allocation(pas::SequentialAllocation, mu::Real) model, S = pas.model, pas.S - Θ, β, Π, G, Uc, Ucc, Un, Unn = - model.Θ, model.β, model.Π, model.G, + Theta, beta, Pi, G, Uc, Ucc, Un, Unn = + model.Theta, model.beta, model.Pi, model.G, model.Uc, model.Ucc, model.Un, model.Unn function FOC!(out, z::Vector) c = z[1:S] n = z[S+1:2S] - Ξ = z[2S+1:end] - out[1:S] = Uc.(c, n) - μ * (Ucc.(c, n) .* c + Uc.(c, n)) - Ξ # FOC c - out[S+1:2S] = Un.(c, n) - μ * (Unn(c, n) .* n .+ Un.(c, n)) + Θ .* Ξ # FOC n - out[2S+1:end] = Θ .* n - c .- G # resource constraint + Xi = z[2S+1:end] + out[1:S] = Uc.(c, n) - mu * (Ucc.(c, n) .* c + Uc.(c, n)) - Xi # FOC c + out[S+1:2S] = Un.(c, n) - mu * (Unn(c, n) .* n .+ Un.(c, n)) + Theta .* Xi # FOC n + out[2S+1:end] = Theta .* n - c .- G # resource constraint return out end # Find the root of the FOC @@ -497,35 +492,35 @@ function time1_allocation(pas::SequentialAllocation, μ::Real) error("Could not find LS allocation.") end z = res.zero - c, n, Ξ = z[1:S], z[S+1:2S], z[2S+1:end] + c, n, Xi = z[1:S], z[S+1:2S], z[2S+1:end] # Now compute x I = Uc(c, n) .* c + Un(c, n) .* n - x = \(LinearAlgebra.I - β * model.Π, I) - return c, n, x, Ξ + x = \(LinearAlgebra.I - beta * model.Pi, I) + return c, n, x, Xi end function time0_allocation(pas::SequentialAllocation, B_::AbstractFloat, s_0::Integer) model = pas.model - Π, Θ, G, β = model.Π, model.Θ, model.G, model.β + Pi, Theta, G, beta = model.Pi, model.Theta, model.G, model.beta Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn # First order conditions of planner's problem function FOC!(out, z) - μ, c, n, Ξ = z[1], z[2], z[3], z[4] - xprime = time1_allocation(pas, μ)[3] + mu, c, n, Xi = z[1], z[2], z[3], z[4] + xprime = time1_allocation(pas, mu)[3] out .= vcat( - Uc(c, n) .* (c - B_) + Un(c, n) .* n + β * dot(Π[s_0, :], xprime), - Uc(c, n) .- μ * (Ucc(c, n) .* (c - B_) + Uc(c, n)) .- Ξ, - Un(c, n) .- μ * (Unn(c, n) .* n + Un(c, n)) + Θ[s_0] .* Ξ, - (Θ .* n .- c .- G)[s_0] + Uc(c, n) .* (c - B_) + Un(c, n) .* n + beta * dot(Pi[s_0, :], xprime), + Uc(c, n) .- mu * (Ucc(c, n) .* (c - B_) + Uc(c, n)) .- Xi, + Un(c, n) .- mu * (Unn(c, n) .* n + Un(c, n)) + Theta[s_0] .* Xi, + (Theta .* n .- c .- G)[s_0] ) end # Find root - res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.ΞFB[s_0]]) + res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.XiFB[s_0]]) if res.f_converged == false error("Could not find time 0 LS allocation.") end @@ -533,18 +528,18 @@ function time0_allocation(pas::SequentialAllocation, end -function time1_value(pas::SequentialAllocation, μ::Real) +function time1_value(pas::SequentialAllocation, mu::Real) model = pas.model - c, n, x, Ξ = time1_allocation(pas, μ) + c, n, x, Xi = time1_allocation(pas, mu) U_val = model.U.(c, n) - V = \(LinearAlgebra.I - model.β*model.Π, U_val) + V = \(LinearAlgebra.I - model.beta*model.Pi, U_val) return c, n, x, V end -function Τ(model::Model, c::Union{Real,Vector}, n::Union{Real,Vector}) +function Omega(model::Model, c::Union{Real,Vector}, n::Union{Real,Vector}) Uc, Un = model.Uc.(c, n), model.Un.(c, n) - return 1. .+ Un./(model.Θ .* Uc) + return 1. .+ Un./(model.Theta .* Uc) end @@ -554,7 +549,7 @@ function simulate(pas::SequentialAllocation, sHist::Union{Vector, Nothing}=nothing) model = pas.model - Π, β, Uc = model.Π, model.β, model.Uc + Pi, beta, Uc = model.Pi, model.beta, model.Uc if isnothing(sHist) sHist = QuantEcon.simulate(pas.mc, T, init=s_0) @@ -562,26 +557,26 @@ function simulate(pas::SequentialAllocation, cHist = zeros(T) nHist = zeros(T) Bhist = zeros(T) - ΤHist = zeros(T) - μHist = zeros(T) + OmegaHist = zeros(T) + muHist = zeros(T) RHist = zeros(T-1) # time 0 - μ, cHist[1], nHist[1], _ = time0_allocation(pas, B_, s_0) - ΤHist[1] = Τ(pas.model, cHist[1], nHist[1])[s_0] + mu, cHist[1], nHist[1], _ = time0_allocation(pas, B_, s_0) + OmegaHist[1] = Omega(pas.model, cHist[1], nHist[1])[s_0] Bhist[1] = B_ - μHist[1] = μ + muHist[1] = mu # time 1 onward for t in 2:T - c, n, x, Ξ = time1_allocation(pas,μ) + c, n, x, Xi = time1_allocation(pas,mu) u_c = Uc(c,n) s = sHist[t] - ΤHist[t] = Τ(pas.model, c, n)[s] - Eu_c = dot(Π[sHist[t-1],:], u_c) + OmegaHist[t] = Omega(pas.model, c, n)[s] + Eu_c = dot(Pi[sHist[t-1],:], u_c) cHist[t], nHist[t], Bhist[t] = c[s], n[s], x[s] / u_c[s] - RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (β * Eu_c) - μHist[t] = μ + RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (beta * Eu_c) + muHist[t] = mu end - return cHist, nHist, Bhist, ΤHist, sHist, μHist, RHist + return cHist, nHist, Bhist, OmegaHist, sHist, muHist, RHist end @@ -603,7 +598,7 @@ end function BellmanEquation(model::Model, xgrid::AbstractVector, policies0::Vector) - S = size(model.Π, 1) # number of states + S = size(model.Pi, 1) # number of states xbar = [minimum(xgrid), maximum(xgrid)] time_0 = false cf, nf, xprimef = policies0 @@ -618,19 +613,19 @@ function get_policies_time1(T::BellmanEquation, i_x::Integer, x::AbstractFloat, s::Integer, Vf::AbstractArray) model, S = T.model, T.S - β, Θ, G, Π = model.β, model.Θ, model.G, model.Π + beta, Theta, G, Pi = model.beta, model.Theta, model.G, model.Pi U, Uc, Un = model.U, model.Uc, model.Un function objf(z::Vector, grad) c, xprime = z[1], z[2:end] n=c+G[s] Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S] - return -(U(c, n) + β * dot(Π[s, :], Vprime)) + return -(U(c, n) + beta * dot(Pi[s, :], Vprime)) end function cons(z::Vector, grad) c, xprime = z[1], z[2:end] n=c+G[s] - return x - Uc(c, n) * c - Un(c, n) * n - β * dot(Π[s, :], xprime) + return x - Uc(c, n) * c - Un(c, n) * n - beta * dot(Pi[s, :], xprime) end lb = vcat(0, T.xbar[1] * ones(S)) ub = vcat(1 - G[s], T.xbar[2] * ones(S)) @@ -657,18 +652,18 @@ end function get_policies_time0(T::BellmanEquation, B_::AbstractFloat, s0::Integer, Vf::Array) model, S = T.model, T.S - β, Θ, G, Π = model.β, model.Θ, model.G, model.Π + beta, Theta, G, Pi = model.beta, model.Theta, model.G, model.Pi U, Uc, Un = model.U, model.Uc, model.Un function objf(z, grad) c, xprime = z[1], z[2:end] n = c+G[s0] Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S] - return -(U(c, n) + β * dot(Π[s0, :], Vprime)) + return -(U(c, n) + beta * dot(Pi[s0, :], Vprime)) end function cons(z::Vector, grad) c, xprime = z[1], z[2:end] n = c + G[s0] - return -Uc(c, n) * (c - B_) - Un(c, n) * n - β * dot(Π[s0, :], xprime) + return -Uc(c, n) * (c - B_) - Un(c, n) * n - beta * dot(Pi[s0, :], xprime) end lb = vcat(0, T.xbar[1] * ones(S)) ub = vcat(1-G[s0], T.xbar[2] * ones(S)) @@ -962,6 +957,7 @@ The recursive formulation is implemented as follows ```{code-cell} julia +# Interpolations.jl doesn't support irregular grids for splines using DataInterpolations mutable struct BellmanEquation_Recursive{TP <: Model, TI <: Integer, TR <: Real} @@ -976,74 +972,76 @@ mutable struct BellmanEquation_Recursive{TP <: Model, TI <: Integer, TR <: Real} zFB::Vector{Vector{TR}} end - struct RecursiveAllocation{TP <: Model, - TI <: Integer, - TVg <: AbstractVector, - TT <: Tuple} + TI <: Integer, + TVg <: AbstractVector, + TT <: Tuple} model::TP mc::MarkovChain S::TI T::BellmanEquation_Recursive - μgrid::TVg + mugrid::TVg xgrid::TVg Vf::Array policies::TT end - -function RecursiveAllocation(model::Model, μgrid::AbstractArray) +function RecursiveAllocation(model::Model, mugrid::AbstractArray) G = model.G - S = size(model.Π, 1) # number of states - mc = MarkovChain(model.Π) + S = size(model.Pi, 1) # number of states + mc = MarkovChain(model.Pi) # now find the first best allocation - Vf, policies, T, xgrid = solve_time1_bellman(model, μgrid) + Vf, policies, T, xgrid = solve_time1_bellman(model, mugrid) T.time_0 = true # Bellman equation now solves time 0 problem - return RecursiveAllocation(model, mc, S, T, μgrid, xgrid, Vf, policies) + return RecursiveAllocation(model, mc, S, T, mugrid, xgrid, Vf, policies) end - -function solve_time1_bellman(model::Model{TR}, μgrid::AbstractArray) where {TR <: Real} - Π = model.Π - S = size(model.Π, 1) +function solve_time1_bellman(model::Model{TR}, + mugrid::AbstractArray) where {TR <: Real} + Pi = model.Pi + S = size(model.Pi, 1) # First get initial fit from lucas stockey solution. # Need to change things to be ex_ante PP = SequentialAllocation(model) function incomplete_allocation(PP::SequentialAllocation, - μ_::AbstractFloat, - s_::Integer) - c, n, x, V = time1_value(PP, μ_) - return c, n, dot(Π[s_, :], x), dot(Π[s_, :], V) + mu_::AbstractFloat, + s_::Integer) + c, n, x, V = time1_value(PP, mu_) + return c, n, dot(Pi[s_, :], x), dot(Pi[s_, :], V) end cf = Array{Function}(undef, S, S) nf = Array{Function}(undef, S, S) xprimef = Array{Function}(undef, S, S) Vf = Vector{Function}(undef, S) - xgrid = Array{TR}(undef, S, length(μgrid)) + xgrid = Array{TR}(undef, S, length(mugrid)) for s_ in 1:S - c = Array{TR}(undef, length(μgrid), S) - n = Array{TR}(undef, length(μgrid), S) - x = Array{TR}(undef, length(μgrid)) - V = Array{TR}(undef, length(μgrid)) - for (i_μ, μ) in enumerate(μgrid) - c[i_μ, :], n[i_μ, :], x[i_μ], V[i_μ] = - incomplete_allocation(PP, μ, s_) + c = Array{TR}(undef, length(mugrid), S) + n = Array{TR}(undef, length(mugrid), S) + x = Array{TR}(undef, length(mugrid)) + V = Array{TR}(undef, length(mugrid)) + for (i_mu, mu) in enumerate(mugrid) + c[i_mu, :], n[i_mu, :], x[i_mu], V[i_mu] = incomplete_allocation(PP, + mu, + s_) end xprimes = repeat(x, 1, S) xgrid[s_, :] = x - for sprime = 1:S - splc = CubicSpline(c[:, sprime][end:-1:1], x[end:-1:1];extrapolate=true) - spln = CubicSpline(n[:, sprime][end:-1:1], x[end:-1:1];extrapolate=true) - splx = CubicSpline(xprimes[:, sprime][end:-1:1], x[end:-1:1];extrapolate=true) + for sprime in 1:S + splc = CubicSpline(c[:, sprime][end:-1:1], x[end:-1:1]; + extrapolate = true) + spln = CubicSpline(n[:, sprime][end:-1:1], x[end:-1:1]; + extrapolate = true) + splx = CubicSpline(xprimes[:, sprime][end:-1:1], x[end:-1:1]; + extrapolate = true) cf[s_, sprime] = y -> splc(y) nf[s_, sprime] = y -> spln(y) xprimef[s_, sprime] = y -> splx(y) end - splV = CubicSpline(V[end:-1:1], x[end:-1:1];extrapolate=true) + splV = CubicSpline(V[end:-1:1], x[end:-1:1]; extrapolate = true) Vf[s_] = y -> splV(y) end @@ -1051,7 +1049,7 @@ function solve_time1_bellman(model::Model{TR}, μgrid::AbstractArray) where {TR # Create xgrid xbar = [maximum(minimum(xgrid)), minimum(maximum(xgrid))] - xgrid = range(xbar[1], xbar[2], length = length(μgrid)) + xgrid = range(xbar[1], xbar[2], length = length(mugrid)) # Now iterate on Bellman equation T = BellmanEquation_Recursive(model, xgrid, policies) @@ -1061,9 +1059,11 @@ function solve_time1_bellman(model::Model{TR}, μgrid::AbstractArray) where {TR Vfnew, policies = fit_policy_function(T, PF, xgrid) diff = 0.0 - for s=1:S - diff = max(diff, maximum(abs, (Vf[s].(xgrid) - Vfnew[s].(xgrid)) ./ - Vf[s].(xgrid))) + for s in 1:S + diff = max(diff, + maximum(abs, + (Vf[s].(xgrid) - Vfnew[s].(xgrid)) ./ + Vf[s].(xgrid))) end println("diff = $diff") @@ -1073,10 +1073,11 @@ function solve_time1_bellman(model::Model{TR}, μgrid::AbstractArray) where {TR return Vf, policies, T, xgrid end - function fit_policy_function(T::BellmanEquation_Recursive, PF::Function, - xgrid::AbstractVector{TF}) where {TF <: AbstractFloat} + xgrid::AbstractVector{TF}) where { + TF <: + AbstractFloat} S = T.S # preallocation PFvec = Array{TF}(undef, 4S + 1, length(xgrid)) @@ -1090,9 +1091,9 @@ function fit_policy_function(T::BellmanEquation_Recursive, for (i_x, x) in enumerate(xgrid) PFvec[:, i_x] = PF(i_x, x, s_) end - splV = CubicSpline(PFvec[1,:], xgrid) + splV = CubicSpline(PFvec[1, :], xgrid) Vf[s_] = y -> splV(y) - for sprime=1:S + for sprime in 1:S splc = CubicSpline(PFvec[1 + sprime, :], xgrid) spln = CubicSpline(PFvec[1 + S + sprime, :], xgrid) splxprime = CubicSpline(PFvec[1 + 2S + sprime, :], xgrid) @@ -1107,18 +1108,16 @@ function fit_policy_function(T::BellmanEquation_Recursive, return Vf, policies end - function Tau(pab::RecursiveAllocation, - c::AbstractArray, - n::AbstractArray) + c::AbstractArray, + n::AbstractArray) model = pab.model Uc, Un = model.Uc(c, n), model.Un(c, n) - return 1. .+ Un ./ (model.Θ .* Uc) + return 1.0 .+ Un ./ (model.Theta .* Uc) end Tau(pab::RecursiveAllocation, c::Real, n::Real) = Tau(pab, [c], [n]) - function time0_allocation(pab::RecursiveAllocation, B_::Real, s0::Integer) T, Vf = pab.T, pab.Vf xbar = T.xbar @@ -1128,12 +1127,14 @@ function time0_allocation(pab::RecursiveAllocation, B_::Real, s0::Integer) return c0, n0, xprime0, T0 end - function simulate(pab::RecursiveAllocation, B_::TF, s_0::Integer, T::Integer, - sHist::Vector=simulate(pab.mc, T, init=s_0)) where {TF <: AbstractFloat} + sHist::Vector = simulate(pab.mc, T, init = s_0)) where { + TF <: + AbstractFloat + } model, mc, Vf, S = pab.model, pab.mc, pab.Vf, pab.S - Π, Uc = model.Π, model.Uc + Pi, Uc = model.Pi, model.Uc cf, nf, xprimef, TTf = pab.policies cHist = Array{TF}(undef, T) @@ -1142,45 +1143,45 @@ function simulate(pab::RecursiveAllocation, xHist = Array{TF}(undef, T) TauHist = Array{TF}(undef, T) THist = Array{TF}(undef, T) - μHist = Array{TF}(undef, T) + muHist = Array{TF}(undef, T) #time0 - cHist[1], nHist[1], xHist[1], THist[1] = time0_allocation(pab, B_, s_0) + cHist[1], nHist[1], xHist[1], THist[1] = time0_allocation(pab, B_, s_0) TauHist[1] = Tau(pab, cHist[1], nHist[1])[s_0] Bhist[1] = B_ - μHist[1] = Vf[s_0](xHist[1]) + muHist[1] = Vf[s_0](xHist[1]) #time 1 onward for t in 2:T - s_, x, s = sHist[t-1], xHist[t-1], sHist[t] + s_, x, s = sHist[t - 1], xHist[t - 1], sHist[t] c = Array{TF}(undef, S) n = Array{TF}(undef, S) xprime = Array{TF}(undef, S) TT = Array{TF}(undef, S) - for sprime=1:S - c[sprime], n[sprime], xprime[sprime], TT[sprime] = - cf[s_, sprime](x), nf[s_, sprime](x), - xprimef[s_, sprime](x), TTf[s_, sprime](x) + for sprime in 1:S + c[sprime], n[sprime], xprime[sprime], TT[sprime] = cf[s_, sprime](x), + nf[s_, sprime](x), + xprimef[s_, + sprime](x), + TTf[s_, sprime](x) end Tau_val = Tau(pab, c, n)[s] u_c = Uc(c, n) - Eu_c = dot(Π[s_, :], u_c) + Eu_c = dot(Pi[s_, :], u_c) - μHist[t] = Vf[s](xprime[s]) + muHist[t] = Vf[s](xprime[s]) - cHist[t], nHist[t], Bhist[t], TauHist[t] = c[s], n[s], x/Eu_c, Tau_val + cHist[t], nHist[t], Bhist[t], TauHist[t] = c[s], n[s], x / Eu_c, Tau_val xHist[t], THist[t] = xprime[s], TT[s] end - return cHist, nHist, Bhist, xHist, TauHist, THist, μHist, sHist + return cHist, nHist, Bhist, xHist, TauHist, THist, muHist, sHist end - function BellmanEquation_Recursive(model::Model{TF}, xgrid::AbstractVector{TF}, policies0::Array) where {TF <: AbstractFloat} - - S = size(model.Π, 1) # number of states + S = size(model.Pi, 1) # number of states xbar = [minimum(xgrid), maximum(xgrid)] time_0 = false z0 = Array{Array}(undef, length(xgrid), S) @@ -1190,18 +1191,18 @@ function BellmanEquation_Recursive(model::Model{TF}, cs = Array{TF}(undef, S) ns = Array{TF}(undef, S) xprimes = Array{TF}(undef, S) - for j = 1:S - cs[j], ns[j], xprimes[j] = cf[s, j](x), nf[s, j](x), xprimef[s, j](x) + for j in 1:S + cs[j], ns[j], xprimes[j] = cf[s, j](x), nf[s, j](x), + xprimef[s, j](x) end z0[i_x, s] = vcat(cs, ns, xprimes, zeros(S)) end end cFB, nFB, IFB, xFB, zFB = find_first_best(model, S, 2) - return BellmanEquation_Recursive(model, S, xbar, time_0, z0, cFB, nFB, xFB, zFB) + return BellmanEquation_Recursive(model, S, xbar, time_0, z0, cFB, nFB, xFB, + zFB) end - - function get_policies_time1(T::BellmanEquation_Recursive, i_x::Integer, x::Real, @@ -1209,44 +1210,45 @@ function get_policies_time1(T::BellmanEquation_Recursive, Vf::AbstractArray{Function}, xbar::AbstractVector) model, S = T.model, T.S - β, Θ, G, Π = model.β, model.Θ, model.G, model.Π - U,Uc,Un = model.U, model.Uc, model.Un + beta, Theta, G, Pi = model.beta, model.Theta, model.G, model.Pi + U, Uc, Un = model.U, model.Uc, model.Un - S_possible = sum(Π[s_, :].>0) - sprimei_possible = findall(Π[s_, :].>0) + S_possible = sum(Pi[s_, :] .> 0) + sprimei_possible = findall(Pi[s_, :] .> 0) function objf(z, grad) - c, xprime = z[1:S_possible], z[S_possible+1:2S_possible] - n = (c .+ G[sprimei_possible]) ./ Θ[sprimei_possible] + c, xprime = z[1:S_possible], z[(S_possible + 1):(2S_possible)] + n = (c .+ G[sprimei_possible]) ./ Theta[sprimei_possible] Vprime = [Vf[sprimei_possible[si]](xprime[si]) for si in 1:S_possible] - return -dot(Π[s_, sprimei_possible], U.(c, n) + β * Vprime) + return -dot(Pi[s_, sprimei_possible], U.(c, n) + beta * Vprime) end function cons(out, z, grad) - c, xprime, TT = - z[1:S_possible], z[S_possible + 1:2S_possible], z[2S_possible + 1:3S_possible] - n = (c .+ G[sprimei_possible]) ./ Θ[sprimei_possible] + c, xprime, TT = z[1:S_possible], z[(S_possible + 1):(2S_possible)], + z[(2S_possible + 1):(3S_possible)] + n = (c .+ G[sprimei_possible]) ./ Theta[sprimei_possible] u_c = Uc.(c, n) - Eu_c = dot(Π[s_, sprimei_possible], u_c) - out .= x * u_c/Eu_c - u_c .* (c - TT) - Un(c, n) .* n - β * xprime + Eu_c = dot(Pi[s_, sprimei_possible], u_c) + out .= x * u_c / Eu_c - u_c .* (c - TT) - Un(c, n) .* n - beta * xprime end function cons_no_trans(out, z, grad) - c, xprime = z[1:S_possible], z[S_possible + 1:2S_possible] - n = (c .+ G[sprimei_possible]) ./ Θ[sprimei_possible] + c, xprime = z[1:S_possible], z[(S_possible + 1):(2S_possible)] + n = (c .+ G[sprimei_possible]) ./ Theta[sprimei_possible] u_c = Uc.(c, n) - Eu_c = dot(Π[s_, sprimei_possible], u_c) - out .= x * u_c / Eu_c - u_c .* c - Un(c, n) .* n - β * xprime + Eu_c = dot(Pi[s_, sprimei_possible], u_c) + out .= x * u_c / Eu_c - u_c .* c - Un(c, n) .* n - beta * xprime end if model.transfers == true - lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1], zeros(S_possible)) + lb = vcat(zeros(S_possible), ones(S_possible) * xbar[1], + zeros(S_possible)) if model.n_less_than_one == true ub = vcat(ones(S_possible) - G[sprimei_possible], - ones(S_possible) * xbar[2], ones(S_possible)) + ones(S_possible) * xbar[2], ones(S_possible)) else ub = vcat(100 * ones(S_possible), - ones(S_possible) * xbar[2], - 100 * ones(S_possible)) + ones(S_possible) * xbar[2], + 100 * ones(S_possible)) end init = vcat(T.z0[i_x, s_][sprimei_possible], T.z0[i_x, s_][2S .+ sprimei_possible], @@ -1254,9 +1256,10 @@ function get_policies_time1(T::BellmanEquation_Recursive, opt = Opt(:LN_COBYLA, 3S_possible) equality_constraint!(opt, cons, zeros(S_possible)) else - lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1]) + lb = vcat(zeros(S_possible), ones(S_possible) * xbar[1]) if model.n_less_than_one == true - ub = vcat(ones(S_possible)-G[sprimei_possible], ones(S_possible)*xbar[2]) + ub = vcat(ones(S_possible) - G[sprimei_possible], + ones(S_possible) * xbar[2]) else ub = vcat(ones(S_possible), ones(S_possible) * xbar[2]) end @@ -1279,15 +1282,16 @@ function get_policies_time1(T::BellmanEquation_Recursive, (minf, minx, ret) = NLopt.optimize(opt, init) if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED && - ret != :FTOL_REACHED && ret != :MAXTIME_REACHED + ret != :FTOL_REACHED && ret != :MAXTIME_REACHED error("optimization failed: ret = $ret") end T.z0[i_x, s_][sprimei_possible] = minx[1:S_possible] - T.z0[i_x, s_][S .+ sprimei_possible] = minx[1:S_possible] .+ G[sprimei_possible] - T.z0[i_x, s_][2S .+ sprimei_possible] = minx[S_possible .+ 1:2S_possible] + T.z0[i_x, s_][S .+ sprimei_possible] = minx[1:S_possible] .+ + G[sprimei_possible] + T.z0[i_x, s_][2S .+ sprimei_possible] = minx[(S_possible .+ 1):(2S_possible)] if model.transfers == true - T.z0[i_x, s_][3S .+ sprimei_possible] = minx[2S_possible + 1:3S_possible] + T.z0[i_x, s_][3S .+ sprimei_possible] = minx[(2S_possible + 1):(3S_possible)] else T.z0[i_x, s_][3S .+ sprimei_possible] = zeros(S) end @@ -1295,26 +1299,25 @@ function get_policies_time1(T::BellmanEquation_Recursive, return vcat(-minf, T.z0[i_x, s_]) end - function get_policies_time0(T::BellmanEquation_Recursive, B_::Real, s0::Integer, Vf::AbstractArray{Function}, xbar::AbstractVector) model = T.model - β, Θ, G = model.β, model.Θ, model.G + beta, Theta, G = model.beta, model.Theta, model.G U, Uc, Un = model.U, model.Uc, model.Un function objf(z, grad) c, xprime = z[1], z[2] - n = (c + G[s0]) / Θ[s0] - return -(U(c, n) + β * Vf[s0](xprime)) + n = (c + G[s0]) / Theta[s0] + return -(U(c, n) + beta * Vf[s0](xprime)) end - function cons(z,grad) + function cons(z, grad) c, xprime, TT = z[1], z[2], z[3] - n = (c + G[s0]) / Θ[s0] - return -Uc(c, n) * (c - B_ - TT) - Un(c, n) * n - β * xprime + n = (c + G[s0]) / Theta[s0] + return -Uc(c, n) * (c - B_ - TT) - Un(c, n) * n - beta * xprime end cons_no_trans(z, grad) = cons(vcat(z, 0), grad) @@ -1326,13 +1329,13 @@ function get_policies_time0(T::BellmanEquation_Recursive, ub = [100.0, xbar[2], 100.0] end init = vcat(T.zFB[s0][1], T.zFB[s0][3], T.zFB[s0][4]) - init = [0.95124922, -1.15926816, 0.0] + init = [0.95124922, -1.15926816, 0.0] opt = Opt(:LN_COBYLA, 3) equality_constraint!(opt, cons) else lb = [0.0, xbar[1]] if model.n_less_than_one == true - ub = [1-G[s0], xbar[2]] + ub = [1 - G[s0], xbar[2]] else ub = [100, xbar[2]] end @@ -1344,7 +1347,6 @@ function get_policies_time0(T::BellmanEquation_Recursive, init[init .> ub] = ub[init .> ub] init[init .< lb] = lb[init .< lb] - min_objective!(opt, objf) lower_bounds!(opt, lb) upper_bounds!(opt, ub) @@ -1354,14 +1356,14 @@ function get_policies_time0(T::BellmanEquation_Recursive, (minf, minx, ret) = NLopt.optimize(opt, init) if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED && - ret != :FTOL_REACHED - error("optimization failed: ret = $ret") + ret != :FTOL_REACHED + error("optimization failed: ret = $ret") end if model.transfers == true - return -minf, minx[1], minx[1]+G[s0], minx[2], minx[3] + return -minf, minx[1], minx[1] + G[s0], minx[2], minx[3] else - return -minf, minx[1], minx[1]+G[s0], minx[2], 0 + return -minf, minx[1], minx[1] + G[s0], minx[2], 0 end end ``` @@ -1429,31 +1431,30 @@ We assume the same utility parameters as in the {doc}`Lucas-Stokey economy <../d This utility function is implemented in the following constructor ```{code-cell} julia -function crra_utility(; - β = 0.9, - σ = 2.0, - γ = 2.0, - Π = 0.5 * ones(2, 2), - G = [0.1, 0.2], - Θ = ones(Float64, 2), - transfers = false - ) +function CRRAModel(; + beta = 0.9, + sigma = 2.0, + gamma = 2.0, + Pi = 0.5 * ones(2, 2), + G = [0.1, 0.2], + Theta = ones(Float64, 2), + transfers = false) function U(c, n) - if σ == 1.0 + if sigma == 1.0 U = log(c) else - U = (c.^(1.0 - σ) - 1.0) / (1.0 - σ) + U = (c .^ (1.0 - sigma) - 1.0) / (1.0 - sigma) end - return U - n.^(1 + γ) / (1 + γ) + return U - n .^ (1 + gamma) / (1 + gamma) end # Derivatives of utility function - Uc(c,n) = c.^(-σ) - Ucc(c,n) = -σ * c.^(-σ - 1.0) - Un(c,n) = -n.^γ - Unn(c,n) = -γ * n.^(γ - 1.0) + Uc(c, n) = c .^ (-sigma) + Ucc(c, n) = -sigma * c .^ (-sigma - 1.0) + Un(c, n) = -n .^ gamma + Unn(c, n) = -gamma * n .^ (gamma - 1.0) n_less_than_one = false - return Model(β, Π, G, Θ, transfers, - U, Uc, Ucc, Un, Unn, n_less_than_one) + return Model(beta, Pi, G, Theta, transfers, + U, Uc, Ucc, Un, Unn, n_less_than_one) end ``` @@ -1468,31 +1469,31 @@ Paths with circles are histories in which there is peace, while those with triangle denote war. ```{code-cell} julia -time_example = crra_utility(G=[0.1, 0.1, 0.1, 0.2, 0.1, 0.1], - Θ = ones(6)) # Θ can in principle be random +time_example = CRRAModel(;G = [0.1, 0.1, 0.1, 0.2, 0.1, 0.1], + Theta = ones(6)) # Theta can in principle be random -time_example.Π = [ 0.0 1.0 0.0 0.0 0.0 0.0; +time_example.Pi = [0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.5 0.5 0.0; 0.0 0.0 0.0 0.0 0.0 1.0; 0.0 0.0 0.0 0.0 0.0 1.0; 0.0 0.0 0.0 0.0 0.0 1.0] -# Initialize μgrid for value function iteration -μgrid = range(-0.7, 0.01, length = 200) +# Initialize mugrid for value function iteration +mugrid = range(-0.7, 0.01, length = 200) time_example.transfers = true # Government can use transfers time_sequential = SequentialAllocation(time_example) # Solve sequential problem -time_bellman = RecursiveAllocation(time_example, μgrid) +time_bellman = RecursiveAllocation(time_example, mugrid) sHist_h = [1, 2, 3, 4, 6, 6, 6] sHist_l = [1, 2, 3, 5, 6, 6, 6] -sim_seq_h = simulate(time_sequential, 1., 1, 7, sHist_h) -sim_bel_h = simulate(time_bellman, 1., 1, 7, sHist_h) -sim_seq_l = simulate(time_sequential, 1., 1, 7, sHist_l) -sim_bel_l = simulate(time_bellman, 1., 1, 7, sHist_l) +sim_seq_h = simulate(time_sequential, 1.0, 1, 7, sHist_h) +sim_bel_h = simulate(time_bellman, 1.0, 1, 7, sHist_h) +sim_seq_l = simulate(time_sequential, 1.0, 1, 7, sHist_l) +sim_bel_l = simulate(time_bellman, 1.0, 1, 7, sHist_l) using Plots @@ -1500,24 +1501,28 @@ titles = hcat("Consumption", "Labor Supply", "Government Debt", "Tax Rate", "Government Spending", "Output") sim_seq_l_plot = hcat(sim_seq_l[1:3]..., sim_seq_l[4], time_example.G[sHist_l], - time_example.Θ[sHist_l] .* sim_seq_l[2]) + time_example.Theta[sHist_l] .* sim_seq_l[2]) sim_bel_l_plot = hcat(sim_bel_l[1:3]..., sim_bel_l[5], time_example.G[sHist_l], - time_example.Θ[sHist_l] .* sim_bel_l[2]) + time_example.Theta[sHist_l] .* sim_bel_l[2]) sim_seq_h_plot = hcat(sim_seq_h[1:3]..., sim_seq_h[4], time_example.G[sHist_h], - time_example.Θ[sHist_h] .* sim_seq_h[2]) + time_example.Theta[sHist_h] .* sim_seq_h[2]) sim_bel_h_plot = hcat(sim_bel_h[1:3]..., sim_bel_h[5], time_example.G[sHist_h], - time_example.Θ[sHist_h] .* sim_bel_h[2]) -p = plot(size = (920, 750), layout =(3, 2), - xaxis=(0:6), grid=false, titlefont=Plots.font("sans-serif", 10)) + time_example.Theta[sHist_h] .* sim_bel_h[2]) +p = plot(size = (920, 750), layout = (3, 2), + xaxis = (0:6), grid = false, titlefont = Plots.font("sans-serif", 10)) plot!(p, title = titles) -for i=1:6 - plot!(p[i], 0:6, sim_seq_l_plot[:, i], marker=:circle, color=:black, lab="") - plot!(p[i], 0:6, sim_bel_l_plot[:, i], marker=:circle, color=:red, lab="") - plot!(p[i], 0:6, sim_seq_h_plot[:, i], marker=:utriangle, color=:black, lab="") - plot!(p[i], 0:6, sim_bel_h_plot[:, i], marker=:utriangle, color=:red, lab="") +for i in 1:6 + plot!(p[i], 0:6, sim_seq_l_plot[:, i], marker = :circle, color = :black, + lab = "") + plot!(p[i], 0:6, sim_bel_l_plot[:, i], marker = :circle, color = :red, + lab = "") + plot!(p[i], 0:6, sim_seq_h_plot[:, i], marker = :utriangle, color = :black, + lab = "") + plot!(p[i], 0:6, sim_bel_h_plot[:, i], marker = :utriangle, color = :red, + lab = "") end p ``` @@ -1574,21 +1579,21 @@ $$ In accordance, we will re-define our utility function ```{code-cell} julia -function log_utility(;β = 0.9, - ψ = 0.69, - Π = 0.5 * ones(2, 2), - G = [0.1, 0.2], - Θ = ones(2), - transfers = false) +function log_utility(; beta = 0.9, + psi = 0.69, + Pi = 0.5 * ones(2, 2), + G = [0.1, 0.2], + Theta = ones(2), + transfers = false) # Derivatives of utility function - U(c,n) = log(c) + ψ * log(1 - n) - Uc(c,n) = 1 ./ c - Ucc(c,n) = -c.^(-2.0) - Un(c,n) = -ψ ./ (1.0 .- n) - Unn(c,n) = -ψ ./ (1.0 .- n).^2.0 + U(c, n) = log(c) + psi * log(1 - n) + Uc(c, n) = 1 ./ c + Ucc(c, n) = -c .^ (-2.0) + Un(c, n) = -psi ./ (1.0 .- n) + Unn(c, n) = -psi ./ (1.0 .- n) .^ 2.0 n_less_than_one = true - return Model(β, Π, G, Θ, transfers, - U, Uc, Ucc, Un, Unn, n_less_than_one) + return Model(beta, Pi, G, Theta, transfers, + U, Uc, Ucc, Un, Unn, n_less_than_one) end ``` @@ -1604,7 +1609,7 @@ log_example = log_utility() log_example.transfers = true # Government can use transfers log_sequential = SequentialAllocation(log_example) # Solve sequential problem -log_bellman = RecursiveAllocation(log_example, μgrid) # Solve recursive problem +log_bellman = RecursiveAllocation(log_example, mugrid) # Solve recursive problem T = 20 sHist = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1] @@ -1614,22 +1619,26 @@ sim_seq = simulate(log_sequential, 0.5, 1, T, sHist) sim_bel = simulate(log_bellman, 0.5, 1, T, sHist) sim_seq_plot = hcat(sim_seq[1:3]..., - sim_seq[4], log_example.G[sHist], log_example.Θ[sHist] .* sim_seq[2]) + sim_seq[4], log_example.G[sHist], + log_example.Theta[sHist] .* sim_seq[2]) sim_bel_plot = hcat(sim_bel[1:3]..., - sim_bel[5], log_example.G[sHist], log_example.Θ[sHist] .* sim_bel[2]) + sim_bel[5], log_example.G[sHist], + log_example.Theta[sHist] .* sim_bel[2]) #plot policies p = plot(size = (920, 750), layout = grid(3, 2), - xaxis=(0:T), grid=false, titlefont=Plots.font("sans-serif", 10)) + xaxis = (0:T), grid = false, titlefont = Plots.font("sans-serif", 10)) labels = fill(("", ""), 6) labels[3] = ("Complete Market", "Incomplete Market") plot!(p, title = titles) -for i = vcat(collect(1:4), 6) - plot!(p[i], sim_seq_plot[:, i], marker=:circle, color=:black, lab=labels[i][1]) - plot!(p[i], sim_bel_plot[:, i], marker=:utriangle, color=:blue, lab=labels[i][2], - legend=:bottomright) +for i in vcat(collect(1:4), 6) + plot!(p[i], sim_seq_plot[:, i], marker = :circle, color = :black, + lab = labels[i][1]) + plot!(p[i], sim_bel_plot[:, i], marker = :utriangle, color = :blue, + lab = labels[i][2], + legend = :bottomright) end -plot!(p[5], sim_seq_plot[:, 5], marker=:circle, color=:blue, lab="") +plot!(p[5], sim_seq_plot[:, 5], marker = :circle, color = :blue, lab = "") ``` When the government experiences a prolonged period of peace, it is able to reduce @@ -1653,38 +1662,29 @@ Random.seed!(42) ```{code-cell} julia T_long = 200 sim_seq_long = simulate(log_sequential, 0.5, 1, T_long) -sHist_long = sim_seq_long[end-2] +sHist_long = sim_seq_long[end - 2] sim_bel_long = simulate(log_bellman, 0.5, 1, T_long, sHist_long) sim_seq_long_plot = hcat(sim_seq_long[1:4]..., - log_example.G[sHist_long], log_example.Θ[sHist_long] .* sim_seq_long[2]) + log_example.G[sHist_long], + log_example.Theta[sHist_long] .* sim_seq_long[2]) sim_bel_long_plot = hcat(sim_bel_long[1:3]..., sim_bel_long[5], - log_example.G[sHist_long], log_example.Θ[sHist_long] .* sim_bel_long[2]) + log_example.G[sHist_long], + log_example.Theta[sHist_long] .* sim_bel_long[2]) -p = plot(size = (920, 750), layout = (3, 2), xaxis=(0:50:T_long), grid=false, - titlefont=Plots.font("sans-serif", 10)) +p = plot(size = (920, 750), layout = (3, 2), xaxis = (0:50:T_long), + grid = false, + titlefont = Plots.font("sans-serif", 10)) plot!(p, title = titles) -for i = 1:6 - plot!(p[i], sim_seq_long_plot[:, i], color=:black, linestyle=:solid, lab=labels[i][1]) - plot!(p[i], sim_bel_long_plot[:, i], color=:blue, linestyle=:dot, lab=labels[i][2], - legend=:bottomright) +for i in 1:6 + plot!(p[i], sim_seq_long_plot[:, i], color = :black, linestyle = :solid, + lab = labels[i][1]) + plot!(p[i], sim_bel_long_plot[:, i], color = :blue, linestyle = :dot, + lab = labels[i][2], + legend = :bottomright) end p ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - # @test sim_seq_long_plot[50, 3] ≈ 0.3951985593686047 - # @test sim_bel_long_plot[50, 3] ≈ 0.05684753244006188 atol = 1e-2 - # @test sim_seq_long_plot[100, 4] ≈ 0.340233842670859 - # @test sim_bel_long_plot[100, 4] ≈ 0.2093423366870517 atol = 1e-3 - # @test sim_seq_long_plot[200, 2] ≈ 0.5839693539786998 - # @test sim_bel_long_plot[200, 2] ≈ 0.6324036099550768 atol = 1e-3 -end -``` - [^fn_a]: In an allocation that solves the Ramsey problem and that levies distorting taxes on labor, why would the government ever want to hand revenues back to the private sector? It would not in an economy with state-contingent debt, since diff --git a/lectures/dynamic_programming_squared/dyn_stack.md b/lectures/dynamic_programming_squared/dyn_stack.md index 6a0f08f8..4527a9d8 100644 --- a/lectures/dynamic_programming_squared/dyn_stack.md +++ b/lectures/dynamic_programming_squared/dyn_stack.md @@ -914,7 +914,8 @@ dynamic program as outlined above ```{code-cell} julia -using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Statistics, Parameters, Random +using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Statistics, + Random ``` @@ -929,16 +930,12 @@ We define named tuples and default values for the model and solver settings, and instantiate one copy of each ```{code-cell} julia -model = @with_kw (a0 = 10, - a1 = 2, - β = 0.96, - γ = 120., - n = 300) +function model(; a0 = 10, a1 = 2, beta = 0.96, gamma = 120.0, n = 300) + return (; a0, a1, beta, gamma, n) +end # things like tolerances, etc. -settings = @with_kw (tol0 = 1e-8, - tol1 = 1e-16, - tol2 = 1e-2) +settings(; tol0 = 1e-8, tol1 = 1e-16, tol2 = 1e-2) = (; tol0, tol1, tol2) defaultModel = model(); defaultSettings = settings(); @@ -947,21 +944,26 @@ defaultSettings = settings(); Now we can compute the actual policy using the LQ routine from QuantEcon.jl ```{code-cell} julia -@unpack a0, a1, β, γ, n = defaultModel -@unpack tol0, tol1, tol2 = defaultSettings +(; a0, a1, beta, gamma, n) = defaultModel +(; tol0, tol1, tol2) = defaultSettings -βs = [β^x for x = 0:n-1] +betas = [beta^x for x in 0:(n - 1)] Alhs = I + zeros(4, 4); -Alhs[4, :] = [β * a0 / (2 * γ), -β * a1 / (2 * γ), -β * a1 / γ, β] # Euler equation coefficients +Alhs[4, :] = [ + beta * a0 / (2 * gamma), + -beta * a1 / (2 * gamma), + -beta * a1 / gamma, + beta, +] # Euler equation coefficients Arhs = I + zeros(4, 4); -Arhs[3, 4] = 1.; +Arhs[3, 4] = 1.0; Alhsinv = inv(Alhs); A = Alhsinv * Arhs; -B = Alhsinv * [0, 1, 0, 0,]; +B = Alhsinv * [0, 1, 0, 0]; R = [0 -a0/2 0 0; -a0/2 a1 a1/2 0; 0 a1/2 0 0; 0 0 0 0]; -Q = γ; -lq = QuantEcon.LQ(Q, R, A, B, bet=β); +Q = gamma; +lq = QuantEcon.LQ(Q, R, A, B, bet = beta); P, F, d = stationary_values(lq) P22 = P[4:end, 4:end]; @@ -970,17 +972,17 @@ P22inv = inv(P22); H_0_0 = -P22inv * P21; # simulate forward -π_leader = zeros(n); +pi_leader = zeros(n); z0 = [1, 1, 1]; x0 = H_0_0 * z0; y0 = vcat(z0, x0); Random.seed!(1) # for reproducibility yt, ut = compute_sequence(lq, y0, n); -π_matrix = R + F' * Q * F; +pi_matrix = R + F' * Q * F; for t in 1:n - π_leader[t] = -(yt[:, t]' * π_matrix * yt[:, t]); + pi_leader[t] = -(yt[:, t]' * pi_matrix * yt[:, t]) end println("Computed policy for Stackelberg leader: $F") @@ -990,7 +992,7 @@ println("Computed policy for Stackelberg leader: $F") --- tags: [remove-cell] --- -#test F ≈ [-1.580044538772657, 0.29461312747031404, 0.6748093760774972, 6.539705936147515]' +@test F ≈ [-1.580044538772657, 0.29461312747031404, 0.6748093760774972, 6.539705936147515]' ``` ### Implied Time Series for Price and Quantities @@ -1001,12 +1003,12 @@ The following code plots the price and quantities q_leader = yt[2, 1:end]; q_follower = yt[3, 1:end]; q = q_leader + q_follower; -p = a0 .- a1*q; +p = a0 .- a1 * q; -plot(1:n+1, [q_leader, q_follower, p], - title = "Output and Prices, Stackelberg Duopoly", - labels = ["leader output" "follower output" "price"], - xlabel = L"t") +plot(1:(n + 1), [q_leader, q_follower, p], + title = "Output and Prices, Stackelberg Duopoly", + labels = ["leader output" "follower output" "price"], + xlabel = L"t") ``` ### Value of Stackelberg Leader @@ -1017,7 +1019,7 @@ We'll compute it two ways (they give identical answers -- just a check on coding and thinking) ```{code-cell} julia -v_leader_forward = sum(βs .* π_leader); +v_leader_forward = sum(betas .* pi_leader); v_leader_direct = -yt[:, 1]' * P * yt[:, 1]; println("v_leader_forward (forward sim) is $v_leader_forward") @@ -1028,13 +1030,13 @@ println("v_leader_direct is $v_leader_direct") --- tags: [remove-cell] --- -#test v_leader_forward ≈ 150.0316212532547 -#test v_leader_direct ≈ 150.03237147548967 +@test v_leader_forward ≈ 150.0316212532547 +@test v_leader_direct ≈ 150.03237147548967 ``` ```{code-cell} julia # manually check whether P is an approximate fixed point -P_next = (R + F' * Q * F + β * (A - B * F)' * P * (A - B * F)); +P_next = (R + F' * Q * F + beta * (A - B * F)' * P * (A - B * F)); all(P - P_next .< tol0) ``` @@ -1042,14 +1044,14 @@ all(P - P_next .< tol0) --- tags: [remove-cell] --- -#test all(P - P_next .< tol0) +@test all(P - P_next .< tol0) ``` ```{code-cell} julia # manually checks whether two different ways of computing the # value function give approximately the same answer v_expanded = -((y0' * R * y0 + ut[:, 1]' * Q * ut[:, 1] + - β * (y0' * (A - B * F)' * P * (A - B * F) * y0))); + beta * (y0' * (A - B * F)' * P * (A - B * F) * y0))); (v_leader_direct - v_expanded < tol0)[1, 1] ``` @@ -1057,7 +1059,7 @@ v_expanded = -((y0' * R * y0 + ut[:, 1]' * Q * ut[:, 1] + --- tags: [remove-cell] --- -#test (v_leader_direct - v_expanded < tol0)[1, 1] +@test (v_leader_direct - v_expanded < tol0)[1, 1] ``` ## Exhibiting Time Inconsistency of Stackelberg Plan @@ -1085,11 +1087,13 @@ for t in 1:n vt_reset_leader[t] = -yt_reset[:, t]' * P * yt_reset[:, t] end -p1 = plot(1:n+1, [(-F * yt)', (-F * yt_reset)'], labels = ["Stackelberg Leader" L"Continuation Leader at $t$"], - title = "Leader Control Variable", xlabel = L"t"); -p2 = plot(1:n+1, [yt[4, :], yt_reset[4, :]], title = "Follower Control Variable", xlabel = L"t", legend = false); +p1 = plot(1:(n + 1), [(-F * yt)', (-F * yt_reset)'], + labels = ["Stackelberg Leader" L"Continuation Leader at $t$"], + title = "Leader Control Variable", xlabel = L"t"); +p2 = plot(1:(n + 1), [yt[4, :], yt_reset[4, :]], + title = "Follower Control Variable", xlabel = L"t", legend = false); p3 = plot(1:n, [vt_leader, vt_reset_leader], legend = false, - xlabel = L"t", title = "Leader Value Function"); + xlabel = L"t", title = "Leader Value Function"); plot(p1, p2, p3, layout = (3, 1), size = (800, 600)) ``` @@ -1102,14 +1106,14 @@ We check that the recursive **Big** $K$ **, little** $k$ formulation of the foll $\vec q_1$ that we computed when we solved the Stackelberg problem ```{code-cell} julia -Ã = I + zeros(5, 5); -Ã[1:4, 1:4] .= A - B * F; -R̃ = [0 0 0 0 -a0/2; 0 0 0 0 a1/2; 0 0 0 0 0; 0 0 0 0 0; -a0/2 a1/2 0 0 a1]; -Q̃ = Q; -B̃ = [0, 0, 0, 0, 1]; - -lq_tilde = QuantEcon.LQ(Q̃, R̃, Ã, B̃, bet=β); -P̃, F̃, d̃ = stationary_values(lq_tilde); +A_tilde = I + zeros(5, 5); +A_tilde[1:4, 1:4] .= A - B * F; +R_tilde = [0 0 0 0 -a0/2; 0 0 0 0 a1/2; 0 0 0 0 0; 0 0 0 0 0; -a0/2 a1/2 0 0 a1]; +Q_tilde = Q; +B_tilde = [0, 0, 0, 0, 1]; + +lq_tilde = QuantEcon.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, bet = beta); +P_tilde, F_tilde, d_tilde = stationary_values(lq_tilde); y0_tilde = vcat(y0, y0[3]); yt_tilde = compute_sequence(lq_tilde, y0_tilde, n)[1]; ``` @@ -1117,7 +1121,7 @@ yt_tilde = compute_sequence(lq_tilde, y0_tilde, n)[1]; ```{code-cell} julia # checks that the recursive formulation of the follower's problem gives # the same solution as the original Stackelberg problem -plot(1:n+1, [yt_tilde[5, :], yt_tilde[3, :]], labels = [L"\tilde{q}" L"q"]) +plot(1:(n + 1), [yt_tilde[5, :], yt_tilde[3, :]], labels = [L"\tilde{q}" L"q"]) ``` Note: Variables with `_tilde` are obtained from solving the follower's @@ -1132,7 +1136,7 @@ max(abs(yt_tilde[5] - yt_tilde[3])) --- tags: [remove-cell] --- -#test max(abs(yt_tilde[5] - yt_tilde[3])) ≈ 0. atol = 1e-15 +@test max(abs(yt_tilde[5] - yt_tilde[3])) ≈ 0. atol = 1e-15 ``` ```{code-cell} julia @@ -1144,7 +1148,7 @@ yt[:, 1][end] - (yt_tilde[:, 2] - yt_tilde[:, 1])[end] < tol0 --- tags: [remove-cell] --- -#test yt[:, 1][end] - (yt_tilde[:, 2] - yt_tilde[:, 1])[end] < tol0 +@test yt[:, 1][end] - (yt_tilde[:, 2] - yt_tilde[:, 1])[end] < tol0 ``` ### Explanation of Alignment @@ -1159,7 +1163,7 @@ Can you spot what features of $\tilde F$ imply this? Hint: remember the components of $X_t$ ```{code-cell} julia -F̃ # policy function in the follower's problem +F_tilde # policy function in the follower's problem ``` ```{code-cell} julia @@ -1167,43 +1171,43 @@ P # value function in the Stackelberg problem ``` ```{code-cell} julia -P̃ # value function in the follower's problem +P_tilde # value function in the follower's problem ``` ```{code-cell} julia # manually check that P is an approximate fixed point -all((P - ((R + F' * Q * F) + β * (A - B * F)' * P * (A - B * F)) .< tol0)) +all((P - ((R + F' * Q * F) + beta * (A - B * F)' * P * (A - B * F)) .< tol0)) ``` ```{code-cell} julia --- tags: [remove-cell] --- -#test all((P - ((R + F' * Q * F) + β * (A - B * F)' * P * (A - B * F)) .< tol0)) +@test all((P - ((R + F' * Q * F) + beta * (A - B * F)' * P * (A - B * F)) .< tol0)) ``` ```{code-cell} julia # compute `P_guess` using `F_tilde_star` -F̃_star = -[0, 0, 0, 1, 0]'; +F_tilde_star = -[0, 0, 0, 1, 0]'; P_guess = zeros(5, 5); for i in 1:1000 - P_guess = ((R̃ + F̃_star' * Q̃ * F̃_star) + - β * (Ã - B̃ * F̃_star)' * P_guess - * (Ã - B̃ * F̃_star)); + P_guess = ((R_tilde + F_tilde_star' * Q_tilde * F_tilde_star) + + beta * (A_tilde - B_tilde * F_tilde_star)' * P_guess + * (A_tilde - B_tilde * F_tilde_star)) end ``` ```{code-cell} julia # value function in the follower's problem --(y0_tilde' * P̃ * y0_tilde)[1, 1] +-(y0_tilde' * P_tilde * y0_tilde)[1, 1] ``` ```{code-cell} julia --- tags: [remove-cell] --- -@test -(y0_tilde' * P̃ * y0_tilde)[1, 1] ≈ 112.65590740578173 +@test -(y0_tilde' * P_tilde * y0_tilde)[1, 1] ≈ 112.65590740578173 ``` ```{code-cell} julia @@ -1220,8 +1224,8 @@ tags: [remove-cell] ```{code-cell} julia # c policy using policy iteration algorithm -F_iter = (β * inv(Q + β * B̃' * P_guess * B̃) - * B̃' * P_guess * Ã); +F_iter = (beta * inv(Q + beta * B_tilde' * P_guess * B_tilde) + * B_tilde' * P_guess * A_tilde); P_iter = zeros(5, 5); dist_vec = zeros(5, 5); @@ -1229,27 +1233,30 @@ for i in 1:100 # compute P_iter dist_vec = similar(P_iter) for j in 1:1000 - P_iter = (R̃ + F_iter' * Q * F_iter) + β * - (Ã - B̃ * F_iter)' * P_iter * - (Ã - B̃ * F_iter); + P_iter = (R_tilde + F_iter' * Q * F_iter) + + beta * + (A_tilde - B_tilde * F_iter)' * P_iter * + (A_tilde - B_tilde * F_iter) # update F_iter - F_iter = β * inv(Q + β * B̃' * P_iter * B̃) * - B̃' * P_iter * Ã; + F_iter = beta * inv(Q + beta * B_tilde' * P_iter * B_tilde) * + B_tilde' * P_iter * A_tilde - dist_vec = P_iter - ((R̃ + F_iter' * Q * F_iter) + - β * (Ã - B̃ * F_iter)' * P_iter * - (Ã - B̃ * F_iter)); + dist_vec = P_iter - ((R_tilde + F_iter' * Q * F_iter) + + beta * (A_tilde - B_tilde * F_iter)' * P_iter * + (A_tilde - B_tilde * F_iter)) end end if maximum(abs.(dist_vec)) < 1e-8 - dist_vec2 = F_iter - (β * inv(Q + β * B̃' * P_iter * B̃) * B̃' * P_iter * Ã) - if maximum(abs.(dist_vec2)) < 1e-8 - @show F_iter - else - println("The policy didn't converge: try increasing the number of outer loop iterations") - end + dist_vec2 = F_iter - + (beta * inv(Q + beta * B_tilde' * P_iter * B_tilde) * B_tilde' * + P_iter * A_tilde) + if maximum(abs.(dist_vec2)) < 1e-8 + @show F_iter + else + println("The policy didn't converge: try increasing the number of outer loop iterations") + end else println("The policy didn't converge: try increasing the number of inner loop iterations") end @@ -1259,22 +1266,23 @@ end yt_tilde_star = zeros(n, 5); yt_tilde_star[1, :] = y0_tilde; -for t in 1:n-1 - yt_tilde_star[t+1, :] = (Ã - B̃ * F̃_star) * yt_tilde_star[t, :]; +for t in 1:(n - 1) + yt_tilde_star[t + 1, :] = (A_tilde - B_tilde * F_tilde_star) * + yt_tilde_star[t, :] end plot([yt_tilde_star[:, 5], yt_tilde[3, :]], labels = [L"\tilde{q}" L"q"]) ``` ```{code-cell} julia -maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:end-1])) +maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:(end - 1)])) ``` ```{code-cell} julia --- tags: [remove-cell] --- -#test maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:end-1])) < 1e-15 +@test maximum(abs.(yt_tilde_star[:, 5] - yt_tilde[3, 1:end-1])) < 1e-15 ``` ## Markov Perfect Equilibrium @@ -1317,21 +1325,21 @@ B1 = [0, 0, 1]; B2 = [0, 1, 0]; R1 = [0 0 -a0/2; 0 0 a1/2; -a0/2 a1/2 a1]; R2 = [0 -a0/2 0; -a0/2 a1 a1/2; 0 a1/2 0]; -Q1 = Q2 = γ; -S1 = S2 = W1 = W2 = M1 = M2 = 0.; +Q1 = Q2 = gamma; +S1 = S2 = W1 = W2 = M1 = M2 = 0.0; # solve using nnash from QE F1, F2, P1, P2 = nnash(A, B1, B2, R1, R2, Q1, Q2, - S1, S2, W1, W2, M1, M2, - beta = β, - tol = tol1); + S1, S2, W1, W2, M1, M2, + beta = beta, + tol = tol1); # simulate forward AF = A - B1 * F1 - B2 * F2; z = zeros(3, n); z[:, 1] .= 1; -for t in 1:n-1 - z[:, t+1] = AF * z[:, t] +for t in 1:(n - 1) + z[:, t + 1] = AF * z[:, t] end println("Policy for F1 is $F1") @@ -1342,8 +1350,8 @@ println("Policy for F2 is $F2") --- tags: [remove-cell] --- -#test round(F1[1], digits = 4) == -0.227 -#test round(F2[2], digits = 4) == 0.0945 +@test round(F1[1], digits = 4) == -0.227 +@test round(F2[2], digits = 4) == 0.0945 ``` ```{code-cell} julia @@ -1351,7 +1359,8 @@ q1 = z[2, :]; q2 = z[3, :]; q = q1 + q2; # total output, MPE p = a0 .- a1 * q; # total price, MPE -plot([q, p], labels = ["total ouput" "total price"], title = "Output and prices, duopoly MPE", xlabel = L"t") +plot([q, p], labels = ["total ouput" "total price"], + title = "Output and prices, duopoly MPE", xlabel = L"t") ``` ```{code-cell} julia @@ -1370,11 +1379,11 @@ tags: [remove-cell] # compute values u1 = -F1 * z; u2 = -F2 * z; -π_1 = (p .* q1)' - γ * u1.^2; -π_2 = (p .* q2)' - γ * u2.^2; +pi_1 = (p .* q1)' - gamma * u1 .^ 2; +pi_2 = (p .* q2)' - gamma * u2 .^ 2; -v1_forward = π_1 * βs; -v2_forward = π_2 * βs; +v1_forward = pi_1 * betas; +v2_forward = pi_2 * betas; v1_direct = -z[:, 1]' * P1 * z[:, 1]; v2_direct = -z[:, 1]' * P2 * z[:, 1]; @@ -1387,16 +1396,16 @@ println("Firm 2: Direct is $v2_direct, Forward is $(v2_forward[1])"); --- tags: [remove-cell] --- -#test round(v1_direct, digits = 4) == 133.3296 -#test round(v2_direct, digits = 4) == 133.3296 -#test round(v1_forward[1], digits = 4) == 133.3303 -#test round(v2_forward[1], digits = 4) == 133.3303 +@test round(v1_direct, digits = 3) == 133.329 +@test round(v2_direct, digits = 3) == 133.329 +@test round(v1_forward[1], digits = 3) == 133.330 +@test round(v2_forward[1], digits = 3) == 133.330 ``` ```{code-cell} julia # sanity check -Λ_1 = A - B2 * F2; -lq1 = QuantEcon.LQ(Q1, R1, Λ_1, B1, bet = β); +Lambda_1 = A - B2 * F2; +lq1 = QuantEcon.LQ(Q1, R1, Lambda_1, B1, bet = beta); P1_ih, F1_ih, d = stationary_values(lq1); v2_direct_alt = -z[:, 1]' * P1_ih * z[:, 1] + d; @@ -1410,13 +1419,15 @@ vt_MPE = zeros(n); vt_follower = zeros(n); for t in 1:n - vt_MPE[t] = -z[:, t]' * P1 * z[:, t]; - vt_follower[t] = -yt_tilde[:, t]' * P̃ * yt_tilde[:, t]; + vt_MPE[t] = -z[:, t]' * P1 * z[:, t] + vt_follower[t] = -yt_tilde[:, t]' * P_tilde * yt_tilde[:, t] end -plot([vt_MPE, vt_leader, vt_follower], labels = ["MPE" "Stackelberg leader" "Stackelberg follower"], - title = "MPE vs Stackelberg Values", - xlabel = L"t") +plot([vt_MPE, vt_leader, vt_follower], + labels = ["MPE" "Stackelberg leader" "Stackelberg follower"], + title = "MPE vs Stackelberg Values", + xlabel = L"t", + legend = :outertopright) ``` ```{code-cell} julia @@ -1428,7 +1439,7 @@ println("vt_MPE(y0) = $(vt_MPE[1])"); ```{code-cell} julia # total difference in value b/t Stackelberg and MPE -vt_leader[1] + vt_follower[1] - 2*vt_MPE[1] +vt_leader[1] + vt_follower[1] - 2 * vt_MPE[1] ``` ```{code-cell} julia diff --git a/lectures/dynamic_programming_squared/lqramsey.md b/lectures/dynamic_programming_squared/lqramsey.md index 9aef62cb..84a22866 100644 --- a/lectures/dynamic_programming_squared/lqramsey.md +++ b/lectures/dynamic_programming_squared/lqramsey.md @@ -582,8 +582,7 @@ using Test ``` ```{code-cell} julia -using LaTeXStrings, QuantEcon, Plots, LinearAlgebra, Parameters - +using LaTeXStrings, QuantEcon, Plots, LinearAlgebra abstract type AbstractStochProcess end @@ -592,14 +591,13 @@ struct ContStochProcess{TF <: AbstractFloat} <: AbstractStochProcess C::Matrix{TF} end - struct DiscreteStochProcess{TF <: AbstractFloat} <: AbstractStochProcess P::Matrix{TF} x_vals::Matrix{TF} end struct Economy{TF <: AbstractFloat, SP <: AbstractStochProcess} - β::TF + beta::TF Sg::Matrix{TF} Sd::Matrix{TF} Sb::Matrix{TF} @@ -613,65 +611,62 @@ function compute_exog_sequences(econ, x) g, d, b, s = [dropdims(S * x, dims = 1) for S in (Sg, Sd, Sb, Ss)] #= solve for Lagrange multiplier in the govt budget constraint - In fact we solve for ν = λ / (1 + 2*λ). Here ν is the - solution to a quadratic equation a(ν^2 - ν) + b = 0 where + In fact we solve for nu = lambda / (1 + 2*lambda). Here nu is the + solution to a quadratic equation a(nu^2 - nu) + b = 0 where a and b are expected discounted sums of quadratic forms of the state. =# Sm = Sb - Sd - Ss return g, d, b, s, Sm end - -function compute_allocation(econ, Sm, ν, x, b) - Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss - - # solve for the allocation given ν and x - Sc = 0.5 .* (Sb + Sd - Sg - ν .* Sm) - Sl = 0.5 .* (Sb - Sd + Sg - ν .* Sm) +function compute_allocation(econ, Sm, nu, x, b) + (; Sg, Sd, Sb, Ss) = econ + # solve for the allocation given nu and x + Sc = 0.5 .* (Sb + Sd - Sg - nu .* Sm) + Sl = 0.5 .* (Sb - Sd + Sg - nu .* Sm) c = dropdims(Sc * x, dims = 1) l = dropdims(Sl * x, dims = 1) p = dropdims((Sb - Sc) * x, dims = 1) # Price without normalization - τ = 1 .- l ./ (b .- c) - rvn = l .* τ + tau = 1 .- l ./ (b .- c) + rvn = l .* tau - return Sc, Sl, c, l, p, τ, rvn + return Sc, Sl, c, l, p, tau, rvn end - -function compute_ν(a0, b0) +function compute_nu(a0, b0) disc = a0^2 - 4a0 * b0 - if disc ≥ 0 - ν = 0.5 *(a0 - sqrt(disc)) / a0 + if disc >= 0 + nu = 0.5 * (a0 - sqrt(disc)) / a0 else println("There is no Ramsey equilibrium for these parameters.") error("Government spending (economy.g) too low") end # Test that the Lagrange multiplier has the right sign - if ν * (0.5 - ν) < 0 + if nu * (0.5 - nu) < 0 print("Negative multiplier on the government budget constraint.") error("Government spending (economy.g) too low") end - return ν + return nu end - -function compute_Π(B, R, rvn, g, ξ) - π = B[2:end] - R[1:end-1] .* B[1:end-1] - rvn[1:end-1] + g[1:end-1] - Π = cumsum(π .* ξ) - return π, Π +function compute_Pi(B, R, rvn, g, xi) + pi = B[2:end] - R[1:(end - 1)] .* B[1:(end - 1)] - rvn[1:(end - 1)] + + g[1:(end - 1)] + Pi = cumsum(pi .* xi) + return pi, Pi end - -function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, T) +function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, + T) # simplify notation - @unpack β, Sg, Sd, Sb, Ss = econ - @unpack P, x_vals = econ.proc + (; beta, Sg, Sd, Sb, Ss) = econ + (; P, x_vals) = econ.proc mc = MarkovChain(P) - state = simulate(mc, T, init=1) + state = simulate(mc, T, init = 1) x = x_vals[:, state] # Compute exogenous sequence @@ -679,38 +674,36 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:DiscreteStochProcess}, T # compute a0, b0 ns = size(P, 1) - F = I - β.*P - a0 = (F \ ((Sm * x_vals)'.^2))[1] ./ 2 - H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss)*x_vals) + F = I - beta .* P + a0 = (F \ ((Sm * x_vals)' .^ 2))[1] ./ 2 + H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss) * x_vals) b0 = (F \ H')[1] ./ 2 # compute lagrange multiplier - ν = compute_ν(a0, b0) + nu = compute_nu(a0, b0) - # Solve for the allocation given ν and x - Sc, Sl, c, l, p, τ, rvn = compute_allocation(econ, Sm, ν, x, b) + # Solve for the allocation given nu and x + Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b) # compute remaining variables - H = ((Sb - Sc) * x_vals) .* ((Sl - Sg) * x_vals) - (Sl * x_vals).^2 + H = ((Sb - Sc) * x_vals) .* ((Sl - Sg) * x_vals) - (Sl * x_vals) .^ 2 temp = dropdims(F * H', dims = 2) B = temp[state] ./ p H = dropdims(P[state, :] * ((Sb - Sc) * x_vals)', dims = 2) - R = p ./ (β .* H) - temp = dropdims(P[state, :] *((Sb - Sc) * x_vals)', dims = 2) - ξ = p[2:end] ./ temp[1:end-1] + R = p ./ (beta .* H) + temp = dropdims(P[state, :] * ((Sb - Sc) * x_vals)', dims = 2) + xi = p[2:end] ./ temp[1:(end - 1)] - # compute π - π, Π = compute_Π(B, R, rvn, g, ξ) + # compute pi + pi, Pi = compute_Pi(B, R, rvn, g, xi) - return (g = g, d = d, b = b, s = s, c = c, - l = l, p = p, τ = τ, rvn = rvn, B = B, - R = R, π = π, Π = Π, ξ = ξ) + return (; g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi) end function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T) # simplify notation - @unpack β, Sg, Sd, Sb, Ss = econ - @unpack A, C = econ.proc + (; beta, Sg, Sd, Sb, Ss) = econ + (; A, C) = econ.proc # generate an initial condition x0 satisfying x0 = A x0 nx, nx = size(A) @@ -725,7 +718,7 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T) w = randn(nw, T) x[:, 1] = x0 for t in 2:T - x[:, t] = A *x[:, t-1] + C * w[:, t] + x[:, t] = A * x[:, t - 1] + C * w[:, t] end # compute exogenous sequence @@ -733,77 +726,74 @@ function compute_paths(econ::Economy{<:AbstractFloat, <:ContStochProcess}, T) # compute a0 and b0 H = Sm'Sm - a0 = 0.5 * var_quadratic_sum(A, C, H, β, x0) - H = (Sb - Sd + Sg)'*(Sg + Ss) - b0 = 0.5 * var_quadratic_sum(A, C, H, β, x0) + a0 = 0.5 * var_quadratic_sum(A, C, H, beta, x0) + H = (Sb - Sd + Sg)' * (Sg + Ss) + b0 = 0.5 * var_quadratic_sum(A, C, H, beta, x0) # compute lagrange multiplier - ν = compute_ν(a0, b0) + nu = compute_nu(a0, b0) - # solve for the allocation given ν and x - Sc, Sl, c, l, p, τ, rvn = compute_allocation(econ, Sm, ν, x, b) + # solve for the allocation given nu and x + Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b) # compute remaining variables - H = Sl'Sl - (Sb - Sc)' *(Sl - Sg) + H = Sl'Sl - (Sb - Sc)' * (Sl - Sg) L = zeros(T) for t in eachindex(L) - L[t] = var_quadratic_sum(A, C, H, β, x[:, t]) + L[t] = var_quadratic_sum(A, C, H, beta, x[:, t]) end B = L ./ p - Rinv = dropdims(β .* (Sb- Sc)*A*x, dims = 1) ./ p + Rinv = dropdims(beta .* (Sb - Sc) * A * x, dims = 1) ./ p R = 1 ./ Rinv AF1 = (Sb - Sc) * x[:, 2:end] - AF2 = (Sb - Sc) * A * x[:, 1:end-1] - ξ = AF1 ./ AF2 - ξ = dropdims(ξ, dims = 1) + AF2 = (Sb - Sc) * A * x[:, 1:(end - 1)] + xi = AF1 ./ AF2 + xi = dropdims(xi, dims = 1) - # compute π - π, Π = compute_Π(B, R, rvn, g, ξ) + # compute pi + pi, Pi = compute_Pi(B, R, rvn, g, xi) - return (g = g, d = d, b = b, s = s, - c = c, l = l, p = p, τ = τ, - rvn = rvn, B = B, R = R, - π = π, Π = Π, ξ = ξ) + return (; g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi) end function gen_fig_1(path) T = length(path.c) - plt_1 = plot(path.rvn, lw=2, label = L"\tau_t l_t") - plot!(plt_1, path.g, lw=2, label= L"g_t") - plot!(plt_1, path.c, lw=2, label= L"c_t") - plot!(xlabel="Time", grid=true) + plt_1 = plot(path.rvn, lw = 2, label = L"\tau_t l_t") + plot!(plt_1, path.g, lw = 2, label = L"g_t") + plot!(plt_1, path.c, lw = 2, label = L"c_t") + plot!(xlabel = "Time", grid = true) - plt_2 = plot(path.rvn, lw=2, label=L"\tau_t l_t") - plot!(plt_2, path.g, lw=2, label=L"g_t") - plot!(plt_2, path.B[2:end], lw=2, label=L"B_{t+1}") - plot!(xlabel="Time", grid=true) + plt_2 = plot(path.rvn, lw = 2, label = L"\tau_t l_t") + plot!(plt_2, path.g, lw = 2, label = L"g_t") + plot!(plt_2, path.B[2:end], lw = 2, label = L"B_{t+1}") + plot!(xlabel = "Time", grid = true) - plt_3 = plot(path.R, lw=2, label=L"R_{t-1}") - plot!(plt_3, xlabel="Time", grid=true) + plt_3 = plot(path.R, lw = 2, label = L"R_{t-1}") + plot!(plt_3, xlabel = "Time", grid = true) - plt_4 = plot(path.rvn, lw=2, label=L"\tau_t l_t") - plot!(plt_4, path.g, lw=2, label=L"g_t") - plot!(plt_4, path.π, lw=2, label=L"\pi_t") - plot!(plt_4, xlabel="Time", grid=true) + plt_4 = plot(path.rvn, lw = 2, label = L"\tau_t l_t") + plot!(plt_4, path.g, lw = 2, label = L"g_t") + plot!(plt_4, path.pi, lw = 2, label = L"\pi_t") + plot!(plt_4, xlabel = "Time", grid = true) - plot(plt_1, plt_2, plt_3, plt_4, layout=(2,2), size = (800,600)) + plot(plt_1, plt_2, plt_3, plt_4, layout = (2, 2), size = (800, 600)) end function gen_fig_2(path) - T = length(path.c) - paths = [path.ξ, path.Π] + paths = [path.xi, path.Pi] labels = [L"\xi_t", L"\Pi_t"] plt_1 = plot() plt_2 = plot() plots = [plt_1, plt_2] for (plot, path, label) in zip(plots, paths, labels) - plot!(plot, 2:T, path, lw=2, label=label, xlabel="Time", grid=true) + plot!(plot, 2:T, path, lw = 2, label = label, xlabel = "Time", + grid = true) end - plot(plt_1, plt_2, layout=(2,1), size = (600,500)) + plot(plt_1, plt_2, layout = (2, 1), size = (600, 500)) end ``` @@ -847,39 +837,23 @@ using Random Random.seed!(42) # parameters -β = 1 / 1.05 -ρ, mg = .7, .35 -A = [ρ mg*(1 - ρ); 0.0 1.0] -C = [sqrt(1 - ρ^2) * mg / 10 0.0; 0 0] +beta = 1 / 1.05 +rho, mg = 0.7, 0.35 +A = [rho mg*(1 - rho); 0.0 1.0] +C = [sqrt(1 - rho^2) * mg/10 0.0; 0 0] Sg = [1.0 0.0] Sd = [0.0 0.0] Sb = [0 2.135] Ss = [0.0 0.0] proc = ContStochProcess(A, C) -econ = Economy(β, Sg, Sd, Sb, Ss, proc) +econ = Economy(beta, Sg, Sd, Sb, Ss, proc) T = 50 path = compute_paths(econ, T) gen_fig_1(path) ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test path.p[3] ≈ 1.5395294981420302 atol = 1e-3 # randomness check. - #test path.g[31] ≈ 0.366487070014081 atol = 1e-3 # stuff we plot - #test path.c[36] ≈ 0.6291101011610297 atol = 1e-3 - #test path.B[9] ≈ 0.07442403655989423 atol = 1e-3 - #test path.rvn[27] ≈ 0.35013269833342753 atol = 1e-3 - #test path.π[31] ≈ -0.05846215388377568 atol = 1e-3 - #test path.R[43] ≈ 1.0437715852385672 atol = 1e-3 - #test path.ξ[43] ≈ 1.001895202392805 atol = 1e-3 - #test path.Π[43] ≈ -0.4185282208457552 atol = 1e-3 # plot tests -end -``` The legends on the figures indicate the variables being tracked. @@ -911,48 +885,31 @@ Random.seed!(42); ```{code-cell} julia # Parameters -β = 1 / 1.05 +beta = 1 / 1.05 P = [0.8 0.2 0.0 - 0.0 0.5 0.5 - 0.0 0.0 1.0] + 0.0 0.5 0.5 + 0.0 0.0 1.0] # Possible states of the world # Each column is a state of the world. The rows are [g d b s 1] x_vals = [0.5 0.5 0.25; - 0.0 0.0 0.0; - 2.2 2.2 2.2; - 0.0 0.0 0.0; - 1.0 1.0 1.0] + 0.0 0.0 0.0; + 2.2 2.2 2.2; + 0.0 0.0 0.0; + 1.0 1.0 1.0] Sg = [1.0 0.0 0.0 0.0 0.0] Sd = [0.0 1.0 0.0 0.0 0.0] Sb = [0.0 0.0 1.0 0.0 0.0] Ss = [0.0 0.0 0.0 1.0 0.0] proc = DiscreteStochProcess(P, x_vals) -econ = Economy(β, Sg, Sd, Sb, Ss, proc) +econ = Economy(beta, Sg, Sd, Sb, Ss, proc) T = 15 path = compute_paths(econ, T) gen_fig_1(path) ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test path.p[3] ≈ 1.5852129146694405 - #test path.B[13] ≈ 0.003279632025474284 - #@test path.g ≈ [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, - # 0.25, 0.25] - #test path.rvn[7] ≈ 0.3188722725349599 - #test path.c[2] ≈ 0.6147870853305598 - #@test path.R ≈ [1.05, 1.05, 1.05, 1.05, 1.05, 1.0930974212983846, 1.05, 1.05, 1.05, 1.05, - # 1.05, 1.05, 1.05, 1.05, 1.05] - #@test path.ξ ≈ [1.0, 1.0, 1.0, 1.0, 1.0, 0.9589548368586813, 1.0, 1.0, 1.0, 1.0, 1.0, - # 1.0, 1.0, 1.0] -end -``` The call `gen_fig_2(path)` generates @@ -997,41 +954,26 @@ Random.seed!(42); tags: [hide-output] --- # parameters -β = 1 / 1.05 -ρ, mg = .95, .35 -A = [0. 0. 0. ρ mg*(1-ρ); +beta = 1 / 1.05 +rho, mg = .95, .35 +A = [0. 0. 0. rho mg*(1-rho); 1. 0. 0. 0. 0.; 0. 1. 0. 0. 0.; 0. 0. 1. 0. 0.; 0. 0. 0. 0. 1.] C = zeros(5, 5) -C[1, 1] = sqrt(1 - ρ^2) * mg / 8 +C[1, 1] = sqrt(1 - rho^2) * mg / 8 Sg = [1. 0. 0. 0. 0.] Sd = [0. 0. 0. 0. 0.] Sb = [0. 0. 0. 0. 2.135] Ss = [0. 0. 0. 0. 0.] proc = ContStochProcess(A, C) -econ = Economy(β, Sg, Sd, Sb, Ss, proc) +econ = Economy(beta, Sg, Sd, Sb, Ss, proc) T = 50 path = compute_paths(econ, T) ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test path.p[3] ≈ 1.524261187305079 - #test path.B[13] ≈ -0.053219518947408805 - #test path.g[7] ≈ 0.36908804521710115 - #test path.rvn[7] ≈ 0.35146870025913474 - #test path.c[2] ≈ 0.6259521929536346 - #test path.R[5][1] ≈ 1.0501742289013196 - #test path.ξ[10] ≈ 1.002202281639002 -end -``` - ```{code-cell} julia gen_fig_1(path) ``` diff --git a/lectures/dynamic_programming_squared/opt_tax_recur.md b/lectures/dynamic_programming_squared/opt_tax_recur.md index 93571209..b58ab291 100644 --- a/lectures/dynamic_programming_squared/opt_tax_recur.md +++ b/lectures/dynamic_programming_squared/opt_tax_recur.md @@ -726,10 +726,10 @@ import QuantEcon: simulate mutable struct Model{TF <: AbstractFloat, TM <: AbstractMatrix{TF}, TV <: AbstractVector{TF}} - β::TF - Π::TM + beta::TF + Pi::TM G::TV - Θ::TV + Theta::TV transfers::Bool U::Function Uc::Function @@ -747,31 +747,31 @@ import QuantEcon: simulate S::TI cFB::TV nFB::TV - ΞFB::TV + XiFB::TV zFB::TV end function SequentialAllocation(model) - β, Π, G, Θ = model.β, model.Π, model.G, model.Θ - mc = MarkovChain(Π) - S = size(Π, 1) # Number of states + beta, Pi, G, Theta = model.beta, model.Pi, model.G, model.Theta + mc = MarkovChain(Pi) + S = size(Pi, 1) # Number of states # now find the first best allocation - cFB, nFB, ΞFB, zFB = find_first_best(model, S, 1) + cFB, nFB, XiFB, zFB = find_first_best(model, S, 1) - return SequentialAllocation(model, mc, S, cFB, nFB, ΞFB, zFB) + return SequentialAllocation(model, mc, S, cFB, nFB, XiFB, zFB) end function find_first_best(model, S, version) if version != 1 && version != 2 throw(ArgumentError("version must be 1 or 2")) end - β, Θ, Uc, Un, G, Π = - model.β, model.Θ, model.Uc, model.Un, model.G, model.Π + beta, Theta, Uc, Un, G, Pi = + model.beta, model.Theta, model.Uc, model.Un, model.G, model.Pi function res!(out, z) c = z[1:S] n = z[S+1:end] - out[1:S] = Θ .* Uc(c, n) + Un(c, n) - out[S+1:end] = Θ .* n - c - G + out[1:S] = Theta .* Uc(c, n) + Un(c, n) + out[S+1:end] = Theta .* n - c - G end res = nlsolve(res!, 0.5 * ones(2 * S)) @@ -782,31 +782,31 @@ function find_first_best(model, S, version) if version == 1 cFB = res.zero[1:S] nFB = res.zero[S+1:end] - ΞFB = Uc(cFB, nFB) # Multiplier on the resource constraint - zFB = vcat(cFB, nFB, ΞFB) - return cFB, nFB, ΞFB, zFB + XiFB = Uc(cFB, nFB) # Multiplier on the resource constraint + zFB = vcat(cFB, nFB, XiFB) + return cFB, nFB, XiFB, zFB elseif version == 2 cFB = res.zero[1:S] nFB = res.zero[S+1:end] IFB = Uc(cFB, nFB) .* cFB + Un(cFB, nFB) .* nFB - xFB = \(I - β * Π, IFB) + xFB = \(I - beta * Pi, IFB) zFB = [vcat(cFB[s], xFB[s], xFB) for s in 1:S] return cFB, nFB, IFB, xFB, zFB end end -function time1_allocation(pas::SequentialAllocation, μ) +function time1_allocation(pas::SequentialAllocation, mu) model, S = pas.model, pas.S - Θ, β, Π, G, Uc, Ucc, Un, Unn = - model.Θ, model.β, model.Π, model.G, + Theta, beta, Pi, G, Uc, Ucc, Un, Unn = + model.Theta, model.beta, model.Pi, model.G, model.Uc, model.Ucc, model.Un, model.Unn function FOC!(out, z) c = z[1:S] n = z[S+1:2S] - Ξ = z[2S+1:end] - out[1:S] = Uc(c, n) .- μ * (Ucc(c, n) .* c .+ Uc(c, n)) .- Ξ # FOC c - out[S+1:2S] = Un(c, n) .- μ * (Unn(c, n) .* n .+ Un(c, n)) + Θ .* Ξ # FOC n - out[2S+1:end] = Θ .* n - c - G # Resource constraint + Xi = z[2S+1:end] + out[1:S] = Uc(c, n) .- mu * (Ucc(c, n) .* c .+ Uc(c, n)) .- Xi # FOC c + out[S+1:2S] = Un(c, n) .- mu * (Unn(c, n) .* n .+ Un(c, n)) + Theta .* Xi # FOC n + out[2S+1:end] = Theta .* n - c - G # Resource constraint return out end # Find the root of the FOC @@ -815,56 +815,56 @@ function time1_allocation(pas::SequentialAllocation, μ) error("Could not find LS allocation.") end z = res.zero - c, n, Ξ = z[1:S], z[S+1:2S], z[2S+1:end] + c, n, Xi = z[1:S], z[S+1:2S], z[2S+1:end] # Now compute x Inv = Uc(c, n) .* c + Un(c, n) .* n - x = \(I - β * model.Π, Inv) - return c, n, x, Ξ + x = \(I - beta * model.Pi, Inv) + return c, n, x, Xi end function time0_allocation(pas::SequentialAllocation, B_, s_0) model = pas.model - Π, Θ, G, β = model.Π, model.Θ, model.G, model.β + Pi, Theta, G, beta = model.Pi, model.Theta, model.G, model.beta Uc, Ucc, Un, Unn = model.Uc, model.Ucc, model.Un, model.Unn # First order conditions of planner's problem function FOC!(out, z) - μ, c, n, Ξ = z[1], z[2], z[3], z[4] - xprime = time1_allocation(pas, μ)[3] + mu, c, n, Xi = z[1], z[2], z[3], z[4] + xprime = time1_allocation(pas, mu)[3] out .= vcat( - Uc(c, n) .* (c - B_) .+ Un(c, n) .* n + β * dot(Π[s_0, :], xprime), - Uc(c, n) .- μ * (Ucc(c, n) .* (c - B_) + Uc(c, n)) - Ξ, - Un(c, n) .- μ * (Unn(c, n) .* n .+ Un(c, n)) + Θ[s_0] .* Ξ, - (Θ .* n .- c - G)[s_0] + Uc(c, n) .* (c - B_) .+ Un(c, n) .* n + beta * dot(Pi[s_0, :], xprime), + Uc(c, n) .- mu * (Ucc(c, n) .* (c - B_) + Uc(c, n)) - Xi, + Un(c, n) .- mu * (Unn(c, n) .* n .+ Un(c, n)) + Theta[s_0] .* Xi, + (Theta .* n .- c - G)[s_0] ) end # Find root - res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.ΞFB[s_0]]) + res = nlsolve(FOC!, [0.0, pas.cFB[s_0], pas.nFB[s_0], pas.XiFB[s_0]]) if res.f_converged == false error("Could not find time 0 LS allocation.") end return (res.zero...,) end -function time1_value(pas::SequentialAllocation, μ) +function time1_value(pas::SequentialAllocation, mu) model = pas.model - c, n, x, Ξ = time1_allocation(pas, μ) + c, n, x, Xi = time1_allocation(pas, mu) U_val = model.U.(c, n) - V = \(I - model.β*model.Π, U_val) + V = \(I - model.beta*model.Pi, U_val) return c, n, x, V end function Τ(model, c, n) Uc, Un = model.Uc.(c, n), model.Un.(c, n) - return 1 .+ Un ./ (model.Θ .* Uc) + return 1 .+ Un ./ (model.Theta .* Uc) end function simulate(pas::SequentialAllocation, B_, s_0, T, sHist = nothing) model = pas.model - Π, β, Uc = model.Π, model.β, model.Uc + Pi, beta, Uc = model.Pi, model.beta, model.Uc if isnothing(sHist) sHist = QuantEcon.simulate(pas.mc, T, init=s_0) @@ -873,25 +873,25 @@ function simulate(pas::SequentialAllocation, B_, s_0, T, sHist = nothing) nHist = similar(cHist) Bhist = similar(cHist) ΤHist = similar(cHist) - μHist = similar(cHist) + muHist = similar(cHist) RHist = zeros(T-1) # time 0 - μ, cHist[1], nHist[1], _ = time0_allocation(pas, B_, s_0) + mu, cHist[1], nHist[1], _ = time0_allocation(pas, B_, s_0) ΤHist[1] = Τ(pas.model, cHist[1], nHist[1])[s_0] Bhist[1] = B_ - μHist[1] = μ + muHist[1] = mu # time 1 onward for t in 2:T - c, n, x, Ξ = time1_allocation(pas,μ) + c, n, x, Xi = time1_allocation(pas,mu) u_c = Uc(c,n) s = sHist[t] ΤHist[t] = Τ(pas.model, c, n)[s] - Eu_c = dot(Π[sHist[t-1],:], u_c) + Eu_c = dot(Pi[sHist[t-1],:], u_c) cHist[t], nHist[t], Bhist[t] = c[s], n[s], x[s] / u_c[s] - RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (β * Eu_c) - μHist[t] = μ + RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (beta * Eu_c) + muHist[t] = mu end - return cHist, nHist, Bhist, ΤHist, sHist, μHist, RHist + return cHist, nHist, Bhist, ΤHist, sHist, muHist, RHist end mutable struct BellmanEquation{TP <: Model, @@ -911,7 +911,7 @@ end end function BellmanEquation(model, xgrid, policies0) - S = size(model.Π, 1) # Number of states + S = size(model.Pi, 1) # Number of states xbar = collect(extrema(xgrid)) time_0 = false cf, nf, xprimef = policies0 @@ -924,19 +924,19 @@ end function get_policies_time1(T, i_x, x, s, Vf) model, S = T.model, T.S - β, Θ, G, Π = model.β, model.Θ, model.G, model.Π + beta, Theta, G, Pi = model.beta, model.Theta, model.G, model.Pi U, Uc, Un = model.U, model.Uc, model.Un function objf(z, grad) c, xprime = z[1], z[2:end] n = c + G[s] Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S] - return -(U(c, n) + β * dot(Π[s, :], Vprime)) + return -(U(c, n) + beta * dot(Pi[s, :], Vprime)) end function cons(z, grad) c, xprime = z[1], z[2:end] n = c+G[s] - return x - Uc(c, n) * c - Un(c, n) * n - β * dot(Π[s, :], xprime) + return x - Uc(c, n) * c - Un(c, n) * n - beta * dot(Pi[s, :], xprime) end lb = vcat(0, T.xbar[1] * ones(S)) ub = vcat(1 - G[s], T.xbar[2] * ones(S)) @@ -962,18 +962,18 @@ end function get_policies_time0(T, B_, s0, Vf) model, S = T.model, T.S - β, Θ, G, Π = model.β, model.Θ, model.G, model.Π + beta, Theta, G, Pi = model.beta, model.Theta, model.G, model.Pi U, Uc, Un = model.U, model.Uc, model.Un function objf(z, grad) c, xprime = z[1], z[2:end] n = c + G[s0] Vprime = [Vf[sprime](xprime[sprime]) for sprime in 1:S] - return -(U(c, n) + β * dot(Π[s0, :], Vprime)) + return -(U(c, n) + beta * dot(Pi[s0, :], Vprime)) end function cons(z, grad) c, xprime = z[1], z[2:end] n = c + G[s0] - return -Uc(c, n) * (c - B_) - Un(c, n) * n - β * dot(Π[s0, :], xprime) + return -Uc(c, n) * (c - B_) - Un(c, n) * n - beta * dot(Pi[s0, :], xprime) end lb = vcat(0, T.xbar[1] * ones(S)) ub = vcat(1-G[s0], T.xbar[2] * ones(S)) @@ -1295,41 +1295,41 @@ The above steps are implemented in a type called RecursiveAllocation --- tags: [output_scroll] --- - struct RecursiveAllocation{TP <: Model, TI <: Integer, +struct RecursiveAllocation{TP <: Model, TI <: Integer, TVg <: AbstractVector, TVv <: AbstractVector, TVp <: AbstractArray} model::TP mc::MarkovChain S::TI T::BellmanEquation - μgrid::TVg + mugrid::TVg xgrid::TVg Vf::TVv policies::TVp end -function RecursiveAllocation(model, μgrid) - mc = MarkovChain(model.Π) +function RecursiveAllocation(model, mugrid) + mc = MarkovChain(model.Pi) G = model.G - S = size(model.Π, 1) # Number of states + S = size(model.Pi, 1) # Number of states # Now find the first best allocation - Vf, policies, T, xgrid = solve_time1_bellman(model, μgrid) + Vf, policies, T, xgrid = solve_time1_bellman(model, mugrid) T.time_0 = true # Bellman equation now solves time 0 problem - return RecursiveAllocation(model, mc, S, T, μgrid, xgrid, Vf, policies) + return RecursiveAllocation(model, mc, S, T, mugrid, xgrid, Vf, policies) end -function solve_time1_bellman(model, μgrid) - μgrid0 = μgrid - S = size(model.Π, 1) +function solve_time1_bellman(model, mugrid) + mugrid0 = mugrid + S = size(model.Pi, 1) # First get initial fit PP = SequentialAllocation(model) - c = zeros(length(μgrid), 2) + c = zeros(length(mugrid), 2) n = similar(c) x = similar(c) V = similar(c) - for (i, μ) in enumerate(μgrid0) - c[i, :], n[i, :], x[i, :], V[i, :] = time1_value(PP, μ) + for (i, mu) in enumerate(mugrid0) + c[i, :], n[i, :], x[i, :], V[i, :] = time1_value(PP, mu) end Vf = Vector{AbstractInterpolation}(undef, 2) cf = similar(Vf) @@ -1346,7 +1346,7 @@ function solve_time1_bellman(model, μgrid) policies = [cf, nf, xprimef] # Create xgrid xbar = [maximum(minimum(x, dims = 1)), minimum(maximum(x, dims = 1))] - xgrid = range(xbar[1], xbar[2], length = length(μgrid0)) + xgrid = range(xbar[1], xbar[2], length = length(mugrid0)) # Now iterate on bellman equation T = BellmanEquation(model, xgrid, policies) diff = 1.0 @@ -1410,19 +1410,19 @@ end function simulate(pab::RecursiveAllocation, B_, s_0, T, sHist = QuantEcon.simulate(mc, s_0, T)) model, S, policies = pab.model, pab.S, pab.policies - β, Π, Uc = model.β, model.Π, model.Uc + beta, Pi, Uc = model.beta, model.Pi, model.Uc cf, nf, xprimef = policies[1], policies[2], policies[3] cHist = zeros(T) nHist = similar(cHist) Bhist = similar(cHist) ΤHist = similar(cHist) - μHist = similar(cHist) + muHist = similar(cHist) RHist = zeros(T - 1) # time 0 cHist[1], nHist[1], xprime = time0_allocation(pab, B_, s_0) ΤHist[1] = Τ(pab.model, cHist[1], nHist[1])[s_0] Bhist[1] = B_ - μHist[1] = 0.0 + muHist[1] = 0.0 # time 1 onward for t in 2:T s, x = sHist[t], xprime[sHist[t]] @@ -1431,12 +1431,12 @@ function simulate(pab::RecursiveAllocation, B_, s_0, T, xprime = [xprimef[s, sprime](x) for sprime in 1:S] ΤHist[t] = Τ(pab.model, c, n)[s] u_c = Uc(c, n) - Eu_c = dot(Π[sHist[t-1], :], u_c) - μHist[t] = pab.Vf[s](x) - RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (β * Eu_c) + Eu_c = dot(Pi[sHist[t-1], :], u_c) + muHist[t] = pab.Vf[s](x) + RHist[t-1] = Uc(cHist[t-1], nHist[t-1]) / (beta * Eu_c) cHist[t], nHist[t], Bhist[t] = c[s], n, x / u_c[s] end - return cHist, nHist, Bhist, ΤHist, sHist, μHist, RHist + return cHist, nHist, Bhist, ΤHist, sHist, muHist, RHist end ``` @@ -1490,30 +1490,29 @@ This utility function is implemented in the type CRRAutility ```{code-cell} julia function crra_utility(; - β = 0.9, - σ = 2.0, - γ = 2.0, - Π = 0.5 * ones(2, 2), - G = [0.1, 0.2], - Θ = ones(2), - transfers = false - ) + beta = 0.9, + sigma = 2.0, + gamma = 2.0, + Pi = 0.5 * ones(2, 2), + G = [0.1, 0.2], + Theta = ones(2), + transfers = false) function U(c, n) - if σ == 1.0 + if sigma == 1.0 U = log(c) else - U = (c.^(1.0 .- σ) .- 1.0) / (1.0 - σ) + U = (c .^ (1.0 .- sigma) .- 1.0) / (1.0 - sigma) end - return U .- n.^(1 + γ) / (1 + γ) + return U .- n .^ (1 + gamma) / (1 + gamma) end # Derivatives of utility function - Uc(c,n) = c.^(-σ) - Ucc(c,n) = -σ * c.^(-σ - 1.0) - Un(c,n) = -n.^γ - Unn(c,n) = -γ * n.^(γ - 1.0) + Uc(c, n) = c .^ (-sigma) + Ucc(c, n) = -sigma * c .^ (-sigma - 1.0) + Un(c, n) = -n .^ gamma + Unn(c, n) = -gamma * n .^ (gamma - 1.0) n_less_than_one = false - return Model(β, Π, G, Θ, transfers, - U, Uc, Ucc, Un, Unn, n_less_than_one) + return Model(beta, Pi, G, Theta, transfers, + U, Uc, Ucc, Un, Unn, n_less_than_one) end ``` @@ -1528,15 +1527,15 @@ We can now plot the Ramsey tax under both realizations of time $t = 3$ governme using Random Random.seed!(42) # For reproducible results. -M_time_example = crra_utility(G=[0.1, 0.1, 0.1, 0.2, 0.1, 0.1], - Θ=ones(6)) # Θ can in principle be random +M_time_example = crra_utility(G = [0.1, 0.1, 0.1, 0.2, 0.1, 0.1], + Theta = ones(6)) # Theta can in principle be random -M_time_example.Π = [0.0 1.0 0.0 0.0 0.0 0.0; - 0.0 0.0 1.0 0.0 0.0 0.0; - 0.0 0.0 0.0 0.5 0.5 0.0; - 0.0 0.0 0.0 0.0 0.0 1.0; - 0.0 0.0 0.0 0.0 0.0 1.0; - 0.0 0.0 0.0 0.0 0.0 1.0] +M_time_example.Pi = [0.0 1.0 0.0 0.0 0.0 0.0; + 0.0 0.0 1.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.5 0.5 0.0; + 0.0 0.0 0.0 0.0 0.0 1.0; + 0.0 0.0 0.0 0.0 0.0 1.0; + 0.0 0.0 0.0 0.0 0.0 1.0] PP_seq_time = SequentialAllocation(M_time_example) # Solve sequential problem @@ -1556,38 +1555,23 @@ titles = hcat("Consumption", "Output") sim_seq_l_plot = [sim_seq_l[1:4]..., M_time_example.G[sHist_l], - M_time_example.Θ[sHist_l].*sim_seq_l[2]] + M_time_example.Theta[sHist_l] .* sim_seq_l[2]] sim_seq_h_plot = [sim_seq_h[1:4]..., M_time_example.G[sHist_h], - M_time_example.Θ[sHist_h].*sim_seq_h[2]] - + M_time_example.Theta[sHist_h] .* sim_seq_h[2]] #plots = plot(layout=(3,2), size=(800,600)) plots = [plot(), plot(), plot(), plot(), plot(), plot()] -for i = 1:6 - plot!(plots[i], sim_seq_l_plot[i], color=:black, lw=2, - marker=:circle, markersize=2, label="") - plot!(plots[i], sim_seq_h_plot[i], color=:red, lw=2, - marker=:circle, markersize=2, label="") - plot!(plots[i], title=titles[i], grid=true) +for i in 1:6 + plot!(plots[i], sim_seq_l_plot[i], color = :black, lw = 2, + marker = :circle, markersize = 2, label = "") + plot!(plots[i], sim_seq_h_plot[i], color = :red, lw = 2, + marker = :circle, markersize = 2, label = "") + plot!(plots[i], title = titles[i], grid = true) end -plot(plots[1], plots[2], plots[3], plots[4], plots[5], plots[6], layout=(3,2), size=(800,600)) +plot(plots[1], plots[2], plots[3], plots[4], plots[5], plots[6], + layout = (3, 2), size = (800, 600)) ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - @test M_time_example.G[sHist_l] ≈ [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] - #test M_time_example.Θ[sHist_l] .* sim_seq_l[2] ≈ [1.026385289423105, 0.9945696863679917, - # 0.9945696863679917, 0.9945696863679917, - # 0.9945696863679917, 0.9945696863679917, - # 0.9945696863679917] - @test M_time_example.G[sHist_h] ≈ [0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1] - #test sim_seq_l[end] ≈ [1.0361020796451619, 1.111111111111111, 1.052459380877434, - # 1.111111111111111, 1.111111111111111, 1.111111111111111] -end -``` **Tax smoothing** @@ -1628,11 +1612,11 @@ The following plot illustrates how the government lowers the interest rate at time 0 by raising consumption ```{code-cell} julia -plot(sim_seq_l[end], color=:black, lw=2, - marker=:circle, markersize=2, label="") -plot!(sim_seq_h[end], color=:red, lw=2, - marker=:circle, markersize=2, label="") -plot!(title="Gross Interest Rate", grid=true) +plot(sim_seq_l[end], color = :black, lw = 2, + marker = :circle, markersize = 2, label = "") +plot!(sim_seq_h[end], color = :red, lw = 2, + marker = :circle, markersize = 2, label = "") +plot!(title = "Gross Interest Rate", grid = true) ``` ### Government Saving @@ -1701,35 +1685,25 @@ Random.seed!(42); # For reproducible results. ``` ```{code-cell} julia -M2 = crra_utility(G=[0.15], Π=ones(1, 1), Θ=[1.0]) +M2 = crra_utility(G = [0.15], Pi = ones(1, 1), Theta = [1.0]) PP_seq_time0 = SequentialAllocation(M2) # solve sequential problem B_vec = range(-1.5, 1.0, length = 100) taxpolicy = Matrix(hcat([simulate(PP_seq_time0, B_, 1, 2)[4] for B_ in B_vec]...)') -interest_rate = Matrix(hcat([simulate(PP_seq_time0, B_, 1, 3)[end] for B_ in B_vec]...)') +interest_rate = Matrix(hcat([simulate(PP_seq_time0, B_, 1, 3)[end] + for B_ in B_vec]...)') titles = ["Tax Rate" "Gross Interest Rate"] labels = [[L"Time , $t = 0$" L"Time , $t \geq 0$"], ""] -plots = plot(layout=(2,1), size =(700,600)) +plots = plot(layout = (2, 1), size = (700, 600)) for (i, series) in enumerate((taxpolicy, interest_rate)) - plot!(plots[i], B_vec, series, linewidth=2, label=labels[i]) - plot!(plots[i], title=titles[i], grid=true, legend=:topleft) + plot!(plots[i], B_vec, series, linewidth = 2, label = labels[i]) + plot!(plots[i], title = titles[i], grid = true, legend = :topleft) end plot(plots) ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test B_vec[3] ≈ -1.4494949494949494 - #test taxpolicy[2, 2] ≈ 0.0020700125847712414 - #test interest_rate[3, 1] ≈ 1.113064964490116 -end -``` - The figure indicates that if the government enters with positive debt, it sets a tax rate at $t=0$ that is less than all later tax rates. @@ -1781,9 +1755,10 @@ B1_vec = hcat([simulate(PP_seq_time0, B_, 1, 2)[3][2] for B_ in B_vec]...)' # Compute the optimal policy if the government could reset tau1_reset = Matrix(hcat([simulate(PP_seq_time0, B1, 1, 1)[4] for B1 in B1_vec]...)') -plot(B_vec, taxpolicy[:, 2], linewidth=2, label=L"\tau_1") -plot!(B_vec, tau1_reset, linewidth=2, label=L"\tau_1^R") -plot!(title="Tax Rate", xlabel="Initial Government Debt", legend=:topleft, grid=true) +plot(B_vec, taxpolicy[:, 2], linewidth = 2, label = L"\tau_1") +plot!(B_vec, tau1_reset, linewidth = 2, label = L"\tau_1^R") +plot!(title = "Tax Rate", xlabel = "Initial Government Debt", legend = :topleft, + grid = true) ``` The tax rates in the figure are equal for only two values of initial government debt. @@ -1819,21 +1794,21 @@ $$ We will write a new constructor LogUtility to represent this utility function ```{code-cell} julia -function log_utility(;β = 0.9, - ψ = 0.69, - Π = 0.5 * ones(2, 2), - G = [0.1, 0.2], - Θ = ones(2), - transfers = false) +function log_utility(; beta = 0.9, + psi = 0.69, + Pi = 0.5 * ones(2, 2), + G = [0.1, 0.2], + Theta = ones(2), + transfers = false) # Derivatives of utility function - U(c,n) = log(c) + ψ * log(1 - n) - Uc(c,n) = 1 ./ c - Ucc(c,n) = -c.^(-2.0) - Un(c,n) = -ψ ./ (1.0 .- n) - Unn(c,n) = -ψ ./ (1.0 .- n).^2.0 + U(c, n) = log(c) + psi * log(1 - n) + Uc(c, n) = 1 ./ c + Ucc(c, n) = -c .^ (-2.0) + Un(c, n) = -psi ./ (1.0 .- n) + Unn(c, n) = -psi ./ (1.0 .- n) .^ 2.0 n_less_than_one = true - return Model(β, Π, G, Θ, transfers, - U, Uc, Ucc, Un, Unn, n_less_than_one) + return Model(beta, Pi, G, Theta, transfers, + U, Uc, Ucc, Un, Unn, n_less_than_one) end ``` @@ -1853,9 +1828,9 @@ Random.seed!(42); # For reproducible results. ```{code-cell} julia M1 = log_utility() -μ_grid = range(-0.6, 0.0, length = 200) +mu_grid = range(-0.6, 0.0, length = 200) PP_seq = SequentialAllocation(M1) # Solve sequential problem -PP_bel = RecursiveAllocation(M1, μ_grid) # Solve recursive problem +PP_bel = RecursiveAllocation(M1, mu_grid) # Solve recursive problem T = 20 sHist = [1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1] @@ -1865,8 +1840,8 @@ sim_seq = simulate(PP_seq, 0.5, 1, T, sHist) sim_bel = simulate(PP_bel, 0.5, 1, T, sHist) # Plot policies -sim_seq_plot = [sim_seq[1:4]..., M1.G[sHist], M1.Θ[sHist].*sim_seq[2]] -sim_bel_plot = [sim_bel[1:4]..., M1.G[sHist], M1.Θ[sHist].*sim_bel[2]] +sim_seq_plot = [sim_seq[1:4]..., M1.G[sHist], M1.Theta[sHist] .* sim_seq[2]] +sim_bel_plot = [sim_bel[1:4]..., M1.G[sHist], M1.Theta[sHist] .* sim_bel[2]] titles = hcat("Consumption", "Labor Supply", @@ -1874,36 +1849,26 @@ titles = hcat("Consumption", "Tax Rate", "Government Spending", "Output") -labels = [["Sequential", "Recursive"], ["",""], ["",""], ["",""], ["",""], ["",""]] -plots=plot(layout=(3,2), size=(850,780)) - -for i = 1:6 - plot!(plots[i], sim_seq_plot[i], color=:black, lw=2, marker=:circle, - markersize=2, label=labels[i][1]) - plot!(plots[i], sim_bel_plot[i], color=:blue, lw=2, marker=:xcross, - markersize=2, label=labels[i][2]) - plot!(plots[i], title=titles[i], grid=true, legend=:topright) +labels = [ + ["Sequential", "Recursive"], + ["", ""], + ["", ""], + ["", ""], + ["", ""], + ["", ""], +] +plots = plot(layout = (3, 2), size = (850, 780)) + +for i in 1:6 + plot!(plots[i], sim_seq_plot[i], color = :black, lw = 2, marker = :circle, + markersize = 2, label = labels[i][1]) + plot!(plots[i], sim_bel_plot[i], color = :blue, lw = 2, marker = :xcross, + markersize = 2, label = labels[i][2]) + plot!(plots[i], title = titles[i], grid = true, legend = :topright) end plot(plots) ``` -```{code-cell} julia ---- -tags: [remove-cell] ---- -@testset begin - #test sim_seq_plot[1][14] ≈ 0.38396935397869975 - #test sim_seq_plot[2][14] ≈ 0.5839693539786998 - #test sim_seq_plot[3][14] ≈ 0.3951985593686047 - #test sim_seq_plot[4][14] ≈ 0.3631746680706347 - #test sim_seq_plot[5][14] ≈ 0.2 - #test sim_seq_plot[6][14] ≈ 0.5839693539786998 - #test sim_bel_plot[3][5] ≈ 0.5230509296608254 atol = 1e-3 - #test sim_bel_plot[5][7] ≈ 0.1 - #test sim_bel_plot[2][3] ≈ 0.5402933557593538 atol = 1e-3 -end -``` - As should be expected, the recursive and sequential solutions produce almost identical allocations.