diff --git a/examples/basic_example.jl b/examples/basic_example.jl
index 025f5db..e5a494e 100644
--- a/examples/basic_example.jl
+++ b/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")
diff --git a/examples/benchmark_w_matlab.jl b/examples/benchmark_w_matlab.jl
index c0849bb..a44f7b5 100644
--- a/examples/benchmark_w_matlab.jl
+++ b/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")
\ No newline at end of file
diff --git a/examples/change_expectations.jl b/examples/change_expectations.jl
index 397d4e1..f38939e 100644
--- a/examples/change_expectations.jl
+++ b/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
diff --git a/examples/get_simulations.jl b/examples/get_simulations.jl
index 34ecaa5..9b6dc2e 100644
--- a/examples/get_simulations.jl
+++ b/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)
diff --git a/examples/multithreading_speedup.jl b/examples/multithreading_speedup.jl
index 26e7304..a8b96ec 100644
--- a/examples/multithreading_speedup.jl
+++ b/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?
diff --git a/examples/scenario_analysis_via_overload.jl b/examples/scenario_analysis_via_overload.jl
index 91e64cd..9d14206 100644
--- a/examples/scenario_analysis_via_overload.jl
+++ b/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)
diff --git a/examples/scenario_analysis_via_shock.jl b/examples/scenario_analysis_via_shock.jl
index c78240f..4d84da9 100644
--- a/examples/scenario_analysis_via_shock.jl
+++ b/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)
diff --git a/main.jl b/main.jl
index 28e6435..fe078f1 100644
--- a/main.jl
+++ b/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)
diff --git a/src/model_init/init.jl b/src/model_init/init.jl
index 7bd0ef2..557db5c 100644
--- a/src/model_init/init.jl
+++ b/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
diff --git a/src/model_init/init_aggregates.jl b/src/model_init/init_aggregates.jl
index 354fcdc..0332d53 100644
--- a/src/model_init/init_aggregates.jl
+++ b/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
\ No newline at end of file
diff --git a/src/model_init/init_banks.jl b/src/model_init/init_banks.jl
index 5c26136..9669fe4 100644
--- a/src/model_init/init_banks.jl
+++ b/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
\ No newline at end of file
diff --git a/src/model_init/init_firms.jl b/src/model_init/init_firms.jl
index fd52e9d..a2689fc 100644
--- a/src/model_init/init_firms.jl
+++ b/src/model_init/init_firms.jl
@@ -1,4 +1,20 @@
 
