+: Memory efficient
-:Non-standard, Code complexity, Cannot use external arrays

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.

@JuliaRegistrator register()

How can I use this package to solve a finite horizon lifecycle model such as this one by Moll?
His Matlab code is here.

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?

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))

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?

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?

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> 

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)

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?

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

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.

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.

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:

https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/about-citation-files

For now, I use the same formula as in the case of the uniform grid. Check that it is correct.

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.

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!

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?

At some point rely exclusively on names rather than orders. But then I would need namedarrays to specify initial guess with multidimensional states.

Hey would it be possible to add an example solving the NCGM & compare w/ closed form solutions?

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.