Simple macroeconomic model in DynamicalSystems.jl

[josec]

Hi. I am very new to Julia.
I have a very simple macroeconomic model, which I managed to simulate and plot with numpy for Julia.
I tried to simulate it in DynamicalSystems.jl, but not all the equations are of the type: x[t] = f(x[t-1]). This is the model:
“”"
(1) y[t] = m * (a + G[t])
(2) c[t] = a + alpha * y[t-1]
(3) I[t] = beta * (y[t] - y[t-1])
“”"

Parameters: m = 1/(1-alpha), a = 50, alpha = 0.75, beta = 4, G = 100.
Steady state is y = 600.

My question: Is there a way to run such a model in DynamicalSystems? My problem is with equations 1 & 3.

Thanks in advance for your help.

[jlperla]

I think you will need to give the boundary/initial values to make the problem well formed.

If this is a macroeconomic/control problem then it probably has a transversality condition, which is a very particular type of boundary condition that probably isn’t handled in that package. But I know that @Albert_Zevelev was trying to collect packages that could handle those

[josec]

Thanks very much for the reply.
It is not a control problem or a rbc kind of model; it is an old fashioned keynesian model, with no microfoundations.
You are right I am sorry that I forgot the initial conditions: zeros for y, c, I, while G=100.
But my problem is that DynamicalSystems.jl requires the equations to have the form y[t+1] = f(x[t]), but equations 1 and 3 have variables indexed [t+1] on both sides of the “equal” sign.
Thanks again.

[jlperla]

Cool, if it is an initial value problem it is a lot easier.

I am still having trouble seeing this as a dynamical system though. Unless G[t] is an exogenously given function? But if so then I think that you can substitute to get almost everything in closed form to simulate on one dimension and this package is overkill?

If you misspecified this a little but still can’t substitute to get equations withiut the “differential” (ie, if some equations have terms all at the same time) then my suggestion is to convert this to continuous time. With that, you can write things as a differential algebraic equation with a mass matrix and solve it that way with differentialequations.jl

But as I said, I think either something is missing in your setup or else you can solve it with simpler methods

[Albert_Zevelev]

I can look at this tomorrow.

In this example gov spending, G[t] is an exogenously given sequence, which happens to be equal to 100 for all t.

I think it’s very easy to solve w/o any package.
First simulate y[t] given y0
Then: together c[t], I[t] given c0, I0

If you wanna use DynamicalSystems.jl for this, ask @Datseris

[Datseris]

I don’t understand the problem. Make a DiscreteDynamicalSystem with function

function f(u, p, t)
y, c, I = u
α, β, ... = p
yn = m*(a + G[t])
cn = a + α * y
In = β * (yn - y)
return SVector(yn, cn, In)
end

I don’t know what G is. You say its 100, but then how can you access it with an index…?

[josec]

Thank you so much for the input.
In fact G[t] is an exogenous variable: it is zero at the beginning, and becomes 100 in period 15.
y[t] is national income, c[t] is consumption, I[t] is investment, G[t] is government spending, an exogenous variable.
I apologize for the omission of the details, I guess I was just trying to avoid too much detail for the people to read.
Thanks, and again I am sorry for the omissions.

[josec]

@Datseris

function f(u, p, t)
y, c, I = u
α, β, ... = p
yn = m*(a + G[t])
cn = a + α * y
In = β * (yn - y)
return SVector(yn, cn, In)
end

EDIT:
The first equation y[t] = m*(a+G[t]) is actually the steady state solution (I thought I should put that one). But the original model should have been:
y[t] = c[t] + I[t] + G[t]
c[t] = a + alpha*(y[t-1])
I[t] = beta*(y[t] - y[t-1])

Initial value for y, c, I is zero.
G[t] is zero up until period 14. In period 15 it becomes 100.
Parameters: a = 50, alpha = 0.75, beta = 4.

I guess this should have been my first message. I solved this in Python, and in Julia with numpy, but thought perhaps there was a more direct way of doing it with DynamicalSystems (but the indexing in equations 1 and 3 got me confused).

As [jlperla] indicated previously this could easily be transformed into a continous system that can be analyzed with DifferentialEquations.jl, but I just wanted to get a better understanding of discrete systems in Julia.

I really appreciate your attention and help.

