Better initialisation

examples/basic_example.jl

@@ -15,7 +15,7 @@ initial_conditions = Bit.AUSTRIA2010Q1.initial_conditions
 # We can now initialise the model, by specifying in advance the maximum number of epochs.
 
 T = 16
-model = Bit.initialise_model(parameters, initial_conditions, T)
+model = Bit.init_model(parameters, initial_conditions, T)
 
 
 # Note that the it is very simple to inspect the model by typing
@@ -28,7 +28,7 @@ fieldnames(typeof(model.bank))
 
 # We can now define a data tracker, which will store the time series of the model.
 
-data = Bit.initialise_data(model);
+data = Bit.init_data(model);
 
 # We can run now the model for a number of epochs and progressively update the data tracker.
 
@@ -53,11 +53,12 @@ p7 = plot(data.wages, title = "wages", titlefont = 10)
 p8 = plot(data.euribor, title = "euribor", titlefont = 10)
 p9 = plot(data.nominal_gdp ./ data.real_gdp, title = "gdp deflator", titlefont = 10)
 
+
 plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, layout = (3, 3), legend = false)
 
 # To run multiple monte-carlo repetitions in parallel we can use
 
-model = Bit.initialise_model(parameters, initial_conditions, T)
+model = Bit.init_model(parameters, initial_conditions, T)
 data_vector = Bit.run_n_sims(model, 4)
 
 # Note that this will use the number of threads specified when activating the Julia environment.
@@ -96,3 +97,10 @@ p9 = errorline(
     titlefont = 10,
 )
 plot(p1, p2, p3, p4, p5, p6, p7, p8, p9, layout = (3, 3), legend = false)
+
+plot(p1, p4, p5, p3, p8, p9, layout = (3, 2), legend = false, size = (400, 600), dpi = 300, left_margin = 3Plots.mm) 
+
+plot(p1, p4, p5, p3, p8, p9, layout = (2, 3), legend = false, size = (600, 400), dpi = 300)#, left_margin = 3Plots.mm) 
+
+
+savefig("output.png")

examples/benchmark_w_matlab.jl

@@ -11,8 +11,8 @@ T = 12*2
 
 # We will run the model without any output to avoid the overhead of printing the results.
 function run_no_output(;multi_threading = false)
-    model = BeforeIT.initialise_model(parameters, initial_conditions, T)
-    data = BeforeIT.initialise_data(model);
+    model = BeforeIT.init_model(parameters, initial_conditions, T)
+    data = BeforeIT.init_data(model);
 
     for _ in 1:T
         BeforeIT.one_epoch!(model; multi_threading = multi_threading)
@@ -25,7 +25,7 @@ end
 @time run_no_output(;multi_threading = true) 
 
 # time taken by the MATLAB code, computed independently on an Apple M1 chip
-matlab_times = [3.1465, 3.0795, 3.0274, 3.0345, 3.0873]
+matlab_times = [3.1919, 3.2454, 3.1501, 3.1074, 3.1551]
 matlab_time = mean(matlab_times)
 matlab_time_std = std(matlab_times)
 
@@ -50,11 +50,16 @@ julia_time_multi_thread_std = std(julia_times_multi_thread)
 # get the number of threads used
 n_threads = Threads.nthreads()
 
+theme(:default, bg = :white)
 # bar chart of the time taken vs the time taken by the MATLAB code, also plot the stds as error bars
-bar(["MATLAB", "Julia, 1 thread", "Julia, $n_threads threads"], [matlab_time, julia_time_1_thread, julia_time_multi_thread], yerr = [matlab_time_std, julia_time_1_thread_std, julia_time_multi_thread_std], legend = false, dpi=300)
-ylabel!("Time for one simulation (s)")
+# make a white background with no grid
+bar(["MATLAB", "Julia, 1 thread", "Julia, $n_threads threads"], [matlab_time, julia_time_1_thread, julia_time_multi_thread], 
+yerr = [matlab_time_std, julia_time_1_thread_std, julia_time_multi_thread_std],
+legend = false, dpi=300, size=(400, 300), grid = false, ylabel = "Time for one simulation (s)")
 
 # the Julia implementation is faster than the MATLAB implementation, and the multi-threaded version is faster than the single-threaded version.
 
+# increase
+
 # save the image
 savefig("benchmark_w_matlab.png")

examples/change_expectations.jl

@@ -15,7 +15,7 @@ init = Bit.AUSTRIA2010Q1.initial_conditions
 
 Random.seed!(1234)
 T = 40
