Skip to content

Commit

Permalink
update HK and BS
Browse files Browse the repository at this point in the history
  • Loading branch information
matthieugomez committed Jul 24, 2024
1 parent 3f41442 commit 654c0a2
Show file tree
Hide file tree
Showing 5 changed files with 98 additions and 105 deletions.
3 changes: 2 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
160 changes: 84 additions & 76 deletions examples/AssetPricing/BrunnermeirSannikov.jl
Original file line number Diff line number Diff line change
@@ -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)
2 changes: 1 addition & 1 deletion examples/AssetPricing/GarleanuPanageas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
32 changes: 8 additions & 24 deletions examples/AssetPricing/HeKrishnamurthy.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
6 changes: 3 additions & 3 deletions src/finiteschemesolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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 =.+ eps() .<= y .<=.- eps() # only unconstrained ydot is relevant for residual_norm calculation
Expand All @@ -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
Expand Down

0 comments on commit 654c0a2

Please sign in to comment.