SOLVING DELAY DIFFERENTIAL EQUATIONS.
SOME EXAMPLES FROM ECONOMIC DYNAMICS

Henryk Zawadzki

University of Economics in Katowice, Poland

Delay differential equations (DDEs)

Some DDEs (or time - delay differential equations) appeared in XVIII century in works of Leonhard Euler, an eminent Swiss mathematician and physicist, who spent most of his adult life in St. Petersburg, Russia, and in Berlin, Prussia.

Leonhard Euler (1707 – 1783)

Delay differential equations) appeared in economic scientific papers in early thirties of XX century among others in works of R. Frisch, M. Kalecki and (slightly later) in works of J. Tinbergen and R.M. Goodwin.

R. Frisch                  M. Kalecki                J. Tinbergen                   R. M. Goodwin

There are different kinds of DDs. I will focus on just one kind, namely equations of the form

x'(t)=F[x(t),x(t-θ)]

where x = x(t)  is the function we are looking for and  θ > 0 is a constant called a delay. In order to solve such an equation we also need to specify an initial function ϕ which gives a behaviour of  x  prior to time 0 (assuming that we start at t = 0).

x(t)=ϕ(t)  for t∈[-θ,0]

Kalecki’s business cycle model

The model consists of four basic equations

L(t)=I(t-θ),
A(t)=(1/θ)(∫_(t-θ)^t)I(τ)dτ,
K'(t)=L(t)-U,
I(t)=m[(Subscript[C, 1])+A(t)]-n K(t),

where L denotes deliveries of finished capital goods, I denotes investment orders, A is a production of capital goods and K is a capital stock. (Subscript[C, 1]), U, m, n and θ are parameters. In particular, parameter θ >0 is an average lag of investment (gestation lag) - the time interval between the decision to invest and the delivery of the finished capital goods.
Denoting by

J(t)=I(t)-U

the net investment, after some transformations we can obtain a linear delay differential equation (let us call it Kalecki’s equation)

J'(t)=(m/θ)J(t)-((m/θ)+n)J(t-θ)

which describes the dynamics of the variable J, or more exactly - deviations of J from some equilibrium level J* . In order to solve this equation we have also define initial function ϕ describing the behaviour of J for t ∈ [- θ,0].

Numerical solution of Kalecki’s equation

Clear[m,n,θ]

m:=0.05;n:=1;θ:=1;

Let the initial function be of the form

ϕ[t_]:=1+t

