# Modeling dynamic economic systems using ModelingToolkit (MTK)

**URL:** <https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861>\
**Category:** Modelling & Simulations\
**Created:** [December 9, 2021, 9:22pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861 "2021-12-09T21:22:43Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 9, 2021, 9:22pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/1 "2021-12-09T21:22:43Z")

</div>

I am seeking help in specifying a dynamic economic model in MTK. Although I have followed the evolution of MTK since its creation, I have difficulties with “function” as a function can be the whole model, a function of parameters or a function of variables. My aim is to simulate the dynamic model written in MTK to generate synthetic timeseries and then take it across to DataDrivenDiffEq for model discovery.

The nonlinear model consists of a system of three first order ODEs:

\begin{aligned} \frac{\dot{ω}}{ω} & = \Big(\frac{η}{1+η}\Big) \Big[Φ(λ) −a\_l\Big] \ \ (1) \\ \frac{\dot{λ}}{λ} & = \Big[κ(π) C \Big(\frac{1 - ω}{b}\Big)^\frac{1}{η} −δ−a\_l−β - \frac{1}{η(1-ω)} \frac{\dot{ω}}{ω}\Big] \ \ (2) \\ \frac{\dot{d\_f}}{d\_f} & = \Big[r − g (κ(π), ω)\Big] + \frac{κ(π) −(1 − ω)}{d\_f} \ \ (3) \\ \end{aligned}

with four behavioural/variable functions π, Φ(λ) , κ(π) and g(κ(π), ω) of the form:

\begin{aligned} π & : = 1 - ω - r d\_f \\ Φ(λ) & = Φ₁ + Φ₂e^{(Φ₃λ)} \\ κ(π) & = κ₁ + κ₂e^{(κ₃π)} \\ g(κ(π), ω) & = κ(π) C \Big(\frac{1 - ω}{b}\Big)^\frac{1}{η} -δ - \frac{\dot{ω}}{(1 - ω)η} \\ \end{aligned}

the parameters are: t η κ C b δ β a Φ₁ Φ₂ Φ₃ κ₁ κ₂ κ₃ r

and the variables are: ω(t), λ(t), d\_f (t)

This nonlinear formulation is already relatively simple as it is specified in logarithmic derivative terms. Otherwise it is even more complicated and nonlinear. Help in defining and registering the various types of functions would be much appreciated.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 10, 2021, 11:40am UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/2 "2021-12-10T11:40:23Z")

</div>

I don’t understand the question. What’s your current code?

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 10, 2021, 1:53pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/3 "2021-12-10T13:53:53Z")

</div>

using DataDrivenDiffEq  
using LinearAlgebra  
using ModelingToolkit  
using OrdinaryDiffEq  
using Plots  
using Random  
using Symbolics: scalarize  
using Latexify

Random.seed!(1111) # For the noise

@parameters t η κ C b δ β a π Φ₁ Φ₂ Φ₃ κ₁ κ₂ κ₃ r  
@variables ω(t) λ(t) d\_f(t)  
D = Differential(t)

function π(ω, r, d\_f) # Accounting identity  
1 - ω - r\*d\_f  
end  
@register π

function Φ(λ, Φ₁, Φ₂, Φ₃) # Phillip’s curve  
Φ₁ + Φ₂_exp(Φ₃_λ)  
end  
@register Φ(λ)

function κ(π, κ₁, κ₂, κ₃) # Investment function  
κ₁ + κ₂_exp(κ₃_π)  
end  
@register κ(π)

function g(κ, ω) # Growth rate of the economy  
κ_C_((1-ω)/b)^(1/η) - δ - D(ω)/((1 - ω)\*η)  
end  
@ register g(κ, ω)

# Define a nonlinear system

sfc = [D(ω) ~ ω\*(η/(1+η))_(Φ(λ) - a),  
D(λ) ~ λ_(κ(π)_C_((1-ω)/b)^(1/η) - δ -β -a -(1/(η\*(1-ω)))_D(ω)/ω),  
D(d\_f) ~ d\_f_(r - g(κ, ω)) + κ(π) - (1 - ω)]

latexify(sfc)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 16, 2021, 12:34pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/4 "2021-12-16T12:34:27Z")

</div>

That seems like it would just work if you didn’t register the functions? You want to register as little as possible so it can define symbolic derivatives by tracing.

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 16, 2021, 2:06pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/5 "2021-12-16T14:06:13Z")

</div>

For now, I get the right math for the model with this:

@variables t ω(t) λ(t) Φ(t) κ(t) π(t) d\_f(t) g(t) # independent and dependent variables  
@parameters Φ₁ Φ₂ Φ₃ η a\_l κ₁ κ₂ κ₃ δ β r C b # parameters  
D = Differential(t) # define an operator for the differentiation w.r.t. time

