Merge branch 'main' into compathelper/new_version/2023-10-11-01-04-44…

…-475-03462036025

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]

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.

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

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

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
+    )

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

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

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

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

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)

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}

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

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