[jlperla]

@Datseris can dynamical systems handle “difference algebraic equations” or whatever they are called? I think that is what he has written down? Maybe this exact case can be transformed by substitution into a standard first order difference equation, I am not sure if the “algebraic” things are standard algorithms or even should be.

@josec when I think of the dynamical systems package, I think of thinks like chaos, cycles, bifurcation analysis, etc. If those are the things you expect to emerge out of your setup then cool. Not an expert here, but I think that because you are specifying a linear problem, it is impossible to have chaos/cycles/etc. If those can’t occur, this package may be a sledgehammer in this specific case since you can manually iterate forward.

[RJDennis]

I don’t understand why this system cannot be simulated using a for loop.

[Albert_Zevelev]

This one clearly can.
I think he wants to know how solve more general systems of difference equations in Julia, using this simple system as an example

[jlperla]

I think the problem here is that it has an algebraic parts of it. Without those, there wouldn’t be a problem. In continuous time it would be a DAE - of which the SciML ecosystem has plenty of algorithms.

[josec]

Thanks for your help. I tried to reorganize the code you sent in the following way. I would like to plot y,c, and print their values.

function f(u, p, t)
u = y, c, I, G
p = 50, 0.75, 4 # vector with the parameters.
u0 = 0, 0, 0, 0 # the initial conditions.

for loop for the exogenous government spending G (it is zero initially and becomes 100 in period 15).

Total number of periods is 200.

for n = 1:14
G[n] = 0
end
for n = 15:26
G[n] = 100
end
yn = cn + In + Gn
cn = 50 + 0.75 * y
In = 4 * (yn - y)
return SVector(yn, cn, In, Gn)
end

Thanks again for your advice.

Jose A. Cordero

[josec]

Yes, exactly, that’s what I would like to do.
On this path, I would eventually like to solve a system like the one in Samuelson multiplier-accelerator model.

Thank you very much.

[josec]

@jlperla you are right. This is a linear system and definitely cannot generate cycles or chaotic behavior. I am probably trying to to this with the wrong package. I didn’t think about that.

Thanks!!!

[josec]

