matthieugomez / econpdes.jl Goto Github PK
View Code? Open in Web Editor NEWSolve non-linear HJB equations.
License: Other
Solve non-linear HJB equations.
License: Other
This issue is used to trigger TagBot; feel free to unsubscribe.
If you haven't already, you should update your TagBot.yml
to include issue comment triggers.
Please see this post on Discourse for instructions and more details.
If you'd like for me to do this for you, comment TagBot fix
on this issue.
I'll open a PR within a few hours, please be patient!
At some point rely exclusively on names rather than orders. But then I would need namedarrays
to specify initial guess with multidimensional states.
Hi,
I've been using this package for a while and would love to acknowledge it in my work. Do you want to create a citation file so we could properly cite it? Something like below:
The Aiyagari Continuous Time example has not reappeared in the EconPDEs example. AiyagariContinuousTime.jl.
More generally a current and present value Hamiltonian growth model would be nice.
Wachter (2005) calibration should have b=0.011, not 0.011*4 as it is right now in the code. It is already annualized; there is no need to multiply it by 4 (Wachter reports already annualized parameters).
I verified it against published 2005 FRL paper by computing E[r_f]. With b=0.011 it matches the published result.
I'm having trouble solving a simple finite horizon bvp.
A simple example:
V_t + 2*V_x =0
V(0,t) = t
The closed form solution is V(x,t) = .5*(2t-x).
Is there a direct way to solve this using EconPDEs?
In the AchdouHanLasryLionsMoll one asset example, the crucial pdesolve command has disappeared:
m = AchdouHanLasryLionsMollModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0)
Hello:
I'm very new to Julia and I'm not sure if this is my problem or the code has issue because of the Julia version change or etc.
All of you codes in example folder shows following error.
ERROR: LoadError: syntax: invalid assignment location "; c, r, ρ, σ, a, T"
I removed ';' sign in the function code like this
function (m::ArbitrageHoldingCosts)(state::NamedTuple, y::NamedTuple, τ::Number)
(c, r, ρ, σ, a, T) = m
(z) = state
(F, Fz_up, Fz_down, Fzz) = y
but then Julia shows:
ERROR: LoadError: MethodError: no method matching iterate
I really want to learn from your code. Could you please help me out with this issue?
pdesolve(m, state,y0) returns the error "MethodError: no method matching norm(::Array{Float64,3},::Float64)" for DiTella model. I think the problem is that since we have a 3D array in the form of (pA, pB, and p), at some point the FiniteDifferenceScheme tries to calculate norm of this 3D array, which is not acceptable. Any idea how to fix this? Thanks.
@JuliaRegistrator register()
I tried README example w/ Julia v 1.6.2:
julia> using EconPDEs
julia> stategrid = (;μ = range(-0.05, 0.1, length = 500))
(μ = -0.05:0.0003006012024048096:0.1,)
julia> solend = (; V = ones(500))
(V = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],)
julia> function f(state::NamedTuple, sol::NamedTuple)
μbar = 0.018 ; ϑ = 0.00073 ; θμ = 0.252 ; νμ = 0.528 ; ρ = 0.025 ; ψ = 1.5 ; γ = 7.5
(; μ) = state
(; V, Vμ_up, Vμ_down, Vμμ) = sol
# note that pdesolve will directly compute the derivatives of the valuefunction.
# up and down correspond to the upward and downard derivative obtained by first difference
Vμ = (μ <= μbar) ? Vμ_up : Vμ_down
Vt = - (1 - ρ * V + (1 - 1 / ψ) * (μ - 0.5 * γ * ϑ) * V + θμ * (μbar - μ) * Vμ +
0.5 * νμ^2 * ϑ * Vμμ + 0.5 * (1 / ψ - γ) / (1- 1 / ψ) * νμ^2 * ϑ * Vμ^2/V)
(; Vt)
end
ERROR: syntax: invalid assignment location "; μ" around REPL[7]:3
Stacktrace:
[1] top-level scope
@ REPL[7]:1
julia>
Hi, thanks for a great package. I am having trouble with running the pdesolve function. In running your example codes, I keep bumping into this error :
pdesolve(m, grid, y0)
ERROR: MethodError: no method matching pdesolve(::CampbellCochraneModel, ::Plots.#grid, ::DataStructures.OrderedDict{Symbol,Array{Float64,1}})
Closest candidates are:
pdesolve(::Any, ::DataStructures.OrderedDict, ::DataStructures.OrderedDict; is_algebraic, kwargs...) at C:\Users\jhji.julia\v0.6\EconPDEs\src\pdesolve.jl:235
Stacktrace:
[1] eval(::Module, ::Any) at .\boot.jl:235
Above is case when I run Campbell and Cochrane example, but all others return similar. I am using following packages :
IJulia,DataFrames,Plots,TaylorSeries,QuantEcon,Interpolations,Optim,Roots,Distributions,Combinatorics,DataStructures,BandedMatrices, AiyagariContinuousTime,EconPDEs
Thanks in advance!
Edit
I think I figured out the typo in the example. I think you meant to type pdesolve(m, state, y0) in example script.
Since your implementation calculates finite differences in a fully-implicit fashion, I cannot just solve an eigenvalue problem on the generator matrix to find the stationary distribution over the state space. My thinking was that I should solve for the stationary distribution simultaneously with the HJB, since the policy function and its derivatives may appear in the KFE.
However, I am not successful in constructing a simple extension to AchdouHanLasryLionsMollModel
. I am unclear about whether this can be done with EconPDEs at all, because my hand-written FD schemes for KFEs always needed some sort of normalization in every iteration.
using EconPDEs, Distributions
struct AchdouHanLasryLionsMollModel
# income process parameters
κy::Float64
ybar::Float64
σy::Float64
r::Float64
# utility parameters
ρ::Float64
γ::Float64
amin::Float64
amax::Float64
end
function AchdouHanLasryLionsMollModel(;κy = 0.1, ybar = 1.0, σy = 0.07, r = 0.03, ρ = 0.05, γ = 2.0, amin = 0.0, amax = 500.0)
AchdouHanLasryLionsMollModel(κy, ybar, σy, r, ρ, γ, amin, amax)
end
function initialize_state(m::AchdouHanLasryLionsMollModel; yn = 5, an = 50)
κy = m.κy ; ybar = m.ybar ; σy = m.σy ; ρ = m.ρ ; γ = m.γ ; amin = m.amin ; amax = m.amax
distribution = Gamma(2 * κy * ybar / σy^2, σy^2 / (2 * κy))
ymin = quantile(distribution, 0.001)
ymax = quantile(distribution, 0.999)
ys = collect(range(ymin, stop = ymax, length = yn))
as = collect(range(amin, stop = amax, length = an))
OrderedDict(:y => ys, :a => as)
end
function initialize_y(m::AchdouHanLasryLionsMollModel, state)
distribution = Normal((minimum(state[:y]) + maximum(state[:y]))/2, (minimum(state[:y]) + maximum(state[:y]))/8)
g_prep = [pdf(distribution, y) * 1 / (√a+1) for y in state[:y], a in state[:a]]
OrderedDict(:v => [(y + m.r * a)^(1-m.γ)/(1-m.γ)/m.ρ for y in state[:y], a in state[:a]],
:g => g_prep/sum(g_prep))
end
function (m::AchdouHanLasryLionsMollModel)(state, value)
κy = m.κy ; σy = m.σy ; ybar = m.ybar ; r = m.r ; ρ = m.ρ ; γ = m.γ ; amin = m.amin ; amax = m.amax
y, a = state.y, state.a
ymin = minimum(state[:y]); ymax = maximum(state[:y])
v, vy, va, vyy, vya, vaa = value.v, value.vy, value.va, value.vyy, value.vya, value.vaa
g, gy, ga, gyy, gya, gaa = value.g, value.gy, value.ga, value.gyy, value.gya, value.gaa
μy = κy * (ybar - y)
va = max(va, eps())
# There is no second derivative at 0 so just specify first order derivative
c = va^(-1 / γ)
μa = y + r * a - c
if (a ≈ amin) && (μa <= 0.0)
va = (y + r * amin)^(-γ)
c = y + r * amin
μa = 0.0
end
# this branch is unnecessary when individuals dissave at the top (default)
if (a ≈ amax) && (μa >= 0.0)
va = (((ρ - r) / γ + r) * a)^(-γ)
c = ((ρ - r) / γ + r) * a
μa = y + (r - ρ) / γ * a
end
# - ∂μa/∂a * g - μa*∂g/∂a -∂μy/∂y*g -μy*∂g/∂y + 1/2 *σy^2 *∂^2g/∂y^2
gt = - (r - (-1 / γ) * va^(-1 / γ - 1) * vaa) * g - μa * ga + κy * g - μy * gy + 1/2 * σy^2 * gyy
vt = c^(1 - γ) / (1 - γ) + va * μa + vy * μy + 0.5 * vyy * σy^2 - ρ * v
return (vt, gt), (μy, μa), (v = v, c = c, va = va, vy = vy, y = y, a = a, μa = μa, vaa = vaa)
end
m = AchdouHanLasryLionsMollModel()
state = initialize_state(m)
y0 = initialize_y(m, state)
y, result, distance = pdesolve(m, state, y0, iterations = 20)
Basically, I added the distribution g
to the solution object and added its law of motion, although I am not sure about the boundary condition on the a
axis, but I tried different versions. The solution iteration does not converge and the result is odd in any case.
I am working through your GarleanuPanageasModel and I think I found a typo, but I cannot be sure, because I do not know with certainty how your pA
and mcA
maps into the paper's g^A
and n^A
. The former is suggested by the definition of κ
, the latter by the definition of r
.
Anyway, the denominator of equation 9 has the terms concerning A
with a plus sign and the terms with B
with a minus sign. This is an error, unless pAx/pA = - gAx/gA
.
Furthermore, equation A.16 in the appendix defines n^A
which resembles your mcA
almost perfectly. But, the coefficient in the second line agrees with you definition of mcA
, only if γA/(1-γA) = 1
, which is not true.
Is there a documentation that could accompany the example?
Right now, I assume that derivative of value function at boundary is zero. But it does not have to be. The method works for any value of the derivative at the boundary (because it is enough to pinpoint value of function outside the grid).
In other words, I should allow users to input their own conditions about derivative at boundary. Would be useful for investment problem of firm in Bolton, Chen, Wang.
Something like a OrderedDict or a tuple
(Vμ = [0.0, 0.0])
(Vμ = zeros(30), Vσ = zeros(10))
I would like to replace my own HJB solver with this package. In the process I realized that it cannot handle more than two state variables.
One obvious reason is that the differentiate
function has methods for one and two states. Did you also assume 2 states in the rest of the package, or would it be enough to provide another method for differentiate
?
+: Memory efficient
-:Non-standard, Code complexity, Cannot use external arrays
For now, I use the same formula as in the case of the uniform grid. Check that it is correct.
Is there an easy way to simulate the variables after solving?
For example after I solve a simple consumption problem I would like to use the program output to plot for one path of shocks:
income over time
consumption over time
wealth level over time
etc
Consider the simple problem (w/ closed form solution):
I tried coding it up following your investment example:
using EconPDEs
stategrid = OrderedDict(:k => range(0.0, 5.0, length = 1000))
solend = OrderedDict(:V => ones(1000))
function f(state::NamedTuple, sol::NamedTuple)
r = 0.05; δ = 0.10; z = 0.20;
k = state.k
V, Vk_up, Vk_down, Vkk = sol.V, sol.Vk_up, sol.Vk_down, sol.Vkk
#
Vk = Vk_up
iter = 0
@label start
i = Vk - 1.0
μk = i - δ*k
if (iter == 0) & (μk <= 0)
iter += 1
Vk = Vk_down
@goto start
end
Vt = - (z*k - i -0.5*(i^2.0) + μk*Vk - r*V)
(Vt = Vt,)
end
y, residual_norm = pdesolve(f, stategrid, solend)
I also tried:
using EconPDEs
stategrid = OrderedDict(:k => range(0.0, 5.0, length = 1000))
solend = OrderedDict(:V => ones(1000))
function f(state::NamedTuple, sol::NamedTuple)
r = 0.05; δ = 0.10; z = 0.20;
k = state.k
V, Vk_up, Vk_down, Vkk = sol.V, sol.Vk_up, sol.Vk_down, sol.Vkk
#
Vk = Vk_up
iter = 0
@label start
i = Vk - 1.0
μk = i - δ*k
Vk = (μk >= 0) ? Vk_up : Vk_down
Vt = - (z*k - i -0.5*(i^2.0) + μk*Vk - r*V)
(Vt = Vt,), (; V, Vk, Vkk, k)
end
y, residual_norm = pdesolve(f, stategrid, solend)
I'm not sure I understand how solve a model correctly w/ your package.
In both cases I get the same wrong value function not equal to the closed form solution.
The affine, closed-form solution:
How can I solve this simple firm investment problem with your package?
How can I use this package to solve a finite horizon lifecycle model such as this one by Moll?
His Matlab code is here.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.