diff --git a/Project.toml b/Project.toml
index 88be086..42985f9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,11 +5,15 @@ version = "1.0.0-DEV"
 
 [deps]
 Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 
 [compat]
 NLsolve = "4"
+MacroTools = "0.5"
+Crayons = "4"
 julia = "1"
 
 [extras]
diff --git a/README.md b/README.md
index 8c0a8c4..221cae7 100644
--- a/README.md
+++ b/README.md
@@ -7,3 +7,6 @@ This package provides solution and calibration methods for stock-flow consistent
 
 ## Basic usage
 
+## Remarks
+
+The instantiation of a model is not thread-safe.
diff --git a/examples/calibrate.jl b/examples/calibrate.jl
new file mode 100644
index 0000000..ea56f59
--- /dev/null
+++ b/examples/calibrate.jl
@@ -0,0 +1,15 @@
+using Consistent
+using DataFrames
+using ForwardDiff
+
+function calibrate!(data, model, parameters=Dict())
+    x -> model.f!(y, x, ...) + x
+    function loss(data,)
+
+    end
+    ForwardDiff.derivative(() -> loss, input)
+end
+
+model = SIM()
+df = DataFrame(reverse(lags[:, end-20:end-1], dims=2)', model.endogenous_variables)
+df
\ No newline at end of file
diff --git a/examples/combine.jl b/examples/combine.jl
new file mode 100644
index 0000000..25d2641
--- /dev/null
+++ b/examples/combine.jl
@@ -0,0 +1,30 @@
+using Consistent
+
+PC_gdp = @model begin
+    @endogenous Y YD T V C
+    @exogenous r G B_h
+    @parameters α_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]
+    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
+        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
+
+PC_complete = PC_gdp + PC_hh
\ No newline at end of file
diff --git a/examples/solve.jl b/examples/solve.jl
new file mode 100644
index 0000000..f49c105
--- /dev/null
+++ b/examples/solve.jl
@@ -0,0 +1,44 @@
+using Consistent
+using NLsolve
+using DataFrames
+using Gadfly
+using Pipe
+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
+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
+
+# Solve model for 50 periods
+for i in 1:50
+    solution = solve(model, lags, exos, param_values) # solve first period
+    lags = hcat(lags, solution)
+end
+
+# Convert results to DataFrame
+df = DataFrame(lags', model.endogenous_variables)
+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
+    )
\ No newline at end of file
diff --git a/src/CombineModels.jl b/src/CombineModels.jl
new file mode 100644
index 0000000..0d6bf44
--- /dev/null
+++ b/src/CombineModels.jl
@@ -0,0 +1,19 @@
+import Base: +
+
+function +(model1::Model, model2::Model)
+    undetermined = findall(in(model2.endogenous_variables), model1.endogenous_variables)
+    @assert isempty(undetermined) "The endogenous variables $(model1[undetermined]) appear twice"
+    exos1 = filter(x -> !(x in model2.endogenous_variables), model1.exogenous_variables)
+    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)
+end
\ No newline at end of file
diff --git a/src/Consistent.jl b/src/Consistent.jl
index 714c69c..5449d43 100644
--- a/src/Consistent.jl
+++ b/src/Consistent.jl
@@ -7,5 +7,12 @@ include("Model.jl")
 include("Variables.jl")
 include("ConstructResiduals.jl")
 include("Macros.jl")
+include("CombineModels.jl")
 
-end
+include("models/SIM.jl")
+include("models/SIM_stoch.jl")
+include("models/LP.jl")
+include("models/DIS.jl")
+include("models/PC.jl")
+
+end # module Consistent
\ No newline at end of file
diff --git a/src/ConstructResiduals.jl b/src/ConstructResiduals.jl
index 416b42b..b85e83d 100644
--- a/src/ConstructResiduals.jl
+++ b/src/ConstructResiduals.jl
@@ -1,4 +1,4 @@
-function construct_residuals(name, function_body, expr)
+function construct_residuals(name, function_body, args)
     f! = quote
         function f!($(name[1]), $(name[2]), $(name[3]), $(name[4]), $(name[5]))
             nothing
