Solving Boundary Value Differential Equation Problems for Economics

[MainBatDealer]

Hi all,

I am interested in using the BoundaryValueDiffEq part of DifferentialEquations to solve economic growth problems in continuous time. The problem is, I can’t get a relatively simple example to work!

The example I am trying to solve is the Cass-Koopmans-Ramsey model from Aghion and Howitt’s Endogenous Growth Theory. It is a nonlinear system of differential equations

\dot{K} = K^\alpha - \delta K - 1 / \lambda

\dot{C} = -C (\alpha K ^ {\alpha - 1} - \delta - \rho)

With an initial capital stock K(0) = K_0 and an end-of-time condition

\lim_{t\rightarrow \infty} e^{-\rho t} \frac{K}{C} = 0.

\alpha, \delta, and \rho are parameters in the model and represent the productivity of capital, the depreciation rate of capital, and the intertemporal preference of a utility maximizer.

The steady-state for this problem is easy to solve for analytically, but I’m interested in what happens along the transitional path, hence, the need for numerical analysis.

Recognizing that BoundaryValueDiffEq will not appreciate a boundary condition at t = \infty, I made the following transformation:

\tau = t / (1 + t),

which transforms my system to

\partial K / \partial \tau = (K^\alpha - \delta K - 1/\lambda) / (1 - \tau) ^ 2

\partial C/ \partial \tau = -C (\alpha K ^ {\alpha - 1} - \delta - \rho) / (1 - \tau) ^ 2

and my terminal condition to

\lim_{\tau \rightarrow 1} e^{-1/(1 - \tau)} \frac{K}{C} = 0

The terminal condition is still somewhat convoluted, so I can exploit what I know about the steady-state to replace it. I know that the system should reach the steady-state at \tau = 1, so, after solving for the steady-state, I replace my terminal condition with

C(1) = K(1) ^\alpha - \delta K(1)

where

K(1) = \left[\frac{\rho + \delta}{\alpha}\right]^{\frac{1}{\alpha - 1}}

Now that my system is defined to my satisfaction, I turn to the code:

using DifferentialEquations
using Plots

#Parameters
ρ = 0.04
δ = 0.10
α = 0.66
p = [ρ, δ, α]

#Steady-States
K_ss = ((δ + ρ) / α) ^ (1 / (α - 1)) 
C_ss = (K_ss ^ α - δ * K_ss)

#Initial Capital
K0 = 90.0

#Defining the system
function end_sav!(du, u, p, t)
    ρ = p[1]
    δ = p[2]
    α = p[3]

    K = u[1]
    C = u[2]

    du[1] = (K^α - δ * K - C ) / (1 - t) ^ 2
    du[2] = - C * (ρ + δ - α * K ^ (α - 1)) / (1 - t) ^ 2
end

# Define the boundary conditions
function bc!(residual, u, t)

    residual[1] = u[1][1] - K0
    residual[2] = u[2][end] - C_ss
end

# Prevent the solver from trying
# negative values of capital and consumption
function g(u, resid)
    resid[1] = u[1] - 0.001
    resid[2] = u[2] - 0.001
end

end_sav_prob = BVProblem(end_sav!, bc!, [0.90, C_ss], (0.0, 1.), p; callback = GeneralDomain(g))

sol = solve(end_sav_prob)

I’ve tried different algorithms, different starting values, different time scales, varying the parameters, etc. Nothing seems to work. The most common issue is instability. Any suggestions on what’s going on? Additionally, I can’t be the first person to try to use Julia to solve these kinds of problems. Does anyone know of any resources/packages for working with continuous-time models in Julia? QuantEcon has some really good stuff, but as far as I can tell they stick to discrete-time methods.

[ChrisRackauckas]

It could be that the solver that is used is unstable for this. Try some of the ODEInterface.jl BVP solvers and see. BVPs are probably the weakest point of DifferentialEquations.jl right now.

[jlperla]

The basic problem is that these are not just BVPs but rather are a particular case of IVP and BVP found in control-theory. Parts of the system of equations naturally look forward and others backwards. I think they need specialized algorithms (if trivial, such as variations on the shooting method) that haven’t been built into the DifferentialEquations.jl ecosystem. Maybe there is an existing combination but I am not sure of it.
Possibly it wouldn’t take that much work, but would require someone to analyze, implement, and test these sorts of problems.

Part of this is also that economists often work with the recursive version of this problem - which leads to a different sort of differential equation, and one best solved through specialized upwind finite difference methods (e.g. see the neoclassical growth models in codes - HACT - Benjamin Moll - University of Princeton).

