Skip to content

Commit

Permalink
add HK + BS + print more intermediates
Browse files Browse the repository at this point in the history
  • Loading branch information
matthieugomez committed Jul 22, 2024
1 parent 842d341 commit 3f41442
Show file tree
Hide file tree
Showing 11 changed files with 274 additions and 96 deletions.
12 changes: 6 additions & 6 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -88,12 +88,12 @@ where `S(x)` is the value of exercising the option. Notice the traditional "valu

# Examples
The [examples folder](https://github.com/matthieugomez/EconPDEs.jl/tree/master/examples) solves a variety of models:
- *Habit Model* (Campbell Cochrane (1999) and Wachter (2005))
- *Long Run Risk Model* (Bansal Yaron (2004))
- *Disaster Model* (Wachter (2013))
- *Heterogeneous Agent Models* (Garleanu Panageas (2015), Di Tella (2017), Haddad (2015))
- *Consumption with Borrowing Constraint* (Wang Wang Yang (2016), Achdou Han Lasry Lions Moll (2018))
- *Investment with Borrowing Constraint* (Bolton Chen Wang (2009))
- *Habit Model*: Campbell Cochrane (1999) and Wachter (2005)
- *Long Run Risk Model*: Bansal Yaron (2004)
- *Disaster Model*: Wachter (2013)
- *Heterogeneous Agent Models*: He Krishnamurthy (2013), Brunnermeir Sannikov (2013), Garleanu Panageas (2015), Di Tella (2017), Haddad (JMP)
- *Consumption with Borrowing Constraint*: Wang Wang Yang (2016), Achdou Han Lasry Lions Moll (2018)
- *Investment with Borrowing Constraint*: Bolton Chen Wang (2009)

## Installation
The package is registered in the [`General`](https://github.com/JuliaRegistries/General) registry and so can be installed at the REPL with `] add EconPDEs`.
Expand Down
28 changes: 13 additions & 15 deletions examples/AssetPricing/BansalYaron.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,19 +17,6 @@ Base.@kwdef struct BansalYaronModel
ψ::Float64 = 1.5
end

function initialize_stategrid(m::BansalYaronModel; μn = 30, vn = 30)
μdistribution = Normal(m.μbar, sqrt(m.νμ^2 * m.vbar / (2 * m.κμ)))
μs = range(quantile(μdistribution, 0.025), quantile(μdistribution, 0.975), length = μn)
νdistribution = Gamma(2 * m.κv * m.vbar / m.νv^2, m.νv^2 / (2 * m.κv))
vs = range(quantile(νdistribution, 0.025), quantile(νdistribution, 0.975), length = vn)
OrderedDict( => μs, :v => vs)
end

function initialize_y(m::BansalYaronModel, stategrid)
μn = length(stategrid[])
vn = length(stategrid[:v])
OrderedDict(:p => ones(μn, vn))
end

function (m::BansalYaronModel)(state::NamedTuple, y::NamedTuple)
(; μbar, vbar, κμ, νμ, κv, νv, ρ, γ, ψ) = m
Expand Down Expand Up @@ -71,8 +58,19 @@ end

# Bansal Yaron (2004)
m = BansalYaronModel()
stategrid = initialize_stategrid(m)
yend = initialize_y(m, stategrid)

## define state space
μn = 30
vn = 30
μdistribution = Normal(m.μbar, sqrt(m.νμ^2 * m.vbar / (2 * m.κμ)))
μs = range(quantile(μdistribution, 0.025), quantile(μdistribution, 0.975), length = μn)
νdistribution = Gamma(2 * m.κv * m.vbar / m.νv^2, m.νv^2 / (2 * m.κv))
vs = range(quantile(νdistribution, 0.025), quantile(νdistribution, 0.975), length = vn)
stategrid = OrderedDict( => μs, :v => vs)

## define initial guess
yend = OrderedDict(:p => ones(μn, vn))

result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5

Expand Down
120 changes: 120 additions & 0 deletions examples/AssetPricing/BrunnermeirSannikov.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
# 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

Base.@kwdef mutable struct BrunnermeierSannikov
# See Table H Numerical examples
ρE::Float64 = 0.06
ρH::Float64 = 0.06
γ::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
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 / κ)
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)
end
else
# in the second region Ψ = 1
Ψ = 1
# We get q through market clearing equation after plugging i = (p - 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)
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) / κ
Φ = log(q) / κ
σx =* Ψ - x) * σ / (1 -* Ψ - x) * 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))
if (iter == 0) && (μx < 0)
iter += 1
pEx, pHx = pEx_down, pHx_down
@goto start
end
μ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
# 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)
end



