New structure

Project.toml

@@ -5,15 +5,18 @@ version = "1.0.0-DEV"
 
 [deps]
 Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+Pipe = "b98c9c47-44ae-5843-9183-064241ee97a0"
 
 [compat]
-NLsolve = "4"
-MacroTools = "0.5"
 Crayons = "4"
+MacroTools = "0.5"
+NLsolve = "4"
 julia = "1"
 
 [extras]

README.md

@@ -7,7 +7,157 @@ This package provides solution and calibration methods for stock-flow consistent
 
 ## Basic usage
 
+### Model definition
+
+Consider SIM from Godley and Lavoie 2007:
+
+```julia
+# Define parameter values
+params_dict = @parameters begin
+    θ = 0.2
+    α_1 = 0.6
+    α_2 = 0.4
+end
+
+# Define endogenous variables
+endogenous = @variables Y, T, YD, C, H_s, H_h, H
+# Define exogenous variables
+exogenous = @variables G
+
+# Define model
+my_first_model = model(
+    endos = endogenous,
+    exos = exogenous,
+    params = params_dict,
+    equations = @equations begin
+        Y = C + G
+        T = θ * Y
+        YD = Y - T
+        C = α_1 * YD + α_2 * H[-1]
+        H_s + H_s[-1] = G - T
+        H_h + H_h[-1] = YD - C
+        H = H_s + H_s[-1] + H[-1]
+    end
+)
+```
+
+The difference between parameters and exogenous parameters is that the latter *can change over time* which is especially important if they appear in lagged terms. In this case, we need to provide several values for one exogenous variable.
+
+Note also that the model is not aware of any concrete values of parameters or exogenous variables. Instead, data is always supplied externaly to solution/calibration functions. Thus, `params_dict` is just syntactical sugar for `@variables [k for (k, v) in params_dict]`.
+
+Lastly, the specification of endogenous variables is *optional* and might be omitted if is much effort for larger models. However, it enables easier debugging. If not endogenous variables are not specified, the package will assume the symbol farthest on the left hand side of each equation to be endogenous.
+
+### Model solution
+If we want to solve a model we need data on
+1. exogenous variables (and their lags)
+2. lags of endogenous variables
+3. parameters
+
+```julia
+# data on exogenous parameter G
+exos = [20.0][:, :]
+# lagged values of endogenous variables are all 0.0
+lags = fill(0.0, length(my_first_model.endogenous_variables), 1)
+# get raw parameter values
+param_values = map(x -> params_dict[x], my_first_model.parameters)
+```
+
+
+
+```julia
+# Solve model for 59 periods
+for i in 1:59
+    solution = solve(my_first_model, lags, exos, param_values)
+    lags = hcat(lags, solution)
+end
+```
+
+### Data handling and plotting
+`DataFrames` and `Pipe` give us functionality similar to R's `dplyr`; the package ```Gadfly``` is very similar to R's `ggplots2`:
+
+```julia
+using DataFrames
+using Pipe
+using Gadfly
+
+# Convert results to DataFrame
+df = DataFrame(lags', my_first_model.endogenous_variables)
+# Add time column
+df[!, :period] = 1:nrow(df)
+# Select variables, convert to long format, and plot variables
+@pipe df |>
+    select(_, [:Y, :C, :YD, :period]) |>
+    stack(_, Not(:period), variable_name=:variable) |>
+    plot(
+        _,
+        x=:period,
+        y=:value,
+        color=:variable,
+        Geom.line
+    )
+```
+
+An example with actual data can be found here.
+
+## Advanced usage
+
+### Probabilistic models
+
+### Model calibration
+
+
+
+## Syntax
+
+There are plenty different syntax options for defining variables *outside the model function* enabled:
+
+```julia
+# As single variables (slurping)
+endogenous = @variables Y T YD C H_s H_h H
+# As tuple
+endogenous = @variables Y, T, YD, C, H_s, H_h, H
+# As array
+endogenous = @variables [
+    Y,
+    T,
+    YD,
+    C,
+    H_s,
+    H_h,
+    H
+]
+# As block
+endogenous = @variables begin
+    Y
+    T
+    YD
+    C
+    H_s
+    H_h
+    H
+end
+```
+
+Inside the function we need parantheses and can not use whitespace seperation.
+
+## Internals
+
+Internally, we just generate a function `f!` for our model which can be used together with an arbitrary root finding solver:
+
+```julia
+function f!(diff, endos, lags, exos, params)
+    diff[1] = endos[1] - (endos[4] + exos[1, end - 0])
+    diff[2] = endos[2] - params[1] * endos[1]
+    diff[3] = endos[3] - (endos[1] - endos[2])
+    diff[4] = endos[4] - (params[2] * endos[3] + params[3] * lags[7, end - 0])
+    diff[5] = (endos[5] + lags[5, end - 0]) - (exos[1, end - 0] - endos[2])
+    diff[6] = (endos[6] + lags[6, end - 0]) - (endos[3] - endos[4])
+    diff[7] = endos[7] - (endos[5] + lags[5, end - 0] + lags[7, end - 0])
+end
+```
+
 ## Remarks
 