@@ -9,10 +9,8 @@ function construct_residuals(name, function_body, expr)
     for i in eachindex(function_body)
         function_body[i] = :($(name[1])[$i] = $(function_body[i]))
     end
-    body = deepcopy(expr)
-    body.args = function_body
 
     # add function body to function
-    f!.args[2].args[end] = body
+    f!.args[2].args[end] = Expr(:block, function_body...)
     return f!
 end
\ No newline at end of file
diff --git a/src/Loss.jl b/src/Loss.jl
new file mode 100644
index 0000000..66cc0c3
--- /dev/null
+++ b/src/Loss.jl
@@ -0,0 +1,24 @@
+"""
+Root-mean-square error of a time series.
+"""
+rmse(x, y) = sqrt(sum((x .- y).^2) / length(x))
+
+"""
+RMSE standardized by the mean of our target values.
+"""
+standardized_rmse(x, y) = rmse(x, y) / mean(y)
+
+"""
+Mean standardized RMSE for a vector of time series. Makes only sense for strictly positive variables.
+
+Note: Without some kind of stationarity this estimation does not really make sense either,
+but this is often a general problem of stock-flow consistent models.
+"""
+function msrmse(x, y)
+    # loop over columns
+    average = 0.0
+    for i in axes(x, 2)
+        average += standardized_rmse(x[:, i], y[:, i])
+    end
+    return average
+end
\ No newline at end of file
diff --git a/src/Macros.jl b/src/Macros.jl
index 5c31984..9c79726 100644
--- a/src/Macros.jl
+++ b/src/Macros.jl
@@ -4,12 +4,12 @@ Macro to specify the endogenous variables.
 # Example:
     @endogenous Y C
 """
-macro endogenous(input...)
+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.
+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
@@ -19,7 +19,7 @@ macro exogenous(input...)
 end
 
 """
-Macro to specify the parameters.
+Macro to specify the parameters. Parameters typically can not change over time and can be calibrated to fit given data.
 
 # Example:
     @parameters α θ
@@ -28,32 +28,17 @@ macro parameters(input...)
     Consistent.sfc_model.parameters = remove_expr([input...])
 end
 
-"""
-Macro to specify the model equations. Use `begin ... end`.
+function build_f!(args)
+    function_body = deepcopy(args)
 
-# Example:
-    @equations begin
-        Y = G + C
-    end
-"""
-macro equations(input...) # FIXME: refactor
-    ex = remove_blocks(MacroTools.striplines(input...))
-    Consistent.sfc_model.equations = deepcopy(ex.args)
-    function_body = deepcopy(ex.args)
-
-    # we need block input (begin ... end)
-    if (ex.head == :block)
-        for i in eachindex(function_body)
-            # check if we really have a proper equation
-            if (function_body[i].head == :(=))
-                function_body[i] =
-                    :($(function_body[i].args[1]) - $(function_body[i].args[2]))
-            else
-                error("Only equalities are supported.")
-            end
+    for i in eachindex(function_body)
+        # check if we really have a proper equation
+        if (function_body[i].head == :(=))
+            function_body[i] =
+                :($(function_body[i].args[1]) - $(function_body[i].args[2]))
+        else
+            error("Only equalities are supported.")
         end
-    else
-        error("Block input expected")
     end
 
     # construct arrays for different types of variables
@@ -71,7 +56,22 @@ macro equations(input...) # FIXME: refactor
     end
 
     # construct function for residuals of model variables
-    return MacroTools.striplines(:(Consistent.sfc_model.f! = $(construct_residuals(name, function_body, ex))))
+    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)
 end
 
 """
@@ -106,7 +106,7 @@ Equations:
            (7)  H = H_s + H_s[-1] + H[-1]
 ```
 """