5 changes: 1 addition & 4 deletions examples/AssetPricing/CampbellCochrane.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,6 @@ function initialize_stategrid(m::CampbellCochraneModel; sn = 1000)
OrderedDict(:s => vcat(slow[1:(end-1)], shigh[2:end]))
end

function initialize_y(m::CampbellCochraneModel, stategrid)
OrderedDict(:p => ones(length(stategrid[:s])))
end

function (m::CampbellCochraneModel)(state::NamedTuple, y::NamedTuple)
(; μ, σ, γ, ρ, κs, b) = m
Expand Down Expand Up @@ -62,7 +59,7 @@ end

m = CampbellCochraneModel()
stategrid = initialize_stategrid(m)
yend = initialize_y(m, stategrid)
yend = OrderedDict(:p => ones(length(stategrid[:s])))
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5

Expand Down
15 changes: 4 additions & 11 deletions examples/AssetPricing/GarleanuPanageas.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,20 +26,13 @@ Base.@kwdef struct GarleanuPanageasModel
ω::Float64 = 0.92
end

function initialize_stategrid(m::GarleanuPanageasModel; xn = 200)
OrderedDict(:x => range(0.0, 1.0, length = xn))
end

function initialize_y(m::GarleanuPanageasModel, stategrid)
xn = length(stategrid[:x])
OrderedDict(:pA => ones(xn), :pB => ones(xn), :ϕ1 => ones(xn), :ϕ2 => ones(xn))

function (m::GarleanuPanageasModel)(state::NamedTuple, y::NamedTuple)
# pA is wealth / consumption ratio of agent A
# pB is wealth / consumption ratio of agent B
# ϕ1 is value of claim that promises 1 today and then grows at rate μ - δ- δ1
# ϕ2 is value of claim that promises 1 today and then grows at rate μ - δ- δ2
end

function (m::GarleanuPanageasModel)(state::NamedTuple, y::NamedTuple)
(; γA, ψA, γB, ψB, ρ, δ, νA, μ, σ, B1, δ1, B2, δ2, ω) = m
(; x) = state
(; pA, pAx_up, pAx_down, pAxx, pB, pBx_up, pBx_down, pBxx, ϕ1, ϕ1x_up, ϕ1x_down, ϕ1xx, ϕ2, ϕ2x_up, ϕ2x_down, ϕ2xx) = y
Expand Down Expand Up @@ -91,7 +84,7 @@ function (m::GarleanuPanageasModel)(state::NamedTuple, y::NamedTuple)
end

m = GarleanuPanageasModel()
stategrid = initialize_stategrid(m; xn = 1000)
yend = initialize_y(m, stategrid)
stategrid = OrderedDict(:x => range(0.0, 1.0, length = 100))
yend = OrderedDict(:pA => ones(length(stategrid[:x])), :pB => ones(length(stategrid[:x])), :ϕ1 => ones(length(stategrid[:x])), :ϕ2 => ones(length(stategrid[:x])))
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5
44 changes: 16 additions & 28 deletions examples/AssetPricing/Haddad.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,31 +22,6 @@ Base.@kwdef struct HaddadModel
ψ::Float64 = 1.5
end


function initialize_stategrid(m::HaddadModel; μn = 30, vn = 30)
(; μbar, vbar, κμ, νμ, κv, νv, αbar, λ, ρ, γ, ψ) = m

σ = sqrt(νμ^2 * vbar / (2 * κμ))
μmin = quantile(Normal(μbar, σ), 0.001)
μmax = quantile(Normal(μbar, σ), 0.999)
μs = range(μmin, μmax, length = μn)

α = 2 * κv * vbar / νv^2
β = νv^2 / (2 * κv)
vmin = quantile(Gamma(α, β), 0.001)
vmax = quantile(Gamma(α, β), 0.999)
vs = range(vmin, vmax, length = vn)

OrderedDict( => μs, :v => vs)
end

function initialize_y(m::HaddadModel, stategrid)
μn = length(stategrid[])
vn = length(stategrid[:v])
OrderedDict(:p => ones(μn, vn))
end