-- The instantiation of a model is not thread-safe yet.
-- The point I would like to see proven: Julia is a long-term Pareto improvement over R, Matlab et al.
+Currently the model instantiation is not thread-safe.
+
+Feel free to ask questions and report bugs via issues!

examples/combine.jl

@@ -1,30 +1,31 @@
 using Consistent
 
-PC_gdp = @model begin
-    @endogenous Y YD T V C
-    @exogenous r G B_h
-    @parameters α_1 α_2 θ
-    @equations begin
+PC_gdp = model(
+    endos = @variables(Y, YD, T, V, C),
+    exos = @variables(r, G, B_h),
+    params = @variables(α_1, α_2, θ),
+    eqs = @equations begin
         Y = C + G
         YD = Y - T + r[-1] * B_h[-1]
         T = θ * (Y + r[-1] * B_h[-1])
         V = V[-1] + (YD - C)
         C = α_1 * YD + α_2 * V[-1]
     end
-end
+)
 
-PC_hh = @model begin
-    @endogenous H_h B_h B_s H_s B_cb r
-    @exogenous r_exo G V YD T
-    @parameters λ_0 λ_1 λ_2
-    @equations begin
+PC_hh = model(
+    endos = @variables(H_h, B_h, B_s, H_s, B_cb, r),
+    exos = @variables(r_exo, G, V, YD, T),
+    params = @variables(λ_0, λ_1, λ_2),
+    eqs = @equations begin
         H_h = V - B_h
         B_h = (λ_0 + λ_1 * r - λ_2 * (YD / V)) * V
         B_s = (G + r[-1] * B_s[-1]) - (T + r[-1] * B_cb[-1]) + B_s[-1]
         H_s = B_cb - B_cb[-1] + H_s[-1]
         B_cb = B_s - B_h
         r = r_exo
     end
-end
+)
 
+# Note: Since the variables and equations are ordered this operation is not commutative!
 PC_complete = PC_gdp + PC_hh

examples/solve.jl

@@ -1,35 +1,57 @@
 using Consistent
-using NLsolve
 using DataFrames
 using Gadfly
 using Pipe
+using BenchmarkTools
 Gadfly.push_theme(:dark)
 
-function solve(model, lags, exos, params)
-    nlsolve(
-        (x, y) -> model.f!(x, y, lags, exos, params),
-        fill(1.0, length(model.endogenous_variables)),
-        autodiff = :forward,
-    ).zero
+# Define parameter values
+params_dict = @parameters begin
+    θ = 0.2
+    α_1 = 0.6
+    α_2 = 0.4
 end
 
-# Setup SIM
-model = Consistent.SIM() # load predefined SIM model
-params_dict = Consistent.SIM(true) # load default parameters
-exos = [20.0] # this is only G
-exos = exos[:, :] # bring in matrix format
-lags = fill(0.0, length(model.endogenous_variables), 1) # lagged values of endogenous variables are all 0.0
-param_values = map(x -> params_dict[x], model.parameters) # get raw parameter values
-solution = solve(model, lags, exos, param_values) # solve first period
+# Define endogenous variables
+endogenous = @variables Y, T, YD, C, H_s, H_h, H
+# Define exogenous variables
+exogenous = @variables G
 
-# Solve model for 50 periods
-for i in 1:50
-    solution = solve(model, lags, exos, param_values) # solve first period
+# Define model equations
+equations = @equations begin
+    Y = C + G
+    T = θ * Y
+    YD = Y - T
+    C = α_1 * YD + α_2 * H[-1]
+    H_s + H_s[-1] = G - T
+    H_h + H_h[-1] = YD - C
+    H = H_s + H_s[-1] + H[-1]
+end
+
+# Define model
+my_first_model = model(
+    endos = endogenous,
+    exos = exogenous,
+    params = params_dict,
+    eqs = equations
+)
+
+# Data on exogenous parameter G
+exos = [20.0][:, :]
+# Lagged values of endogenous variables are all 0.0
+lags = fill(0.0, length(my_first_model.endogenous_variables), 1)
+# Get raw parameter values
+param_values = map(x -> params_dict[x], my_first_model.parameters)
+
+# Solve model for 59 periods
+for i in 1:59
+    solution = solve(my_first_model, lags, exos, param_values)
     lags = hcat(lags, solution)
 end
 
 # Convert results to DataFrame