-macro model(input...)
+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)
diff --git a/src/Model.jl b/src/Model.jl
index cea7cb8..3c7e4e2 100644
--- a/src/Model.jl
+++ b/src/Model.jl
@@ -2,6 +2,10 @@ using Crayons
 
 """
 Type for a stock-flow consistent model.
+
+The most important part is the automatically generated function `f!` which has the following form:
+    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}
diff --git a/src/Variables.jl b/src/Variables.jl
index 2617285..dba7286 100644
--- a/src/Variables.jl
+++ b/src/Variables.jl
@@ -61,7 +61,7 @@ function create_vars( # FIXME: read properly
                         completed_line.args = [name[1], :($position)]
                     else
                         if args[2] < 0
-                            completed_line.args = [name[2], :($position), :($(-args[2]))]
+                            completed_line.args = [name[2], :($position), Expr(:call, :-, :end, args[2] + 1)]
                         else
                             error("future indices are not allowed!")
                         end
@@ -74,7 +74,7 @@ function create_vars( # FIXME: read properly
                 position = findall(x -> x == args[1], exos)[1]
                 if length(args) == 2
                     if args[2] <= 0
-                        completed_line.args = [name[3], :($position), :($(-args[2] + 1))]
+                        completed_line.args = [name[3], :($position), Expr(:call, :-, :end, args[2])]
                     else
                         error("future indices are not allowed!")
                     end
diff --git a/src/models/BMW.jl b/src/models/BMW.jl
new file mode 100644
index 0000000..beafbe6
--- /dev/null
+++ b/src/models/BMW.jl
@@ -0,0 +1,32 @@
+"""
+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
+end
\ No newline at end of file
diff --git a/examples/models/DIS.jl b/src/models/DIS.jl
similarity index 100%
rename from examples/models/DIS.jl
rename to src/models/DIS.jl
diff --git a/src/models/LP.jl b/src/models/LP.jl
new file mode 100644
index 0000000..0931aa1
--- /dev/null
+++ b/src/models/LP.jl
@@ -0,0 +1,47 @@
+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
+    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
new file mode 100644
index 0000000..7843f35
--- /dev/null
+++ b/src/models/PC.jl
@@ -0,0 +1,18 @@
+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
+end
\ No newline at end of file
diff --git a/examples/models/SIM.jl b/src/models/SIM.jl
similarity index 73%
rename from examples/models/SIM.jl
rename to src/models/SIM.jl
index f3ab32e..9c44310 100644
--- a/examples/models/SIM.jl
+++ b/src/models/SIM.jl
@@ -13,4 +13,10 @@ SIM() = @model begin
     end
 end
 
-params = Dict(:θ => 0.2, :α_1 => 0.6, :α_2 => 0.4)
\ No newline at end of file
+function SIM(bool)
+    if bool
+        Dict(:θ => 0.2, :α_1 => 0.6, :α_2 => 0.4)
+    else
+        SIM()
+    end
+end
\ No newline at end of file
diff --git a/examples/models/SIM_stoch.jl b/src/models/SIM_stoch.jl
similarity index 95%
rename from examples/models/SIM_stoch.jl
rename to src/models/SIM_stoch.jl
index ad4d60c..67684b6 100644
--- a/examples/models/SIM_stoch.jl
+++ b/src/models/SIM_stoch.jl
@@ -1,7 +1,7 @@
 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
+    @parameters θ α_1 α_2 α_3
     @equations begin
         G = G_0 + u_G
         Y = C + G
@@ -19,6 +19,7 @@ params = Dict(
     :θ => 0.2, 
     :α_1 => 0.6,
     :α_2 => 0.4,
+    :α_3 => 0.1
 )
 
 # possible extension: C = f(C_e) + U_e, C_e = f(..., C_e[-1]), where C_e is not observable.