These also don’t fit neatly into DifferentialEquations yet (at least the upwind stuff for more complicated problems) since they require specialized tricks that break the differential operators into a forward and backward part. These control problems have a forwards and a backwards looking component in them, so formalizing the interface in DifferentialEquations isn’t obvious - though certainly possible.

[MainBatDealer]

Hmm. Well, I was hoping to avoid writing a new solver for these kinds of problems. I know that lots of economists have written code and algorithms to solve these, but it’s pretty much all in Matlab. Maybe I’ll get around to porting some of it to Julia but until then I guess I’ll make use of my University license.

[jlperla]

The differences between Matlab vs. julia is irrelevant for a specialized algorithm and simple problems like this if you aren’t designing a package. You could port over the matlab code for the neoclassical growth from Ben Moll’s site almost verbatim in an hour or less if you aren’t trying to improve a fancy interface - and I am sure that Ben would love to host a julia version on his site next to the matlab.

For designing variations using these algorithms for fancier models, you need to code things from scratch anyways so you can choose whatever languabe you find most convenient.

With all that said, I would love to see these things formalized as solvers in differentialequations, but the path to get there isn’t obvious for these non-LQ control-problems.

[davidanthoff]

I worked with a student on that kind of problem a while ago (might actually have been that exact same model, I don’t remember the details), and we ended up solving them via what amounts to a relaxation method and then used JuMP to solve the set of non-linear equations that we got there. That worked really well.

[ptoche]

I have solved these problems the following way: Make the stationary state your initial value plus or minus a small deviation and run time backwards (from 0 to minus something big). This works a lot better because the saddle-path of the system is unstable (or in the language of economists: there is a unique converging path which you will attain only if you select a point that is exactly on the trajectory). To compute the path from below and from above you take deviations of opposite signs. Pick a point like (C_{ss}+\epsilon, K_{ss}+\epsilon) with \epsilon>0 to compute trajectories that start from “above” and \epsilon<0 for trajectories that start from below.

[rveltz]

Have you tried ApproxFun?

[ptoche]

Here’s an implementation of my cryptic comment above. There may be a better way to leverage the functionalities of DifferentialEquations, but it should get you started!

using DifferentialEquations
using Plots

# Parameters in a dictionary 
p = Dict(:A => 0.4, :α => 0.5, :δ => 0.1, :ρ => 0.1, :σ => 0.5)

"""Neoclassical Production Function. To plot, run the following:
   plot(x->F(x,p), 0, 2, color=:black, legend=false)"""
function F(x, p::Dict)
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]
    return A * x^α - δ * x
end

"""The Ramsey-Cass-Koompans system
   K: state variable, with given initial value K0
   C: control variable, free to jump"""
function ramsey!(du,u,p::Dict,t)

    # Parameters passed via dictionary
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]

    # Variables
    K, C = u

    # ODE system
    du[1] = (A * K^α - δ * K - C )
    du[2] = (α * A * K^(α - 1) - δ - ρ) * σ * C

    nothing
end

"""compute the stationary state"""
function ss(p::Dict)

    # Parameters passed via dictionary
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]

    # Steady-State
    Kss = (α * A / (δ + ρ))^(1/(1-α))
    Css = A * Kss^α - δ * Kss

    return Dict(:Kss => Kss, :Css => Css)
end

# The stationary state
SS = ss(p)  # returns 1.0 and 0.3 for reference parameter values
Kss, Css = ss(p)[:Kss], ss(p)[:Css]

# Run time "backwards" | start below the stationary state
tspan = (0.0, -50.0)
tweak = 0.98
u0 = [0.999999 * Kss; 0.999999 * tweak * Css]
ode = ODEProblem(ramsey!, u0, tspan, p)
sol1 = solve(ode, abstol=1e-8, reltol=1e-8)

# Policy function C = ϕ(K)

# initialize plot
plt = plot(title="Policy Function in Ramsey-Cass-Koopmans Model", xlims=(0,2*Kss), ylims=(0,2*Css), xlabel="K(t)", ylabel="C(t)", legend=false)

# plot nullclines
plot!(plt, x->F(x,p), 0, 2, color=:black, linewidth=3, legend=false)
vline!(plt, [Kss], color=:black, linewidth=3, legend=false)
plot!(plt, sol1, xflip=false, vars=(1,2), xlims=(0,2*Kss), ylims=(0,2*Css), color=:blue, linewidth=2)

# Run time "backwards" | start above the stationary state
tspan = (0.0, -100.0)
tweak = 1.01
u0 = [1.00001 * Kss; 1.00001 * tweak * Css]
ode = ODEProblem(ramsey!, u0, tspan, p)
sol2 = solve(ode, abstol=1e-8, reltol=1e-8)
plot!(plt, sol2, xflip=false, vars=(1,2), xlims=(0,2*Kss), ylims=(0,2*Css), color=:blue, linewidth=2)

