diff --git a/Project.toml b/Project.toml index 80fe228..7731225 100644 --- a/Project.toml +++ b/Project.toml @@ -24,7 +24,8 @@ julia = "1.2" [extras] Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" InfinitesimalGenerators = "2fce0c6f-5f0b-5c85-85c9-2ffe1d5ee30d" +Roots ="f2b01f46-fcfa-551c-844a-d8ac1e96c665" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Distributions", "InfinitesimalGenerators"] +test = ["Test", "Distributions", "InfinitesimalGenerators", "Roots"] diff --git a/examples/AssetPricing/BrunnermeirSannikov.jl b/examples/AssetPricing/BrunnermeirSannikov.jl index 001ed48..9dc99cc 100644 --- a/examples/AssetPricing/BrunnermeirSannikov.jl +++ b/examples/AssetPricing/BrunnermeirSannikov.jl @@ -1,120 +1,128 @@ -# Brunnermeier Sannikov "A Macroeconomic Model with a Financial Sector" (2014) -# The idfference with Ditella is that Ψ appears in the last algebratic equation, making it much much harder to solve for it -# using non linear solver -using EconPDEs +using EconPDEs, Roots + +# There are households and experts, with different productivity +# Experts can issue equity but must retain at least a fraction χlow of their equity. This is similar to HK: households can only invest up to w^E * m in intermediaries; with m = 1 / χlow - 1 +# Ψ is the fraction of capital held by experts. Compared to HK, can be lower than one +# Moreover, experts can issue equity at lower required returns than the return on capital (ie. they earn management fee) Base.@kwdef mutable struct BrunnermeierSannikov - # See Table H Numerical examples + # Calibration uses Section 3.6 of the textbook chapter (I think r is typo for ρH) ρE::Float64 = 0.06 - ρH::Float64 = 0.06 + ρH::Float64 = 0.05 γ::Float64 = 2.0 ψ::Float64 = 0.5 δ::Float64 = 0.05 σ::Float64 = 0.1 κ::Float64 = 10.0 a::Float64 = 0.11 - alow::Float64 = 0.1 - χlow::Float64 = 1.0 - firstregion::Bool = true - q::Float64 = 1.0 - qx::Float64 = 0.0 + alow::Float64 = 0.03 + χlow::Float64 = 0.5 + q_old::Float64 = 1.0 + qx_old::Float64 = 0.0 + Ψ_old::Float64 = 1.0 end function (m::BrunnermeierSannikov)(state::NamedTuple, y::NamedTuple, Δx, xmin) (; ρE, ρH, γ, ψ, δ, σ, κ, a, alow, χlow) = m (; x) = state (; pE, pEx_up, pEx_down, pExx, pH, pHx_up, pHx_down, pHxx) = y - # χ is the fraction of risk held by experts. - # ψ is the fraction of (real) capital held by experts - # so χ * ψ / x corresponds to total risk held by experts / their networth - # σE = κE / γ + (1 - γ) / (γ * (ψ - 1)) * σpE - # implies κE = γ * σE - (1 - γ) / ((ψ - 1)) * σpE - # so κE - κH = γ * (σE -σH) - (1 - γ) / ((ψ - 1)) * (σpE -σpH) - # = γ * (χ * ψ / x - (1-χ * ψ) / (1-x) - (1 - γ) / (γ * (ψ - 1)) * (σpE - σpH) - # . Specifically, we suppose that experts must retain at least a fraction χbar ∈ (0, 1] of equity. (i - # drift and volatility of state variable ν - #@show x, p, ψ pE = max(pE, 1e-3) pH = max(pH, 1e-3) - - - χ = max(χlow, x) - # market clear givges - # x / pE + (1-x) / pH = (ψ * a + (1-ψ) * alow - i) / q - # this implies ψ pEx, pHx = pEx_up, pHx_up iter = 0 @label start - q = m.q - if x ≈ xmin - m.firstregion = true - # poin toutside the gri, so zero - Ψ = 0.0 - q = (Ψ * a + (1 - Ψ) * alow + 1 / κ) / ((x / pE + (1 - x) / pH) + 1 / κ) + q_old = m.q_old + qx_old = m.qx_old + Ψ_old = m.Ψ_old + + # First, solve for q, qx, and qxx + if abs(x - xmin) < 1e-9 + # initial condition at lower boundery of x grid (point outside the grid technically) + # We must have that expert consumption approaches zero as xt approaches zero (linked to transversality condition) + Ψ_old = 0.0 + q_old = (alow + 1 / κ) / (1 / pH + 1 / κ) + qx_old = 0.0 end - if m.firstregion - # in the second region (a - alow) / q = (κE - κH) * (σ + σq) - # here q corresponds to value in the previous x grif - i = (q - 1) / κ - Ψ = max(((x / pE + (1 - x) / pH) * q + i - alow) / (a - alow), (1e-10 + x) / χ) - K = q / (a - alow) * χ * σ^2 * (γ / (x * (1 - x)) - (1 - γ) / (ψ - 1) * (pEx / pE - pHx / pH)) - X =χ * Ψ - x - qx = (1 - sqrt(max(eps(), X * K))) / X * q - qx = max(min(qx, 4.0), 0.0) - q = q + qx * Δx - #@assert abs((1 - qx / q * X)^2 - K * X) < 1e-8 - if Ψ > 1.0 - m.firstregion = false - Ψ = 1 - q = (a + 1 / κ) / ((x / pE + (1 - x) / pH) + 1 / κ) - qx = - (a + 1 / κ) / ((x / pE + (1 - x) / pH) + 1 / κ)^2 * (1 / pE - 1 / pH - x / pE^2 * pEx - (1-x) / pH^2 * pHx) + if Ψ_old < 1.0 + # crisis region. + # in this case, we do not have Ψ = 1.0 but we have + # (a - alow) / q = (κE - κH) * (σ + σq) + # note that + out = find_zero(Ψ_old) do Ψ + local q = (Ψ * a + (1 - Ψ) * alow + 1 / κ) / (x / pE + (1 - x) / pH + 1 / κ) + local qx = (q - q_old) / Δx + local αE = χ * Ψ / x + local σx = x * (αE - 1) * σ / (1 - x * (αE - 1) * qx / q) + local σpE = pEx / pE * σx + local σpH = pHx / pH * σx + local σq = qx / q * σx + local σE = αE * (σ + σq) + local αH = (1 - αE * x) / (1 - x) + local σH = αH * (σ + σq) + local κE = γ * (σE - (1 / γ - 1) / (ψ - 1) * σpE) + local κH = γ * (σH - (1 / γ - 1) / (ψ - 1) * σpH) + return (a - alow) / q - χ * (κE - κH) * (σ + σq) + end + Ψ = out[1] + q = (Ψ * a + (1 - Ψ) * alow + 1 / κ) / (x / pE + (1 - x) / pH + 1 / κ) + qx = (q - q_old) / Δx + if Ψ > 1 + Ψ_old = 1.0 end - else - # in the second region Ψ = 1 - Ψ = 1 - # We get q through market clearing equation after plugging i = (p - 1) / κ + end + if Ψ_old ≥ 1.0 + # normal region. + Ψ = 1.0 q = (a + 1 / κ) / ((x / pE + (1 - x) / pH) + 1 / κ) qx = - (a + 1 / κ) / ((x / pE + (1 - x) / pH) + 1 / κ)^2 * (1 / pE - 1 / pH - x / pE^2 * pEx - (1-x) / pH^2 * pHx) - i = (q - 1) / κ - @assert abs((a - i) / q - (x / pE + (1 - x) / pH)) < 1e-8 end - #@show m.firstregion, x, qx - #@assert x isa Float64 - qxx = (qx - m.qx) / Δx - m.q = q - m.qx = qx - q = max(1e-3, q) - i = (q - 1) / κ + qxx = (qx - qx_old) / Δx + m.q_old = q + m.qx_old = qx + m.Ψ_old = Ψ + + + # having solved for q, do the time step + q = max(q, sqrt(eps())) + ι = (q - 1) / κ Φ = log(q) / κ - σx = (χ * Ψ - x) * σ / (1 - (χ * Ψ - x) * qx / q) + αE = χ * Ψ / x + σx = x * (αE - 1) * σ / (1 - x * (αE - 1) * qx / q) σpE = pEx / pE * σx σpH = pHx / pH * σx σq = qx / q * σx - σE = Ψ * χ / x * (σ + σq) - σH = (1 - Ψ * χ) / (1 - x) * (σ + σq) - κE = γ * σE - (1 - γ) / (ψ - 1) * σpE - κH = γ * σH - (1 - γ) / (ψ - 1) * σpH - μx = x * ((a - i) / q - 1 / pE) + σx * (κE - σ - σq) + x * (σ + σq) * (1 - χ) * (κE - κH) - #μx = x * (1 - x) * (κE * σE - κH * σH + 1 / pH - 1 / pE - (σE - σH) * (σ + σq)) + σE = αE * (σ + σq) + αH = (1 - αE * x) / (1 - x) + σH = αH * (σ + σq) + κE = γ * (σE - (1 / γ - 1) / (ψ - 1) * σpE) + κH = γ * (σH - (1 / γ - 1) / (ψ - 1) * σpH) + μx = x * (1 - x) * (κE * σE - κH * σH + 1 / pH - 1 / pE - (σE - σH) * (σ + σq)) + μx2 = x * ((a - ι) / q - 1 / pE) + σx * (κE - σ - σq) + x * (1 - χ) * (κE - κH) * (σ + σq) if (iter == 0) && (μx < 0) + # upwinding iter += 1 pEx, pHx = pEx_down, pHx_down @goto start end + #@assert abs((a - ι) / q - (x / pE + (1 - x) / pH)) < 1e-6 μpE = pEx / pE * μx + 0.5 * pExx / pE * σx^2 μpH = pHx / pH * μx + 0.5 * pHxx / pH * σx^2 μq = qx / q * μx + 0.5 * qxx / q * σx^2 - r = min(max((a - i) / q + Φ - δ + μq + σ * σq - (χ * κE + (1 - χ) * κH) * (σ + σq), -0.1), 0.1) - #if ψ < 1 - # @assert r ≈ (alow - i) / p + Φ - δ + μp + σ * σp - κH * (σ + σp) - #end + r = clamp((a - ι) / q + Φ - δ + μq + σ * σq - (χ * κE + (1 - χ) * κH) * (σ + σq), -3.0, 3.0) # Market Pricing pEt = - pE * (1 / pE - ψ * ρE + (ψ - 1) * (r + κE * σE) + μpE - (ψ - 1) * γ / 2 * σE^2 + (2 - ψ - γ) / (2 * (ψ - 1)) * σpE^2 + (1 - γ) * σpE * σE) pHt = - pH * (1 / pH - ψ * ρH + (ψ - 1) * (r + κH * σH) + μpH - (ψ - 1) * γ / 2 * σH^2 + (2 - ψ - γ) / (2 * (ψ - 1)) * σpH^2 + (1 - γ) * σpH * σH) - (; pEt, pHt), (; x, r, μx, σx, Ψ, κ, κE, κH, σE, σH, σq, pEt, pHt, q, qx, qxx, μq) + d = Ψ * a + (1 - Ψ) * alow - ι + μR = r + (χ * κE + (1 - χ) * κH) * (σ + σq) + (; pEt, pHt), (; x, r, μx, σx, Ψ, κ, κE, κH, σE, σH, σq, pEt, pHt, q, qx, qxx, μq, μx2, Φ, χ, d, μR, pExx, pHxx) end - - +m = BrunnermeierSannikov() +xn = 200 +stategrid = OrderedDict(:x => range(0, 1.0, length = xn+2)[2:(end-1)]) +Δx = step(stategrid[:x]) +xmin = minimum(stategrid[:x]) +yend = OrderedDict(:pE => 9 .* ones(xn), :pH => 10 .* ones(xn)) +y, residual_norm, a = pdesolve((state, y) -> m(state, y, Δx, xmin), stategrid, yend; autodiff = :finite, maxΔ = 1000.0) \ No newline at end of file diff --git a/examples/AssetPricing/GarleanuPanageas.jl b/examples/AssetPricing/GarleanuPanageas.jl index ca40289..54d5540 100644 --- a/examples/AssetPricing/GarleanuPanageas.jl +++ b/examples/AssetPricing/GarleanuPanageas.jl @@ -80,7 +80,7 @@ function (m::GarleanuPanageasModel)(state::NamedTuple, y::NamedTuple) ϕ1t = - ϕ1 * (1 / ϕ1 + μ - δ - δ1 + μϕ1 + σ * σϕ1 - r - κ * (σϕ1 + σ)) ϕ2t = - ϕ2 * (1 / ϕ2 + μ - δ - δ2 + μϕ2 + σ * σϕ2 - r - κ * (σϕ2 + σ)) - return (; pAt, pBt, ϕ1t, ϕ2t), (;r) + return (; pAt, pBt, ϕ1t, ϕ2t), (;r, κ) end m = GarleanuPanageasModel() diff --git a/examples/AssetPricing/HeKrishnamurthy.jl b/examples/AssetPricing/HeKrishnamurthy.jl index 487805b..40c8a04 100644 --- a/examples/AssetPricing/HeKrishnamurthy.jl +++ b/examples/AssetPricing/HeKrishnamurthy.jl @@ -20,8 +20,7 @@ end # Denote αI ratio of risky asset holdings to total capital w + H # households have log utility. Assume that they die continuously so that we can treat financial wealth as the thigns to maximize in t + dt (instead of total wealth), so that we can treat labor income easily. -# We make assumptions so that the household sector chooses to keep a minimum of λ w^h in short-term deb issued by intermediary. Formally, we assume that a fraction λ can only ever invest in riskless bond while remaining can invest in intermediation market. but both are poooled at the end of the period. - +# We make assumptions so that the household sector chooses to keep a minimum of λ w^h in short-term deb issued by intermediary. Formally, we assume that a fraction λ can only ever invest in riskless bond while remaining can invest in intermediation market # Denote pS the wealth-to-consumption ratio of specialists function (m::HeKrishnamurthy)(state::NamedTuple, y::NamedTuple) @@ -38,18 +37,11 @@ function (m::HeKrishnamurthy)(state::NamedTuple, y::NamedTuple) pxx = - (1 + l) / (x / pS + (1 - x) * ρ)^2 * (- 1 / pS^2 * pSx - 1 / pS^2 * pSx - x / pS^2 * pSxx + 2 * x / pS^3 * pSx^2) + 2 * (1 + l) / (x / pS + (1 - x) * ρ)^3 * (1 / pS - ρ - x / pS^2 * pSx)^2 @assert x / pS + (1 - x) * ρ ≈ (1 + l) / p - # market clearing for equity gives ((1-x) * (1-λ) * αH + x) * αI = 1 - # If unconstrained, assumption is that αI = 1.0 - αH = 1 + αI = 1 / (1 - λ * (1 - x)) - # total equity of households cannot be higher than m * networth of specialists - if (1 - x) * (1 - λ) * αH > x * m - # constrained case - αH = x * m / ((1 - x) * (1 - λ)) + if x * (1 + m) * αI < 1 αI = 1 / (x * (1 + m)) end - @assert ((1-x) * (1-λ) * αH + x) * αI ≈ 1.0 - # we have the usual feeback loop between asset prices and wealth # σx = x * (αI - 1) * (σ + σp) σx = x * (αI - 1) * σ / (1 - x * (αI - 1) * px / p) @@ -58,8 +50,7 @@ function (m::HeKrishnamurthy)(state::NamedTuple, y::NamedTuple) @assert σx ≈ x * (αI - 1) * (σ + σp) σR = σ + σp κ = γ * (αI * σR - σpS) - # as usual, μx = x * (1 - x) * (growth rate of wealth of specialists minus growth rate of wealth of households) - μx = x * (1 - x) * ((αI - (1 - λ) * αI * αH) * κ * σR - 1 / pS + ρ - l / ((1 - x) * p)) + μx = x * (1 - x) * ((αI - (1 - αI * x) / (1-x)) * κ * σR - 1 / pS + ρ - (αI - (1 - αI * x) / (1 - x)) * σR^2 - l / ((1 - x) * p)) if (iter == 0) && (μx < 0) iter += 1 pSx = pSx_down @@ -68,24 +59,17 @@ function (m::HeKrishnamurthy)(state::NamedTuple, y::NamedTuple) μpS = pSx / pS * μx + 0.5 * pSxx / pS * σx^2 μp = px / p * μx + 0.5 * pxx / p * σx^2 r = 1 / p + g + μp + σ * σp - κ * σR - @assert x * (r + αI * κ * σR - 1 / pS) + (1 - x) * (r + (1 - λ) * αI * αH * κ * σR - ρ + l / ((1-x) * p)) ≈ g + μp + σ * σp σCS = κ / γ μCS = (r - ρ) / γ + (1 + 1 / γ) / (2 * γ) * κ^2 pSt = - (1 / pS + μCS + μpS + σCS * σpS - r - κ * (σCS + σpS)) μR = r + κ * σR - (; pSt), (; pS, p, r, κ, σp, μR, σR, αI, μx, σx, x) + σI = αI * σR + (; pSt), (; pS, p, r, κ, σp, μR, σR, αI, μx, σx, x, px, pxx, σI) end - - m = HeKrishnamurthy() -xn = 1000 +xn = 100 stategrid = OrderedDict(:x => range(0, 1, length = xn+2)[2:(end-1)].^1.5) yend = OrderedDict(:pS => ones(length(stategrid[:x]))) y, residual_norm, a = pdesolve(m, stategrid, yend) -@assert residual_norm <= 1e-5 - - - - - +@assert residual_norm <= 1e-5 \ No newline at end of file diff --git a/src/finiteschemesolve.jl b/src/finiteschemesolve.jl index 0cdc2f2..20fb09b 100644 --- a/src/finiteschemesolve.jl +++ b/src/finiteschemesolve.jl @@ -43,7 +43,7 @@ function implicit_timestep(G!, ypost, Δ; is_algebraic = fill(false, size(ypost) end # Solve for steady state -function finiteschemesolve(G!, y0; Δ = 1.0, is_algebraic = fill(false, size(y0)...), iterations = 100, inner_iterations = 10, verbose = true, inner_verbose = false, method = :newton, autodiff = :forward, maxdist = sqrt(eps()), scale = 10.0, J0c = (nothing, nothing), minΔ = 1e-9, y̲ = fill(-Inf, length(y0)), ȳ = fill(Inf, length(y0)), reformulation = :smooth, maxΔ = Inf, autoscale = true, kwargs...) +function finiteschemesolve(G!, y0; Δ = 1.0, is_algebraic = fill(false, size(y0)...), iterations = 100, inner_iterations = 10, verbose = true, inner_verbose = false, method = :newton, autodiff = :forward, maxdist = sqrt(eps()), innerdist = sqrt(eps()), scale = 10.0, J0c = (nothing, nothing), minΔ = 1e-9, y̲ = fill(-Inf, length(y0)), ȳ = fill(Inf, length(y0)), reformulation = :smooth, maxΔ = Inf, autoscale = true, kwargs...) ypost = y0 ydot = zero(y0) G!(ydot, ypost) @@ -65,7 +65,7 @@ function finiteschemesolve(G!, y0; Δ = 1.0, is_algebraic = fill(false, size(y0) end while (iter < iterations) & (Δ >= minΔ) & (residual_norm > maxdist) iter += 1 - y, nlresidual_norm = implicit_timestep(G!, ypost, Δ; is_algebraic = is_algebraic, verbose = inner_verbose, iterations = inner_iterations, method = method, autodiff = autodiff, maxdist = maxdist, J0c = J0c, y̲ = y̲, ȳ = ȳ, reformulation = reformulation, kwargs...) + y, nlresidual_norm = implicit_timestep(G!, ypost, Δ; is_algebraic = is_algebraic, verbose = inner_verbose, iterations = inner_iterations, method = method, autodiff = autodiff, maxdist = innerdist, J0c = J0c, y̲ = y̲, ȳ = ȳ, reformulation = reformulation, kwargs...) G!(ydot, y) if any(y̲ .!= -Inf) || any(ȳ .!= Inf) mask = y̲ .+ eps() .<= y .<= ȳ .- eps() # only unconstrained ydot is relevant for residual_norm calculation @@ -74,7 +74,7 @@ function finiteschemesolve(G!, y0; Δ = 1.0, is_algebraic = fill(false, size(y0) residual_norm, oldresidual_norm = norm(ydot) / length(ydot), residual_norm end residual_norm = isnan(residual_norm) ? Inf : residual_norm - if nlresidual_norm <= maxdist + if nlresidual_norm <= innerdist # if the implicit time step is correctly solved if verbose @printf "%4d %8.4e %8.4e\n" iter Δ residual_norm