Hi every one, I am really new to Julia and I would like to evaluate a dynamic system like this one:
Y(t) = C(t) + I(t) + Z
C(t) = c*Y(t)
I(t) = h*Y(t)
Z = 100; c=0.7; h = 0.2; T = 100
There is a simple way (without for loops) to evaluate every simple point like t=1:T? If so, does the order of the equation matter?
This could be setup as a discrete system and iterated through one of the fancy packages out there, but why the aversion to loops? The loop code will probably be the most clear for this exact problem.
You can use dot syntax to evaluate a function for many arguments f.(1:T) (which is more compact than loops, but not any faster — loops are fast too in Julia). But you have to write a function f first that actually solves your 3x3 system of equations in (Y,C,I).
But maybe I’m misunderstanding you. If you’re thinking more of an iterated functions system, then the easiest syntax is probably a loop.
(It would be clearer to understand you if you actually posted working code; the code you posted is an infinite-recursion loop and hence does not work.)
Note that putting your parameters c and h as non-const global variables will slow things down by a lot in Julia.
Thank you both for the quick responses.
What I am afraid about loops is both performance and that I am not sure if I need to specify the equations in a proper order (this is may not be the case with a package).
I do not have any working code with Julia (started this week 
, yey! ). The example that I wrote was based on a python package that I use for simulate my models and I intend to do a similar thing in Julia.
This is the package that I mentioned: pysolve3
Loops are fast in Julia.
I don’t understand the role of t in your equations.
Should the left-hand side be e.g. Y(t+1)?
t only indicates that this variable is evaluated at period t.
In this (really) short example, it Y(t) and not Y(t+1). Analytically, what I have is:
Y(t) = Z/(1-c-h)
and Z describes Y’s dynamics and could be
Z(t+1) = (1+g)*Z(t)
in which g is an exogenous variable.
The model that I have in mind is much bigger than this one and I can provide it to make myself clearer
(Thank you all for the support 
)
What I have in mind is something similar to this:
using Calculus
using Plots
# Functions declaration
# Cobb-Douglas production function
# Default beta = 1-alpha gives constant return to scale (CRS)
CobbDouglas(K, AL, alpha, beta = 1-alpha) = K^alpha * AL^beta
# Production function to be used in model
Production(K, A, L) = CobbDouglas(K, A*L, 0.7)
# Capital in (t+1) as function of Capital in (t), Investment in (t)
# and depreciation rate delta
CapitalChange(K,Inv, delta) = (1-delta)*K + Inv
# Investment (equal to savings) as function of Production in (t) and savings rate s
Investment(Y, s) = s * Y
# Labour in (t+1) as function of Labour in (t) and growth rate of population n
LabourChange(L, n) = (1+n)*L
# Technology in (t+1) as function of Technology in (t) and growth rate of technology g
TechnologyChange(A, g) = (1+g)*A
function ExtractParameters(parameters_dict)
# Takes a Dictionary and creates variables with names corresponding to its key names
# and values corresponding to its values
expr = quote end
for (k, v) in parameters_dict
push!(expr.args, :($(Symbol(k)) = $v))
end
eval(expr)
return
end
function SolveFirmsProblem(ProductionFunction, K, A, L)
production_function_gradient = Calculus.gradient(x -> ProductionFunction(x[1], x[2], x[3]), [K, A, L])
W = production_function_gradient[3] # W = dF/dL
Rk = production_function_gradient[1] # Rk = dF/dK
W, Rk
end
mutable struct model_parameters
s::Float64
n::Float64
g::Float64
delta::Float64
end
parameters = model_parameters(0.04, 0.02, 0.01, 0.01)
function Solow(K0, A0, L0, parameters, T)
# In t = 0:
Y0 = Production(K0, A0, L0)
Inv0 = Investment(Y0, parameters.s)
C0 = Y0 - Inv0 # Goods market clearing condition
W0, Rk0 = SolveFirmsProblem(Production, K0, A0, L0)
# Initializing sequences of allocations and prices for times t=1..T
for x = [:C, :Y, :K, :Inv, :Rk, :W, :A, :L]
@eval $x = zeros(T)
end
t = 1
K[t] = CapitalChange(K0, Inv0, parameters.delta)
A[t] = TechnologyChange(A0, parameters.g)
L[t] = LabourChange(L0, parameters.n)
Y[t] = Production(K[t], A[t], L[t])
Inv[t] = Investment(Y[t], parameters.s)
C[t] = Y[t] - Inv[t] # Goods market clearing condition
W[t], Rk[t] = SolveFirmsProblem(Production, K[t], A[t], L[t])
for t = 2:T
K[t] = CapitalChange(K[t-1], Inv[t-1], parameters.delta)
A[t] = TechnologyChange(A[t-1], parameters.g)
L[t] = LabourChange(L[t-1], parameters.n)
Y[t] = Production(K[t], A[t], L[t])
Inv[t] = Investment(Y[t], parameters.s)
C[t] = Y[t] - Inv[t] # Goods market clearing condition
W[t], Rk[t] = SolveFirmsProblem(Production, K[t], A[t], L[t])
end
C, Y, K, Inv, Rk, W, A, L
end
K0 = 100
A0 = 1.1
L0 = 300
T = 100
C, Y, K, Inv, Rk, W, A, L = Solow(K0, A0, L0, parameters, T)
plot(C, color="blue")
In this case, what I am worried about is the following chunk:
for t = 2:T
K[t] = CapitalChange(K[t-1], Inv[t-1], parameters.delta)
A[t] = TechnologyChange(A[t-1], parameters.g)
L[t] = LabourChange(L[t-1], parameters.n)
Y[t] = Production(K[t], A[t], L[t])
Inv[t] = Investment(Y[t], parameters.s)
C[t] = Y[t] - Inv[t] # Goods market clearing condition
W[t], Rk[t] = SolveFirmsProblem(Production, K[t], A[t], L[t])
end
In this case, it looks that the equation order matter. I would like to write a system of equation without thinking about this and let a function (or a package) handle that.
The model described in your original post is not a “dynamic system” but rather a system of (linear) equations.
As you explain, one way it could be made dynamic: Z(t+1) = (1+g)*Z(t)
This package: GitHub - RJDennis/SolveDSGE.jl: A Julia package to solve, simulate, and analyze nonlinear DSGE models.
uses powerful Projection & Perturbation methods to solve dynamic, non-linear systems, including yours.
For fun let’s try to work through your example:
Z(0)=100; #initial condition
Z(t) = (1+g)*Z(t-1) #Dynamics (law of motion)
# system of 3 eq in 3 var (Y, C, I) given Z(t) & parameters
Y(t) = C(t) + I(t) + Z(t)
C(t) = c*Y(t)
I(t) = h*Y(t)
Analytical solution:
Z_t = Z_0 (1+g)^t, Y_t =\frac{Z_t}{1-c-h}, C_t =c\times Y_t , I_t =h\times Y_t
One way to simulate in Julia:
(note you should really const the parameters or put loop inside function…)
#1: create function to solve system @ time t, given param & Zt
function SolveModel(Zt, c, h)
Yt= Zt/(1-c-h)
Ct= c*Yt
It= h*Yt
return Yt, Ct, It
end
#2: Parameters, initial cond
c=0.7; h = 0.2; T = 100; g=0.01; #parameters
Z0=100; #initial condition
Z=[]; Y=[]; C=[]; I=[]; #store results
#3: simulate, loop over time
for t = 0:T
Zt = Z0*(1+g)^t
(Yt, Ct, It) = SolveModel(Zt, c, h)
push!(Z,Zt); push!(Y,Yt); push!(C,Ct); push!(I,It);
end
#4: plot time-series
using Plots
plot(0:T, [Z,Y,C,I], label = ["Z" "Y" "C" "I"], legend=:topleft, xlabel="t")
Before ppl complain, I’ll put the loop inside a function:
function SolveModel(Zt, c, h)
Yt= Zt/(1-c-h)
Ct= c*Yt
It= h*Yt
return Yt, Ct, It
end
function SimulateModel(c, h, T, g, Z0, Z, Y, C, I)
for t = 0:T
Zt = Z0*(1+g)^t
(Yt, Ct, It) = SolveModel(Zt, c, h)
push!(Z,Zt); push!(Y,Yt); push!(C,Ct); push!(I,It);
end
return Z, Y, C, I
end
#
c=0.7; h = 0.2; T = 100; g=0.01; #parameters
Z0=100; #initial condition
Z=[]; Y=[]; C=[]; I=[]; #store results
(Z, Y, C, I) = SimulateModel(c, h, T, g, Z0, Z, Y, C, I)
Wow
Thank you a lot for the example code!
This is exactly what I need!
There is a way to not worry about equation’s order? In other words, does the order that I specify Z (or any other variable) changes the results? I am saying that because I intend to work with larger models in which I may not know the best order between the variables’ definition.
Once again, thank you a lot!
Do you mean the order of the equations being solved?
For example:
f(x)=1
g(x)=2
#versus
g(x)=2
f(x)=1
In general the order shouldn’t matter, but that all depends on the underlying solver you are using to solve that system of equations. Some solvers are slower when you write a system of equations in a certain order.
Based on your code, I mean something like this:
function SolveModel(Zt, c, h)
Yt= Zt/(1-c-h)
Ct= c*Yt
It= h*Yt
return Yt, Ct, It
end
versus (or any other sequence)
function SolveModel(Zt, c, h)
Ct= c*Yt
It= h*Yt
Yt= Zt/(1-c-h)
return Yt, Ct, It
end
In other words, does julia solve algebric loops or I need to take care of the equations especification order?
Both ways should work.
Try and see yourself
Unfortunately it does not work 
When I run the second code chunk, it returns this error
UndefVarError: Yt not defined