# add decoration
scatter!((1.,0.3), label=["Stationary State" "."], color=:red)

savefig(plt, "ramsey-stable-arm.png")

[donk_fish]

If you’re unaware, there’s https://github.com/matthieugomez/EconPDEs.jl that might be more suitable

[donk_fish]

I’ve experienced around 2-3x speedup relative to matlab translating some of the codes on the Moll website, but I guess it’s arguable how relevant that really is.

But it could be relevant if you need to recalculate some model a massive amount of times, e.g. during estimation.

A bit off-topic maybe, sorry about that.

[MainBatDealer]

Wow! This is very comprehensive, thank you! I was about to try putting together an example of running time backward, but this is great!

I do have one quick question, though. I understand that A is a productivity parameter, but what is \sigma?

Based on the other comments I’ve gotten here, I may have to change up my strategy as my model gets more complicated (I’m planning on adding another state variable and several more control variables), but I’ll definitely try this first.

[ptoche]

@MainBatDealer

A is a productivity parameter, but what is σ ?

\sigma is the elasticity of intertemporal substitution. With time-separable preferences and isoelastic preferences, it happens to be equal to the inverse of the coefficient of relative risk aversion (CRRA). That is, \sigma is the inverse of: \frac{-C U^{\prime\prime}(C)}{U^{\prime}(C)}. These are constant for isoelastic utility functions. In your original code, you had the limit case \sigma=1, aka as “logarithmic utility”, that is your utility function (aka felicity function because it is the utility that accrues at one instant, not over the entire planning horizon):

U(C) = \frac{C^{1-1/\sigma}-1}{1-1/\sigma} \longrightarrow \log(C) \quad \text{as} \quad \sigma \rightarrow 1

Sometimes \sigma is denoted \gamma, sometimes \theta (notation’s a mess - Edit: Oh! wikipedia has \eta for the inverse). In models with risk, it’s common to focus on CRRA (=1/\sigma), but in deterministic models, it’s more common to define \sigma as the parameter of interest. It’s common to omit the minus one in the numerator. You’ll even occasionally see the extra lazy notation U(C)=C^{\phi}/\phi.

More along the lines of what your original intent was is this 20 year-old Matlab code I put together when I was doing my thesis. In words, it’s like a shooting algorithm. Here’s another piece of old Matlab code for a model with a habit stock. These may not run on modern Matlab anymore, but you can get the main idea and then write it in Julia!

[ptoche]

If your models are “representative agent” models that aren’t too different from the standard two-sector Romer/Rebelo/Aghion-Howitt/etc. models, then the approaches I have outlined here will work: My old Matlab code implements a shooting algorithm which will be very easy to adapt and works in most cases (you will need to tweak initial conditions and it will be harder to deal with problems with large curvature, i.e. high CRRA/ low \sigma). The Julia code I posted above with time running backwards leverages the fact that there’s a unique stable path and all the other paths are exploding. Both approaches are suitable to simulate one-agent problems (including central planner problems). They can handle a certain amount of risk like Brownian motion or simple Markov switches. However, once you move to multiple-agents, you have to give up the non-linearity of individual decisions to focus on the aggregation of multiple decisions. That is what Benjamin Moll’s HACT code mentioned above does. In other words, for multi-agent models you’ll need to change your approach entirely. And that’s where you really need Julia and all the Julia efficiency tricks, like keeping your parameters is structs, dealing with matrices intelligently, etc. So to sum up: for one-agent models with non-linearity the methods I outlined work well; for multi-agents you will have to sacrifice non-linearity at the individual level and focus on aggregation/distribution issues: for that follow the Moll-type approach (check the Aiyagari model). For the state of the art, look at Matthieu Gomez’s EconPDEs code: His repo has solved examples with multiple stocks. Update us on what you manage to put together!

[ptoche]

To understand the local dynamics of the system, it’s useful to draw a “fieldplot” (also known as “phase diagram” among economists). Here’s some code:

using DifferentialEquations
using Plots
using GeometryTypes # to make a grid

# Parameters in a dictionary 
p = Dict(:A => 0.4, :α => 0.5, :δ => 0.1, :ρ => 0.1, :σ => 0.5)

"""Neoclassical Production Function. To plot, run the following:
   plot(x->F(x,p), 0, 2, color=:black, legend=false)"""
function F(x, p::Dict)
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]
    return A * x^α - δ * x
end

"""The Ramsey-Cass-Koompans system
   K: state variable, with given initial value K0
   C: control variable, free to jump"""