sfc = [ Φ ~ Φ₁ + Φ₂_exp(Φ₃_λ),  
κ ~ κ₁ + κ₂_exp(κ₃_π),  
π ~ (1 - ω - r_d\_f),  
g ~ κ_C\*((1-ω)/b)^(1/η) - δ - D(ω)/((1 - ω)_η),  
D(ω) ~ ω_((η/(1+η))_(Φ - a\_l)),  
D(λ) ~ λ_((κ_C_((1-ω)/b)^(1/η) - δ -β -a\_l -(1/(η\*(1-ω)))_D(ω)/ω)),  
D(d\_f) ~ d\_f_(r - g) + κ - (1 - ω)]

latexify(sfc)

The next step is to see how “NonlinearSystem” and “NonlinearProblem” will accept this sfc model. If it solves as it is I do not need functions at all and this reads as textbook presentation.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 16, 2021, 3:02pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/6 "2021-12-16T15:02:42Z")

</div>

That’s not a nonlinear system. It’s an ODESystem representing a DAE, which you can solve with `DAEProblem`’s constructor.

---

<div class="post-metadata">

**Author:** ![ptoche](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ptoche/32/23554_2.png) [@ptoche](https://discourse.julialang.org/u/ptoche)\
**Post date:** [December 16, 2021, 4:16pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/7 "2021-12-16T16:16:12Z")

</div>

As this is economics, isn’t this a system derived from an optimization problem with state, control and co-state variable? And the derivatives are with respect to time? And you have initial values and a terminal/limit condition? And do you know the nature of this system, e.g. saddle-point stable?

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 16, 2021, 9:01pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/8 "2021-12-16T21:01:29Z")

</div>

I am just lost in the process of setting up the DAEProblem. The docs and test of the Robertson DAE problem is too cryptic for me. Any other DAEProblem examples using ModelingToolkit anywhere?

This is where I stand:  
@variables t ω(t) λ(t) Φ(t) κ(t) π(t) d\_f(t) g(t) # independent and dependent variables  
@parameters Φ₁ Φ₂ Φ₃ η a\_l κ₁ κ₂ κ₃ δ β r C b # parameters  
D = Differential(t) # define an operator for the differentiation w.r.t. time

sfc = [ res1 = Φ₁ + Φ₂_exp(Φ₃_λ) - Φ,  
res2 = κ₁ + κ₂_exp(κ₃_π) - κ,  
res3 = (1 - ω - r_d\_f) - π,  
res4 = κ_C\*((1-ω)/b)^(1/η) - δ - D(ω)/((1 - ω)_η) - g,  
res5 = ω_((η/(1+η))_(Φ - a\_l)) - D(ω),  
res6 = λ_((κ_C_((1-ω)/b)^(1/η) - δ -β -a\_l -(1/(η\*(1-ω)))_D(ω)/ω)) - D(λ),  
res7 = d\_f_(r - g) + κ - (1 - ω) - D(d\_f)]

latexify(sfc)

and then  
@named keen = NonlinearSystem(sfc, t, [ω, λ, d\_f, Φ, κ, π, g], [Φ₁, Φ₂, Φ₃, η, a\_l, κ₁, κ₂, κ₃, δ, β, r, C, b])

u0 = [ω =\> 0.75,  
λ =\> 0.75,  
d\_f =\> 0.80]  
du0 = []

param = [Φ₁ =\> -0.01, Φ₂ =\> , Φ₃ =\> ,  
η =\> , a\_l =\> 2.0 ,  
κ₁ =\> , κ₂ =\> , κ₃ =\> ,  
δ =\> 1, β =\> 1, r =\> 4.0, C =\> 0.333, b =\> 0.135]

prob = DAEProblem(keen,du0,u0,(0.0,50.0), param)  
sol = solve(prob,IDA())

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 16, 2021, 9:07pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/9 "2021-12-16T21:07:54Z")

</div>

Not yet. For now it is just a stock-flow consistent model of three nonlinear ODEs that can exhibit chaotic behaviour for some values of parameters and initial values. Will get fancier when this basic setup works.

---

<div class="post-metadata">

**Author:** ![tlorans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlorans/32/19849_2.png) [@tlorans](https://discourse.julialang.org/u/tlorans)\
**Post date:** [December 16, 2021, 9:17pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/10 "2021-12-16T21:17:41Z")

</div>

Nice to see a Steve Keen fan here. Are you considering to develop an SFC package in Julia using MLT or just exploring?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 16, 2021, 10:25pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/11 "2021-12-16T22:25:25Z")

</div>

