3rd order pullback improvements

README.md

@@ -236,11 +236,13 @@ The package contains the following models in the `models` folder:
 - [Ascari and Sbordone (2014)](https://www.aeaweb.org/articles?id=10.1257/jel.52.3.679) `Ascari_Sbordone_2014.jl`
 - [Backus, Kehoe, and Kydland (1992)](https://www.jstor.org/stable/2138686) `Backus_Kehoe_Kydland_1992.jl`
 - [Baxter and King (1993)](https://www.jstor.org/stable/2117521) `Baxter_King_1993.jl`
+- [Bernanke, Gertler, and Gilchrist (1999)](https://doi.org/10.1016/S1574-0048(99)10034-X) `Bernanke_Gertler_Gilchrist_1999.jl`
 - [Caldara et al. (2012)](https://www.sciencedirect.com/science/article/abs/pii/S1094202511000433) `Caldara_et_al_2012.jl`
 - [Gali (2015)](https://press.princeton.edu/books/hardcover/9780691164786/monetary-policy-inflation-and-the-business-cycle) - Chapter 3 `Gali_2015_chapter_3_nonlinear.jl`
 - [Gali and Monacelli (2005)](https://crei.cat/wp-content/uploads/users/pages/roes8739.pdf) - CPI inflation-based Taylor rule `Gali_Monacelli_2005_CITR.jl`
 - [Federal Reserve Board U.S. model (FRB/US) - LINVER (2024))](https://www.federalreserve.gov/econres/us-models-about.htm) `FRBUS.jl`
 - [Gerali, Neri, Sessa, and Signoretti (2010)](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1538-4616.2010.00331.x) `GNSS_2010.jl`
+- [Gertler and Karadi (2011)](https://doi.org/10.1016/j.jmoneco.2010.10.004) `Gertler_Karadi_2011.jl`
 - [Ghironi and Melitz (2005)](https://faculty.washington.edu/ghiro/GhiroMeliQJE0805.pdf) `Ghironi_Melitz_2005.jl`
 - [Ireland (2004)](http://irelandp.com/pubs/tshocksnk.pdf) `Ireland_2004.jl`
 - [Jermann and Quadrini (2012)](https://www.aeaweb.org/articles?id=10.1257/aer.102.1.238) - RBC `JQ_2012_RBC.jl`

docs/src/index.md

@@ -38,11 +38,13 @@ The package contains the following models in the `models` folder:
 - [Ascari and Sbordone (2014)](https://www.aeaweb.org/articles?id=10.1257/jel.52.3.679) `Ascari_Sbordone_2014.jl`
 - [Backus, Kehoe, and Kydland (1992)](https://www.jstor.org/stable/2138686) `Backus_Kehoe_Kydland_1992.jl`
 - [Baxter and King (1993)](https://www.jstor.org/stable/2117521) `Baxter_King_1993.jl`
+- [Bernanke, Gertler, and Gilchrist (1999)](https://doi.org/10.1016/S1574-0048(99)10034-X) `Bernanke_Gertler_Gilchrist_1999.jl`
 - [Caldara et al. (2012)](https://www.sciencedirect.com/science/article/abs/pii/S1094202511000433) `Caldara_et_al_2012.jl`
 - [Gali (2015)](https://press.princeton.edu/books/hardcover/9780691164786/monetary-policy-inflation-and-the-business-cycle) - Chapter 3 `Gali_2015_chapter_3_nonlinear.jl`
 - [Gali and Monacelli (2005)](https://crei.cat/wp-content/uploads/users/pages/roes8739.pdf) - CPI inflation-based Taylor rule `Gali_Monacelli_2005_CITR.jl`
 - [Federal Reserve Board U.S. model (FRB/US) - LINVER (2024))](https://www.federalreserve.gov/econres/us-models-about.htm) `FRBUS.jl`
 - [Gerali, Neri, Sessa, and Signoretti (2010)](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1538-4616.2010.00331.x) `GNSS_2010.jl`
+- [Gertler and Karadi (2011)](https://doi.org/10.1016/j.jmoneco.2010.10.004) `Gertler_Karadi_2011.jl`
 - [Ghironi and Melitz (2005)](https://faculty.washington.edu/ghiro/GhiroMeliQJE0805.pdf) `Ghironi_Melitz_2005.jl`
 - [Ireland (2004)](http://irelandp.com/pubs/tshocksnk.pdf) `Ireland_2004.jl`
 - [Jermann and Quadrini (2012)](https://www.aeaweb.org/articles?id=10.1257/aer.102.1.238) - RBC `JQ_2012_RBC.jl`

docs/src/unfinished_docs/todo.md

@@ -52,7 +52,6 @@
 - [ ] implement forwarddiff for find_shocks
 - [ ] redo inversion filter 1st order rrule based on the higher order ones. the accumulated matmul might not be necessary at all
 - [ ] inversion filter: use subset of observables and states when propagating states (see kalman filter)
-- [x] start filter from initial values provided by user
 - [ ] higher order estimation should start from mean not the stochastic steady state as the mean is the most likely starting point
 - [ ] large models will need functions to be compiled individually as done for higher order; when tackling that, also separate steady state related equations from the steady state, so that speed issue is addresses due to replacing parameters with the steady state equations from the parameter block; also creat non allocating (residuals) steady state function
 - [ ] check tols throughout. adopt max(abs,rel*norm) tols
@@ -73,9 +72,6 @@
 - [ ] do proper testing of ss solver with random set of params, equal across configs
 - [ ] load create parts of derivatives later and not directly after parameters block
 - [ ] fix model estimate plot. data not above estimate (should be red but is blue)
-- [x] implement higher order (pruned) variance decomposition
-- [x] add `marginal_contribution` (Shapley) option to `get_variance_decomposition` for pruned higher-order solutions
-- [x] add `marginal_contribution` (Shapley) option to `get_shock_decomposition`/`plot_shock_decomposition` for pruned higher-order solutions (allocates the `Nonlinearities` term across shocks)
 - [ ] try slicesampler instead of pigeons
 - [ ] speed up sensitivity by caching matrix inversion from implicit diff with LRUcache
 - [ ] fix this inference errors for large functions. they are slow. fix derivatives in general.
@@ -103,9 +99,7 @@
 - [ ] functions to reverse state_update (input: previous shock and current state, output previous state), find shocks corresponding to bringing one state to the next
 - [ ] cover nested case: min(50,a+b+max(c,10))
 - [ ] add balanced growth path handling
-- [x] autocorr and corr with derivatives. return 3d array
 - [ ] add pydsge and econpizza to overview
-- [ ] add for loop parser in @parameters
 - [ ] implement more multi country models
 - [ ] speed benchmarking (focus on ImplicitDiff part)
 - [ ] for cond forecasting allow less shocks than conditions with a warning. should be svd then
@@ -129,6 +123,12 @@
 - [ ] figure out combinations for inputs (parameters and variables in different formats for get_irf for example)
 - [ ] weed out SS solver and saved objects
 
+- [x] start filter from initial values provided by user
+- [x] implement higher order (pruned) variance decomposition
+- [x] add `marginal_contribution` (Shapley) option to `get_variance_decomposition` for pruned higher-order solutions
+- [x] add `marginal_contribution` (Shapley) option to `get_shock_decomposition`/`plot_shock_decomposition` for pruned higher-order solutions (allocates the `Nonlinearities` term across shocks)
+- [x] autocorr and corr with derivatives. return 3d array
+- [x] add for loop parser in @parameters
 - [x] add filter free estimation to one of the estimation tests (think about filter free filtering instead of smoothing)
 - [x] check dispatch doctor coverage in new functions
 - [x] dont use variable, parameter, function names starting with _

ext/MooncakeExt.jl

@@ -187,6 +187,7 @@ function increment_nested_initial_state_rdata!(
     cr_tangent,
 )
     cr_tangent isa ChainRulesCore.AbstractZero && return NoRData()
+    cr_tangent = ChainRulesCore.unthunk(cr_tangent)
     initial_state_fdata = Mooncake.tangent(initial_state_cd)
     @inbounds for i in eachindex(cr_tangent)
         component_tangent = cr_tangent[i]
@@ -450,4 +451,4 @@ function Mooncake.rrule!!(
     return CoDual(nothing, Mooncake.NoFData()), pb!!
 end
 
-end  # module MooncakeExt
+end  # module MooncakeExt

ext/StatsPlotsExt.jl

@@ -2,7 +2,7 @@ module StatsPlotsExt
 
 using MacroModelling
 
-import MacroModelling: ParameterType, ℳ, Symbol_input, String_input, Tolerances, NsssTolerances, SolverTolerances, merge_calculation_options, MODEL®, DATA®, PARAMETERS®, ALGORITHM®, FILTER®, VARIABLES®, SMOOTH®, SHOW_PLOTS®, SAVE_PLOTS®, SAVE_PLOTS_NAME®, SAVE_PLOTS_FORMAT®, SAVE_PLOTS_PATH®, PLOTS_PER_PAGE®, MAX_ELEMENTS_PER_LEGENDS_ROW®, EXTRA_LEGEND_SPACE®, PLOT_ATTRIBUTES®, QME®, SYLVESTER®, LYAPUNOV®, TOLERANCES®, VERBOSE®, DATA_IN_LEVELS®, PERIODS®, SHOCKS®, SHOCK_SIZE®, NEGATIVE_SHOCK®, GENERALISED_IRF®, GENERALISED_IRF_WARMUP_ITERATIONS®, CONDITIONS_IN_LEVELS®, GENERALISED_IRF_DRAWS®, INITIAL_STATE®, IGNORE_OBC®, CONDITIONS®, SHOCK_CONDITIONS®, LEVELS®, LABEL®, RENAME_DICTIONARY®, STEADY_STATE_FUNCTION®, parse_shocks_input_to_index, parse_variables_input_to_index, replace_indices, replace_indices_special, filter_data_with_model, get_relevant_steady_states, replace_indices_in_symbol, parse_algorithm_to_state_update, girf, decompose_name, obc_objective_optim_fun, obc_constraint_optim_fun, compute_irf_responses, process_ignore_obc_flag, adjust_generalised_irf_flag, process_shocks_input, normalize_filtering_options, normalize_presample_periods, trim_informative_sample, adjust_initial_state, SteadyStateFunctionType
+import MacroModelling: ParameterType, ℳ, Symbol_input, String_input, Tolerances, NsssTolerances, SolverTolerances, merge_calculation_options, MODEL®, DATA®, PARAMETERS®, ALGORITHM®, FILTER®, VARIABLES®, SMOOTH®, SHOW_PLOTS®, SAVE_PLOTS®, SAVE_PLOTS_NAME®, SAVE_PLOTS_FORMAT®, SAVE_PLOTS_PATH®, PLOTS_PER_PAGE®, MAX_ELEMENTS_PER_LEGENDS_ROW®, EXTRA_LEGEND_SPACE®, PLOT_ATTRIBUTES®, QME®, SYLVESTER®, LYAPUNOV®, TOLERANCES®, VERBOSE®, DATA_IN_LEVELS®, PERIODS®, SHOCKS®, SHOCK_SIZE®, NEGATIVE_SHOCK®, GENERALISED_IRF®, GENERALISED_IRF_WARMUP_ITERATIONS®, CONDITIONS_IN_LEVELS®, GENERALISED_IRF_DRAWS®, INITIAL_STATE®, IGNORE_OBC®, CONDITIONS®, SHOCK_CONDITIONS®, LEVELS®, LABEL®, RENAME_DICTIONARY®, STEADY_STATE_FUNCTION®, parse_shocks_input_to_index, parse_variables_input_to_index, replace_indices, replace_indices_special, filter_data_with_model, get_relevant_steady_states, replace_indices_in_symbol, parse_algorithm_to_state_update, girf, decompose_name, obc_objective_optim_fun, obc_constraint_optim_fun, compute_irf_responses, process_ignore_obc_flag, adjust_generalised_irf_flag, process_shocks_input, normalize_filtering_options, normalize_presample_periods, trim_informative_sample, adjust_initial_state, SteadyStateFunctionType, normalize_superscript
 import MacroModelling: DEFAULT_CACHING, DEFAULT_USE_WORKSPACES, DEFAULT_ALGORITHM, DEFAULT_FILTER_SELECTOR, DEFAULT_WARMUP_ITERATIONS, DEFAULT_VARIABLES_EXCLUDING_OBC, DEFAULT_SHOCK_SELECTION, DEFAULT_PRESAMPLE_PERIODS, DEFAULT_DATA_IN_LEVELS, DEFAULT_SHOCK_DECOMPOSITION_SELECTOR, DEFAULT_SMOOTH_SELECTOR, DEFAULT_LABEL, DEFAULT_SHOW_PLOTS, DEFAULT_SAVE_PLOTS, DEFAULT_SAVE_PLOTS_FORMAT, DEFAULT_SAVE_PLOTS_PATH, DEFAULT_PLOTS_PER_PAGE_SMALL, DEFAULT_TRANSPARENCY, DEFAULT_MAX_ELEMENTS_PER_LEGEND_ROW, DEFAULT_EXTRA_LEGEND_SPACE, DEFAULT_VERBOSE, DEFAULT_QME_ALGORITHM, DEFAULT_SYLVESTER_SELECTOR, DEFAULT_SYLVESTER_THRESHOLD, DEFAULT_LARGE_SYLVESTER_ALGORITHM, DEFAULT_SYLVESTER_ALGORITHM, DEFAULT_LYAPUNOV_ALGORITHM, DEFAULT_PLOT_ATTRIBUTES, DEFAULT_ARGS_AND_KWARGS_NAMES, DEFAULT_PLOTS_PER_PAGE_LARGE, DEFAULT_SHOCKS_EXCLUDING_OBC, DEFAULT_VARIABLES_EXCLUDING_AUX_AND_OBC, DEFAULT_PERIODS, DEFAULT_SHOCK_SIZE, DEFAULT_NEGATIVE_SHOCK, DEFAULT_GENERALISED_IRF, DEFAULT_GENERALISED_IRF_WARMUP, DEFAULT_GENERALISED_IRF_DRAWS, DEFAULT_INITIAL_STATE, DEFAULT_IGNORE_OBC, DEFAULT_PLOT_TYPE, DEFAULT_CONDITIONS_IN_LEVELS, DEFAULT_SIGMA_RANGE, DEFAULT_FONT_SIZE, DEFAULT_VARIABLE_SELECTION, DEFAULT_FORECAST_PERIODS
 import DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF, TYPEDSIGNATURES, TYPEDFIELDS
 import LaTeXStrings
@@ -116,44 +116,6 @@ function apply_custom_name(symbol::R, custom_names::AbstractDict{S, T})::R where
     return symbol
 end
 
-function normalize_superscript(x::Symbol)
-    return normalize_superscript(string(x))
-end
-
-function normalize_superscript(x::AbstractString)
-    sub_map = Dict(
-        '₀' => '0', '₁' => '1', '₂' => '2', '₃' => '3', '₄' => '4',
-        '₅' => '5', '₆' => '6', '₇' => '7', '₈' => '8', '₉' => '9',
-        '₊' => '+', '₋' => '-', '₌' => '=', '₍' => '(', '₎' => ')',
-        'ₐ' => 'a', 'ₑ' => 'e', 'ₕ' => 'h', 'ᵢ' => 'i', 'ⱼ' => 'j',
-        'ₖ' => 'k', 'ₗ' => 'l', 'ₘ' => 'm', 'ₙ' => 'n', 'ₒ' => 'o',
-        'ₚ' => 'p', 'ᵣ' => 'r', 'ₛ' => 's', 'ₜ' => 't', 'ᵤ' => 'u',
-        'ᵥ' => 'v', 'ₓ' => 'x'
-    )
-    super_map = Dict(
-        '⁰' => '0', '¹' => '1', '²' => '2', '³' => '3', '⁴' => '4',
-        '⁵' => '5', '⁶' => '6', '⁷' => '7', '⁸' => '8', '⁹' => '9',
-        '⁺' => '+', '⁻' => '-', '⁼' => '=', '⁽' => '(', '⁾' => ')',
-        'ᵃ' => 'a', 'ᵇ' => 'b', 'ᶜ' => 'c', 'ᵈ' => 'd', 'ᵉ' => 'e',
-        'ᶠ' => 'f', 'ᵍ' => 'g', 'ʰ' => 'h', 'ᶦ' => 'i', 'ʲ' => 'j',
-        'ᵏ' => 'k', 'ˡ' => 'l', 'ᵐ' => 'm', 'ⁿ' => 'n', 'ᵒ' => 'o',
-        'ᵖ' => 'p', 'ʳ' => 'r', 'ˢ' => 's', 'ᵗ' => 't', 'ᵘ' => 'u',
-        'ᵛ' => 'v', 'ʷ' => 'w', 'ˣ' => 'x', 'ʸ' => 'y', 'ᶻ' => 'z'
-    )
-
-    buf = IOBuffer()
-    for c in x
-        if haskey(sub_map, c)
-            write(buf, sub_map[c])
-        elseif haskey(super_map, c)
-            write(buf, super_map[c])
-        else
-            write(buf, c)
-        end
-    end
-    return String(take!(buf))
-end
-
 function solver_tol_to_dict(st::SolverTolerances)
     return Dict{Symbol,Any}(
         :atol => st.atol,
@@ -2233,6 +2195,7 @@ function plot_irf(𝓂::ℳ;
             shock_sort_perm = sortperm(shock_names_display, by = normalize_superscript)
             shock_idx = shock_idx[shock_sort_perm]
             shock_names_display = shock_names_display[shock_sort_perm]
+            Y = Y[:, :, shock_sort_perm]
         end
     else
         shock_names_display = ["shock_matrix"]
@@ -2869,6 +2832,7 @@ function plot_irf!(𝓂::ℳ;
             shock_sort_perm = sortperm(shock_names_display, by = normalize_superscript)
             shock_idx = shock_idx[shock_sort_perm]
             shock_names_display = shock_names_display[shock_sort_perm]
+            Y = Y[:, :, shock_sort_perm]
         end
     else
         shock_names_display = ["shock_matrix"]

models/Aguiar_Gopinath_2007.jl

@@ -1,3 +1,6 @@
+# Aguiar, M. and Gopinath, G. (2007). Emerging market business cycles: The cycle is the trend.
+# Journal of Political Economy, 115(1), 69–102. https://doi.org/10.1086/511283
+
 @model Aguiar_Gopinath_2007 begin
 	y[0] = (exp(g[0]) * l[0]) ^ alpha * exp(z[0]) * k[-1] ^ (1 - alpha)
 

models/Ascari_Sbordone_2014.jl

@@ -1,3 +1,6 @@
+# Ascari, G. and Sbordone, A.M. (2014). The macroeconomics of trend inflation.
+# Journal of Economic Literature, 52(3), 679–739. https://doi.org/10.1257/jel.52.3.679
+
 @model Ascari_Sbordone_2014 begin
 	1 / y[0] ^ sigma = beta * (1 + i[0]) / (pi[1] * y[1] ^ sigma)
 

models/Backus_Kehoe_Kydland_1992.jl

@@ -1,44 +1,52 @@
+# Backus, D.K., Kehoe, P.J., and Kydland, F.E. (1992). International real business cycles.
+# Journal of Political Economy, 100(4), 745–775. https://www.jstor.org/stable/2138686
+
 @model Backus_Kehoe_Kydland_1992 begin
+    # To scale to N countries, add/remove symbols in every [H, F] list below.
+    # Update the @parameters block rho definitions and steady-state calibration accordingly.
     for co in [H, F]
-        Y{co}[0] = ((LAMBDA{co}[0] * K{co}[-4]^theta{co} * N{co}[0]^(1-theta{co}))^(-nu{co}) + sigma{co} * Z{co}[-1]^(-nu{co}))^(-1/nu{co})
+        Y{co}[0] = ((LAMBDA{co}[0] * K{co}[-4]^theta{co} * N{co}[0]^(1 - theta{co}))^(-nu{co}) + sigma{co} * Z{co}[-1]^(-nu{co}))^(-1 / nu{co})
 
-        K{co}[0] = (1-delta{co})*K{co}[-1] + S{co}[0]
+        K{co}[0] = (1 - delta{co}) * K{co}[-1] + S{co}[0]
 
-        X{co}[0] = for lag in (-4+1):0 phi{co} * S{co}[lag] end
+        X{co}[0] = for lag in (-4 + 1):0 phi{co} * S{co}[lag] end
 
-        A{co}[0] = (1-eta{co}) * A{co}[-1] + N{co}[0]
+        A{co}[0] = (1 - eta{co}) * A{co}[-1] + N{co}[0]
 
-        L{co}[0] = 1 - alpha{co} * N{co}[0] - (1-alpha{co})*eta{co} * A{co}[-1]
+        L{co}[0] = 1 - alpha{co} * N{co}[0] - (1 - alpha{co}) * eta{co} * A{co}[-1]
 
-        U{co}[0] = (C{co}[0]^mu{co}*L{co}[0]^(1-mu{co}))^gamma{co}
+        U{co}[0] = (C{co}[0]^mu{co} * L{co}[0]^(1 - mu{co}))^gamma{co}
 
-        psi{co} * mu{co} / C{co}[0]*U{co}[0] = LGM[0]
+        psi{co} * mu{co} / C{co}[0] * U{co}[0] = LGM[0]
 
-        psi{co} * (1-mu{co}) / L{co}[0] * U{co}[0] * (-alpha{co}) = - LGM[0] * (1-theta{co}) / N{co}[0] * (LAMBDA{co}[0] * K{co}[-4]^theta{co}*N{co}[0]^(1-theta{co}))^(-nu{co})*Y{co}[0]^(1+nu{co})
+        psi{co} * (1 - mu{co}) / L{co}[0] * U{co}[0] * (-alpha{co}) = - LGM[0] * (1 - theta{co}) / N{co}[0] * (LAMBDA{co}[0] * K{co}[-4]^theta{co} * N{co}[0]^(1 - theta{co}))^(-nu{co}) * Y{co}[0]^(1 + nu{co})
 
-        for lag in 0:(4-1)  
-            beta{co}^lag * LGM[lag]*phi{co}
+        for lag in 0:(4 - 1)  
+            beta{co}^lag * LGM[lag] * phi{co}
         end +
         for lag in 1:4
-            -beta{co}^lag * LGM[lag] * phi{co} * (1-delta{co})
-        end = beta{co}^4 * LGM[+4] * theta{co} / K{co}[0] * (LAMBDA{co}[+4] * K{co}[0]^theta{co} * N{co}[+4]^(1-theta{co})) ^ (-nu{co})* Y{co}[+4]^(1+nu{co})
+            -beta{co}^lag * LGM[lag] * phi{co} * (1 - delta{co})
+        end = beta{co}^4 * LGM[+4] * theta{co} / K{co}[0] * (LAMBDA{co}[+4] * K{co}[0]^theta{co} * N{co}[+4]^(1 - theta{co}))^(-nu{co}) * Y{co}[+4]^(1 + nu{co})
 
-        LGM[0] = beta{co} * LGM[+1] * (1+sigma{co} * Z{co}[0]^(-nu{co}-1)*Y{co}[+1]^(1+nu{co}))
+        LGM[0] = beta{co} * LGM[+1] * (1 + sigma{co} * Z{co}[0]^(-nu{co} - 1) * Y{co}[+1]^(1 + nu{co}))
 
-        NX{co}[0] = (Y{co}[0] - (C{co}[0] + X{co}[0] + Z{co}[0] - Z{co}[-1]))/Y{co}[0]
+        NX{co}[0] = (Y{co}[0] - (C{co}[0] + X{co}[0] + Z{co}[0] - Z{co}[-1])) / Y{co}[0]
     end
 
-    (LAMBDA{H}[0]-1) = rho{H}{H}*(LAMBDA{H}[-1]-1) + rho{H}{F}*(LAMBDA{F}[-1]-1) + Z_E{H} * E{H}[x]
-
-    (LAMBDA{F}[0]-1) = rho{F}{F}*(LAMBDA{F}[-1]-1) + rho{F}{H}*(LAMBDA{H}[-1]-1) + Z_E{F} * E{F}[x]
+    # Shock process: each country's lambda depends on its own lag and spillovers from all others.
+    # The inner accumulator sums over co2 != co, handling any number of countries.
+    for co in [H, F]
+        (LAMBDA{co}[0] - 1) = rho{co}{co} * (LAMBDA{co}[-1] - 1) + for co2 in [H, F] if co2 != co rho{co}{co2} * (LAMBDA{co2}[-1] - 1) end end + Z_E{co} * E{co}[x]
+    end
 
-    for co in [H,F] C{co}[0] + X{co}[0] + Z{co}[0] - Z{co}[-1] end = for co in [H,F] Y{co}[0] end
+    # World resource constraint: sum of expenditures = sum of outputs
+    for co in [H, F] C{co}[0] + X{co}[0] + Z{co}[0] - Z{co}[-1] end = for co in [H, F] Y{co}[0] end
 end
 
 @parameters Backus_Kehoe_Kydland_1992 begin
     K_ss = 11
     K[ss] = K_ss | beta
-    
+
     mu      =    0.34
     gamma   =    -1.0
     alpha   =    1
@@ -51,9 +59,15 @@ end
     psi     =    0.5
 
     Z_E = 0.00852
-
-    rho{H}{H} = 0.906
-    rho{F}{F} = rho{H}{H}
-    rho{H}{F} = 0.088
-    rho{F}{H} = rho{H}{F}
+
+    # To scale to N countries, update the country list in the loops below
+    for co1 in [H, F]
+        for co2 in [H, F]
+            if co1 == co2
+                rho{co1}{co2} = 0.906
+            else
+                rho{co1}{co2} = 0.088 / 1 # divide by N - 1 for the model to be stable with more countries 
+            end
+        end
+    end
 end

models/Baxter_King_1993.jl

@@ -1,3 +1,6 @@
+# Baxter, M. and King, R.G. (1993). Fiscal policy in general equilibrium.
+# American Economic Review, 83(3), 315–334. https://www.jstor.org/stable/2117521
+
 @model Baxter_King_1993 begin
     uc[0] = c[0]^-1