@Datseris I apologize for the careless use of the pound sign within markdown (didn`t realize this was in markdown format).

[Albert_Zevelev]

System:
Y_t = C_{t} + I_t + G_t
C_t = a + \alpha Y_{t-1}
I_t = \beta (Y_{t} - Y_{t-1})
G_t = 0\times 1\{t\leq 14\} + 100\times 1\{t\geq 15\}
Parameters: a=50, \alpha =.75, \beta =4
Initial Condition: Y_{0}

System of 4 equations in 4 variables Y_{t}, C_{t}, I_{t}, G_{t}
A solution gives each variable as a function of t.
G_{t} already given (exogenously).

Solution 1 (Closed form)
Plug the 2nd/3rd equation into the 1st:
Y_t = (a + \alpha Y_{t-1}) + (\beta Y_{t} - \beta Y_{t-1}) + G_t
\Leftrightarrow (1-\beta)Y_t = a + (\alpha -\beta) Y_{t-1} + G_t
\Leftrightarrow Y_t = \frac{a}{(1-\beta)} + \frac{(\alpha -\beta)}{(1-\beta)} Y_{t-1} + \frac{1}{(1-\beta)}G_t , Y_0 given is a 1st order difference equation, solve for Y_{t}.
Next, plug Y_{t} into the 2nd/3rd equations to get C_{t}, I_{t}
Hopefully Symbolics.jl will one day solve systems of difference equations @ChrisRackauckas.
For now Mathematica:

Solution 2 (numerical in Julia)

using Plots;
#
a=50.0; α=.75; β=4.0; Y0=0.0;
G(t)= (t <= 14) ? 0.0 : 100.0 
Y_rec(Ym, t; α=α, β=β, a=a) = (a/(1-β)) + ((α-β)/(1-β))*Ym + (1/(1-β))*G(t)
I_rec(Y, Ym; β=β) = β*(Y-Ym)
C_rec(Ym; a=a, α=α) = a + α*Ym
# Simulate
Δt = 1.0; T = 250; time = 0.0:Δt:T
G_sim,    Y_sim,    I_sim,    C_sim    = [zeros(length(time),1) for i in 1:4]
Y_sim[1] = Y0 
#
for i in 1:(length(time)-1)
    G_sim[i+1] =      G(i)
    Y_sim[i+1] =      Y_rec(Y_sim[i], i)
    I_sim[i+1] =      I_rec(Y_sim[i+1], Y_sim[i])
    C_sim[i+1] =      C_rec(Y_sim[i])
end
#
ix = 2:21 #(100) # time=0, to time=23
p1 = plot(size=(400,400),legend=:bottomleft)
plot!(time[ix], G_sim[ix], lab = "G")
plot!(time[ix], Y_sim[ix], lab = "Y")
plot!(time[ix], I_sim[ix], lab = "I")
plot!(time[ix], C_sim[ix], lab = "C")

If we change one parameter \beta =0.4 instead of \beta=4, things look more as expected:

As we lengthen the time horizon, variables approach Steady State:

@josec Based on the simulation, I believe there was a typo somewhere in your description of the problem (or parameter value, or initial value). Can you please post the original source?

[josec]

@Albert_Zevelev: I appreciate very much your help; extremely useful.
Just three comments:

• This example, I made it up, so I am sorry that I cant provide the source because there isnt one.
• The steady state solution (Y_rec = Ym, in your notation) is:
  Y = (1/1-alpha) * (a + G). This means that the value of beta should not affect the steady state solution. I thus wonder why, with the code, changing the value of beta changes the results.
• In the simulation on the first graph, Y, C, and I seem to move down to zero, while the closed solution indicates that, in steady state, Y=600, C=500 and I should be zero from the investment equation.

Thanks again for your help, and I apologize for asking so many questions.

[josec]

@Albert_Zevelev: Hi again. I made some slight changes and got the correct steady state results, even with beta =4. This is what I did:

1. Eliminated the equation with Y_rec(Ym, t; α=α, β=β, a=a) = …
2. After the equations for I_rec and C_rec I included this equation: Y_rec(t)= C_rec + I_rec + G(t)
3. The new equation for Y_rec is an equilibrium condition. The original equation for Y_rec provided the “short-run” equilibrium for Y.
4. After the above change, the plot provided the correct steady state values for Y, C, I, G; approaching the steady state with more periods.

Thanks again for your help!!!

[Albert_Zevelev]

Before I forget.

We start w/ a system of 3 DAEs in 3 variables Y_{t}, C_{t}, I_{t}. G_{t}=100\times 1\{t\geq 15\} is exogenous.
Y_t = C_{t} + I_t + G_t
C_t = a + \alpha Y_{t-1}
I_t = \beta (Y_{t} - Y_{t-1})
Parameters: (a, \alpha, \beta)
Initial Condition: Y_{0}

We reduce this system of 3 DAEs into 3 1^{st} order non-autonomous difference equations:
Y_t =F_{1}(Y_{t-1}, t)= \frac{a}{(1-\beta)} + \frac{(\alpha -\beta)}{(1-\beta)} Y_{t-1} + \frac{1}{(1-\beta)}G_t , Y_0 given
C_t =F_{2}(Y_{t-1}, t) = a + \alpha Y_{t-1}
I_t =F_{3}(Y_{t-1}, t)= \beta (F_{1}(Y_{t-1}, t) - Y_{t-1})
Note: the reduced system only depends on Y_{t-1} so we only need an initial condition for Y_t.

SS:
Y_{ss} = \frac{a}{(1-\beta)} + \frac{(\alpha -\beta)}{(1-\beta)} Y_{ss} + \frac{100}{(1-\beta)} \Leftrightarrow \frac{(1-\alpha)}{(1-\beta)} Y_{ss} = \frac{a}{(1-\beta)} + \frac{100}{(1-\beta)} \Leftrightarrow Y_{ss} = \frac{a}{(1-\alpha)} + \frac{100}{(1-\alpha)}
C_{ss} = a + \alpha Y_{ss} = a + \frac{a \alpha}{(1-\alpha)} + \frac{100 \alpha}{(1-\alpha)} =\frac{a}{(1-\alpha)} + \frac{100 \alpha}{(1-\alpha)}
I_{ss} = \beta (Y_{ss} - Y_{ss}) = 0
G_{ss} = 100
Note: the system acts weird for \beta > \alpha.