function ramsey!(du,u,p::Dict,t)

    # Parameters passed via dictionary
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]

    # Variables
    K, C = u

    # ODE system
    du[1] = (A * K^α - δ * K - C )
    du[2] = (α * A * K^(α - 1) - δ - ρ) * σ * C

end

"""compute the stationary state"""
function ss(p::Dict)

    # Parameters passed via dictionary
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]

    # Steady-State
    Kss = (α * A / (δ + ρ))^(1/(1-α))
    Css = A * Kss^α - δ * Kss

    return Dict(:Kss => Kss, :Css => Css)
end

# event handler to prevent variables going out of bounds
function set_callback()
    c1(u, t, integrator) = u[1]  # stay in positive range
    c2(u, t, integrator) = u[2]  # stay in positive range
    c3(u, t, integrator) = 2*Kss - u[1]  # don't stray too far from stationary state
    c4(u, t, integrator) = 2*Css - u[2]  # don't stray too far from stationary state
    function terminate_affect!(integrator)
        terminate!(integrator)
    end
    # set up affect upon callback
    cb1 = ContinuousCallback(c1, terminate_affect!)
    cb2 = ContinuousCallback(c2, terminate_affect!)
    cb3 = ContinuousCallback(c3, terminate_affect!)
    cb4 = ContinuousCallback(c4, terminate_affect!)
    cb = CallbackSet(cb1, cb2, cb3, cb4)
    return cb
end

"""make a rectangular grid"""
make_grid(ranges::NTuple{N, <: AbstractRange}) where N = Point.(Iterators.product(ranges...))

"""fieldplot of the unstable saddle-point system
   prints a dot for every initial position where the solver fails"""
function fieldplot(f, U0, tspan, p)
    # 1. 
    # unpack parameter values
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]
    # 2.
    # compute stationary state
    Kss, Css = ss(p)[:Kss], ss(p)[:Css]
    # 3. 
    # Plot the nullclines
    plt = plot(title="Fieldplot of Ramsey-Cass-Koopmans Model", xlims=(0,2), ylims=(0,0.6), xlabel="K(t)", ylabel="C(t)", legend=false)
    plot!(plt, x->F(x,p), 0, 2, color=:black, linewidth=3, legend=false)
    vline!(plt, [Kss], color=:black, linewidth=3, legend=false)
    # 4. 
    # Set integration limits: stop integrating variables turn negative or become too large
    # conditions are C < 0, K < 0, C > 2*Css, K > 2*Kss
    cb = set_callback() 
    # 5.
    # draw trajectories that start from pair of values stored in U0
    # initialize the problem (may or may not be needed...)
    ode = ODEProblem(f, [0.5, 0.2], (0.0, 10.0), p)
    sol = solve(ode, callback=cb)
    """convenience function to solve the ODE system and  plot a trajectory for given initial value"""
    function pathplot(f, u0, tspan, p)
        # set up system
        ode = ODEProblem(f, u0, tspan, p)
        # solve system
        try
            sol = solve(ode, callback=cb, abstol=1e-8, reltol=1e-8)
        catch 
            print(".")
        end
        # add to existing plot
        plot!(plt, sol, vars=(1,2), xlims=(0,2), ylims=(0,0.6), xaxis="K(t)", yaxis="C(t)", xflip=false, yflip=false)
        return sol
    end
    # loop over all initial conditions stored in U0
    for (i, z) in enumerate(U0)
        # reset initial conditions
        u0 = [z[1]; z[2]]
        # plot forward path
        sol = pathplot(f, u0, tspan, p)
        # reverse time
        tspan2 = (tspan[1], -tspan[2])
        # plot backward path
        sol = pathplot(f, u0, tspan2, p)
        # save intermediate result for debugging
        # savefig(plt, "ramsey-fieldplot-debug-"*string(i)*".png")
    end
    # 6. 
    # return the plot
    return plt
end

# The stationary state
SS = ss(p)  # returns 1.0 and 0.3 for reference parameter values
Kss, Css = ss(p)[:Kss], ss(p)[:Css]

# make a grid of initial conditions
U0 = make_grid((0:(Kss/10):2*Kss, 0:(Css/10):2*Css))

# make a fieldplot from initial conditions
plt = fieldplot(ramsey!, U0, (0.0, 1.0), p)
display(plt)
savefig(plt, "ramsey-fieldplot.png")

This produces a bunch of

└ @ SciMLBase ~/.julia/packages/SciMLBase/BOcum/src/integrator_interface.jl:351 ....┌ Warning: Instability detected. Aborting

