diff --git a/Project.toml b/Project.toml
index 42985f9..9afdbe9 100644
--- a/Project.toml
+++ b/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]
diff --git a/README.md b/README.md
index 2058569..0f23fc6 100644
--- a/README.md
+++ b/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!
diff --git a/examples/combine.jl b/examples/combine.jl
index 25d2641..54dc51b 100644
--- a/examples/combine.jl
+++ b/examples/combine.jl
@@ -1,23 +1,23 @@
 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]
@@ -25,6 +25,7 @@ PC_hh = @model begin
         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
\ No newline at end of file
diff --git a/examples/solve.jl b/examples/solve.jl
index f49c105..048abe0 100644
--- a/examples/solve.jl
+++ b/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 |>
diff --git a/examples/solve_SIM_stoch.jl b/examples/solve_SIM_stoch.jl
new file mode 100644
index 0000000..8d718d2
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/src/CombineModels.jl b/src/CombineModels.jl
index 0d6bf44..71e32f0 100644
--- a/src/CombineModels.jl
+++ b/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
\ No newline at end of file
diff --git a/src/Consistent.jl b/src/Consistent.jl
index 5449d43..d758400 100644
--- a/src/Consistent.jl
+++ b/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
\ No newline at end of file
diff --git a/src/Helpers.jl b/src/Helpers.jl
index 0b5fc41..b7bd9ad 100644
--- a/src/Helpers.jl
+++ b/src/Helpers.jl
@@ -7,6 +7,9 @@ Thus, we decompose expressions like :(a in b) in individual Symbols.
 function remove_expr(x::Expr)
     if x.head == :call
         return [x.args[2], x.args[1], x.args[3]]
+    # TODO: test whether this might have unintentional side effects
+    elseif x.head == :block
+        return x.args
     else
         error("Can not handle $x")
     end
@@ -61,11 +64,43 @@ function find_symbols(line::Expr)
     found::Set{Symbol} = Set([])
     args = line.args
     for i in eachindex(args)
-        if typeof(args[i]) == Symbol
+        if args[i] isa Symbol
             push!(found, args[i])
-        elseif typeof(args[i]) == Expr
+        elseif args[i] isa Expr
             union!(found, find_symbols(args[i]))
         end
     end
     return setdiff(found, math_operators)
+end
+
+"""
+Get the symbol farthest to the left in an Expr.
+"""
+function left_symbol(line::Expr) # TODO: why do we not need to check heads?
+    args = line.args
+    for i in eachindex(args)
+        if args[i] isa Symbol && !(args[i] in math_operators)
+            return args[i]
+        elseif args[i] isa Expr
+            found = left_symbol(args[i])
+            if found isa Symbol
+                return found
+            end
+        end
+    end
+    return nothing
+end
+
+"""
+Handle input variables in
+- array form
+- coma-seperated form
+- whitespace seperated form
+"""
+function handle_input(input)
+    if (length(input) == 1) && isa(input[1], Expr) && (input[1].head in (:vect, :tuple))
+        input[1].args
+    else
+        input
+    end
 end
\ No newline at end of file
diff --git a/src/Macros.jl b/src/Macros.jl
index 9c79726..350c99f 100644
--- a/src/Macros.jl
+++ b/src/Macros.jl
@@ -1,36 +1,9 @@
-"""
-Macro to specify the endogenous variables.
-
-# Example:
-    @endogenous Y C
-"""
-macro endogenous(input...) # TODO: handle coma-seperated input; assert len endos = nr of equations
-    Consistent.sfc_model.endogenous_variables = remove_expr([input...])
-end
-
-"""
-Macro to specify the exogenous variables. Exogenous variables can change over time and are assumed to be readily available from data source.
-
-# Example:
-    @exogenous G X
-"""
-macro exogenous(input...)
-    Consistent.sfc_model.exogenous_variables = remove_expr([input...])
-end
-
-"""
-Macro to specify the parameters. Parameters typically can not change over time and can be calibrated to fit given data.
-
-# Example:
-    @parameters α θ
-"""
-macro parameters(input...)
-    Consistent.sfc_model.parameters = remove_expr([input...])
-end
+function build_f!(endos, exos, params, args)
+    endos = endos.variables
+    exos = exos.variables
+    params = params.variables
 
-function build_f!(args)
     function_body = deepcopy(args)