function (m::HaddadModel)(state::NamedTuple, y::NamedTuple)
(; μbar, vbar, κμ, νμ, κv, νv, αbar, λ, ρ, γ, ψ) = m
(; μ, v) = state
Expand Down Expand Up @@ -78,13 +53,26 @@ function (m::HaddadModel)(state::NamedTuple, y::NamedTuple)
# Interest rate r
#r = ρ + μc / ψ - (1 - 1 / ψ) / (2 * γ) * κ2 - (1 / ψ - γ) / (2 * γ * (ψ - 1)) * σp2 - (1 - γ) / (γ * ψ) * (κμ * σμ + κv * σv) + 1 / ψ * (κ_Zc * σc + κμ * σμ + κv * σv)
#out = p * (1 / p + μc + μp - r - κ_Zc * σc - κ_Zμ * σp_Zμ - κ_Zv * σp_Zv)

pt = - p * (1 / p - ρ + (1 - 1 / ψ) * (μc - 0.5 * γ * σc^2 * (1 - (αstar - 1)^2)) + μp + (0.5 * (1 / ψ - γ) / (1 - 1 / ψ) + 0.5 * γ * (1 - 1 / ψ) * (αstar - 1)^2) * σp2)
return (; pt)
end

m = HaddadModel()
stategrid = initialize_stategrid(m)
yend = initialize_y(m, stategrid)
# initialize grid of μ and v
σ = sqrt(m.νμ^2 * m.vbar / (2 * m.κμ))
μmin = quantile(Normal(m.μbar, σ), 0.001)
μmax = quantile(Normal(m.μbar, σ), 0.999)
μs = range(μmin, μmax, length = 30)
α = 2 * m.κv * m.vbar / m.νv^2
β = m.νv^2 / (2 * m.κv)
vmin = quantile(Gamma(α, β), 0.001)
vmax = quantile(Gamma(α, β), 0.999)
vs = range(vmin, vmax, length = 30)
stategrid = OrderedDict( => μs, :v => vs)

# guess for price-dividend ratio
yend = OrderedDict(:p => ones(length(stategrid[]), length(stategrid[:v])))

# solve
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5
91 changes: 91 additions & 0 deletions examples/AssetPricing/HeKrishnamurthy.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
# He Krishnamurthy "Intermediary Asset Pricing" AER (2013)

using EconPDEs

Base.@kwdef mutable struct HeKrishnamurthy
# See Table 2—Parameters and Targets
m::Float64 = 4 # intermediation multiplier
λ::Float64 = 0.6 # Debt/assets of intermediary sector
g::Float64 = 0.02 # dividend growth
σ::Float64 = 0.09 # dividend volatility
ρ::Float64 = 0.04 # time discount rate
γ::Float64 = 2.0 # risk aversion specialist
l::Float64 = 1.84 # labor income ratio
end

# Description of the model
# Households cannot invest directly in the risky asset market: can only invest in bonds and intermediaries. Households unwilling to invest more than m * w_t in the equity of intermediary, where w_t is networth intermediary (constant on outside equity financing)
# Specialists have the knowledge to invest in the risky assets, and can invest in the risky asset on behalf of the households
# Specialists accept household funds and allocate between risky and riskless. networth of specialist is w and wealth allocated by households is H.
# 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.


# Denote pS the wealth-to-consumption ratio of specialists
function (m::HeKrishnamurthy)(state::NamedTuple, y::NamedTuple)
(; m, λ, g, σ, ρ, γ, l) = m
(; pS, pSx_up, pSx_down, pSxx) = y
(; x) = state
iter = 0
pSx = pSx_up
@label start
# market clearing for goods gives x / pS + (1 - x) * ρ = (1 + l) / p
# This allows to obtain p (and its derivatives) in terms of pS (and its derivatives)
p = (1 + l) / (x / pS + (1 - x) * ρ)
px = - (1 + l) / (x / pS + (1 - x) * ρ)^2 * (1 / pS - ρ - x / pS^2 * pSx)
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 - λ))
α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)
σpS = pSx / pS * σx
σp = px / p * σx
@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))
if (iter == 0) && (μx < 0)
iter += 1
pSx = pSx_down
@goto start
end
μ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)
end



m = HeKrishnamurthy()
xn = 1000
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





Loading

0 comments on commit 3f41442

Please sign in to comment.