NetInvestment= NDSolve[{J'[t] ==(m/θ)*J[t]-((m/θ)+n)*J[t-θ], J[t /; t≤ 0]==ϕ[t]},J,{t,-θ,15}]

Plot1=Plot[Evaluate[J[t] /. First[NetInvestment]],{t,-θ,15}, PlotRange->All,AxesLabel->{"t","J(t)"}]

Evaluate[J[t]/.First[NetInvestment]/.{t->5}]

Clear[m];Clear[n];Clear[θ];Clear[ϕ];

Manipulate[
Module[{NetInvestment = NDSolve[{J'[t]==(m/θ)*J[t]-((m/θ)+n)*J[t-θ],J[t/;t≤0]==ϕ},J,{t,-θ,15}]},
Plot[Evaluate[J[t] /. First[NetInvestment]],{t,-θ,15}]],{ϕ,{1,t+1,Exp[t], Cos[t]}},{m,0,1,0.05},{n,0,5,1},{θ,0.5,5,0.5}]

Exact “step by step” solution

Let us copy,paste and slightly modify some programs coming from one of Mathematica tutorials.

DiscontinuityTree[t0_,  Tend_, delays_] :=
Module[{dt, next, ord},
ord[t_] := Infinity;
ord[t0] = 0;
next[b_, order_, del_] := Map[dt[b, #, order , del]&, del];
dt[t_ , {d_, nq_},order_, del_] :=Module[{b= t + d},
If[b ≤ Tend,
o  = order + Boole[!nq];
ord[b] = Min[ord[b], o];
Sow[{t->b, d}];
next[b,o, del]]];
rules = Reap[next[t0,0,delays]][[2,1]];
rules = Tally[rules][[All,1]];
f[J_?NumericQ] := {J,ord[J]};
f[a_->b_]:= f[a]->f[b];
rules[[All,1]] = Map[f,rules[[All,1]]];
rules]

IntegrateSmooth[rhs_, history_, delayvars_, pfun_, dvars_, {t_, t0_, t1_}] := Module[{delayvals, dvt, tau, hrule,dvrule, dvrules, oderhs, ode, init, sol},
dvt[tau_] = Map[#[tau]&,dvars];
hrule[pos_] := Thread[dvars-> Map[Function[Evaluate[{t}],#]&, history[[pos]]]];
dvrule[(dv_)[z_]]:= Module[{delay,pos},
delay = t - z;
pos = pfun[t0-delay];
dv[z] -> (dv[z] /. hrule[pos])];
dvrules = Map[dvrule, delayvars];
oderhs = rhs /. dvrules;
ode = Thread[D[dvt[t],t] == oderhs];
init =Thread[dvt[t0]==(dvt[t0] /.hrule[-1])];
sol = DSolve[{ode, init}, dvars, t];
If[Head[sol]===DSolve || Length[sol]==0,
Message[DDESteps::stuck, ode, init];
Throw[$Failed]];
dvt[t] /. First[sol]
];
DDESteps::stuck="DSolve was not able to find a solution for `1` with initial conditions `2`.";

DDESteps[rhsin_, phin_, dvarsin_,{t_, tinit_, tend_}] := Module[{rhs = Listify[rhsin],phi = Listify[phin], dvars = Listify[dvarsin],history, delayvars, delays, dtree,intervals, p,pfun, next,pieces,hfuns},
history = {phi};
delayvars = Cases[rhs, (v:((dv_[z_] | Derivative[1][dv_][z_]) /; (MemberQ[dvars, dv] && UnsameQ[z, t])))->{v,t - z}, Infinity];
{delayvars, delays} = Map[Union, Transpose[delayvars]];
dtree = DiscontinuityTree[tinit, tend,Map[{#, True}&, delays]];
dtree = Union[Flatten[{tinit, tend, dtree[[All, 1,2,1]]}]];
dtree = SortBy[dtree, N];
intervals = Partition[dtree, 2,1];
p= 2;
pfun = Join[{{1,t<tinit}},Apply[Function[{p++,#1≤ t<#2}], intervals,{1}]];
pfun = Function[Evaluate[{t}], Evaluate[Piecewise[pfun, p]]];
Catch[Do[
next = IntegrateSmooth[rhs, history, delayvars, pfun, dvars, Prepend[interval, t]];
history = Append[history, next],
{interval, intervals}]];
pieces =Flatten[{t < tinit, Apply[(#1 ≤ t<#2)&, Drop[intervals, -1],{1}], Apply[(#1≤t≤#2)&, Last[intervals]]}];
pieces = Take[pieces, Length[history]];
hfuns = Map[Function[Evaluate[{t}],Evaluate[Piecewise[Transpose[{#, pieces}], Indeterminate]]]&,Transpose[history]];
Thread[dvars->hfuns]
];
Listify[J_List] := J;
Listify[J_] := {J};

sol = DDESteps[0.05*J[t]-1.05*J[t-1],1+t,J,{t,0,5}]

Evaluate[J[t]/.sol/.{t->5}]

Plot[Evaluate[J[t] /. sol],{t,-1,5}, AxesLabel->{"t","J[t]"}]

The neoclassical growth model of Solow and Swan

R.M. Solow                                    T.W. Swan

The fundamental equation of the Solow-Swan growth model has a form

k'(t)=s f[k(t)]-(ν+λ)k(t)

where k = K/L is the capital-labor ratio, s ∈ (0, 1) is the saving rate, ν > 0  is the population grow rate,  λ > 0 is the depreciacion rate of physical capital and  f is a neoclassical production function (i.e. continuous, increasing and strictly concave in variable k). When f is neoclassical and some additional conditions (Ikeda conditions) are fullfilled, the dynamics of Solow-Swan model (6) is very simple. For every initial conditon

k(0)=(Subscript[k, 0]) >0

it has a solution k = k(t) which tends monotonically to some stable equilibrium (time independent solution)

k(t)≡(SuperStar[k])>0

When f is the Cobb-Douglas function

f(k)=A (k^α)

where A > 0 is the level of technology and α ∈ (0, 1) is a constant, we have a solution

k(t)=(([((s A)/(ν+λ))+c (e^(-(1- α)(ν + λ)t))])^(1/(1-α)))

where c is a constant depending on initial condition (7). The above solution tends to

k(t)≡(SuperStar[k])=(((((s A)/(ν+λ))))^(1/(1-α)))

A:=1;α:=0.5;ν:=0.2;λ:=0.2;s:=0.5;c:=-0.25;

Plot[{(((((s A)/(ν+λ))+c Exp[-(1-α) (ν+λ) t]))^(1/(1-α))),(((((s A)/(ν+λ))))^(1/(1-α)))},{t,0,100},AxesLabel->{"t","k(t)"},PlotRange->{0,2}]

The above figure shows two solutions of equation (6): a solution with initial condition k(0) = 1 (blue line) and stationary solution k(t) = k*=1.5625.

When the rate of change of k at the moment t is a function of k at t - θ, we obtain the following nonlinear delay differential equation (delay Solow-Swan model)

k'(t)=s f[k(t-θ)]-(ν+λ)k(t-θ)

The dynamics of this model is more complex and depends on the values of parameters (first of all the value of delay θ) and initial function. In addition to trajectories similar to those shown at the picture above, we can also observe cycles and (as proved in [Krawiec, Szydłowski (2009) even stable cycles.

Numerical solutions of delay Solow-Swan model

Let us assume, that the production function is the Cobb-Douglas function given by formula (9).

Clear[A,α,ν,λ,s,θ]

f[k_]:=A*(k^α)

A:=1;α:=0.5;ν:=0.2;λ:=0.2;s:=0.5;θ:=5

Capital= NDSolve[{k'[t] ==s*f[k[t-θ]]-(ν+λ)*k[t-θ], k[t /; t≤ 0]==Exp[t]},k,{t,-θ,100}]

Plot[Evaluate[k[t] /. First[Capital]],{t,-θ,100}, PlotRange->All,AxesLabel->{"t","k(t)"}]

Manipulate[
Module[{Capital = NDSolve[{k'[t] ==s*f[k[t-θ]]-(ν+λ)*k[t-θ],k[t/;t≤0]==ϕ},k,{t,-θ,100}]},
Plot[Evaluate[k[t] /. First[Capital]],{t,-θ,100},PlotRange->All]],{ϕ,{1,Exp[t]}},{α,0,1,0.05},{s,0,1,0.05},{ν,0,1,0.05},{λ,0,1,0.05},{θ,0.05,10,0.5}]

References

Barro R.J., Sala-i-Martin X. Economic growth. The MIT Press, Cambridge, Massachusetts,1995.

Gandolfo G. Economic Dynamics. Springer-Verlag, Berlin, 1996.

Krawiec A, Szydłowski M. Investment delay and economic system dynamics. International Workshop, Delayed Complex Systems, 5-9 October 2009, Dresden, Germany.