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
With an initial capital stock K(0) = K_0 and an end-of-time condition
\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.
BoundaryValueDiffEq will not appreciate a boundary condition at t = \infty, I made the following transformation:
which transforms my system to
and my terminal condition to
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
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 δ = p α = p K = u C = u du = (K^α - δ * K - C ) / (1 - t) ^ 2 du = - C * (ρ + δ - α * K ^ (α - 1)) / (1 - t) ^ 2 end # Define the boundary conditions function bc!(residual, u, t) residual = u - K0 residual = u[end] - C_ss end # Prevent the solver from trying # negative values of capital and consumption function g(u, resid) resid = u - 0.001 resid = u - 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.