but the desired field plot still comes out. You can tell from this that only one trajectory will reach the stationary state (the trajectory drawn earlier in blue).

(as I can’t post more than 3 unanswered replies, I’m editing this one)

Or “out of the box” with CairoMakie

using CairoMakie

"""Ramsey model equations as points in (K,C) space"""
function ramseyEq(K, C, p)
    A, α, δ, ρ, σ = p[:A], p[:α], p[:δ], p[:ρ], p[:σ]
    dK = A * K^α - δ * K - C
    dC = (α * A * K^(α - 1) - δ - ρ) * σ * C
    return Point(dK, dC)
end
let
    fig = Figure(fontsize = 20)
    ax = Axis(fig, xlabel="K", ylabel="C", title="Streamplot of Ramsey-Cass-Koopmans Model", backgroundcolor = :black)
    streamplot!(ax, (x,y)->ramseyEq(x,y,p), 0..2*Kss, 0..2*Css, colormap = Reverse(:plasma),
        gridsize= (32,32), arrow_size = 10)
    fig[1, 1] = ax
    display(fig)
    save("ramsey-streamplot.png", fig)
end

[MainBatDealer]

Just a quick update: ApproxFun works really well for this problem. Here’s the code I used to solve it. It’s very similar to @ptoche’s response using DifferentialEquations.

using ApproxFun
using Plots

maxt = 100 # Time periods

x=Fun(identity, 0..maxt) # Let x() be an empty Chebychev fit on the range
                         # [0, 100]

#Parameters
ρ = 0.10
δ = 0.10
α = 0.5
A = 0.4
σ = 0.5
p = [ρ, δ, α, A, σ]


Y(K) = A *  K ^ α - δ * K # output 

#Steady-States
K_ss = ((δ + ρ) / (A * α)) ^ (1 / (α - 1)) 
C_ss = (A * K_ss ^ α - δ * K_ss)

#Starting from above steady-state
ramseyhigh = (K,C) -> [K' - (A * K^α - δ * K - C );
                        C' + C * σ * (ρ + δ - α * A * K ^ (α - 1));
                        K(0) - K_ss * 10.0;
                        K(maxt) - K_ss * 1.00001]

#Starting from below steady-state
ramseylow = (K,C) -> [K' - (A * K^α - δ * K - C );
                        C' + C * σ * (ρ + δ - α * A * K ^ (α - 1));
                        K(0) - K_ss * 0.01;
                        K(maxt) - K_ss * 0.99999]              

# Choosing a starting function for capital
# and consumption.
# The approximation would probably be 
# faster If I had a better starting function.
K10 = one(x) * K_ss
C10 = one(x) * C_ss

# Run the function approximation
Kh,Ch = newton(ramseyhigh, [K10,C10], maxiterations=50, tolerance = tolerance=1E-10)
Kl, Cl = newton(ramseylow, [K10,C10], maxiterations=50, tolerance = tolerance=1E-10)


#Recreating the plot from ptoche's post:

plt = plot(title="Policy Function in Ramsey-Cass-Koopmans Model", xlims=(0,2), ylims=(0,0.6), xlabel="K(t)", ylabel="C(t)", legend=false)
plot!(plt, x -> Y(x), 0, 2, color="black", linewidth=3, legend=false)
vline!(plt, [K_ss], color="black", linewidth=3, legend=false)
plot!(plt, Kl, Cl, xlims=(0,2), ylims=(0,0.6), color=:blue, linewidth=2)
scatter!((K_ss,C_ss), label=["Stationary State" "."], color=:red)
plot!(plt, Kh, Ch, xlims=(0,2), ylims=(0,0.6), color=:blue, linewidth=2)

savefig(plt, "ramsey-approx_fun.png")

[MainBatDealer]

I went ahead and marked the answers by @ptoche as solutions, since they do indeed solve my original question. I would also mark my answer using ApproxFun as a solution, but I’m not sure what the etiquette is about marking your own answers as solutions.

Based on the discussion here, there’s definitely more to be said about solving these kinds of problems, especially as more state/control variables are added to the mix. I’m especially interested in seeing the pros/cons of using ApproxFun vs. EconPDEs after reformulating into a value function problem.

[ptoche]

Feel free to change the “solution”, I don’t think there’s much to worry about since no “points” are awarded. It’s more about guiding first-time readers to a potentially quick answer, if I understand correctly, and to “close” the thread, as it were. I personally learned most from the ApproxFun solution since I’d never heard of it and had no idea it could solve such problems!

[Freya_the_Goddess]

Hi @ptoche and @MainBatDealer ,
can we calculate the time to solution and print it at the plot? So we can compare between ApproxFun and DiffEqs