Add quadrature algorithm for stochastic steady state calculation

Project.toml

@@ -29,6 +29,7 @@ NLopt = "76087f3c-5699-56af-9a33-bf431cd00edd"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
+QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
@@ -89,6 +90,7 @@ Pigeons = "0.3, 0.4"
 PrecompileTools = "1"
 Preferences = "1"
 PythonCall = "0.9"
+QuasiMonteCarlo = "0.3"
 REPL = "1"
 Random = "1"
 RecursiveFactorization = "0.2"

src/MacroModelling.jl

@@ -5,7 +5,7 @@ import DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF, TYPEDSIGNATURES, TYPEDF
 # import StatsFuns: normcdf
 import ThreadedSparseArrays
 using PrecompileTools
-import SpecialFunctions: erfcinv, erfc # can't use constants because of SymPy (e.g. sqrt2)
+import SpecialFunctions: erfcinv, erfc, erfinv # can't use constants because of SymPy (e.g. sqrt2)
 import SpecialFunctions
 import SymPyPythonCall as SPyPyC
 import PythonCall
@@ -48,6 +48,7 @@ import Krylov
 import Krylov: GmresWorkspace, DqgmresWorkspace, BicgstabWorkspace
 import LinearOperators
 import DataStructures: CircularBuffer
+import QuasiMonteCarlo
 import MacroTools: unblock, postwalk, prewalk, @capture, flatten
 
 # import SpeedMapping: speedmapping
@@ -107,15 +108,16 @@ using DispatchDoctor
 
 # Imports
 include("default_options.jl")
-include("common_docstrings.jl")
 include("options_and_caches.jl")
+include("common_docstrings.jl")
 include("structures.jl")
 include("macros.jl")
 include("get_functions.jl")
 include("dynare.jl")
 include("inspect.jl")
 include("moments.jl")
 include("perturbation.jl")
+include("quadrature_sss.jl")
 
 include("./algorithms/sylvester.jl")
 include("./algorithms/lyapunov.jl")

src/get_functions.jl

@@ -1487,21 +1487,33 @@ function get_steady_state(𝓂::ℳ;
     end
 
     if stochastic
-        solve!(𝓂, 
-                opts = opts, 
-                dynamics = true, 
-                algorithm = algorithm, 
-                silent = silent, 
-                obc = length(𝓂.obc_violation_equations) > 0)
-
-        if  algorithm == :third_order
-            SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.third_order.stochastic_steady_state
-        elseif  algorithm == :pruned_third_order
-            SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.pruned_third_order.stochastic_steady_state
-        elseif  algorithm == :pruned_second_order
-            SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.pruned_second_order.stochastic_steady_state
+        if algorithm == :quadrature
+            # For quadrature algorithm, compute SSS directly without calling solve!
+            # The quadrature algorithm handles the solution internally
+            sss, converged, _, _, _, _, _, _ = calculate_quadrature_stochastic_steady_state(𝓂.parameter_values, 𝓂, opts = opts)
+
+            if converged
+                SS[1:length(𝓂.var)] = sss[1:length(𝓂.var)]
+            else
+                @warn "Quadrature stochastic steady state calculation did not converge. Using NSSS instead."
+            end
         else
-            SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.second_order.stochastic_steady_state#[indexin(sort(union(𝓂.var,𝓂.exo_present)),sort(union(𝓂.var,𝓂.aux,𝓂.exo_present)))]
+            solve!(𝓂, 
+                    opts = opts, 
+                    dynamics = true, 
+                    algorithm = algorithm, 
+                    silent = silent, 
+                    obc = length(𝓂.obc_violation_equations) > 0)
+
+            if  algorithm == :third_order
+                SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.third_order.stochastic_steady_state
+            elseif  algorithm == :pruned_third_order
+                SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.pruned_third_order.stochastic_steady_state
+            elseif  algorithm == :pruned_second_order
+                SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.pruned_second_order.stochastic_steady_state
+            else
+                SS[1:length(𝓂.var)] = 𝓂.solution.perturbation.second_order.stochastic_steady_state#[indexin(sort(union(𝓂.var,𝓂.exo_present)),sort(union(𝓂.var,𝓂.aux,𝓂.exo_present)))]
+            end
         end
     end
 
@@ -1534,7 +1546,16 @@ function get_steady_state(𝓂::ℳ;
 
     if derivatives 
         if stochastic
-                if algorithm == :third_order
+                if algorithm == :quadrature
+                    # Calculate derivatives for quadrature algorithm using automatic differentiation
+                    dSSS = 𝒟.jacobian(x->begin 
+                                SSS = calculate_quadrature_stochastic_steady_state(x, 𝓂, opts = opts)
+                                return [collect(SSS[1])[var_idx]...,collect(SSS[3])[calib_idx]...]
+                            end, backend, 𝓂.parameter_values)[:,param_idx]
+
+                    return KeyedArray(hcat(SS[[var_idx...,calib_idx...]], dSSS);  Variables_and_calibrated_parameters = axis1, Steady_state_and_∂steady_state∂parameter = axis2)
+
+                elseif algorithm == :third_order
 
                     # dSSS = 𝒜.jacobian(𝒷(), x->begin 
                     #             SSS = SSS_third_order_parameter_derivatives(x, param_idx, 𝓂, verbose = verbose)

src/macros.jl

@@ -94,6 +94,7 @@ macro model(𝓂,ex...)
     SS_solve_func = x->x
     # SS_calib_func = x->x
     SS_check_func = x->x
+    quadrature_dynamic_func = x->x
     ∂SS_equations_∂parameters = (zeros(0,0), x->x) # ([], SparseMatrixCSC{Float64, Int64}(ℒ.I, 0, 0))
     ∂SS_equations_∂SS_and_pars = (zeros(0,0), x->x)  # ([], Int[], zeros(1,1))
     SS_dependencies = nothing