-
     for i in eachindex(function_body)
         # check if we really have a proper equation
         if (function_body[i].head == :(=))
@@ -43,9 +16,6 @@ function build_f!(args)
 
     # construct arrays for different types of variables
     found = vars(function_body) # all found variables
-    endos = Consistent.sfc_model.endogenous_variables # endogenous variables
-    exos = Consistent.sfc_model.exogenous_variables # exogenous variables
-    params = Consistent.sfc_model.parameters # parameters
     variables = union(Set(endos), Set(exos), Set(params)) # all variables
 
     # check for unused variables from specification
@@ -56,22 +26,7 @@ function build_f!(args)
     end
 
     # construct function for residuals of model variables
-    return MacroTools.striplines(:(Consistent.sfc_model.f! = $(construct_residuals(name, function_body, args))))
-end
-
-"""
-Macro to specify the model equations. Use `begin ... end`.
-
-# Example:
-    @equations begin
-        Y = G + C
-    end
-"""
-macro equations(input...) # FIXME: refactor
-    ex = remove_blocks(MacroTools.striplines(input...))
-    @assert (ex.head == :block) "Block input expected" # we need block input (begin ... end)
-    Consistent.sfc_model.equations = deepcopy(ex.args)
-    return build_f!(ex.args)
+    return MacroTools.striplines(:(Consistent.f! = $(construct_residuals(name, function_body, args))))
 end
 
 """
@@ -79,11 +34,11 @@ Macro to build a stock-flow consistent model.
 
 # Example:
 ```julia-repl
-julia> @model begin
-    @endogenous Y T YD C H_s H_h H
-    @exogenous G
-    @parameters θ α_1 α_2
-    @equations begin
+julia> model(
+    endos = @variables(Y, T, YD, C, H_s, H_h, H),
+    exos = @variables(G),
+    params = @variables(θ, α_1, α_2),
+    eqs = @equations begin
         Y = C + G
         T = θ * Y
         YD = Y - T
@@ -92,22 +47,48 @@ julia> @model begin
         H_h + H_h[-1] = YD - C
         H = H_s + H_s[-1] + H[-1]
     end
-end
+)
+Stock-flow consistent model
 Endogenous Variables: [:Y, :T, :YD, :C, :H_s, :H_h, :H]
 Exogenous Variables:  [:G]
 Parameters:           [:θ, :α_1, :α_2]
 Equations:            
-           (1)  Y = C + G
-           (2)  T = θ * Y
-           (3)  YD = Y - T
-           (4)  C = α_1 * YD + α_2 * H[-1]
-           (5)  H_s + H_s[-1] = G - T
-           (6)  H_h + H_h[-1] = YD - C
-           (7)  H = H_s + H_s[-1] + H[-1]
+                      (1)  Y = C + G
+                      (2)  T = θ * Y
+                      (3)  YD = Y - T
+                      (4)  C = α_1 * YD + α_2 * H[-1]
+                      (5)  H_s + H_s[-1] = G - T
+                      (6)  H_h + H_h[-1] = YD - C
+                      (7)  H = H_s + H_s[-1] + H[-1]
 ```
 """
-macro model(input...) # FIXME: missing specification (e.g. endo YD) not always working
-    global Consistent.sfc_model = Model()
-    eval(input...)
-    return deepcopy(Consistent.sfc_model)
+function model(;
+    endos=nothing::Union{Variables, Nothing},
+    exos=Variables(),
+    params=Variables()::Union{Variables, OrderedDict},
+    eqs,
+    verbose=false
+)
+    if params isa OrderedDict # FIXME: use promotion
+        parameters = Variables(params)
+    else # FIXME
+        parameters = params
+    end
+
+    if isnothing(endos)
+        endos = Variables(left_symbol.(eqs.exprs))
+    end
+
+    if verbose
+        println(build_f!(endos, exos, parameters, eqs.exprs))
+    end
+
+    eval(build_f!(endos, exos, parameters, eqs.exprs))
+    return Model(
+        endos,
+        exos,
+        parameters,
+        eqs,
+        deepcopy(Consistent.f!)
+    )
 end
