diff --git a/README.md b/README.md index 73d0efe6c..a679384a9 100644 --- a/README.md +++ b/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` diff --git a/docs/src/index.md b/docs/src/index.md index b7efe3592..2c0fa65d9 100644 --- a/docs/src/index.md +++ b/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` diff --git a/docs/src/unfinished_docs/todo.md b/docs/src/unfinished_docs/todo.md index 173396b25..1b0d2aa11 100644 --- a/docs/src/unfinished_docs/todo.md +++ b/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 _ diff --git a/ext/MooncakeExt.jl b/ext/MooncakeExt.jl index ebb08dbfa..e1a4912b7 100644 --- a/ext/MooncakeExt.jl +++ b/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 \ No newline at end of file diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 702cc0567..981bcf626 100644 --- a/ext/StatsPlotsExt.jl +++ b/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"] diff --git a/models/Aguiar_Gopinath_2007.jl b/models/Aguiar_Gopinath_2007.jl index 04bab5e88..e6477f236 100644 --- a/models/Aguiar_Gopinath_2007.jl +++ b/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) diff --git a/models/Ascari_Sbordone_2014.jl b/models/Ascari_Sbordone_2014.jl index c246ab450..a73543548 100644 --- a/models/Ascari_Sbordone_2014.jl +++ b/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) diff --git a/models/Backus_Kehoe_Kydland_1992.jl b/models/Backus_Kehoe_Kydland_1992.jl index c6eeb4fae..4870e90ba 100644 --- a/models/Backus_Kehoe_Kydland_1992.jl +++ b/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 \ No newline at end of file diff --git a/models/Baxter_King_1993.jl b/models/Baxter_King_1993.jl index 31fbd8d26..51f310ab6 100644 --- a/models/Baxter_King_1993.jl +++ b/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 diff --git a/models/Bernanke_Gertler_Gilchrist_1999.jl b/models/Bernanke_Gertler_Gilchrist_1999.jl new file mode 100644 index 000000000..3c7d51513 --- /dev/null +++ b/models/Bernanke_Gertler_Gilchrist_1999.jl @@ -0,0 +1,128 @@ +# Bernanke, B., Gertler, M., and Gilchrist, S. (1999). The financial accelerator in a +# quantitative business cycle framework. In J. B. Taylor and M. Woodford (Eds.), +# Handbook of Macroeconomics (Vol. 1, pp. 1341–1393). Elsevier. +# https://doi.org/10.1016/S1574-0048(99)10034-X + +@model Bernanke_Gertler_Gilchrist_1999 begin + # [B.3] Household Euler equation: intertemporal optimality for deposits. + C[0]^(-1) = β * C[1]^(-1) * R[0] / π[1] + + # [B.4] Household intratemporal condition: consumption-leisure tradeoff. + ζ * C[0] / (1 - H[0]) = W[0] + + # [4.6] Labor aggregation with entrepreneurial labor share Ω (Hᵉ normalized to 1). + L[0] = H[0]^(1 - Ω) * Hᵉ^Ω + + # [4.1] Final goods technology: Cobb-Douglas production. + Y[0] = A[0] * K[-1]^α * L[0]^(1 - α) + + # [4.11] Household labor demand: wage equals marginal product. + W[0] = (1 - α) * (1 - Ω) * Y[0] / (X[0] * H[0]) + + # [4.12] Entrepreneurial labor income from marginal product. + Wᵉ[0] = (1 - α) * Ω * Y[0] / X[0] + + # [4.4] Gross return on capital (Rᵏ = R^k in the paper). + Rᵏ[0] = (α * Y[0] / (X[0] * K[-1]) + (1 - δ) * Q[0]) / Q[-1] + + # [4.2] Capital accumulation with convex investment adjustment costs. + K[0] = (1 - δ) * K[-1] + (I[0] / K[-1] - ψⁱ / 2 * (I[0] / K[-1] - δ)^2) * K[-1] + + # [4.3] Tobin's Q implied by the marginal adjustment cost of investment. + Q[0] = 1 / (1 - ψⁱ * (I[0] / K[-1] - δ)) + + # [B.8] Aggregate resource constraint including entrepreneurial consumption and monitoring costs. + Y[0] = C[0] + Cᵉ[0] + I[0] + Gov[0] + monitoring_cost[0] + + # [4.25] Monetary policy rule used in the paper's quantitative model. + R[0] = (πˢˢ / β)^(1 - ρʳ) * R[-1]^ρʳ * (π[0] / πˢˢ)^((1 - ρʳ) * ϕ_pi) * (Y[0] / Y[ss])^((1 - ρʳ) * ϕʸ) * exp(σʳ * εʳ[x]) + + # [B.12] Calvo price index. + 1.0 = θ * π[0]^(ϵᵖ - 1) + (1 - θ) * πstar[0]^(1 - ϵᵖ) + + # [B.11] Optimal reset price from the Calvo price-setting FOC. + πstar[0] = ϵᵖ / (ϵᵖ - 1) * price_aux_1[0] / price_aux_2[0] + + # [B.11a] Calvo numerator recursion. + price_aux_1[0] = Y[0] / X[0] / C[0] + β * θ * π[1]^ϵᵖ * price_aux_1[1] + + # [B.11b] Calvo denominator recursion. + price_aux_2[0] = Y[0] / C[0] + β * θ * π[1]^(ϵᵖ - 1) * price_aux_2[1] + + # [4.27] Aggregate technology process. + log(A[0]) = ρᵃ * log(A[-1]) + σᵃ * εᵃ[x] + + # [4.26] Government spending level (level form from the log-linearized process). + Gov[0] = Govʸ * Y[0] * exp(gov_gap[0]) + + # [4.26] Government spending AR(1) disturbance. + gov_gap[0] = ρᵍ * gov_gap[-1] + σᵍ * εᵍ[x] + + # [B.5] Money demand equation. + M_real[0] = χ * C[0] * R[0] / (R[0] - 1) + + # [4.8] Entrepreneurial equity (net of lender repayment and monitoring). + V[0] = (1 - Γ[0]) * Rᵏ[0] * Q[-1] * K[-1] + + # [4.8] Entrepreneurial consumption: dying entrepreneurs consume their equity. + Cᵉ[0] = (1 - γᵉ) * V[0] + + # [4.7] Entrepreneurial net worth: retained equity plus labor income. + N[0] = γᵉ * V[0] + Wᵉ[0] + + # [3.2] Balance sheet identity defining leverage. + leverage[0] * N[0] = Q[0] * K[0] + + # [A.2] Lognormal distribution: default probability F(ω̄). + F[0] = normcdf((log(ω̄[0]) + σ_ω^2 / 2) / σ_ω) + + # [A.2] Lognormal distribution: partial expectation G(ω̄). + G[0] = normcdf((log(ω̄[0]) - σ_ω^2 / 2) / σ_ω) + + # [A.2] Contract object Γ(ω̄) = G(ω̄) + ω̄(1 - F(ω̄)). + Γ[0] = G[0] + ω̄[0] * (1 - F[0]) + + # [3.5] Lender participation condition (inverted for external finance premium s). + s[0] = (1 - 1 / leverage[0]) / (Γ[0] - μ * G[0]) + + # [A.1]–[A.2] Optimal contract FOC linking leverage and the default threshold. + leverage[0] - 1 = (1 - F[0]) * (Γ[0] - μ * G[0]) / ((1 - F[0] - μ * ω̄[0] * normpdf((log(ω̄[0]) + σ_ω^2 / 2) / σ_ω) / σ_ω) * (1 - Γ[0])) + + # [4.5] External finance premium wedge (ex-post). + Rᵏ[1] = s[0] * R[0] / π[1] + + # [B.8] Aggregate monitoring costs paid in default states. + monitoring_cost[0] = μ * G[0] * Rᵏ[0] * Q[-1] * K[-1] + + # Observable/log-output convenience variable. + log_y[0] = log(Y[0]) + + # Annualized external finance premium convenience variable. + premium_ann[0] = 400 * log(s[0]) +end + +@parameters Bernanke_Gertler_Gilchrist_1999 begin + β = 0.99 + α = 0.35 + δ = 0.025 + Ω = 0.015384615384615385 + Hᵉ = 1.0 + ζ = 2.0 + ψⁱ = 4.0 + μ = 0.12 + σ_ω = 0.28 + γᵉ = 0.9728 + πˢˢ = 1.0 + θ = 0.75 + ϵᵖ = 11.0 + ρʳ = 0.9 + ϕ_pi = 1.1 + ϕʸ = 0.0 + ρᵃ = 0.95 + σᵃ = 0.01 + Govʸ = 0.20 + ρᵍ = 0.95 + σᵍ = 0.01 + σʳ = 0.0025 + χ = 0.0099 +end diff --git a/models/Caldara_et_al_2012.jl b/models/Caldara_et_al_2012.jl index a93eae445..e38ad9ac7 100644 --- a/models/Caldara_et_al_2012.jl +++ b/models/Caldara_et_al_2012.jl @@ -1,3 +1,7 @@ +# Caldara, D., Fernández-Villaverde, J., Rubio-Ramírez, J.F., and Yao, W. (2012). Computing +# DSGE models with recursive preferences and stochastic volatility. +# Review of Economic Dynamics, 15(2), 188–206. https://doi.org/10.1016/j.red.2011.10.001 + @model Caldara_et_al_2012 begin V[0] = ((1 - β) * (c[0] ^ ν * (1 - l[0]) ^ (1 - ν)) ^ (1 - 1 / ψ) + β * V[1] ^ (1 - 1 / ψ)) ^ (1 / (1 - 1 / ψ)) diff --git a/models/FS2000.jl b/models/FS2000.jl index 3e4d86e78..da5d89f60 100644 --- a/models/FS2000.jl +++ b/models/FS2000.jl @@ -1,3 +1,6 @@ +# Schorfheide, F. (2000). Loss function-based evaluation of DSGE models. +# Journal of Applied Econometrics, 15(6), 645–670. https://doi.org/10.1002/jae.582 + @model FS2000 begin dA[0] = exp(gam + z_e_a * e_a[x]) diff --git a/models/Fernandez_Villaverde_2010.jl b/models/Fernandez_Villaverde_2010.jl new file mode 100644 index 000000000..2824ea5b6 --- /dev/null +++ b/models/Fernandez_Villaverde_2010.jl @@ -0,0 +1,137 @@ +# Fernández-Villaverde, J. and Rubio-Ramírez, J.F. (2006). A Baseline DSGE Model. +# Manuscript, University of Pennsylvania. +# +# Fernández-Villaverde, J. (2010). The econometrics of DSGE models. +# SERIEs, 1(1–2), 3–49. https://doi.org/10.1007/s13209-009-0001-0 + +@model Fernandez_Villaverde_2010 begin + d[0] / (c[0] - h * c[-1] / μᶻ[0]) - h * β * d[1] / (c[1] * μᶻ[1] - h * c[0]) = λ[0] + + λ[0] = β * λ[1] / μᶻ[1] / π[1] * R[0] + + r[0] = γ₁ + γ₂ * (u[0] - 1) + + q[0] = β * λ[1] / λ[0] / μᶻ[1] / μᴵ[1] * ((1 - δ) * q[1] + r[1] * u[1] - (γ₁ * (u[1] - 1) + γ₂ / 2 * (u[1] - 1)^2)) + + 1 = q[0] * (1 - (κ / 2 * (x[0] / x[-1] * μᶻ[0] - Λˣ)^2) - (κ * (x[0] / x[-1] * μᶻ[0] - Λˣ) * x[0] / x[-1] * μᶻ[0])) + β * q[1] * λ[1] / λ[0] / μᶻ[1] * κ * (x[1] / x[0] * μᶻ[1] - Λˣ) * (x[1] / x[0] * μᶻ[1])^2 + + f[0] = (η - 1) / η * w_star[0]^(1 - η) * λ[0] * w[0]^η * lᵈ[0] + β * θʷ * (π[0]^χʷ / π[1])^(1 - η) * (w_star[1] / w_star[0] * μᶻ[1])^(η - 1) * f[1] + + f[0] = ψ * d[0] * φ[0] * π_w_star[0]^(-η * (1 + γ)) * lᵈ[0]^(1 + γ) + β * θʷ * (π[0]^χʷ / π[1])^(-η * (1 + γ)) * (w_star[1] / w_star[0] * μᶻ[1])^(η * (1 + γ)) * f[1] + + g₁[0] = λ[0] * mc[0] * yᵈ[0] + β * θᵖ * (π[0]^χ / π[1])^(-ε) * g₁[1] + + g₂[0] = λ[0] * π_star[0] * yᵈ[0] + β * θᵖ * (π[0]^χ / π[1])^(1 - ε) * π_star[0] / π_star[1] * g₂[1] + + ε * g₁[0] = (ε - 1) * g₂[0] + + u[0] * k[-1] / lᵈ[0] = α / (1 - α) * w[0] / r[0] * μᶻ[0] * μᴵ[0] + + mc[0] = (1 / (1 - α))^(1 - α) * (1 / α)^α * w[0]^(1 - α) * r[0]^α + + 1 = θʷ * (π[-1]^χʷ / π[0])^(1 - η) * (w[-1] / w[0] / μᶻ[0])^(1 - η) + (1 - θʷ) * π_w_star[0]^(1 - η) + + 1 = θᵖ * (π[-1]^χ / π[0])^(1 - ε) + (1 - θᵖ) * π_star[0]^(1 - ε) + + R[0] / R̄ = (R[-1] / R̄)^γᴿ * ((π[0] / π̄)^γ_π * ((yᵈ[0] / yᵈ[-1] * μᶻ[0]) / exp(Λʸᵈ))^γʸ)^(1 - γᴿ) * exp(σᵐ * ϵᵐ[x]) + + yᵈ[0] = c[0] + x[0] + (γ₁ * (u[0] - 1) + γ₂ / 2 * (u[0] - 1)^2) * k[-1] / μᶻ[0] / μᴵ[0] + + yᵈ[0] = (μᴬ[0] / μᶻ[0] * (u[0] * k[-1])^α * lᵈ[0]^(1 - α) - Φ) / vᵖ[0] + + l[0] = vʷ[0] * lᵈ[0] + + vᵖ[0] = θᵖ * (π[-1]^χ / π[0])^(-ε) * vᵖ[-1] + (1 - θᵖ) * π_star[0]^(-ε) + + vʷ[0] = θʷ * (w[-1] / w[0] / μᶻ[0] * π[-1]^χʷ / π[0])^(-η) * vʷ[-1] + (1 - θʷ) * π_w_star[0]^(-η) + + k[0] = (1 - δ) * k[-1] / (μᶻ[0] * μᴵ[0]) + (1 - κ / 2 * (x[0] / x[-1] * μᶻ[0] - Λˣ)^2) * x[0] + + F[0] = yᵈ[0] - 1 / (1 - α) * w[0] * lᵈ[0] + + π_w_star[0] = w_star[0] / w[0] + + log(d[0]) = ρᵈ * log(d[-1]) + σᵈ * ϵᵈ[x] + + log(φ[0]) = ρᵠ * log(φ[-1]) + σᵠ * ϵᵠ[x] + + log(μᴵ[0]) = Λ_μ + σ_μ * ϵ_μ_I[x] + + log(μᴬ[0]) = Λᴬ + σᴬ * ϵᴬ[x] + + μᶻ[0] = μᴬ[0]^(1 / (1 - α)) * μᴵ[0]^(α / (1 - α)) +end + +@parameters Fernandez_Villaverde_2010 begin + # Steady-state growth rates (computed from drift parameters) + Λʸᵈ = (Λᴬ + α * Λ_μ) / (1 - α) + + Λˣ = exp(Λʸᵈ) + + γ₁ = exp(Λʸᵈ) * exp(Λ_μ) / β - (1 - δ) + + R̄ = 1 + (π̄ * exp(Λʸᵈ) / β - 1) + + # Structural parameters + h = 0.97 + + β = 0.998 + + δ = 0.025 + + κ = 9.51 + + η = 10 + + ε = 10 + + ψ = 8.92 + + γ = 1.17 + + χʷ = 0.62 + + χ = 0.63 + + θᵖ = 0.82 + + θʷ = 0.68 + + α = 0.21 + + Φ = 0 + + γ₂ = 0.001 + + # Monetary policy rule + γᴿ = 0.77 + + γʸ = 0.19 + + γ_π = 1.29 + + π̄ = 1.01 + + # Shock persistence + ρᵈ = 0.12 + + ρᵠ = 0.93 + + # Shock drift parameters + Λ_μ = 0.0034 + + Λᴬ = 0.0028 + + # Shock standard deviations (median posterior estimates, FV 2010 Table 1) + σᵈ = exp(-1.51) + + σᵠ = exp(-2.36) + + σ_μ = exp(-5.43) + + σᴬ = exp(-3.97) + + σᵐ = exp(-5.85) + + u > 0 +end diff --git a/models/GNSS_2010.jl b/models/GNSS_2010.jl index a71337a38..da30f1c00 100644 --- a/models/GNSS_2010.jl +++ b/models/GNSS_2010.jl @@ -1,3 +1,7 @@ +# Gerali, A., Neri, S., Sessa, L., and Signoretti, F.M. (2010). Credit and banking in a DSGE +# model of the euro area. Journal of Money, Credit and Banking, 42(s1), 107–141. +# https://doi.org/10.1111/j.1538-4616.2010.00331.x + @model GNSS_2010 begin (1 - a_i) * exp(ee_z[0]) * (c_p[0] - a_i * c_p[-1]) ^ (-1) = lam_p[0] diff --git a/models/Gali_2015_chapter_3_nonlinear.jl b/models/Gali_2015_chapter_3_nonlinear.jl index db71d8a1f..89bb982c0 100644 --- a/models/Gali_2015_chapter_3_nonlinear.jl +++ b/models/Gali_2015_chapter_3_nonlinear.jl @@ -1,3 +1,6 @@ +# Gali, J. (2015). Monetary Policy, Inflation, and the Business Cycle: An Introduction to the +# New Keynesian Framework and Its Applications (2nd ed.). Princeton University Press. + @model Gali_2015_chapter_3_nonlinear begin W_real[0] = C[0] ^ σ * N[0] ^ φ diff --git a/models/Gali_2015_chapter_3_obc.jl b/models/Gali_2015_chapter_3_obc.jl index 09212af7a..68555fc88 100644 --- a/models/Gali_2015_chapter_3_obc.jl +++ b/models/Gali_2015_chapter_3_obc.jl @@ -1,3 +1,6 @@ +# Gali, J. (2015). Monetary Policy, Inflation, and the Business Cycle: An Introduction to the +# New Keynesian Framework and Its Applications (2nd ed.). Princeton University Press. + @model Gali_2015_chapter_3_obc begin W_real[0] = C[0] ^ σ * N[0] ^ φ diff --git a/models/Gali_Monacelli_2005_CITR.jl b/models/Gali_Monacelli_2005_CITR.jl index 28f810a09..8c75a5b81 100644 --- a/models/Gali_Monacelli_2005_CITR.jl +++ b/models/Gali_Monacelli_2005_CITR.jl @@ -1,3 +1,7 @@ +# Gali, J. and Monacelli, T. (2005). Monetary policy and exchange rate volatility in a small +# open economy. Review of Economic Studies, 72(3), 707–734. +# https://doi.org/10.1111/j.1467-937X.2005.00349.x + @model Gali_Monacelli_2005_CITR begin x[0] = x[1] - σᵃ ^ (-1) * (r[0] - pih[1] - rnat[0]) diff --git a/models/Gertler_Karadi_2011.jl b/models/Gertler_Karadi_2011.jl new file mode 100644 index 000000000..28279b561 --- /dev/null +++ b/models/Gertler_Karadi_2011.jl @@ -0,0 +1,152 @@ +# Gertler, M. and Karadi, P. (2011). A model of unconventional monetary policy. +# Journal of Monetary Economics, 58(1), 17–34. https://doi.org/10.1016/j.jmoneco.2010.10.004 + +@model Gertler_Karadi_2011 begin + # [4] Marginal utility of consumption with external habit. + ϱ[0] = (C[0] - h * C[-1])^(-σ) - β * h * (C[1] - h * C[0])^(-σ) + + # [4'] Household Euler equation for the riskless asset. + β * R[0] * Λ[1] = 1.0 + + # [--] Stochastic discount factor. + Λ[0] = ϱ[0] / ϱ[-1] + + # [3] Labor market equilibrium. + χ * L[0]^φ_l = ϱ[0] * Pₘ[0] * (1 - α) * Y[0] / L[0] + + # [11] Marginal value of bankers' capital (ν_t). + ν[0] = (1 - θ) * β * Λ[1] * (Rᵏ[1] - R[0]) + β * Λ[1] * θ * x[1] * ν[1] + + # [11] Marginal value of bankers' net worth (η_t). + η[0] = (1 - θ) + β * Λ[1] * θ * z[1] * η[1] + + # [13] Incentive-constraint leverage relation (φ_t = η_t/(λ - ν_t)). + φ[0] = η[0] / (λ - ν[0]) + + # [14] Gross growth rate of bankers' capital (z_{t-1,t} = N_t/N_{t-1}). + z[0] = (Rᵏ[0] - R[-1]) * φ[-1] + R[-1] + + # [--] Gross growth rate of assets (x_{t-1,t} = Q_t S_t/(Q_{t-1} S_{t-1})). + x[0] = φ[0] / φ[-1] * z[0] + + # [15] Aggregate intermediary balance sheet. + Q[0] * K[0] = φ[0] * N[0] + + # [16] Aggregate bankers' net worth (survivors + entrants). + N[0] = Nᵉ[0] + Nⁿ[0] + + # [17] Existing bankers' net worth accumulation. + Nᵉ[0] = θ * z[0] * N[-1] * exp(-σ_Ne * ε_Ne[x]) + + # [18] Entering bankers' net worth. + Nⁿ[0] = ω * Q[0] * ξ[0] * K[-1] + + # [25] Gross return on capital. + Rᵏ[0] = (Pₘ[0] * α * Ym[0] / K[-1] + ξ[0] * (Q[0] - δ_rate[0])) / Q[-1] + + # [22] Intermediate goods production. + Ym[0] = A[0] * (ξ[0] * U[0] * K[-1])^α * L[0]^(1 - α) + + # [27] Capital producer's optimal investment condition (Tobin's Q). + Q[0] = 1 + η_i / 2 * ((In[0] + I[ss]) / (In[-1] + I[ss]) - 1)^2 + + η_i * ((In[0] + I[ss]) / (In[-1] + I[ss]) - 1) * (In[0] + I[ss]) / (In[-1] + I[ss]) - + β * Λ[1] * η_i * ((In[1] + I[ss]) / (In[0] + I[ss]) - 1) * ((In[1] + I[ss]) / (In[0] + I[ss]))^2 + + # [--] Utilization-dependent depreciation. + δ_rate[0] = δ + b_u / (1 + ζ_u) * (U[0]^(1 + ζ_u) - 1) + + # [23] Optimal capacity utilization. + Pₘ[0] * α * Ym[0] / U[0] = b_u * U[0]^ζ_u * ξ[0] * K[-1] + + # [26] Net investment. + In[0] = I[0] - δ_rate[0] * ξ[0] * K[-1] + + # [35] Capital accumulation. + K[0] = ξ[0] * K[-1] + In[0] + + # [--] Government consumption. + G[0] = g_y * Y[ss] * exp(ĝ[0]) + + # [34] Aggregate resource constraint. + Y[0] = C[0] + G[0] + I[0] + η_i / 2 * ((In[0] + I[ss]) / (In[-1] + I[ss]) - 1)^2 * (In[0] + I[ss]) + + # [--] Wholesale and retail output relation. + Ym[0] = Y[0] * D[0] + + # [--] Calvo price dispersion. + D[0] = γ_p * D[-1] * π[-1]^(-γ_p_index * ϵ_p) * π[0]^ϵ_p + + (1 - γ_p) * ((1 - γ_p * π[-1]^(γ_p_index * (1 - ϵ_p)) * π[0]^(ϵ_p - 1)) / (1 - γ_p))^(-ϵ_p / (1 - ϵ_p)) + + # [--] Markup definition. + X[0] = 1 / Pₘ[0] + + # [31] Optimal reset-price numerator recursion (F_t). + F[0] = Y[0] * Pₘ[0] + β * γ_p * Λ[1] * π[1]^ϵ_p * π[0]^(-ϵ_p * γ_p_index) * F[1] + + # [31] Optimal reset-price denominator recursion (Z_t). + Z[0] = Y[0] + β * γ_p * Λ[1] * π[1]^(ϵ_p - 1) * π[0]^(γ_p_index * (1 - ϵ_p)) * Z[1] + + # [31]–[32] Optimal reset-price inflation. + π_star[0] = ϵ_p / (ϵ_p - 1) * F[0] / Z[0] * π[0] + + # [33] Aggregate price index. + π[0]^(1 - ϵ_p) = γ_p * π[-1]^(γ_p_index * (1 - ϵ_p)) + (1 - γ_p) * π_star[0]^(1 - ϵ_p) + + # [38] Fisher equation. + i[0] = R[0] * π[1] + + # [37] Interest-rate rule (Taylor rule with smoothing). + i[0] = i[-1]^ρ_i * ((1 / β) * π[0]^κ_pi * (X[0] / (ϵ_p / (ϵ_p - 1)))^κ_y)^(1 - ρ_i) * exp(σ_i * ε_i[x]) + + # [--] TFP process. + A[0] = A[-1]^ρ_A * exp(-σ_A * ε_A[x]) + + # [--] Capital quality process. + ξ[0] = ξ[-1]^ρ_ξ * exp(-σ_ξ * ε_ξ[x]) + + # [--] Government spending process. + ĝ[0] = ρ_g * ĝ[-1] - σ_g * ε_g[x] + + # [--] Effective capital convenience variable. + Keff[0] = ξ[0] * K[-1] + + # [--] Wage convenience variable. + w[0] = Pₘ[0] * (1 - α) * Y[0] / L[0] + + # [--] Marginal value product of capital convenience variable. + VMPK[0] = Pₘ[0] * α * Y[0] / (ξ[0] * K[-1]) + + # [--] Premium convenience variable (Rᵏ_{t+1}/R_{t+1}). + prem[0] = Rᵏ[1] / R[0] +end + +@parameters Gertler_Karadi_2011 begin + β = 0.99 + σ = 1.0 + h = 0.815 + χ = 3.410808502195193 + φ_l = 0.276 + ζ_u = 7.2 + θ = 0.97155955 + α = 0.33 + δ = 0.025 + g_y = 0.2 + η_i = 1.728 + ϵ_p = 4.167 + γ_p = 0.779 + γ_p_index = 0.241 + ρ_i = 0.8 + κ_pi = 1.5 + κ_y = -0.5 / 4 + ρ_ξ = 0.66 + σ_ξ = 0.05 + ρ_A = 0.95 + σ_A = 0.01 + ρ_g = 0.95 + σ_g = 0.01 + σ_Ne = 0.01 + σ_i = 0.01 + λ = 0.38149498593241726 + ω = 0.0022277804179292593 + b_u = 0.037601010101010155 +end diff --git a/models/Ghironi_Melitz_2005.jl b/models/Ghironi_Melitz_2005.jl index a66372eb1..34ebc38c1 100644 --- a/models/Ghironi_Melitz_2005.jl +++ b/models/Ghironi_Melitz_2005.jl @@ -1,3 +1,7 @@ +# Ghironi, F. and Melitz, M.J. (2005). International trade and macroeconomic dynamics with +# heterogeneous firms. Quarterly Journal of Economics, 120(3), 865–915. +# https://doi.org/10.1093/qje/120.3.865 + @model Ghironi_Melitz_2005 begin 1 = Nd[0] * ρ̃d[0] ^ (1 - θ) + Nx̄[0] * ρ̃x̄[0] ^ (1 - θ) diff --git a/models/Guerrieri_Iacoviello_2017.jl b/models/Guerrieri_Iacoviello_2017.jl index f8b5412eb..be3f92fa8 100644 --- a/models/Guerrieri_Iacoviello_2017.jl +++ b/models/Guerrieri_Iacoviello_2017.jl @@ -1,3 +1,7 @@ +# Guerrieri, L. and Iacoviello, M. (2017). Collateral constraints and macroeconomic +# asymmetries. Journal of Monetary Economics, 90, 28–49. +# https://doi.org/10.1016/j.jmoneco.2017.06.004 + @model Guerrieri_Iacoviello_2017 begin c[0] + c1[0] + ik[0] = y[0] diff --git a/models/Iacoviello_2005_linear.jl b/models/Iacoviello_2005_linear.jl index 55d405e70..2f15f6ffc 100644 --- a/models/Iacoviello_2005_linear.jl +++ b/models/Iacoviello_2005_linear.jl @@ -1,3 +1,7 @@ +# Iacoviello, M. (2005). House prices, borrowing constraints, and monetary policy in the +# business cycle. American Economic Review, 95(3), 739–764. +# https://doi.org/10.1257/0002828054201477 + @model Iacoviello_2005_linear begin Ŷ[0] = ctoY * ĉ[0] + c1toY * c1̂[0] + c2toY * c2̂[0] + ItoY * Î[0] diff --git a/models/Ireland_2004.jl b/models/Ireland_2004.jl index 5740a53de..77f9c25e3 100644 --- a/models/Ireland_2004.jl +++ b/models/Ireland_2004.jl @@ -1,3 +1,6 @@ +# Ireland, P.N. (2004). Technology shocks in the New Keynesian model. +# Review of Economics and Statistics, 86(4), 923–936. https://doi.org/10.1162/0034653043125185 + @model Ireland_2004 begin a[0] = ρᵃ * a[-1] + σᵃ * ϵᵃ[x] diff --git a/models/JQ_2012_RBC.jl b/models/JQ_2012_RBC.jl index 0c9b29f88..ff5e4b8e5 100644 --- a/models/JQ_2012_RBC.jl +++ b/models/JQ_2012_RBC.jl @@ -1,3 +1,6 @@ +# Jermann, U. and Quadrini, V. (2012). Macroeconomic effects of financial shocks. +# American Economic Review, 102(1), 238–271. https://doi.org/10.1257/aer.102.1.238 + @model JQ_2012_RBC begin w[0] / c[0] ^ σ = α / (1 - n[0]) diff --git a/models/NAWM_EAUS_2008.jl b/models/NAWM_EAUS_2008.jl index 6881ca613..4db37f41e 100644 --- a/models/NAWM_EAUS_2008.jl +++ b/models/NAWM_EAUS_2008.jl @@ -1,3 +1,7 @@ +# European Central Bank. (2008). The New Area-Wide Model of the euro area: A micro-founded +# open-economy model for forecasting and policy analysis. ECB Working Paper No. 944. +# https://www.ecb.europa.eu/pub/pdf/scpwps/ecbwp944.pdf + @model NAWM_EAUS_2008 begin EA_R[0] ^ 4 - 1 = EA_PHIRR * (EA_R[-1] ^ 4 - 1) + (1 - EA_PHIRR) * (EA_RRSTAR ^ 4 * EA_PI4TARGET - 1 + EA_PHIRPI * (EA_PIC4[0] - EA_PI4TARGET)) + EA_PHIRGY * (EA_Y[0] / EA_Y[-1] - 1) + σ_EA_R * EA_EPSR[x] diff --git a/models/QUEST3_2009.jl b/models/QUEST3_2009.jl index eb941fed9..c5f01b730 100644 --- a/models/QUEST3_2009.jl +++ b/models/QUEST3_2009.jl @@ -1,3 +1,7 @@ +# Ratto, M., Roeger, W., and in 't Veld, J. (2009). QUEST III: An estimated open-economy +# DSGE model of the euro area with fiscal and monetary policy. +# Economic Modelling, 26(1), 222–233. https://doi.org/10.1016/j.econmod.2008.06.014 + @model QUEST3_2009 begin interest[0] = ((1 + E_INOM[0]) ^ 4 - interestq_exog ^ 4) / interestq_exog ^ 4 diff --git a/models/RBC_baseline.jl b/models/RBC_baseline.jl index d5fcaeec8..18c949f74 100644 --- a/models/RBC_baseline.jl +++ b/models/RBC_baseline.jl @@ -1,3 +1,7 @@ +# King, R.G., Plosser, C.I., and Rebelo, S.T. (1988). Production, growth and business cycles: +# I. The basic neoclassical model. Journal of Monetary Economics, 21(2–3), 195–232. +# https://doi.org/10.1016/0304-3932(88)90030-X + using MacroModelling @model RBC_baseline begin diff --git a/models/SGU_2003_debt_premium.jl b/models/SGU_2003_debt_premium.jl index 2cbc6b8df..3844eaf39 100644 --- a/models/SGU_2003_debt_premium.jl +++ b/models/SGU_2003_debt_premium.jl @@ -1,3 +1,7 @@ +# Schmitt-Grohé, S. and Uribe, M. (2003). Closing small open economy models. +# Journal of International Economics, 61(1), 163–185. +# https://doi.org/10.1016/S0022-1996(02)00056-9 + @model SGU_2003_debt_premium begin d[0] = (1 + r[-1]) * d[-1] - y[0] + c[0] + i[0] + ϕ / 2 * (k[0] - k[-1]) ^ 2 diff --git a/models/Smets_Wouters_2003.jl b/models/Smets_Wouters_2003.jl index 55751a9ce..f25c7d3b4 100644 --- a/models/Smets_Wouters_2003.jl +++ b/models/Smets_Wouters_2003.jl @@ -1,3 +1,7 @@ +# Smets, F. and Wouters, R. (2003). An estimated dynamic stochastic general equilibrium model +# of the euro area. Journal of the European Economic Association, 1(5), 1123–1175. +# https://doi.org/10.1162/154247603770383415 + @model Smets_Wouters_2003 begin -q[0] + beta * ((1 - tau) * q[1] + epsilon_b[1] * (r_k[1] * z[1] - psi^-1 * r_k[ss] * (-1 + exp(psi * (-1 + z[1])))) * (C[1] - h * C[0])^(-sigma_c)) diff --git a/models/Smets_Wouters_2003_obc.jl b/models/Smets_Wouters_2003_obc.jl index 6ecf5ea07..636152dfb 100644 --- a/models/Smets_Wouters_2003_obc.jl +++ b/models/Smets_Wouters_2003_obc.jl @@ -1,3 +1,9 @@ +# Smets, F. and Wouters, R. (2003). An estimated dynamic stochastic general equilibrium model +# of the euro area. Journal of the European Economic Association, 1(5), 1123–1175. +# https://doi.org/10.1162/154247603770383415 +# +# This variant includes an occasionally binding constraint (OBC) on the interest rate. + @model Smets_Wouters_2003_obc begin -q[0] + beta * ((1 - tau) * q[1] + epsilon_b[1] * (r_k[1] * z[1] - psi^-1 * r_k[ss] * (-1 + exp(psi * (-1 + z[1])))) * (C[1] - h * C[0])^(-sigma_c)) diff --git a/models/Smets_Wouters_2007.jl b/models/Smets_Wouters_2007.jl index e2979d985..6e83b1fdb 100644 --- a/models/Smets_Wouters_2007.jl +++ b/models/Smets_Wouters_2007.jl @@ -1,3 +1,7 @@ +# Smets, F. and Wouters, R. (2007). Shocks and frictions in US business cycles: A Bayesian +# DSGE approach. American Economic Review, 97(3), 586–606. +# https://doi.org/10.1257/aer.97.3.586 + @model Smets_Wouters_2007 begin y[0] = c[0] + inve[0] + y[ss] * gy[0] + afunc[0] * kp[-1] / cgamma diff --git a/models/Smets_Wouters_2007_linear.jl b/models/Smets_Wouters_2007_linear.jl index 40eb7e5f0..0b1bbd2c2 100644 --- a/models/Smets_Wouters_2007_linear.jl +++ b/models/Smets_Wouters_2007_linear.jl @@ -1,3 +1,9 @@ +# Smets, F. and Wouters, R. (2007). Shocks and frictions in US business cycles: A Bayesian +# DSGE approach. American Economic Review, 97(3), 586–606. +# https://doi.org/10.1257/aer.97.3.586 +# +# Linearised variant of the Smets and Wouters (2007) model. + @model Smets_Wouters_2007_linear begin a[0] = calfa * rkf[0] + (1 - calfa) * wf[0] diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index f10ccf6b0..40eb9bb1c 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -386,6 +386,43 @@ Base.show(io::IO, 𝓂::ℳ) = println(io, # ForwardDiffExt extends this for ForwardDiff.Dual numbers. primal(x::Real) = x +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 normalize_filtering_options(filter::Symbol, smooth::Bool, @@ -862,6 +899,10 @@ end if count(x -> abs(x) > tol, A) / length(A) < density_threshold && length(A) > min_length # Use dense_to_sparse to avoid Julia 1.12 SparseArrays bug in SparseMatrixCSC(::Matrix) a = dense_to_sparse(A, tol) + + # a = sparse(A) + # droptol!(a, tol) + if multithreaded return ThreadedSparseArrays.ThreadedSparseMatrixCSC(a) else @@ -2315,7 +2356,16 @@ end @warn "Invalid `shocks` argument. Provide a Symbol, Tuple, Vector, Matrix, or one of the documented selectors such as `:all`." shock_idx = Int64[] end - return shock_idx + + if shock_idx isa Integer + return shock_idx + else + shock_idx = unique(collect(shock_idx)) + if length(shock_idx) > 1 + sort!(shock_idx, by = i -> normalize_superscript(string(T.exo[i]))) + end + return shock_idx + end end diff --git a/src/algorithms/lyapunov.jl b/src/algorithms/lyapunov.jl index aad87f7b2..35415b77b 100644 --- a/src/algorithms/lyapunov.jl +++ b/src/algorithms/lyapunov.jl @@ -1,6 +1,6 @@ @stable default_mode = "disable" begin -# Available algorithms: +# Available algorithms: # :doubling - fast and precise # :bartels_stewart - fast for small matrices and precise, dense matrices only # :bicgstab - less precise @@ -113,7 +113,7 @@ end denom = max(ℒ.norm(guess), ℒ.norm(C)) reached_tol = denom == 0 ? 0.0 : ℒ.norm(res_buf) / denom if reached_tol < initial_guess_acceptance_tol - if verbose println("Lyapunov equation - initial guess achieves relative tol of $reached_tol (initial guess tol: $initial_guess_acceptance_tol)") end + if verbose println("Lyapunov equation - initial guess achieves relative tol of $reached_tol (initial guess tol: $initial_guess_acceptance_tol) [A: $(typeof(A)), C: $(typeof(C))]") end return choose_matrix_format(guess), true end end @@ -133,7 +133,7 @@ end verbose = verbose) if deflation_solved if verbose - println("Lyapunov equation - solved via Schur deflation (unit roots pre-detected)") + println("Lyapunov equation - solved via Schur deflation (unit roots pre-detected) [A: $(typeof(A_dense)), C: $(typeof(C_dense))]") end return X_deflated, true end @@ -143,7 +143,7 @@ end X, i, reached_tol = solve_lyapunov_equation(A, C, Val(lyapunov_algorithm), workspace; tol = tol) # timer = timer) if verbose - println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: $lyapunov_algorithm") + println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: $lyapunov_algorithm [A: $(typeof(A)), C: $(typeof(C))]") end if reached_tol > acceptance_tol && lyapunov_algorithm ≠ :doubling @@ -152,7 +152,7 @@ end X, i, reached_tol = solve_lyapunov_equation(A, C, Val(:doubling), workspace; tol = tol) # timer = timer) if verbose - println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: doubling") + println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: doubling [A: $(typeof(A)), C: $(typeof(C))]") end end @@ -162,7 +162,7 @@ end X, i, reached_tol = solve_lyapunov_equation(A, C, Val(:bicgstab), workspace; tol = tol) # timer = timer) if verbose - println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: bicgstab") + println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: bicgstab [A: $(typeof(A)), C: $(typeof(C))]") end end @@ -174,7 +174,7 @@ end X, i, reached_tol = solve_lyapunov_equation(A, C, Val(:bartels_stewart), workspace; tol = tol) # timer = timer) if verbose - println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: bartels_stewart") + println("Lyapunov equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: bartels_stewart [A: $(typeof(A)), C: $(typeof(C))]") end end @@ -191,7 +191,7 @@ end X = X_deflated reached_tol = zero(T) # signal success if verbose - println("Lyapunov equation - solved via Schur deflation (unit-root subspace set to NaN)") + println("Lyapunov equation - solved via Schur deflation (unit-root subspace set to NaN) [A: $(typeof(A_dense)), C: $(typeof(C_dense))]") end end end @@ -917,7 +917,7 @@ function solve_lyapunov_schur_deflation(A::DenseMatrix{T}, if sub_tol > tol.acceptance_tol if verbose - println("Schur deflation: stable sub-block Lyapunov failed (tol=$sub_tol)") + println("Schur deflation: stable sub-block Lyapunov failed (tol=$sub_tol) [T_ss: $(typeof(T_ss)), C_ss: $(typeof(C_ss))]") end return Matrix{T}(undef, 0, 0), false end @@ -946,7 +946,8 @@ function solve_lyapunov_schur_deflation(A::DenseMatrix{T}, if verbose println("Schur deflation: $n_unstable unstable eigenvalue(s), ", - "$n_stable stable, $(count(unit_root_vars)) variable(s) set to NaN") + "$n_stable stable, $(count(unit_root_vars)) variable(s) set to NaN ", + "[T_ss: $(typeof(T_ss)), C_ss: $(typeof(C_ss))]") end return Σ, true diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl index 99fb88763..827e5ff3c 100644 --- a/src/algorithms/sylvester.jl +++ b/src/algorithms/sylvester.jl @@ -1,6 +1,6 @@ @stable default_mode = "disable" begin -# Available algorithms: +# Available algorithms: # :doubling - fast, expensive part: B^2 # :bartels_stewart - fast, dense matrices only # :bicgstab - fastest for large problems, might not reach desired precision, warm start not always helpful @@ -89,7 +89,7 @@ reached_tol = denom == 0 ? 0.0 : ℒ.norm(res_buf) / denom if reached_tol < initial_guess_acceptance_tol - if verbose println("Sylvester equation - $guess_name achieves relative tol of $reached_tol (initial guess tol: $initial_guess_acceptance_tol)") end + if verbose println("Sylvester equation - $guess_name achieves relative tol of $reached_tol (initial guess tol: $initial_guess_acceptance_tol) [A: $(size(a)) $(typeof(a)), B: $(size(b)) $(typeof(b)), C: $(size(c)) $(typeof(c))]") end return choose_matrix_format(guess), true end @@ -106,7 +106,7 @@ verbose = verbose) if verbose && i != 0 - println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: $sylvester_algorithm") + println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: $sylvester_algorithm [A: $(size(a)) $(typeof(a)), B: $(size(b)) $(typeof(b)), C: $(size(c)) $(typeof(c))]") end if (!isfinite(reached_tol) || !(reached_tol < acceptance_tol)) && sylvester_algorithm ≠ :doubling @@ -125,7 +125,7 @@ verbose = verbose) if verbose# && i != 0 - println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: doubling") + println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: doubling [A: $(size(aa)) $(typeof(aa)), B: $(size(b)) $(typeof(b)), C: $(size(cc)) $(typeof(cc))]") end end @@ -148,7 +148,7 @@ verbose = verbose) if verbose && i != 0 - println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: bartels_stewart") + println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: bartels_stewart [A: $(size(aa)) $(typeof(aa)), B: $(size(bb)) $(typeof(bb)), C: $(size(cc)) $(typeof(cc))]") end end @@ -173,7 +173,7 @@ end if verbose# && i != 0 - println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $Reached_tol; algorithm: dqgmres (refinement of previous solution)") + println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $Reached_tol; algorithm: dqgmres (refinement of previous solution) [A: $(size(aa)) $(typeof(aa)), B: $(size(b)) $(typeof(b)), C: $(size(cc)) $(typeof(cc))]") end end @@ -193,7 +193,7 @@ verbose = verbose) if verbose# && i != 0 - println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: gmres") + println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $reached_tol; algorithm: gmres [A: $(size(aa)) $(typeof(aa)), B: $(size(b)) $(typeof(b)), C: $(size(cc)) $(typeof(cc))]") end end @@ -205,9 +205,9 @@ cc = 𝕊ℂ.𝐂¹ copyto!(cc, C) - X, i, Reached_tol = solve_sylvester_equation(aa, b, cc, + X, i, Reached_tol = solve_sylvester_equation(aa, b, cc, Val(:dqgmres), 𝕊ℂ, - initial_guess = x, + initial_guess = x, preconditioner = preconditioner, tol = tol, # timer = timer, @@ -218,7 +218,7 @@ end if verbose# && i != 0 - println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $Reached_tol; algorithm: dqgmres (refinement of previous solution)") + println("Sylvester equation - converged to tol $acceptance_tol: $(reached_tol < acceptance_tol); iterations: $i; reached tol: $Reached_tol; algorithm: dqgmres (refinement of previous solution) [A: $(size(aa)) $(typeof(aa)), B: $(size(b)) $(typeof(b)), C: $(size(cc)) $(typeof(cc))]") end end diff --git a/src/parser/equation_processing.jl b/src/parser/equation_processing.jl index 13c793e1b..b74cdf871 100644 --- a/src/parser/equation_processing.jl +++ b/src/parser/equation_processing.jl @@ -257,6 +257,7 @@ function process_model_equations(model_block_in::Expr, max_obc_horizon::Int, pre bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], eps()), min(bounds[x.args[2]][2], 1e12)) : (eps(), 1e12) x end : + !(x.args[2] isa Expr) ? x : x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin @@ -305,6 +306,7 @@ function process_model_equations(model_block_in::Expr, max_obc_horizon::Int, pre bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], eps()), min(bounds[x.args[2]][2], 1e12)) : (eps(), 1e12) x end : + !(x.args[2] isa Expr) ? x : x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # nonnegative variables begin @@ -393,6 +395,7 @@ function process_model_equations(model_block_in::Expr, max_obc_horizon::Int, pre bounds[x.args[2]] = haskey(bounds, x.args[2]) ? (max(bounds[x.args[2]][1], -1e12), min(bounds[x.args[2]][2], 600)) : (-1e12, 600) x end : + !(x.args[2] isa Expr) ? x : # exp of a constant (e.g. after shock → 0 substitution) x.args[2].head == :ref ? x.args[2].args[1] isa Symbol ? # have exp terms bound so they dont go to Inf begin @@ -849,7 +852,10 @@ function process_parameter_definitions(parameter_block_in::Expr, pmm::post_model par_defined_more_than_once = Set() bounded_vars = [] - parameter_definitions = replace_indices(parameter_block_in) + parameter_definitions = parse_for_loops(parameter_block_in) + parameter_definitions = resolve_if_expr(parameter_definitions) + parameter_definitions = remove_nothing(parameter_definitions) + parameter_definitions = replace_indices(parameter_definitions) # parse parameter inputs # label all variables parameters and exogenous variables and timings across all equations @@ -910,7 +916,7 @@ function process_parameter_definitions(parameter_block_in::Expr, pmm::post_model end : x : x : - x.args[2].head == :block ? + x.args[2] isa Expr && x.args[2].head == :block ? x.args[1].args[1] == :| ? x : x.args[2].args[2].args[1] == :| ? @@ -923,7 +929,7 @@ function process_parameter_definitions(parameter_block_in::Expr, pmm::post_model @warn "Invalid parameter input ignored: " * repr(x) x end : - x.args[2].head == :call ? + x.args[2] isa Expr && x.args[2].head == :call ? x.args[1].args[1] == :| ? x : begin # this is calibration by targeting SS values (conditional parameter at the end) @@ -978,7 +984,7 @@ function process_parameter_definitions(parameter_block_in::Expr, pmm::post_model x : begin diffed = intersect(setdiff([x], ss_tmp), get_symbols(cal_eq)) - if !isempty(diffed) + if !isempty(diffed) && diffed[1] ∉ SYMPYWORKSPACE_RESERVED_NAMES push!(par_tmp,diffed[1]) end end : @@ -1034,7 +1040,7 @@ function process_parameter_definitions(parameter_block_in::Expr, pmm::post_model x : begin diffed = setdiff([x],ss_tmp) - if !isempty(diffed) + if !isempty(diffed) && diffed[1] ∉ SYMPYWORKSPACE_RESERVED_NAMES push!(par_tmp,diffed[1]) end end : diff --git a/src/rrules.jl b/src/rrules.jl index 1e2aedaf0..b6673f27a 100644 --- a/src/rrules.jl +++ b/src/rrules.jl @@ -6563,8 +6563,8 @@ function mul_fill_kron_adjoint!(∂A::AbstractMatrix{R}, ∂B::AbstractMatrix{R}, M1::AbstractMatrix, M2::AbstractMatrix, - A::AbstractMatrix{TA}, - B::AbstractMatrix{TB}; + A::DenseMatrix{TA}, + B::DenseMatrix{TB}; tol::Real = 0.0, block::AbstractMatrix{R} = Matrix{R}(undef, size(M1, 1), 0)) where {R <: Real, TA <: Real, TB <: Real} n2, m2 = size(A) @@ -6696,73 +6696,211 @@ end # Mixed-sparsity variant: A is sparse, B is dense -function mul_fill_kron_adjoint!(∂A::AbstractMatrix{R}, +function mul_fill_kron_adjoint!( + ∂A::AbstractMatrix{R}, ∂B::AbstractMatrix{R}, M1::AbstractMatrix, M2::AbstractMatrix, - A::SparseMatrixCSC{TA, Int}, - B::AbstractMatrix{TB}; + A::SparseMatrixCSC{TA, Ti}, + B::DenseMatrix{TB}; tol::Real = 0.0, - block::AbstractMatrix{R} = Matrix{R}(undef, size(M1, 1), 0)) where {R <: Real, TA <: Real, TB <: Real} + block::AbstractMatrix{R} = Matrix{R}(undef, size(M1, 1), 0), +) where {R <: Real, TA <: Real, TB <: Real, Ti <: Integer} + n2, m2 = size(A) n1, m1 = size(B) - @assert size(M1, 1) == n1 * n2 - @assert size(M2, 2) == m1 * m2 - @assert size(M1, 2) == size(M2, 1) + size(M1, 1) == n1 * n2 || throw(DimensionMismatch("size(M1, 1) must equal size(B, 1) * size(A, 1)")) + size(M2, 2) == m1 * m2 || throw(DimensionMismatch("size(M2, 2) must equal size(B, 2) * size(A, 2)")) + size(M1, 2) == size(M2, 1) || throw(DimensionMismatch("inner dimensions of M1 and M2 do not match")) nrows = n1 * n2 - - if size(block, 1) == nrows && size(block, 2) >= m2 - blk = view(block, :, 1:m2) + tol_R = R(tol) + + # Transpose M2 so that M2t[(j - 1) * m2 + l, s] is contiguous in l. + M2t = Matrix{R}(transpose(M2)) + + # Row compressed representation of M1. + # Use conditional handling instead of try/catch to obtain non‑zero indices and values. + # If `M1` is already a sparse matrix we can call `findnz` directly; otherwise we + # convert it to a sparse representation first. This satisfies the request to + # avoid `try … catch` while preserving the original behaviour. + if typeof(M1) <: AbstractSparseMatrix + M1_row_ids, M1_col_ids, M1_values = findnz(M1) else - blk = Matrix{R}(undef, nrows, m2) + M1_row_ids, M1_col_ids, M1_values = findnz(sparse(M1)) end - A_colptr = A.colptr - A_rowval = A.rowval - A_nzval = A.nzval + M1_row_counts = zeros(Int, nrows) + + @inbounds for row in M1_row_ids + M1_row_counts[row] += 1 + end + + n_supported_rows = 0 + + @inbounds for row in 1:nrows + if M1_row_counts[row] != 0 + n_supported_rows += 1 + end + end + + row_to_slot = zeros(Int, nrows) + supported_rows = Vector{Int}(undef, n_supported_rows) + + slot = 0 + + @inbounds for row in 1:nrows + if M1_row_counts[row] != 0 + slot += 1 + row_to_slot[row] = slot + supported_rows[slot] = row + end + end + + M1_row_ptr = Vector{Int}(undef, n_supported_rows + 1) + M1_row_ptr[1] = 1 + + @inbounds for slot in 1:n_supported_rows + M1_row_ptr[slot + 1] = M1_row_ptr[slot] + M1_row_counts[supported_rows[slot]] + end + + M1_nz_col = Vector{Int}(undef, length(M1_values)) + M1_nz_val = Vector{R}(undef, length(M1_values)) + next_ptr = copy(M1_row_ptr) + + @inbounds for index in eachindex(M1_values) + row = M1_row_ids[index] + slot = row_to_slot[row] + position = next_ptr[slot] + next_ptr[slot] += 1 + + M1_nz_col[position] = M1_col_ids[index] + M1_nz_val[position] = R(M1_values[index]) + end + + M1_row_k = Vector{Int}(undef, n_supported_rows) + M1_row_i = Vector{Int}(undef, n_supported_rows) + + @inbounds for slot in 1:n_supported_rows + row = supported_rows[slot] + i = (row - 1) ÷ n2 + 1 + k = row - (i - 1) * n2 + + M1_row_k[slot] = k + M1_row_i[slot] = i + end + + # Row compressed representation of A. + A_row_ids, A_col_ids, A_values = findnz(A) + + A_row_counts = zeros(Int, n2) + + @inbounds for row in A_row_ids + A_row_counts[row] += 1 + end + + A_row_ptr = Vector{Int}(undef, n2 + 1) + A_row_ptr[1] = 1 + + @inbounds for row in 1:n2 + A_row_ptr[row + 1] = A_row_ptr[row] + A_row_counts[row] + end + + A_nz_col = Vector{Int}(undef, length(A_values)) + A_nz_val = Vector{R}(undef, length(A_values)) + next_A_ptr = copy(A_row_ptr) + + @inbounds for index in eachindex(A_values) + row = A_row_ids[index] + position = next_A_ptr[row] + next_A_ptr[row] += 1 + + A_nz_col[position] = A_col_ids[index] + A_nz_val[position] = R(A_values[index]) + end + + work = Vector{R}(undef, m2) @inbounds for j in 1:m1 - col_start = (j - 1) * m2 + 1 - col_end = j * m2 - ℒ.mul!(blk, M1, view(M2, :, col_start:col_end)) + column_offset = (j - 1) * m2 - re_blk = reshape(blk, n2, n1, m2) + for row_slot in 1:n_supported_rows + k = M1_row_k[row_slot] + i = M1_row_i[row_slot] - # ∂A[k,l] += B[i,j] * re_blk[k,i,l] — B is dense, use iszero guard - for l in 1:m2 - for i in 1:n1 - bij = B[i, j] - if abs(bij) > tol - for k in 1:n2 - ∂A[k, l] += bij * re_blk[k, i, l] - end + bij = R(B[i, j]) + has_B_contribution = abs(bij) > tol_R + + A_row_start = A_row_ptr[k] + A_row_stop = A_row_ptr[k + 1] - 1 + has_A_contribution = A_row_start <= A_row_stop + + if !(has_B_contribution || has_A_contribution) + continue + end + + M1_row_start = M1_row_ptr[row_slot] + M1_row_stop = M1_row_ptr[row_slot + 1] - 1 + + if has_B_contribution + fill!(work, zero(R)) + + for M1_index in M1_row_start:M1_row_stop + M1_value = M1_nz_val[M1_index] + M1_col = M1_nz_col[M1_index] + + @simd for l in 1:m2 + work[l] += M1_value * M2t[column_offset + l, M1_col] end end + + for l in 1:m2 + ∂A[k, l] += bij * work[l] end - # ∂B[i,j] += Σ_{k,l} A[k,l] * re_blk[k,i,l] — A is sparse, use nzrange - for i in 1:n1 + if has_A_contribution acc = zero(R) - for l in 1:m2 - for aidx in A_colptr[l]:(A_colptr[l + 1] - 1) - k = A_rowval[aidx] - acc += R(A_nzval[aidx]) * re_blk[k, i, l] + + for A_index in A_row_start:A_row_stop + l = A_nz_col[A_index] + acc += A_nz_val[A_index] * work[l] + end + + ∂B[i, j] += acc end + + else + acc = zero(R) + + for M1_index in M1_row_start:M1_row_stop + M1_value = M1_nz_val[M1_index] + M1_col = M1_nz_col[M1_index] + + local_acc = zero(R) + + for A_index in A_row_start:A_row_stop + l = A_nz_col[A_index] + local_acc += A_nz_val[A_index] * M2t[column_offset + l, M1_col] + end + + acc += M1_value * local_acc end + ∂B[i, j] += acc end end end + return ∂A, ∂B +end # Mixed-sparsity variant: A is dense, B is sparse function mul_fill_kron_adjoint!(∂A::AbstractMatrix{R}, ∂B::AbstractMatrix{R}, M1::AbstractMatrix, M2::AbstractMatrix, - A::AbstractMatrix{TA}, + A::DenseMatrix{TA}, B::SparseMatrixCSC{TB, Int}; tol::Real = 0.0, block::AbstractMatrix{R} = Matrix{R}(undef, size(M1, 1), 0)) where {R <: Real, TA <: Real, TB <: Real} @@ -7683,7 +7821,123 @@ end # Helper: adjoint of compressed_kron³(X) w.r.t. X. # Forward: out[row,col] = (aii*(ajj*akk + ajk*akj) + aij*(aji*akk + ajk*aki) + aik*(aji*akj + ajj*aki)) / divisor # where row ↔ (i1≥j1≥k1) and col ↔ (i2≥j2≥k2) and a_pq = X[p,q]. -function compressed_kron³_pullback!(∂X::AbstractMatrix{T}, ∂Y::AbstractMatrix{T}, X::AbstractMatrix{T}; tol::Real = 0.0) where T <: Real +function compressed_kron³_pullback!( + ∂X::AbstractMatrix{T}, + ∂Y::AbstractSparseMatrix{T, Ti}, + X::AbstractMatrix{T}; + tol::Real = 0.0, +) where {T <: Real, Ti <: Integer} + + Xd = X isa DenseMatrix ? X : collect(X) + + n_x_rows, n_x_cols = size(Xd) + + compressed_rows = n_x_rows * (n_x_rows + 1) * (n_x_rows + 2) ÷ 6 + compressed_cols = n_x_cols * (n_x_cols + 1) * (n_x_cols + 2) ÷ 6 + + size(∂Y, 1) == compressed_rows || throw(DimensionMismatch("size(∂Y, 1) is inconsistent with size(X, 1)")) + size(∂Y, 2) == compressed_cols || throw(DimensionMismatch("size(∂Y, 2) is inconsistent with size(X, 2)")) + + tol_T = T(tol) + + row_values = SparseArrays.rowvals(∂Y) + non_zero_values = nonzeros(∂Y) + + row_first = Vector{Int}(undef, compressed_rows) + row_second = Vector{Int}(undef, compressed_rows) + row_third = Vector{Int}(undef, compressed_rows) + row_scale = Vector{T}(undef, compressed_rows) + + @inbounds for i1 in 1:n_x_rows + for j1 in 1:i1 + for k1 in 1:j1 + row = (i1 - 1) * i1 * (i1 + 1) ÷ 6 + (j1 - 1) * j1 ÷ 2 + k1 + + row_first[row] = i1 + row_second[row] = j1 + row_third[row] = k1 + + if i1 == j1 + divisor = j1 == k1 ? 6 : 2 + else + divisor = (j1 == k1 || i1 == k1) ? 2 : 1 + end + + row_scale[row] = inv(T(divisor)) + end + end + end + + col_first = Vector{Int}(undef, compressed_cols) + col_second = Vector{Int}(undef, compressed_cols) + col_third = Vector{Int}(undef, compressed_cols) + + @inbounds for i2 in 1:n_x_cols + for j2 in 1:i2 + for k2 in 1:j2 + col = (i2 - 1) * i2 * (i2 + 1) ÷ 6 + (j2 - 1) * j2 ÷ 2 + k2 + + col_first[col] = i2 + col_second[col] = j2 + col_third[col] = k2 + end + end + end + + @inbounds for col in 1:compressed_cols + i2 = col_first[col] + j2 = col_second[col] + k2 = col_third[col] + + for stored_index in SparseArrays.nzrange(∂Y, col) + row = Int(row_values[stored_index]) + g = non_zero_values[stored_index] + + abs(g) <= tol_T && continue + + i1 = row_first[row] + j1 = row_second[row] + k1 = row_third[row] + + g_d = g * row_scale[row] + + aii = Xd[i1, i2] + aij = Xd[i1, j2] + aik = Xd[i1, k2] + + aji = Xd[j1, i2] + ajj = Xd[j1, j2] + ajk = Xd[j1, k2] + + aki = Xd[k1, i2] + akj = Xd[k1, j2] + akk = Xd[k1, k2] + + q_i2 = aij * ajk + aik * ajj + q_j2 = aik * aji + aii * ajk + q_k2 = aii * ajj + aij * aji + + ∂X[i1, i2] += g_d * (ajj * akk + ajk * akj) + ∂X[i1, j2] += g_d * (aji * akk + ajk * aki) + ∂X[i1, k2] += g_d * (aji * akj + ajj * aki) + + ∂X[j1, i2] += g_d * (aij * akk + aik * akj) + ∂X[j1, j2] += g_d * (aii * akk + aik * aki) + ∂X[j1, k2] += g_d * (aij * aki + aii * akj) + + ∂X[k1, i2] += g_d * q_i2 + ∂X[k1, j2] += g_d * q_j2 + ∂X[k1, k2] += g_d * q_k2 + end + end + + return ∂X +end + +function compressed_kron³_pullback!(∂X::AbstractMatrix{T}, + ∂Y::DenseMatrix{T}, + X::AbstractMatrix{T}; + tol::Real = 0.0) where T <: Real Xd = X isa DenseMatrix ? X : collect(X) n_rows, n_cols = size(Xd) # Unlike the forward pass, the pullback must iterate over ALL row/column @@ -7746,52 +8000,170 @@ end # Equivalent to: # compressed_kron³_pullback!(∂X, M1 * M2, X) # but avoids the n_compressed³ × n_compressed³ allocation. -function mul_compressed_kron³_pullback!(∂X::AbstractMatrix{T}, +function mul_compressed_kron³_pullback!( + ∂X::AbstractMatrix, M1::AbstractMatrix, M2::AbstractMatrix, - X::AbstractMatrix{T}; - tol::Real = 0.0) where T <: Real + X::AbstractMatrix; + tol::Real = 0.0, +) + Xd = X isa DenseMatrix ? X : collect(X) - n_rows, n_cols = size(Xd) - g_col = Vector{T}(undef, size(M1, 1)) + n_x_rows, n_x_cols = size(Xd) - for i2 in 1:n_cols, j2 in 1:i2 - for k2 in 1:j2 - col = (i2 - 1) * i2 * (i2 + 1) ÷ 6 + (j2 - 1) * j2 ÷ 2 + k2 + compressed_rows = n_x_rows * (n_x_rows + 1) * (n_x_rows + 2) ÷ 6 + compressed_cols = n_x_cols * (n_x_cols + 1) * (n_x_cols + 2) ÷ 6 - # Compute g_col = M1 * M2[:, col] lazily for this triple - ℒ.mul!(g_col, M1, view(M2, :, col)) + size(M1, 1) == compressed_rows || throw(DimensionMismatch("size(M1, 1) is inconsistent with size(X, 1)")) + size(M2, 2) == compressed_cols || throw(DimensionMismatch("size(M2, 2) is inconsistent with size(X, 2)")) + size(M1, 2) == size(M2, 1) || throw(DimensionMismatch("inner dimensions of M1 and M2 do not match")) - for i1 in 1:n_rows - # Hoist i1-dependent reads - @inbounds aii = Xd[i1, i2]; aij = Xd[i1, j2]; aik = Xd[i1, k2] + T = promote_type(eltype(∂X), eltype(M1), eltype(M2), eltype(Xd)) + tol_T = T(tol) + + # Avoid try/catch: directly handle dense vs sparse `M1`. + if typeof(M1) <: AbstractSparseMatrix + row_ids, col_ids, values = findnz(M1) + else + row_ids, col_ids, values = findnz(sparse(M1)) + end + + row_counts = zeros(Int, compressed_rows) + + @inbounds for row in row_ids + row_counts[row] += 1 + end + + n_supported_rows = 0 + + @inbounds for row in 1:compressed_rows + if row_counts[row] != 0 + n_supported_rows += 1 + end + end + + row_to_slot = zeros(Int, compressed_rows) + supported_rows = Vector{Int}(undef, n_supported_rows) + + slot = 0 + + @inbounds for row in 1:compressed_rows + if row_counts[row] != 0 + slot += 1 + row_to_slot[row] = slot + supported_rows[slot] = row + end + end + + row_ptr = Vector{Int}(undef, n_supported_rows + 1) + row_ptr[1] = 1 + + @inbounds for slot in 1:n_supported_rows + row_ptr[slot + 1] = row_ptr[slot] + row_counts[supported_rows[slot]] + end + + nz_col = Vector{Int}(undef, length(values)) + nz_val = Vector{T}(undef, length(values)) + next_ptr = copy(row_ptr) + + @inbounds for index in eachindex(values) + row = row_ids[index] + slot = row_to_slot[row] + position = next_ptr[slot] + next_ptr[slot] += 1 + + nz_col[position] = col_ids[index] + nz_val[position] = values[index] + end + + row_first = Vector{Int}(undef, n_supported_rows) + row_second = Vector{Int}(undef, n_supported_rows) + row_third = Vector{Int}(undef, n_supported_rows) + row_scale = Vector{T}(undef, n_supported_rows) + + @inbounds for i1 in 1:n_x_rows for j1 in 1:i1 - # Hoist j1-dependent reads - @inbounds aji = Xd[j1, i2]; ajj = Xd[j1, j2]; ajk = Xd[j1, k2] - # Precompute sub-expressions for ∂X[k1, ...] updates - q_i2 = aij * ajk + aik * ajj - q_j2 = aik * aji + aii * ajk - q_k2 = aii * ajj + aij * aji - @inbounds for k1 in 1:j1 + for k1 in 1:j1 row = (i1 - 1) * i1 * (i1 + 1) ÷ 6 + (j1 - 1) * j1 ÷ 2 + k1 - g = g_col[row] - abs(g) <= tol && continue + slot = row_to_slot[row] + + if slot != 0 + row_first[slot] = i1 + row_second[slot] = j1 + row_third[slot] = k1 + if i1 == j1 - divisor = (j1 == k1) ? 6 : 2 + divisor = j1 == k1 ? 6 : 2 else divisor = (j1 == k1 || i1 == k1) ? 2 : 1 end - g_d = g / divisor + + row_scale[slot] = inv(T(divisor)) + end + end + end + end + + work = zeros(T, compressed_cols) + + @inbounds for row_slot in 1:n_supported_rows + fill!(work, zero(T)) + + row_start = row_ptr[row_slot] + row_stop = row_ptr[row_slot + 1] - 1 + + for nz_index in row_start:row_stop + m1_value = nz_val[nz_index] + m1_col = nz_col[nz_index] + + @simd for col in 1:compressed_cols + work[col] += m1_value * M2[m1_col, col] + end + end + + i1 = row_first[row_slot] + j1 = row_second[row_slot] + k1 = row_third[row_slot] + scale = row_scale[row_slot] + + col = 0 + + for i2 in 1:n_x_cols + for j2 in 1:i2 + for k2 in 1:j2 + col += 1 + + g = work[col] + + abs(g) <= tol_T && continue + + g_d = g * scale + + aii = Xd[i1, i2] + aij = Xd[i1, j2] + aik = Xd[i1, k2] + + aji = Xd[j1, i2] + ajj = Xd[j1, j2] + ajk = Xd[j1, k2] + aki = Xd[k1, i2] akj = Xd[k1, j2] akk = Xd[k1, k2] + + q_i2 = aij * ajk + aik * ajj + q_j2 = aik * aji + aii * ajk + q_k2 = aii * ajj + aij * aji + ∂X[i1, i2] += g_d * (ajj * akk + ajk * akj) ∂X[i1, j2] += g_d * (aji * akk + ajk * aki) ∂X[i1, k2] += g_d * (aji * akj + ajj * aki) + ∂X[j1, i2] += g_d * (aij * akk + aik * akj) ∂X[j1, j2] += g_d * (aii * akk + aik * aki) ∂X[j1, k2] += g_d * (aij * aki + aii * akj) + ∂X[k1, i2] += g_d * q_i2 ∂X[k1, j2] += g_d * q_j2 ∂X[k1, k2] += g_d * q_k2 @@ -7799,7 +8171,8 @@ function mul_compressed_kron³_pullback!(∂X::AbstractMatrix{T}, end end end - end + + return ∂X end # ===================================================================================== @@ -8103,7 +8476,8 @@ function rrule(::typeof(calculate_third_order_solution), # ∂B_from_sylv = 𝐒₃_stable' * A' * ∂C_adj — reuse ∂𝐗₃ as temp ℒ.mul!(∂𝐗₃, At, ∂C_adj) ℒ.mul!(∂B_from_sylv, 𝐒₃_stable', ∂𝐗₃) - # ∂B_from_sylv = sparse(𝐒₃_stable' * ∂𝐗₃) + + ∂B_from_sylv = choose_matrix_format(∂B_from_sylv) # ∂𝐗₃ = spinv' * ∂C_adj (overwrite temp with real value) # ℒ.mul!(∂𝐗₃, sxpinv', ∂C_adj) ∂𝐗₃ = choose_matrix_format(spinv' * ∂C_adj) @@ -18415,6 +18789,43 @@ function contract_filter_free_me_std_cotangent(d_me_std::AbstractMatrix, user_me return out end +function expand_filter_free_matrix_cotangent(d_red::AbstractMatrix, + needed::AbstractVector{Int}, + n_vars_full::Int, + n_cols::Int) + d_full = zeros(eltype(d_red), n_vars_full, n_cols) + @inbounds d_full[needed, :] .= d_red + return d_full +end + +function expand_filter_free_state_cotangent(d_red::AbstractVector, + needed::AbstractVector{Int}, + n_vars_full::Int) + d_full = zeros(eltype(d_red), n_vars_full) + @inbounds d_full[needed] .= d_red + return d_full +end + +function expand_filter_free_state_cotangent(d_red::AbstractVector{<:AbstractVector}, + needed::AbstractVector{Int}, + n_vars_full::Int) + d_full = [zeros(eltype(d_red[i]), n_vars_full) for i in eachindex(d_red)] + @inbounds for i in eachindex(d_red) + d_full[i][needed] .= d_red[i] + end + return d_full +end + +function scatter_filter_free_ss_and_pars_cotangent(d_SS_obs::AbstractVector, + obs_indices::AbstractVector{Int}, + n_ss_and_pars::Int) + d_SS_and_pars = zeros(eltype(d_SS_obs), n_ss_and_pars) + @inbounds for k in eachindex(obs_indices) + d_SS_and_pars[obs_indices[k]] += d_SS_obs[k] + end + return d_SS_and_pars +end + # Visible-period filter-free pullbacks run on the reduced `needed` row slice of # the policy matrices. They propagate the state cotangent backward across only # the retained visible sample, accumulate reduced-matrix cotangents, and return @@ -19162,35 +19573,38 @@ function rrule(::typeof(get_loglikelihood), end d_shocks_full = expand_filter_free_shock_cotangent(d_shocks, shocks, visible_cols, n_warm) d_me_std_full = contract_filter_free_me_std_cotangent(expand_filter_free_me_std_cotangent(d_me_std, measurement_error_std, period_range), user_me_std) - d_𝐒₁_full_cot = zeros(eltype(d_𝐒₁_red), nVars_full, ncols₁) - @inbounds d_𝐒₁_full_cot[needed, :] .= d_𝐒₁_red - d_SS_and_pars = zeros(eltype(d_SS_obs), length(SS_and_pars)) - @inbounds for k in eachindex(obs_indices) - d_SS_and_pars[obs_indices[k]] += d_SS_obs[k] - end - # backprop through initial_state_deviations = initial_state - SS_and_pars[1:nVars] - # AND compute d_initial_state (positional tangent for the user-facing - # `initial_state` argument). - d_initial_state = NoTangent() + d_params = @thunk begin + d_𝐒₁_full_cot = expand_filter_free_matrix_cotangent(d_𝐒₁_red, needed, nVars_full, ncols₁) + d_SS_and_pars = scatter_filter_free_ss_and_pars_cotangent(d_SS_obs, obs_indices, length(SS_and_pars)) if has_override && !(d_state_red isa Union{NoTangent, AbstractZero}) - d_state_full_init = zeros(eltype(d_state_red), nVars_full_for_init) - @inbounds d_state_full_init[needed] .= d_state_red + d_state_full = expand_filter_free_state_cotangent(d_state_red, needed, nVars_full_for_init) + if initial_state_is_levels + @views d_SS_and_pars[1:nVars_full_for_init] .-= d_state_full + end + end + ss_pb((NoTangent(), d_SS_and_pars, d_𝐒₁_full_cot, NoTangent()))[3] + end + d_initial_state = if has_override + @thunk begin + d_𝐒₁_full_cot = expand_filter_free_matrix_cotangent(d_𝐒₁_red, needed, nVars_full, ncols₁) + d_SS_and_pars = scatter_filter_free_ss_and_pars_cotangent(d_SS_obs, obs_indices, length(SS_and_pars)) + d_state_full_init = expand_filter_free_state_cotangent(d_state_red, needed, nVars_full_for_init) if initial_state_is_levels @views d_SS_and_pars[1:nVars_full_for_init] .-= d_state_full_init - d_initial_state = d_state_full_init + initial_state_pullback_tangent(initial_state, d_state_full_init) else d_is_vec = Vector{Vector{eltype(d_state_red)}}(undef, n_overridden_components) d_is_vec[1] = d_state_full_init for k in 2:n_overridden_components d_is_vec[k] = zeros(eltype(d_state_red), length(initial_state[k])) end - d_initial_state = d_is_vec + initial_state_pullback_tangent(initial_state, d_is_vec) end end - # first_order ss rrule expects bare 𝐒₁ cotangent and ignores Δstate - ss_grads = ss_pb((NoTangent(), d_SS_and_pars, d_𝐒₁_full_cot, NoTangent())) - d_params = ss_grads[3] - return NoTangent(), NoTangent(), NoTangent(), d_params, d_shocks_full, d_me_std_full, initial_state_pullback_tangent(initial_state, d_initial_state) + else + NoTangent() + end + return NoTangent(), NoTangent(), NoTangent(), d_params, d_shocks_full, d_me_std_full, d_initial_state end return isfinite(llh) ? (llh, pullback) : on_failure @@ -19590,4 +20004,4 @@ function rrule(::typeof(get_loglikelihood), return (t[1], t[2], t[3], t[4], t[5], t[6]) end return y, pb_short -end +end \ No newline at end of file diff --git a/src/steady_state/solver_parameters.jl b/src/steady_state/solver_parameters.jl index cde260f68..e4ccdf355 100644 --- a/src/steady_state/solver_parameters.jl +++ b/src/steady_state/solver_parameters.jl @@ -28,5 +28,7 @@ const DEFAULT_SOLVER_PARAMETERS = [ solver_parameters(4.1784912636092235, 1.8166012668623566, 0.5168801279930487, 78.18194336881028, 2.139580134601701, 0.4617967010780055, 33.95219683424897, 17.315839925955242, 2.220446049250313e-16, 12.287343174930065, 2.220446049250313e-16, 6.185479065850274, 88.3014875814592, 36.31304631280673, 5.262437586106421, 2.220446049250313e-16, 2.220446049250313e-16, 6.347784900438273, 0.7130503478600859, 0.6594888633818169, 1, 0.0, 2), - solver_parameters(75.41767502114854, 48.44868207231484, 66.70557675476336, 17.001461658038423, 2.715293366195093, 1.2745896783633328, 37.70721763849395, 50.739875243093444, 66.72525429469775, 54.137861579508154, 12.078847271504216, 19.723433260864525, 2.2648933923720476, 54.735557478829, 98.94783074858547, 17.863177309960086, 85.97559329517274, 64.79678330684743, 26.59975637589043, 24.72319192940016, 1, 0.0, 2) + solver_parameters(75.41767502114854, 48.44868207231484, 66.70557675476336, 17.001461658038423, 2.715293366195093, 1.2745896783633328, 37.70721763849395, 50.739875243093444, 66.72525429469775, 54.137861579508154, 12.078847271504216, 19.723433260864525, 2.2648933923720476, 54.735557478829, 98.94783074858547, 17.863177309960086, 85.97559329517274, 64.79678330684743, 26.59975637589043, 24.72319192940016, 1, 0.0, 2), + + solver_parameters(10.327106226758554, 10.327106226758554, 1.0821485512186284, 52.71196149219814, 2.53193716831508, 1.7586052095638707, 50.434347492548646, 1.7535573088570338, 1.51581318888478, 0.943191530494055, 1.8556283860336538, 6.8221667685490015, 72.18579387599354, 7.7659847231914725, 1.6450426805205853, 1.7020036851998528, 2.352764735403362, 0.24785068495263674, 5.649433169080131, 0.9245175406565802, 1, 0.0, 2) ]