-model = Bit.initialise_model(par, init, T)
+model = Bit.init_model(par, init, T)
 data = Bit.run_one_sim!(model)
 
 # Now we can experiment with changing expectations of the agents in the model.
@@ -30,7 +30,7 @@ end
 # run the model again, with the same seed
 
 Random.seed!(1234)
-model = Bit.initialise_model(par, init, T)
+model = Bit.init_model(par, init, T)
 data_back = Bit.run_one_sim!(model)
 
 # plot the results, comparing the two cases as different lines

examples/get_simulations.jl

@@ -12,7 +12,7 @@ for year in 2010:2019
         parameters = load("data/italy/parameters/" * string(year) * "Q" * string(quarter) * ".jld2")
         initial_conditions = load("data/italy/initial_conditions/" * string(year) * "Q" * string(quarter) * ".jld2")
         T = 12
-        model = BeforeIT.initialise_model(parameters, initial_conditions, T)
+        model = BeforeIT.init_model(parameters, initial_conditions, T)
         n_sims = 4
         data_vector = BeforeIT.run_n_sims(model, n_sims)
         save("data/italy/simulations/" * string(year) * "Q" * string(quarter) * ".jld2", "data_vector", data_vector)

examples/multithreading_speedup.jl

@@ -13,7 +13,7 @@ using FileIO, Plots, StatsPlots
 parameters = Bit.ITALY2010Q1.parameters
 initial_conditions = Bit.ITALY2010Q1.initial_conditions
 T = 50
-model = Bit.initialise_model(parameters, initial_conditions, T);
+model = Bit.init_model(parameters, initial_conditions, T);
 
 
 # The model is in scale 1:2000, so it has around 30,000 households
@@ -31,7 +31,7 @@ println(Threads.nthreads())
 
 @time data = Bit.run_one_sim!(model; multi_threading = false);
 
-model = Bit.initialise_model(parameters, initial_conditions, T);
+model = Bit.init_model(parameters, initial_conditions, T);
 @time data = Bit.run_one_sim!(model; multi_threading = true);
 
 # Is the speedup in line to what we would expect?

examples/scenario_analysis_via_overload.jl

@@ -12,8 +12,8 @@ initial_conditions = Bit.AUSTRIA2010Q1.initial_conditions
 
 # initialise the model and the data collector
 T = 20
-model = Bit.initialise_model(parameters, initial_conditions, T);
-data = Bit.initialise_data(model);
+model = Bit.init_model(parameters, initial_conditions, T);
+data = Bit.init_data(model);
 
 # Simulate the model for T quarters
 data_vec_baseline = Bit.run_n_sims(model, 4)

examples/scenario_analysis_via_shock.jl

@@ -12,7 +12,7 @@ initial_conditions = Bit.AUSTRIA2010Q1.initial_conditions
 
 # initialise the model and the data collector
 T = 20
-model = Bit.initialise_model(parameters, initial_conditions, T);
+model = Bit.init_model(parameters, initial_conditions, T);
 
 # Simulate the model for T quarters
 data_vec_baseline = Bit.run_n_sims(model, 4)

main.jl

@@ -5,8 +5,8 @@ parameters = BeforeIT.AUSTRIA2010Q1.parameters
 initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
 
 T = 20
-model = BeforeIT.initialise_model(parameters, initial_conditions, T)
-data = BeforeIT.initialise_data(model)
+model = BeforeIT.init_model(parameters, initial_conditions, T)
+data = BeforeIT.init_data(model)
 
 for t in 1:T
     println("Epoch: ", t)

src/model_init/init.jl

@@ -3,7 +3,7 @@ recursive_namedtuple(x::Any) = x
 recursive_namedtuple(d::Dict) = MutableNamedTuple(; Dict(k => recursive_namedtuple(v) for (k, v) in d)...)
 
 """
-    initialise_model(parameters, initial_conditions, T, typeInt = Int64, typeFloat = Float64)
+    init_model(parameters, initial_conditions, T, typeInt = Int64, typeFloat = Float64)
 
 Initializes the model with given parameters and initial conditions.
 
@@ -18,54 +18,72 @@ Returns:
 - model::Model: The initialized model.
 
 """
