Solving Boundary Value Differential Equation Problems for Economics

Specific Domains Finance and Economics
,
10

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

ramsey-stable-arm

4 Likes