See the DAE benchmarks for examples, like [Chemical Akzo Nobel Differential-Algebraic Equation (DAE) Work-Precision Diagrams](https://benchmarks.sciml.ai/html/DAE/ChemicalAkzoNobel.html)

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 17, 2021, 8:15am UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/12 "2021-12-17T08:15:58Z")

</div>

Great!

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 23, 2021, 1:42pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/13 "2021-12-23T13:42:15Z")

</div>

Let me give the solution to the MTK-DAE problem for the interested readers. This could be simplified further with a D ln x (t) operator as mentioned in [https://github.com/SciML/ModelingToolkit.jl/issues/1344](https://github.com/SciML/ModelingToolkit.jl/issues/1344)

using DataDrivenDiffEq  
using LinearAlgebra  
using ModelingToolkit  
using DifferentialEquations: solve  
using NonlinearSolve  
using OrdinaryDiffEq  
using Plots: plot  
using Random  
using Symbolics: scalarize  
using Latexify

ModelingToolkit.@parameters begin  
Φ₁=0.04/(1. - 0.04^2)  
Φ₂=0.04^3/(1. - 0.04^2)  
η=500  
a₁=0.02  
κ₁=0.  
κ₂=0.05  
κ₃=1.75  
δ=0.05  
β=0.01  
r=0.03  
C=0.33  
b=0.135  
end

@variables begin  
t  
ω(t) = 0.75  
λ(t) = 0.90  
d₁(t) = 0.5  
π(t) = 0.14  
Y(t) = 100.0  
Φ(t) # = 0.95  
κ(t) # = 0.05  
g(t) # = 0.05  
end

D = Differential(t)

Φ = -Φ₁ + Φ₂/((1. - λ)^2)  
κ = κ₁ + κ₂_exp(κ₃_π)  
g = κ_C_((1 - ω)/b)^(1/η) - δ - ω/((1 - ω)_(1 + η))_(Φ - a₁)

sfc = [  
D(ω) ~ ω\*((η/(1 + η))_(Φ - a₁))  
D(λ) ~ λ_((κ_C_((1 - ω)/b)^(1/η) - δ -β -a₁ -(1/((1 + η)_(1 - ω)))_(Φ - a₁)))  
D(d₁) ~ d₁\*(r - g) + κ - (1 - ω)  
D(Y) ~ Y_g  
0. ~ π - 1 + ω + r_d₁  
]

latexify(sfc)

ModelingToolkit.@named keen = ModelingToolkit.ODESystem(sfc)

tspan = (0.0, 200.0)  
mmprob = ODEProblem(keen, , tspan)  
sol = solve(mmprob, Rodas4(),abstol=1/10^14,reltol=1/10^14);

du = mmprob.f(mmprob.u0,mmprob.p,0.0)  
du0 = D.(states(keen)) .=\> du  
daeprob = DAEProblem(keen,du0,,tspan)  
ref\_sol = solve(daeprob,IDA(),abstol=1/10^14,reltol=1/10^14);

#probs = [mmprob,daeprob]  
#refs = [ref\_sol,ref\_sol];

---

<div class="post-metadata">

**Author:** ![finmod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/finmod/32/26597_2.png) [@finmod](https://discourse.julialang.org/u/finmod)\
**Post date:** [December 23, 2021, 2:03pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/14 "2021-12-23T14:03:20Z")

</div>

The idea is to replicate in Julia using MTK the SFC models in Giraud and Grasselli models ([https://ms.mcmaster.ca/~grasselli/GiraudGrasselli2021.pdf](https://ms.mcmaster.ca/~grasselli/GiraudGrasselli2021.pdf)) and in Bastidas et al ([https://www.parisschoolofeconomics.eu/docs/fabre-adrien/minskyan-classical-growth-cycles--bastidas,-fabre,-mcisaac-2018.pdf](https://www.parisschoolofeconomics.eu/docs/fabre-adrien/minskyan-classical-growth-cycles--bastidas,-fabre,-mcisaac-2018.pdf)).

My roadmap is to show that continuous time dynamical systems in economics that have been around since the 1970s can be discovered and identified from the dataset using DataDrivenDiffEq. This is a two step process: 1) SFC model producing synthetic data =\> data discovery of SFC model and 2) use real data for the same variables to discover which model is compatible with it.

---

<div class="post-metadata">

**Author:** ![tlorans](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlorans/32/19849_2.png) [@tlorans](https://discourse.julialang.org/u/tlorans)\
**Post date:** [December 23, 2021, 9:38pm UTC](https://discourse.julialang.org/t/modeling-dynamic-economic-systems-using-modelingtoolkit-mtk/72861/15 "2021-12-23T21:38:56Z")

</div>

Sounds nice.  
Gaël Giraud is a friend’s thesis supervisor, let me know if you want to be in touch.