+"""
+    init_firms(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
+
+Initialize firms with given parameters and initial conditions.
+
+# Arguments
+- `parameters`: The parameters for initializing the firms.
+- `initial_conditions`: The initial conditions for the firms.
+- `typeInt`: (optional) The integer type to be used. Default is `Int64`.
+- `typeFloat`: (optional) The floating-point type to be used. Default is `Float64`.
+
+# Returns
+- firms::Firms: The initialized firms.
+- firms_args::Tuple: The arguments used to initialize the firms.
+
+"""
 function init_firms(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
 
     # unpacking useful parameters
@@ -97,51 +113,22 @@ function init_firms(parameters, initial_conditions; typeInt = Int64, typeFloat =
     K_h = K_H * Y_h # TODO: K_h[(H_W + H_inact + 1):(H_W + H_inact + I)]
     D_h = D_H * Y_h # TODO: D_h[(H_W + H_inact + 1):(H_W + H_inact + I)]
 
-    firms = Firms(
-        G_i,
-        alpha_bar_i,
-        beta_i,
-        kappa_i,
-        w_i,
-        w_bar_i,
-        delta_i,
-        tau_Y_i,
-        tau_K_i,
-        N_i,
-        Y_i,
-        Q_i,
-        Q_d_i,
-        P_i,
-        S_i,
-        K_i,
-        M_i,
-        L_i,
-        pi_bar_i,
-        D_i,
-        Pi_i,
-        V_i,
-        I_i,
-        E_i,
-        P_bar_i,
-        P_CF_i,
-        DS_i,
-        DM_i,
-        zeros(typeFloat, I), # DL_i
-        zeros(typeFloat, I), # DL_d_i
-        zeros(typeFloat, I), # K_e_i
-        zeros(typeFloat, I), # L_e_i
-        zeros(typeFloat, I), # Q_s_i
-        zeros(typeFloat, I), # I_d_i
-        zeros(typeFloat, I), # DM_d_i
-        zeros(typeInt, I),   # N_d_i
-        zeros(typeFloat, I), # Pi_e_i
-        Y_h,
-        C_d_h,
-        I_d_h,
-        C_h,
-        I_h,
-        K_h,
-        D_h,
-        )
-    return firms
+    # additional tracking variables initialised to zero
+    DL_i = zeros(typeFloat, I)
+    DL_d_i = zeros(typeFloat, I)
+    K_e_i = zeros(typeFloat, I)
+    L_e_i = zeros(typeFloat, I)
+    Q_s_i = zeros(typeFloat, I)
+    I_d_i = zeros(typeFloat, I)
+    DM_d_i = zeros(typeFloat, I)
+    N_d_i = zeros(typeInt, I)
+    Pi_e_i = zeros(typeFloat, I)
+
+
+    firms_args = (G_i, alpha_bar_i, beta_i, kappa_i, w_i, w_bar_i, delta_i, tau_Y_i, tau_K_i, N_i, Y_i, Q_i, Q_d_i, 
+                      P_i, S_i, K_i, M_i, L_i, pi_bar_i, D_i, Pi_i, V_i, I_i, E_i, P_bar_i, P_CF_i, DS_i, DM_i, DL_i, 
+                      DL_d_i, K_e_i, L_e_i, Q_s_i, I_d_i, DM_d_i, N_d_i, Pi_e_i, Y_h, C_d_h, I_d_h, C_h, I_h, K_h, D_h)
+
+    firms = Firms(firms_args...)
+    return firms, firms_args
 end
\ No newline at end of file
diff --git a/src/model_init/init_government.jl b/src/model_init/init_government.jl
index 941dd89..a422299 100644
--- a/src/model_init/init_government.jl
+++ b/src/model_init/init_government.jl
@@ -1,5 +1,21 @@
 
 
+"""
+    init_government(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
+
+Initialize the government agent.
+
+# Arguments
+- `parameters`: The parameters.
+- `initial_conditions`: The initial conditions.
+- `typeInt`: The integer type to be used (default: `Int64`).
+- `typeFloat`: The floating-point type to be used (default: `Float64`).
+
+# Returns
+- The initialized government model.
+- The arguments used to initialize the government model.
+
+"""
 function init_government(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
     alpha_G = parameters["alpha_G"]
     beta_G = parameters["beta_G"]
@@ -16,8 +32,10 @@ function init_government(parameters, initial_conditions; typeInt = Int64, typeFl
     C_j = zero(typeFloat)
     P_j = zero(typeFloat)
     Y_G = zero(typeFloat)
-    government =
-        Government(alpha_G, beta_G, sigma_G, Y_G, C_G[T_prime], L_G, sb_inact, sb_other, C_d_j, C_j, P_j)
+
+    gov_args = (alpha_G, beta_G, sigma_G, Y_G, C_G[T_prime], L_G, sb_inact, sb_other, C_d_j, C_j, P_j)
+
+    government = Government(gov_args...)
         
-    return government
+    return government, gov_args
 end
\ No newline at end of file
diff --git a/src/model_init/init_rest_of_the_world.jl b/src/model_init/init_rest_of_the_world.jl
index c996721..ee2e125 100644
--- a/src/model_init/init_rest_of_the_world.jl
+++ b/src/model_init/init_rest_of_the_world.jl
@@ -1,5 +1,21 @@
 
 
+"""
+    init_rotw(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
+
+Initialize the rest of the world (rotw) agent.
+
+# Arguments
+- `parameters`: The parameters.
+- `initial_conditions`: The initial conditions.
+- `typeInt`: The integer type to be used (default: `Int64`).
+- `typeFloat`: The floating-point type to be used (default: `Float64`).
+
+# Returns
+- rotw::RestOfTheWorld: The initialized rest of the world agent.
+- rotw_args::Tuple: The arguments used to initialize the rest of the world agent.
+
+"""
 function init_rotw(parameters, initial_conditions; typeInt = Int64, typeFloat = Float64)
     L = typeInt(parameters["L"])
     G = typeInt(parameters["G"])
@@ -34,7 +50,8 @@ function init_rotw(parameters, initial_conditions; typeInt = Int64, typeFloat =
     Q_m = Vector{typeFloat}(zeros(G))
     Q_d_m = Vector{typeFloat}(zeros(G))
     P_m = Vector{typeFloat}(zeros(G))
-    rotw = RestOfTheWorld(
+
+    rotw_args = (
         alpha_E,
         beta_E,
         sigma_E,
@@ -61,5 +78,8 @@ function init_rotw(parameters, initial_conditions; typeInt = Int64, typeFloat =
         P_m,
         P_l,
     )
-    return rotw
+
+    rotw = RestOfTheWorld(rotw_args...)
+
+    return rotw, rotw_args
 end
diff --git a/src/model_init/init_workers.jl b/src/model_init/init_workers.jl
index 9cc3de3..7118800 100644
--- a/src/model_init/init_workers.jl
+++ b/src/model_init/init_workers.jl
@@ -1,5 +1,25 @@
 
 
+"""
+    init_workers(parameters, initial_conditions, firms; typeInt = Int64, typeFloat = Float64)
+
+Initialize the workers for the given parameters, initial conditions, and firms.
+
+# Arguments
+- `parameters`: The parameters for the initialization.
+- `initial_conditions`: The initial conditions for the initialization.
+- `firms`: The already initialized firms.
+- `typeInt`: (optional) The type for integer values. Default is `Int64`.
+- `typeFloat`: (optional) The type for floating-point values. Default is `Float64`.
+
+# Returns
+- The initialized active workers.
+- The initialized inactive workers.
+- The updated firm vacancies V_i_new, this is needed to update the firms.
+- The arguments used to initialize the active workers.
+- The arguments used to initialize the inactive workers.
+
+"""
 function init_workers(parameters, initial_conditions, firms; typeInt = Int64, typeFloat = Float64)
  
     H_act = typeInt(parameters["H_act"])
@@ -15,11 +35,10 @@ function init_workers(parameters, initial_conditions, firms; typeInt = Int64, ty
     D_H = initial_conditions["D_H"]
     K_H = initial_conditions["K_H"]
 
-    
     H_W = H_act - I - 1
     P_bar_HH = one(typeFloat)
-    w_h = zeros(typeFloat, 1, H_W)
-    O_h = zeros(typeInt, 1, H_W)
+    w_h = zeros(typeFloat, H_W)
+    O_h = zeros(typeInt, H_W)
 
     h = one(typeInt)
 
@@ -37,7 +56,7 @@ function init_workers(parameters, initial_conditions, firms; typeInt = Int64, ty
 
     w_h[O_h .== 0] .= w_UB / theta_UB
     
-    Y_h = zeros(typeFloat, 1, H_W)
+    Y_h = zeros(typeFloat, H_W)
 
     for h in 1:H_W
         if O_h[h] != 0
@@ -55,8 +74,11 @@ function init_workers(parameters, initial_conditions, firms; typeInt = Int64, ty
     I_d_h = zeros(typeFloat, length(ids))
     C_h = zeros(typeFloat, length(ids))
     I_h = zeros(typeFloat, length(ids))
+
+    
     # active workers (both employed and unemployed)
-    workers_act = Workers(Y_h[1:H_W], D_h[1:H_W], K_h[1:H_W], w_h[1:H_W], O_h[1:H_W], C_d_h, I_d_h, C_h, I_h)
+    w_act_args = (Y_h, D_h, K_h, w_h, O_h, C_d_h, I_d_h, C_h, I_h)
+    workers_act = Workers(w_act_args...)
     
     # inactive workers
     ids = Vector{typeInt}((I + H_W + 1):(I + H_W + H_inact))
@@ -75,17 +97,9 @@ function init_workers(parameters, initial_conditions, firms; typeInt = Int64, ty
     I_d_h = zeros(typeFloat, length(ids))
     C_h = zeros(typeFloat, length(ids))
     I_h = zeros(typeFloat, length(ids))
-    workers_inact = Workers(
-        Y_h,
-        D_h,
-        K_h,
-        w_h_inact,
-        O_h_inact,
-        C_d_h,
-        I_d_h,
-        C_h,
-        I_h,
-    )
-
-    return workers_act, workers_inact, V_i_new
+
+    w_inact_args = (Y_h, D_h, K_h, w_h_inact, O_h_inact, C_d_h, I_d_h, C_h, I_h)
+    workers_inact = Workers(w_inact_args...)
+
+    return workers_act, workers_inact, V_i_new, w_act_args, w_inact_args
 end
diff --git a/src/one_simulation.jl b/src/one_simulation.jl
index 29c6a89..fa65608 100644
--- a/src/one_simulation.jl
+++ b/src/one_simulation.jl
@@ -11,7 +11,7 @@ The simulation runs for a number of epochs specified by `model.prop.T`.
 - `data::Data`: The data collected during the simulation.
 
 # Details
-The function initializes the data using `BeforeIT.initialise_data(model)`, then iteratively updates the model and data
+The function initializes the data using `BeforeIT.init_data(model)`, then iteratively updates the model and data
 for each epoch using `BeforeIT.one_epoch!(model)` and `BeforeIT.update_data!(data, model)` respectively.
 
 # Example
@@ -21,7 +21,7 @@ data = run_one_sim!(model)
 """
 function run_one_sim!(model; multi_threading = false, shock = NoShock())
 
-    data = BeforeIT.initialise_data(model)
+    data = BeforeIT.init_data(model)
 
     T = model.prop.T
 
diff --git a/src/utils/data.jl b/src/utils/data.jl
index 6946712..0d60d33 100644
--- a/src/utils/data.jl
+++ b/src/utils/data.jl
@@ -32,7 +32,7 @@ end
 """
 Initialise the data arrays
 """
-function initialise_data(m)
+function init_data(m)
     p = m.prop
     T = p.T
     d = Data([zeros(T + 1) for _ in 1:25]..., zeros(T + 1, p.G), zeros(T + 1, p.G))
diff --git a/test/accounting_identities.jl b/test/accounting_identities.jl
index f1c9033..b7312db 100644
--- a/test/accounting_identities.jl
+++ b/test/accounting_identities.jl
@@ -8,8 +8,8 @@ using Test
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
 
     T = 1
-    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
         BeforeIT.one_epoch!(model; multi_threading = false)
diff --git a/test/deterministic/deterministic_ouput_t1_t5.jl b/test/deterministic/deterministic_ouput_t1_t5.jl
index f76f96a..4d76339 100644
--- a/test/deterministic/deterministic_ouput_t1_t5.jl
+++ b/test/deterministic/deterministic_ouput_t1_t5.jl
@@ -4,8 +4,8 @@
     T = 1
     parameters = BeforeIT.AUSTRIA2010Q1.parameters
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
-    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)
 
     BeforeIT.one_epoch!(model; multi_threading = false)
     BeforeIT.update_data!(data, model)
@@ -30,8 +30,8 @@
     T = 5
     parameters = BeforeIT.AUSTRIA2010Q1.parameters
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
-    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
         BeforeIT.one_epoch!(model; multi_threading = false)
         BeforeIT.update_data!(data, model)
diff --git a/test/deterministic/initialize_deterministic.jl b/test/deterministic/initialize_deterministic.jl
index 7ab90a9..04330a5 100644
--- a/test/deterministic/initialize_deterministic.jl
+++ b/test/deterministic/initialize_deterministic.jl
@@ -4,7 +4,7 @@
 
     parameters = BeforeIT.AUSTRIA2010Q1.parameters
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
-    model = BeforeIT.initialise_model(parameters, initial_conditions, 1)
+    model = BeforeIT.init_model(parameters, initial_conditions, 1)
 
     properties = model.prop
 
diff --git a/test/deterministic/one_epoch_deterministic.jl b/test/deterministic/one_epoch_deterministic.jl
index 52441a7..12ccd9a 100644
--- a/test/deterministic/one_epoch_deterministic.jl
+++ b/test/deterministic/one_epoch_deterministic.jl
@@ -5,7 +5,7 @@
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
 
     T = 1
-    model = BeforeIT.initialise_model(parameters, initial_conditions, T;)
+    model = BeforeIT.init_model(parameters, initial_conditions, T;)
 
     gov = model.gov # government
     cb = model.cb # central bank
diff --git a/test/deterministic/one_run_deterministic.jl b/test/deterministic/one_run_deterministic.jl
index d891d7e..6f2a090 100644
--- a/test/deterministic/one_run_deterministic.jl
+++ b/test/deterministic/one_run_deterministic.jl
@@ -11,16 +11,16 @@
     initial_conditions1 = deepcopy(initial_conditions)
     initial_conditions2 = deepcopy(initial_conditions)
 
-    model = BeforeIT.initialise_model(parameters1, initial_conditions1, T;)
-    data = BeforeIT.initialise_data(model)
+    model = BeforeIT.init_model(parameters1, initial_conditions1, T;)
+    data = BeforeIT.init_data(model)
     for t in 1:(T - 1)
         BeforeIT.one_epoch!(model; multi_threading = false)
         BeforeIT.update_data!(data, model)
     end
 
 
-    model2 = BeforeIT.initialise_model(parameters2, initial_conditions2, T;)
-    data2 = BeforeIT.initialise_data(model2)
+    model2 = BeforeIT.init_model(parameters2, initial_conditions2, T;)
+    data2 = BeforeIT.init_data(model2)
 
     for t in 1:(T - 1)
         BeforeIT.one_epoch!(model2; multi_threading = false)
diff --git a/test/markets/search_and_matching.jl b/test/markets/search_and_matching.jl
index 4cc1520..8980690 100644
--- a/test/markets/search_and_matching.jl
+++ b/test/markets/search_and_matching.jl
@@ -10,7 +10,7 @@ using Random
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
 
     T = 1
-    model = BeforeIT.initialise_model(parameters, initial_conditions, T;)
+    model = BeforeIT.init_model(parameters, initial_conditions, T;)
 
 
     gov = model.gov # government
diff --git a/test/model_init/initialise_model.jl b/test/model_init/initialise_model.jl
index 5284882..a67c2d4 100644
--- a/test/model_init/initialise_model.jl
+++ b/test/model_init/initialise_model.jl
@@ -5,7 +5,7 @@ dir = @__DIR__
 
 parameters = BeforeIT.STEADY_STATE2010Q1.parameters
 initial_conditions = BeforeIT.STEADY_STATE2010Q1.initial_conditions
-model = BeforeIT.initialise_model(parameters, initial_conditions, 1)
+model = BeforeIT.init_model(parameters, initial_conditions, 1)
 
 properties = model.prop
 
diff --git a/test/monte_carlo_evaluations.jl b/test/monte_carlo_evaluations.jl
index 311ca41..47a25af 100644
--- a/test/monte_carlo_evaluations.jl
+++ b/test/monte_carlo_evaluations.jl
@@ -6,8 +6,8 @@ parameters = matread(joinpath(dir, "../data/austria/parameters/2010Q1.mat"))
 initial_conditions = matread(joinpath(dir, "../data/austria/initial_conditions/2010Q1.mat"))
 
 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)
 
 n_sims = 3
 data_vector = BeforeIT.run_n_sims(model, n_sims)
diff --git a/test/one_epoch_consistency.jl b/test/one_epoch_consistency.jl
index 7825b7a..b61b549 100644
--- a/test/one_epoch_consistency.jl
+++ b/test/one_epoch_consistency.jl
@@ -9,8 +9,8 @@ parameters = matread(joinpath(dir, "../data/steady_state/parameters/2010Q1.mat")
 initial_conditions = matread(joinpath(dir, "../data/steady_state/initial_conditions/2010Q1.mat"))
 
 T = 1
-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)
 
 
 println(BeforeIT.get_accounting_identity_banks(model))
diff --git a/test/steady_state.jl b/test/steady_state.jl
index 3635fea..2662e72 100644
--- a/test/steady_state.jl
+++ b/test/steady_state.jl
@@ -5,8 +5,8 @@ parameters = BeforeIT.STEADY_STATE2010Q1.parameters
 initial_conditions = BeforeIT.STEADY_STATE2010Q1.initial_conditions
 
 T = 5
-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(t)
diff --git a/test/utils/estimations.jl b/test/utils/estimations.jl
index 716a7ac..b8f3223 100644
--- a/test/utils/estimations.jl
+++ b/test/utils/estimations.jl
@@ -6,7 +6,7 @@ using MAT, Random
 
     parameters = BeforeIT.AUSTRIA2010Q1.parameters
     initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
-    model = BeforeIT.initialise_model(parameters, initial_conditions, 1)
+    model = BeforeIT.init_model(parameters, initial_conditions, 1)