Simultanious system equation evaluation

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?

1 Like

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.

1 Like

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.

2 Likes

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 :wink:, 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)?

2 Likes

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 :slight_smile: )

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: https://github.com/RJDennis/SolveDSGE.jl
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")

image

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)
1 Like

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!

1 Like

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?

1 Like

Both ways should work.
Try and see yourself

Unfortunately it does not work :cry:
When I run the second code chunk, it returns this error

UndefVarError: Yt not defined