-function initialise_model(parameters::Dict{String, Any}, initial_conditions::Dict{String, Any}, T, typeInt::DataType = Int64, typeFloat::DataType = Float64)
+function init_model(parameters::Dict{String, Any}, initial_conditions::Dict{String, Any}, T, typeInt::DataType = Int64, typeFloat::DataType = Float64)
 
     # properties
     properties = BeforeIT.init_properties(parameters, T; typeInt = typeInt, typeFloat = typeFloat)
 
     # firms
-    firms = BeforeIT.init_firms(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
+    firms, _ = BeforeIT.init_firms(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
 
     # workers, and update firms vacancies
-    workers_act, workers_inact, V_i_new = BeforeIT.init_workers(parameters, initial_conditions, firms; typeInt = typeInt, typeFloat = typeFloat)
+    workers_act, workers_inact, V_i_new, _, _ = BeforeIT.init_workers(parameters, initial_conditions, firms; typeInt = typeInt, typeFloat = typeFloat)
     firms.V_i = V_i_new
 
     # bank
-    bank = BeforeIT.init_bank(parameters, initial_conditions, firms; typeInt = typeInt, typeFloat = typeFloat)
+    bank, _ = BeforeIT.init_bank(parameters, initial_conditions, firms; typeInt = typeInt, typeFloat = typeFloat)
 
     # central bank
-    central_bank = BeforeIT.init_central_bank(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
+    central_bank, _ = BeforeIT.init_central_bank(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
 
     # government
-    government = BeforeIT.init_government(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
+    government, _ = BeforeIT.init_government(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
 
     # rest of the world
-    rotw = BeforeIT.init_rotw(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
+    rotw, _ = BeforeIT.init_rotw(parameters, initial_conditions; typeInt = typeInt, typeFloat = typeFloat)
 
     # aggregates
-    agg = BeforeIT.init_aggregates(parameters, initial_conditions, T; typeInt = typeInt, typeFloat = typeFloat)
-
-    # obtain total income by summing contributions from firm owners, workers and bank owner
-
-    tot_Y_h = sum(firms.Y_h) + sum(workers_act.Y_h) + sum(workers_inact.Y_h) + bank.Y_h
-
-    # uptade K_h and D_h in all agent types
-    firms.K_h .= firms.K_h / tot_Y_h 
-    firms.D_h .= firms.D_h / tot_Y_h
-    workers_act.K_h .= workers_act.K_h / tot_Y_h
-    workers_act.D_h .= workers_act.D_h / tot_Y_h
-    workers_inact.K_h .= workers_inact.K_h / tot_Y_h
-    workers_inact.D_h .= workers_inact.D_h / tot_Y_h
-    bank.K_h = bank.K_h / tot_Y_h
-    bank.D_h = bank.D_h / tot_Y_h
-
-    # get total deposits and update bank balance sheet
-    tot_D_h = sum(firms.D_h) + sum(workers_act.D_h) + sum(workers_inact.D_h) + bank.D_h
-    bank.D_k += tot_D_h
+    agg, _ = BeforeIT.init_aggregates(parameters, initial_conditions, T; typeInt = typeInt, typeFloat = typeFloat)
 
     # model
     model = Model(workers_act, workers_inact, firms, bank, central_bank, government, rotw, agg, properties)
 
+    # update the model with global quantities (total income, total deposits) obtained from all the agents
+    update_variables_with_totals!(model)
+
     return model
 
 end
+
+"""
+    update_variables_with_totals!(model::Model)
+
+Update the variables in the given `model` with some global quantities obtained from all agents.
+This is the last step in the initialization process and it must be performed after all agents have been initialized.
+
+# Arguments
+- `model::Model`: The model object to update.
+
+# Returns
+- Nothing
+
+"""
+function update_variables_with_totals!(model::Model)
+
+    # obtain total income by summing contributions from firm owners, workers and bank owner
+    tot_Y_h = sum(model.firms.Y_h) + sum(model.w_act.Y_h) + sum(model.w_inact.Y_h) + model.bank.Y_h
+
+    # uptade K_h and D_h in all agent types using total income  
+    model.firms.K_h .= model.firms.K_h / tot_Y_h 
+    model.firms.D_h .= model.firms.D_h / tot_Y_h
+    model.w_act.K_h .= model.w_act.K_h / tot_Y_h
+    model.w_act.D_h .= model.w_act.D_h / tot_Y_h
+    model.w_inact.K_h .= model.w_inact.K_h / tot_Y_h
+    model.w_inact.D_h .= model.w_inact.D_h / tot_Y_h
+    model.bank.K_h = model.bank.K_h / tot_Y_h
+    model.bank.D_h = model.bank.D_h / tot_Y_h
+
+    # get total deposits and update bank balance sheet
+    tot_D_h = sum(model.firms.D_h) + sum(model.w_act.D_h) + sum(model.w_inact.D_h) + model.bank.D_h
+    model.bank.D_k += tot_D_h
+end

src/model_init/init_aggregates.jl

@@ -1,5 +1,22 @@
 
 
+"""
+    init_aggregates(parameters, initial_conditions, T; typeInt = Int64, typeFloat = Float64)
+
+Initialize aggregates for the model.
+
+# Arguments
+- `parameters`: The model parameters.
+- `initial_conditions`: The initial conditions.
+- `T`: The total simulation time.
+- `typeInt`: The integer type to use (default: `Int64`).
+- `typeFloat`: The floating-point type to use (default: `Float64`).
+
+# Returns
+- `agg`: The initialized aggregates.
+- `agg_args`: The arguments used to initialize the aggregates.
+
+"""
 function init_aggregates(parameters, initial_conditions, T; typeInt = Int64, typeFloat = Float64)
 
     Y = initial_conditions["Y"]
@@ -25,7 +42,7 @@ function init_aggregates(parameters, initial_conditions, T; typeInt = Int64, typ
     epsilon_E = zero(typeFloat)
     epsilon_I = zero(typeFloat)
 
-    agg = Aggregates(
+    agg_args = (
         Y,
         pi_,
         P_bar,
@@ -42,5 +59,8 @@ function init_aggregates(parameters, initial_conditions, T; typeInt = Int64, typ
         epsilon_I,
         t,
     )
-    return agg
+
+    agg = Aggregates(agg_args...)
+
+    return agg, agg_args
 end

src/model_init/init_banks.jl

@@ -1,5 +1,22 @@
 
 
+"""
+    init_bank(parameters, initial_conditions, firms; typeInt = Int64, typeFloat = Float64)
+
+Initialize a bank with the given parameters, initial conditions, and firms.
+
+# Arguments
+- `parameters`: The parameters.
+- `initial_conditions`: The initial conditions.
+- `firms`: The already initialized firms.
+- `typeInt`: (optional) The integer type to use. Default is `Int64`.
+- `typeFloat`: (optional) The floating-point type to use. Default is `Float64`.
+
+# Returns
+- bank::Bank: The initialized bank.
+- bank_args::Tuple: The arguments used to initialize the bank.
+
+"""
 function init_bank(parameters, initial_conditions, firms; typeInt = Int64, typeFloat = Float64)
 
     theta_DIV = parameters["theta_DIV"]
@@ -32,12 +49,30 @@ function init_bank(parameters, initial_conditions, firms; typeInt = Int64, typeF
     K_h = K_h
     D_h = D_h
     Pi_e_k = typeFloat(0.0)
-    bank = Bank(E_k, Pi_k, Pi_e_k, D_k, r, Y_h, C_d_h, I_d_h, C_h, I_h, K_h, D_h)
 
-    return bank
+    bank_args = (E_k, Pi_k, Pi_e_k, D_k, r, Y_h, C_d_h, I_d_h, C_h, I_h, K_h, D_h)
+    bank = Bank(bank_args...)
+
+    return bank, bank_args
 end
 
 
+"""
+    init_central_bank(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
+
+Initialize the central bank with the given parameters and initial conditions.
+
+# Arguments
+- `parameters`: The parameters.
+- `initial_conditions`: The initial conditions.
+- `typeInt`: (optional) The integer type to be used. Default is `Int64`.
+- `typeFloat`: (optional) The floating-point type to be used. Default is `Float64`.
+
+# Returns
+- central_bank::CentralBank: The initialized central bank.
+- cb_args::Tuple: The arguments used to initialize the central bank.
+
+"""
 function init_central_bank(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
     r_bar = initial_conditions["r_bar"]
     r_G = parameters["r_G"]
@@ -48,6 +83,8 @@ function init_central_bank(parameters, initial_conditions; typeInt = Int64, type
     xi_gamma = parameters["xi_gamma"]
     E_CB = initial_conditions["E_CB"]
 
-    central_bank = CentralBank(r_bar, r_G, rho, r_star, pi_star, xi_pi, xi_gamma, E_CB)
-return central_bank
+    cb_args = (r_bar, r_G, rho, r_star, pi_star, xi_pi, xi_gamma, E_CB)
+    central_bank = CentralBank(cb_args...)
+
+    return central_bank, cb_args
 end