-df = DataFrame(lags', model.endogenous_variables)
+df = DataFrame(lags', my_first_model.endogenous_variables)
+# Add time column
 df[!, :period] = 1:nrow(df)
 # Select variables, convert to long format, and plot variables
 @pipe df |>

examples/solve_SIM_stoch.jl

@@ -0,0 +1,46 @@
+using Consistent
+using NLsolve
+using Distributions
+using Plots
+using StatsPlots
+
+# SIM_stoch
+model = SIMStoch()
+g_0 = 20.0
+u_G = rand(Normal(0, 1), 1000 * 100)
+u_T = rand(Normal(0, 0.25), 1000 * 100)
+u_C = rand(Normal(0, 0.5), 1000 * 100)
+exos = zeros(4, 1)
+exos[1] = g_0
+initial_guess = fill(1.0, length(model.endogenous_variables))
+lags = fill(0.0, length(model.endogenous_variables))
+lags = lags[:, :]
+param_values = map(x -> params[x], model.parameters)
+a_0 = solve(model, lags, exos, param_values)
+simulation = zeros(101, 1000)
+simulation[1, :] .= a_0[2]
+
+for i in 1:1000
+    a = a_0
+    for j = 1:100
+        exos[2] = u_G[j+(i-1)*100]
+        exos[3] = u_T[j+(i-1)*100]
+        exos[4] = u_C[j+(i-1)*100]
+        a = solve(model, a[:, :], exos, param_values)
+        simulation[j+1, i] = a[2]
+    end
+end
+
+plot(simulation[:, 1])
+violin(transpose(simulation[1:50, :]), linewidth = 0, legend = false)
+# plot(simulation[1:30,1:1000], legend = false, seriestype = :scatter)
+
+# DIS steady state
+initial_guess = fill(1.0, length(sfc_model.endogenous_variables))
+lags = fill(0.1, (length(sfc_model.endogenous_variables), 1)) # quite arbitrary
+param_values = map(x -> get(values, x, nothing), sfc_model.parameters)
+a = solve(initial_guess, lags, Matrix[], param_values)
+for i in 1:1000
+    a = solve(initial_guess, a[:,:], Matrix[], param_values)
+end
+a

src/CombineModels.jl

@@ -7,13 +7,10 @@ function +(model1::Model, model2::Model)
     exos2 = filter(
         x -> !((x in model1.endogenous_variables) || (x in exos1)), model2.exogenous_variables
     )
-    equations = vcat(model1.equations, model2.equations)
-    global Consistent.sfc_model = Model()
-    sfc_model.endogenous_variables = vcat(model1.endogenous_variables, model2.endogenous_variables)
-    sfc_model.exogenous_variables = vcat(exos1, exos2)
-    sfc_model.parameters = vcat(model1.parameters, model2.parameters)
-    sfc_model.math_operators = model1.math_operators # FIXME: does this make sense?
-    sfc_model.equations = equations
-    eval(build_f!(equations))
-    deepcopy(sfc_model)
+    return model(
+        endos=Variables(vcat(model1.endogenous_variables, model2.endogenous_variables)),
+        exos=Variables(vcat(exos1, exos2)),
+        params=Variables(vcat(model1.parameters, model2.parameters)),
+        eqs=Equations(vcat(model1.equations, model2.equations))
+    )
 end

src/Consistent.jl

@@ -1,18 +1,24 @@
 module Consistent
 
-export @model, @endogenous, @exogenous, @parameters, @equations
+export @parameters, @equations, @variables
+export model, solve
 
 include("Helpers.jl")
+include("ModelComponents.jl")
 include("Model.jl")
 include("Variables.jl")
 include("ConstructResiduals.jl")
 include("Macros.jl")
 include("CombineModels.jl")
+include("Solve.jl")
 
+# Godley/Lavoie
 include("models/SIM.jl")
 include("models/SIM_stoch.jl")
 include("models/LP.jl")
 include("models/DIS.jl")
 include("models/PC.jl")
 
+include("models/BMW.jl")
+
 end # module Consistent