\ No newline at end of file
diff --git a/src/Model.jl b/src/Model.jl
index 3c7e4e2..39073d5 100644
--- a/src/Model.jl
+++ b/src/Model.jl
@@ -7,17 +7,14 @@ The most important part is the automatically generated function `f!` which has t
     model.f!(residuals, endos, lags, exos, params)
 Intuitively, we evaluate our function `f(endos, ...)`` (which should equal zero) into residuals.
 """
-mutable struct Model
-    endogenous_variables::Vector{Symbol}
-    exogenous_variables::Vector{Symbol}
-    parameters::Vector{Symbol}
-    math_operators::Set{Symbol}
-    equations::Vector{Expr}
+struct Model
+    endogenous_variables::Variables
+    exogenous_variables::Variables
+    parameters::Variables
+    equations::Equations
     f!
 end
 
-Model() = Model(Symbol[], Symbol[], Symbol[], math_operators, Expr[], x -> nothing)
-
 const math_operators = Set([:+, :-, :*, :/, :÷, :\, :^, :%])
 const name = [:diff, :endos, :lags, :exos, :params]
 
@@ -27,11 +24,13 @@ function Base.show(io::IO, m::Model)
     for i in eachindex(descriptors)
         descriptors[i] = descriptors[i] * ' '^(max_width - length(descriptors[i]))
     end
+    println(io, "Stock-flow consistent model")
     print(io, Crayon(foreground = :green), descriptors[1]); println(io, Crayon(reset=true), m.endogenous_variables)
     print(io, Crayon(foreground = :yellow), descriptors[2]); println(io, Crayon(reset=true), m.exogenous_variables)
     print(io, Crayon(foreground = :blue), descriptors[3]); println(io, Crayon(reset=true), m.parameters)
     print(io, Crayon(foreground = :red), descriptors[4]); print(io, Crayon(reset=true))
     for i in eachindex(m.equations)
-        print(io, "\n", ' '^max_width, "($i)  ", m.equations[i])
+        additional_space = div(length(m.equations), 10) - div(i, 10)
+        print(io, "\n", ' '^(max_width + max(additional_space, 0)), "($i)  ", m.equations[i])
     end
 end
\ No newline at end of file
diff --git a/src/ModelComponents.jl b/src/ModelComponents.jl
new file mode 100644
index 0000000..dfb0af1
--- /dev/null
+++ b/src/ModelComponents.jl
@@ -0,0 +1,68 @@
+using OrderedCollections
+import Base
+
+# Wrapper for vector of symbols
+struct Variables <: AbstractVector{Symbol}
+    variables::Vector{Symbol} # TODO: add documentation for variables
+end
+
+Variables(x::OrderedDict) = Variables([k for (k, v) in x])
+Variables(x::Variables) = x
+Variables() = Variables(Symbol[])
+
+MacroTools.@forward Variables.variables Base.getindex, Base.setindex!, Base.size
+
+macro variables(input...)
+    if (length(input) > 0) && (input[1] isa Expr) && (input[1].head == :block)
+        @assert (length(input) == 1) "Can't handle several blocks"
+        args = input[1].args
+        vars = filter(e -> isa(e, Symbol), args)
+        return Variables(deepcopy(vars))
+    else # convert potential tuple to array
+        return Variables(remove_expr([handle_input(input)...]))
+    end
+end
+
+"""
+Macro to specify the parameters. Parameters typically can not change over time and can be calibrated to fit given data.
+
+Returns an OrderedDict.
+
+# Example:
+    @parameters begin
+        θ = 0.2
+        α_1 = 0.6
+        α_2 = 0.4
+    end
+"""
+macro parameters(block)
+    exprs = block.args
+
+    # Filter for assignments
+    assignments = filter(e -> isa(e, Expr) && e.head == :(=), exprs)
+
+    # Extract variable names and their values as symbols with colons
+    pairs = [Expr(:call, :(=>), :(Symbol($("$(a.args[1])"))), a.args[2]) for a in assignments]
+
+    return esc(Expr(:call, :(Consistent.OrderedDict), pairs...))
+end
+
+struct Equations <: AbstractVector{Expr}
+    exprs::Vector{Expr}
+end
+# FIXME: print as strings
+MacroTools.@forward Equations.exprs Base.getindex, Base.setindex!, Base.size
+
+"""
+Macro to specify the model equations. Use `begin ... end`.
+
+# Example:
+    @equations begin
+        Y = G + C
+    end
+"""
+macro equations(input...)
+    ex = remove_blocks(MacroTools.striplines(input...)) # TODO: better debugging with LineNumberNodes
+    @assert (ex.head == :block) "Block input expected" # we need block input (begin ... end)
+    return Equations(deepcopy(ex.args))
+end
\ No newline at end of file
diff --git a/src/Solve.jl b/src/Solve.jl
new file mode 100644
index 0000000..9f26e8f
--- /dev/null
+++ b/src/Solve.jl
@@ -0,0 +1,9 @@
+using NLsolve
+
+function solve(model, lags, exos, params, initial=fill(1.0, length(model.endogenous_variables)))
+    nlsolve(
+        (x, y) -> model.f!(x, y, lags, exos, params),
+        initial,
+        autodiff = :forward,
+    ).zero
+end
\ No newline at end of file
diff --git a/src/Variables.jl b/src/Variables.jl
index dba7286..25c3f66 100644
--- a/src/Variables.jl
+++ b/src/Variables.jl
@@ -1,7 +1,7 @@
 """
 Find all symbols in an array of expressions.
 """
-function vars(lines::Vector)
+function vars(lines::Vector) # FIXME: consider that some symbols are external functions
     found::Set{Symbol} = Set([])
     for i in eachindex(lines)
         union!(found, find_symbols(lines[i]))
diff --git a/src/models/BMW.jl b/src/models/BMW.jl
index beafbe6..60b3a3d 100644
--- a/src/models/BMW.jl
+++ b/src/models/BMW.jl
@@ -2,31 +2,35 @@
 Note: This is some old model where I don't know the origin anymore.
 It is probably related to some Peter Bofinger textbook.
 """
-BMW() = @model begin
-    @endogenous C_s C_d I_s I_d N_s N_d L_s L_d r_l r_m K K_T DA
-    @exogenous r_exo
-    @parameters α_0 α_1_w α_1_r α_2
-    @equations begin
-        C_s = C_d
-        I_s = I_d
-        N_s = N_d
-        L_s = L_d - L_d[-1] + L_s[-1]
-        Y = C_s + I_s
-        WB_d = Y - r_l[-1] * L_d[-1] - AF
-        AF = delta * K[-1]
-        L_d = I_d - AF + L_d[-1]
-        YD = WB_s + r_m[-1] * M_d[-1]
-        M_h = YD - C_d + M_h[-1]
-        M_s = L_s - L_s[-1] - M_s[-1]
-        r_m = r_l
-        WB_s = W * N_s
-        N_d = Y / pr
-        W = WB_d / N_d
-        C_d = α_0 + α_1_w * WB_s + α_1_r * r_m[-1] * M_h[-1] + α_2 * M_h
-        K = I_d - DA - K[-1]
-        DA = delta * K[-1]
-        K_T = kappa * Y[-1]
-        I_d = gamma * (K_T - K[-1]) + DA
-        r_l = r_exo
-    end
+function BMW()
+    Dict(
+        :model => model(
+            exos = @variables(r_exo),
+            params = @variables(α_0, α_1_w, α_1_r, α_2, δ, κ, γ, pr),
+            eqs = @equations begin
+                C_s = C_d
+                I_s = I_d
+                N_s = N_d
+                L_s = L_d - L_d[-1] + L_s[-1]
+                Y = C_s + I_s
+                WB_d = Y - r_l[-1] * L_d[-1] - AF
+                AF = δ * K[-1]
+                L_d = I_d - AF + L_d[-1]
+                YD = WB_s + r_m[-1] * M_d[-1]
+                M_h = YD - C_d + M_h[-1]
+                M_s = L_s - L_s[-1] - M_s[-1]
+                r_m = r_l
+                WB_s = W * N_s
+                N_d = Y / pr
+                W = WB_d / N_d
+                C_d = α_0 + α_1_w * WB_s + α_1_r * r_m[-1] * M_h[-1] + α_2 * M_h
+                K = I_d - DA - K[-1]
+                DA = δ * K[-1]
+                K_T = κ * Y[-1]
+                I_d = γ * (K_T - K[-1]) + DA
+                r_l = r_exo
+                M_d = M_h
+            end
+        )
+    )
 end
\ No newline at end of file
diff --git a/src/models/DIS.jl b/src/models/DIS.jl
index d7045ab..a8b9391 100644
--- a/src/models/DIS.jl
+++ b/src/models/DIS.jl
@@ -1,48 +1,50 @@
-DIS() = @model begin
-    @endogenous y in_T in_e in s_e s N WB UC IN S p NHUC F L_d L_s M_s r_m F_b YD M_h yd_hs C m_h c yd_e_hs
-    @exogenous
-    @parameters α_0 α_1 α_2 σ_T γ β pr W ϕ r_l add ε
-    @equations begin
-        y = s_e + (in_e[0] - in[-1])
-        in_T = σ_T * s_e
-        in_e[0] = in[-1] + γ * (in_T - in[-1])
-        in = in[-1] + (y - s)
-        s_e = β * s[-1] + (1 - β) * s_e[-1]
-        s = c
-        N = y / (pr)
-        WB = N * W
-        UC = WB / y
-        IN = in * UC
-        S = p * s
-        p = (1 + ϕ) * NHUC
-        NHUC = (1 - σ_T) * UC + σ_T * (1 + r_l) * UC[-1]
-        F = S - WB + (IN[0] - IN[-1]) - r_l * IN[-1]
-        L_d = IN
-        L_s = L_d
-        M_s = L_s
-        r_m = r_l - add
-        F_b = r_l * L_s[-1] - r_m[-1] * M_h[-1]
-        YD = WB + F + F_b + r_m[-1] * M_h[-1]
-        M_h[0] - M_h[-1] = YD - C
-        yd_hs = c + (m_h[0] - m_h[-1])
-        C = c * p
-        m_h = M_h / p
-        c = α_0 + α_1 * yd_e_hs + α_2 * m_h[-1]
-        yd_e_hs = ε * yd_hs[-1] + (1 - ε) * yd_e_hs[-1]
+function DIS()
+    params = @parameters begin
+        α_0 = 15.0
+        α_1 = 0.8
+        α_2 = 0.1
+        σ_T = 0.15
+        γ = 0.25
+        β = 0.75
+        pr = 1.0
+        W = 0.75
+        ϕ = 0.25
+        r_l = 0.025
+        add = 0.02
+        ε = 0.75
     end
-end
-
-params = Dict(
-    :α_0 => 15.0,
-    :α_1 => 0.8,
-    :α_2 => 0.1,
-    :σ_T => 0.15,
-    :γ => 0.25,
-    :β => 0.75,
-    :pr => 1.0,
-    :W => 0.75,
-    :ϕ => 0.25,
-    :r_l => 0.025,
-    :add => 0.02,
-    :ε => 0.75
-)
\ No newline at end of file
+    Dict(
+        :params => params,
+        :model => model(
+            params = params,
+            eqs = @equations begin
+                y = s_e + (in_e[0] - in[-1])
+                in_T = σ_T * s_e
+                in_e[0] = in[-1] + γ * (in_T - in[-1])
+                in = in[-1] + (y - s)
+                s_e = β * s[-1] + (1 - β) * s_e[-1]
+                s = c
+                N = y / (pr)
+                WB = N * W
+                UC = WB / y
+                IN = in * UC
+                S = p * s
+                p = (1 + ϕ) * NHUC
+                NHUC = (1 - σ_T) * UC + σ_T * (1 + r_l) * UC[-1]
+                F = S - WB + (IN[0] - IN[-1]) - r_l * IN[-1]
+                L_d = IN
+                L_s = L_d
+                M_s = L_s
+                r_m = r_l - add
+                F_b = r_l * L_s[-1] - r_m[-1] * M_h[-1]
+                YD = WB + F + F_b + r_m[-1] * M_h[-1]
+                M_h[0] - M_h[-1] = YD - C
+                yd_hs = c + (m_h[0] - m_h[-1])
+                C = c * p
+                m_h = M_h / p
+                c = α_0 + α_1 * yd_e_hs + α_2 * m_h[-1]
+                yd_e_hs = ε * yd_hs[-1] + (1 - ε) * yd_e_hs[-1]
+            end
+        ),
+    )
+end
\ No newline at end of file
diff --git a/src/models/LP.jl b/src/models/LP.jl
index 0931aa1..9a55dab 100644
--- a/src/models/LP.jl
+++ b/src/models/LP.jl
@@ -1,47 +1,53 @@
-LP() = @model begin
-    @endogenous Y YD_r T V CG C V_e H_h H_d B_d BL_d B_h BL_h B_s H_s B_cb BL_s ERr_bL r_bL p_bL_e CG_e YD_r_e r_b p_bL
-    @exogenous G
-    @parameters α_1 α_2 λ_20 λ_22 λ_23 λ_24 λ_30 λ_32 λ_33 λ_34 χ θ r_b_exo p_bL_exo
-    @equations begin
-        Y = C + G
-        YD_r = Y - T + r_b[-1] * B_h[-1] + BL_h[-1]
-        T = θ * (Y + r_b[-1] + BL_h[-1])
-        V = V[-1] + (YD_r - C) + CG
-        CG = (p_bL - p_bL[-1]) * BL_h[-1]
-        C = α_1 * YD_r_e + α_2 * V[-1]
-        V_e = V[-1] + (YD_r_e - C) + CG
-        H_h = V - B_h - p_bL * BL_h
-        H_d = V_e - B_h - p_bL * BL_h
-        B_d = (λ_20 + λ_22 * r_b + λ_23 * ERr_bL + λ_24 * (YD_r_e / V_e)) * V_e
-        BL_d = (λ_30 + λ_32 * r_b + λ_33 * ERr_bL + λ_34 * (YD_r_e / V_e)) * (V_e / p_bL)
-        B_h = B_d
-        BL_h = BL_d
-        B_s = (G + r_b[-1] * B_s[-1] + BL_s[-1]) - (T + r_b[-1] * B_cb[-1]) - ((BL_s - BL_s[-1]) * p_bL) + B_s[-1]
-        H_s = B_cb - B_cb[-1] + H_s[-1]
-        B_cb = B_s - B_h
-        BL_s = BL_h
-        ERr_bL = r_bL + χ * ((p_bL_e - p_bL) / p_bL)
-        r_bL = 1 / p_bL
-        p_bL_e = p_bL
-        CG_e = χ * (p_bL_e - p_bL) * BL_h
-        YD_r_e = YD_r[-1]
-        r_b = r_b_exo
-        p_bL = p_bL_exo
+function LP()
+    params = @parameters begin
+        θ = 0.194
+        α_1 = 0.8
+        α_2 = 0.2
+        λ_20 = 0.442
+        λ_22 = 1.1
+        λ_23 = -1
+        λ_24 = -0.03
+        λ_30 = 0.4
+        λ_32 = -1.
+        λ_33 = 1.1
+        λ_34 = -0.03
+        χ = 0.1
+        r_b_exo = 0.03
+        p_bL_exo = 20.
     end
+    Dict(
+        :params => params,
+        :model => model(
+            exos = @variables(G),
+            params = params,
+            eqs = @equations begin
+                Y = C + G
+                YD_r = Y - T + r_b[-1] * B_h[-1] + BL_h[-1]
+                T = θ * (Y + r_b[-1] + BL_h[-1])
+                V = V[-1] + (YD_r - C) + CG
+                CG = (p_bL - p_bL[-1]) * BL_h[-1]
+                C = α_1 * YD_r_e + α_2 * V[-1]
+                V_e = V[-1] + (YD_r_e - C) + CG
+                H_h = V - B_h - p_bL * BL_h
+                H_d = V_e - B_h - p_bL * BL_h
+                B_d = (λ_20 + λ_22 * r_b + λ_23 * ERr_bL + λ_24 * (YD_r_e / V_e)) * V_e
+                BL_d = (λ_30 + λ_32 * r_b + λ_33 * ERr_bL + λ_34 * (YD_r_e / V_e)) * (V_e / p_bL)
+                B_h = B_d
+                BL_h = BL_d
+                B_s = (G + r_b[-1] * B_s[-1] + BL_s[-1]) - (T + r_b[-1] * B_cb[-1]) - ((BL_s - BL_s[-1]) * p_bL) + B_s[-1]
+                H_s = B_cb - B_cb[-1] + H_s[-1]
+                B_cb = B_s - B_h
+                BL_s = BL_h
+                ERr_bL = r_bL + χ * ((p_bL_e - p_bL) / p_bL)
+                r_bL = 1 / p_bL
+                p_bL_e = p_bL
+                CG_e = χ * (p_bL_e - p_bL) * BL_h
+                YD_r_e = YD_r[-1]
+                r_b = r_b_exo
+                p_bL = p_bL_exo
+            end
+        )
+    )
 end
 
-# theta = 0.194
-# alpha_1 = 0.8
-# alpha_2 = 0.2
-# lambda_20 = 0.442
-# lambda_22 = 1.1
-# lambda_23 = -1
-# lambda_24 = -0.03
-# lambda_30 = 0.4
-# lambda_32 = -1.
-# lambda_33 = 1.1
-# lambda_34 = -0.03
-# chi = 0.1
-# r_b_exo = 0.03
-# p_bL_exo = 20.
 # G = 20.
\ No newline at end of file
diff --git a/src/models/PC.jl b/src/models/PC.jl
index 7843f35..f754405 100644
--- a/src/models/PC.jl
+++ b/src/models/PC.jl
@@ -1,18 +1,22 @@
-PC() = @model begin
-    @endogenous Y YD T V C H_h B_h B_s H_s B_cb r
-    @exogenous r_exo G
-    @parameters α_1 α_2 λ_0 λ_1 λ_2 θ
-    @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]
-        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
+function PC()
+    Dict(
+        :model => model(
+            endos = @variables(Y, YD, T, V, C, H_h, B_h, B_s, H_s, B_cb, r),
+            exos = @variables(r_exo, G),
+            params = @variables(α_1, α_2, λ_0, λ_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]
+                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
\ No newline at end of file
diff --git a/src/models/SIM.jl b/src/models/SIM.jl
index 9c44310..88801de 100644
--- a/src/models/SIM.jl
+++ b/src/models/SIM.jl
@@ -1,22 +1,25 @@
-SIM() = @model begin
-    @endogenous Y T YD C H_s H_h H
-    @exogenous G
-    @parameters θ α_1 α_2
-    @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]
+function SIM()
+    params = @parameters begin
+        θ = 0.2
+        α_1 = 0.6
+        α_2 = 0.4
     end
-end
 
-function SIM(bool)
-    if bool
-        Dict(:θ => 0.2, :α_1 => 0.6, :α_2 => 0.4)
-    else
-        SIM()
-    end
+    Dict(
+        :params => params,
+        :model => model(
+            endos=@variables(Y, T, YD, C, H_s, H_h, H),
+            exos=@variables(G),
+            params=params,
+            eqs=@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
+        )
+    )
 end
\ No newline at end of file
diff --git a/src/models/SIM_stoch.jl b/src/models/SIM_stoch.jl
index 67684b6..af8e700 100644
--- a/src/models/SIM_stoch.jl
+++ b/src/models/SIM_stoch.jl
@@ -1,8 +1,8 @@
-SIMStoch() = @model begin
-    @endogenous G Y T YD C_e C H_s H_h H
-    @exogenous G_0 u_G u_T u_C
-    @parameters θ α_1 α_2 α_3
-    @equations begin
+SIMStoch() = model(
+    endos = @variables(G, Y, T, YD, C_e, C, H_s, H_h, H),
+    exos = @variables(G_0, u_G, u_T, u_C),
+    params = @variables(θ, α_1, α_2, α_3),
+    eqs = @equations begin
         G = G_0 + u_G
         Y = C + G
         T = θ * Y + u_T
@@ -13,7 +13,7 @@ SIMStoch() = @model begin
         H_h + H_h[-1] = YD - C
         H = H_s + H_s[-1] + H[-1]
     end
-end
+)
 
 params = Dict(
     :θ => 0.2, 
diff --git a/test/runtests.jl b/test/runtests.jl
index c85860b..60cda5f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,5 +2,26 @@ using Consistent
 using Test
 
 @testset "Consistent.jl" begin
-    # Write your tests here.
-end
+
+    @testset "Internals" begin
+        @test Consistent.remove_expr([:x :(a in b) :y]) == [:x, :a, :in, :b, :y]
+        test_eqs = quote
+            z = y * (y[-1] + 0.5 * z) * θ + x[-1]
+            y = z[-2] * x * b
+        end
+        @test Consistent.replace_vars(test_eqs.args[[2, 4]], [:z, :y], Symbol[:x], [:θ]) == [
+            :(z = endos[2] * (lags[2, end - 0] + 0.5 * endos[1]) * params[1] + exos[1, end - -1]),
+            :(y = lags[1, end - -1] * exos[1, end - 0] * b)
+        ]
+    end
+
+    @testset "Default models" begin
+        sim = Consistent.SIM()
+        @test sim[:model].exogenous_variables.variables == [:G]
+        Consistent.SIMStoch()
+        Consistent.LP()
+        Consistent.PC()
+        Consistent.DIS()
+        Consistent.BMW()
+    end
+end
\ No newline at end of file