diff --git a/src/filter/inversion.jl b/src/filter/inversion.jl
index dd4420f0..1e15ab6b 100644
--- a/src/filter/inversion.jl
+++ b/src/filter/inversion.jl
@@ -3488,21 +3488,25 @@ end
 
 
 function filter_data_with_model(𝓂::ℳ,
-    data_in_deviations::KeyedArray{Float64},
-    ::Val{:second_order}, # algo
-    ::Val{:inversion}; # filter
-    quadratic_matrix_equation_solver::Symbol = :schur,
-    warmup_iterations::Int = 0,
-    filter_algorithm::Symbol = :LagrangeNewton,
-    smooth::Bool = true,
-    verbose::Bool = false)
+                                data_in_deviations::KeyedArray{Float64},
+                                ::Val{:second_order}, # algo
+                                ::Val{:inversion}; # filter
+                                quadratic_matrix_equation_solver::Symbol = :schur,
+                                sylvester_algorithm::Symbol = :doubling, 
+                                warmup_iterations::Int = 0,
+                                filter_algorithm::Symbol = :LagrangeNewton,
+                                smooth::Bool = true,
+                                verbose::Bool = false)
 
     T = 𝓂.timings
 
     variables = zeros(T.nVars, size(data_in_deviations,2))
     shocks = zeros(T.nExo, size(data_in_deviations,2))
 
-    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, verbose = verbose)
+    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, 
+                            quadratic_matrix_equation_solver = quadratic_matrix_equation_solver,
+                            sylvester_algorithm = sylvester_algorithm,
+                            verbose = verbose)
 
     if !converged || solution_error > 1e-12
         @error "Could not find 2nd order stochastic steady state"
@@ -3708,6 +3712,7 @@ function filter_data_with_model(𝓂::ℳ,
                                 ::Val{:pruned_second_order}, # algo
                                 ::Val{:inversion}; # filter
                                 quadratic_matrix_equation_solver::Symbol = :schur,
+                                sylvester_algorithm::Symbol = :doubling, 
                                 warmup_iterations::Int = 0,
                                 filter_algorithm::Symbol = :LagrangeNewton,
                                 smooth::Bool = true,
@@ -3720,7 +3725,11 @@ function filter_data_with_model(𝓂::ℳ,
 
     observables = get_and_check_observables(𝓂, data_in_deviations)
     
-    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, pruning = true, verbose = verbose)
+    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, 𝐒₁, 𝐒₂ = calculate_second_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, 
+                                quadratic_matrix_equation_solver = quadratic_matrix_equation_solver,
+                                sylvester_algorithm = sylvester_algorithm,
+                                pruning = true, 
+                                verbose = verbose)
 
     if solution_error > 1e-12 || isnan(solution_error)
         @error "No solution for these parameters."
@@ -3993,7 +4002,9 @@ function filter_data_with_model(𝓂::ℳ,
     
     observables = get_and_check_observables(𝓂, data_in_deviations)
 
-    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, verbose = verbose)
+    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, 
+                                quadratic_matrix_equation_solver = quadratic_matrix_equation_solver,
+                                verbose = verbose)
 
     if !converged || solution_error > 1e-12
         @error "Could not find 3rd order stochastic steady state"
@@ -4305,7 +4316,10 @@ function filter_data_with_model(𝓂::ℳ,
     
     observables = get_and_check_observables(𝓂, data_in_deviations)
 
-    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, pruning = true, verbose = verbose)
+    sss, converged, SS_and_pars, solution_error, ∇₁, ∇₂, ∇₃, 𝐒₁, 𝐒₂, 𝐒₃ = calculate_third_order_stochastic_steady_state(𝓂.parameter_values, 𝓂, 
+                    pruning = true, 
+                    quadratic_matrix_equation_solver = quadratic_matrix_equation_solver,
+                    verbose = verbose)
 
     if !converged || solution_error > 1e-12
         @error "Could not find pruned 3rd order stochastic steady state"
diff --git a/src/filter/kalman.jl b/src/filter/kalman.jl
index d51369dc..98194a06 100644
--- a/src/filter/kalman.jl
+++ b/src/filter/kalman.jl
@@ -560,6 +560,7 @@ function filter_data_with_model(𝓂::ℳ,
     ::Val{:first_order}, # algo
     ::Val{:kalman}; # filter
     quadratic_matrix_equation_solver::Symbol = :schur,
+    lyapunov_algorithm::Symbol = :doubling,
     warmup_iterations::Int = 0,
     smooth::Bool = true,
     verbose::Bool = false)
@@ -568,7 +569,10 @@ function filter_data_with_model(𝓂::ℳ,
 
     obs_symbols = obs_axis isa String_input ? obs_axis .|> Meta.parse .|> replace_indices : obs_axis
 
-    filtered_and_smoothed = filter_and_smooth(𝓂, data_in_deviations, obs_symbols; verbose = verbose)
+    filtered_and_smoothed = filter_and_smooth(𝓂, data_in_deviations, obs_symbols; 
+                                                quadratic_matrix_equation_solver = quadratic_matrix_equation_solver,
+                                                lyapunov_algorithm = lyapunov_algorithm,
+                                                verbose = verbose)
 
     variables           = filtered_and_smoothed[smooth ? 1 : 5]
     standard_deviations = filtered_and_smoothed[smooth ? 2 : 6]
@@ -584,6 +588,7 @@ function filter_and_smooth(𝓂::ℳ,
                             data_in_deviations::AbstractArray{Float64}, 
                             observables::Vector{Symbol}; 
                             verbose::Bool = false, 
+                            quadratic_matrix_equation_solver::Symbol = :schur,
                             lyapunov_algorithm::Symbol = :doubling,
                             tol::AbstractFloat = 1e-12)
     # Based on Durbin and Koopman (2012)
@@ -604,7 +609,10 @@ function filter_and_smooth(𝓂::ℳ,
 
 	∇₁ = calculate_jacobian(parameters, SS_and_pars, 𝓂)# |> Matrix
 
-    sol, qme_sol, solved = calculate_first_order_solution(∇₁; T = 𝓂.timings, initial_guess = 𝓂.solution.perturbation.qme_solution, verbose = verbose)
+    sol, qme_sol, solved = calculate_first_order_solution(∇₁; T = 𝓂.timings, 
+                                                            quadratic_matrix_equation_solver = quadratic_matrix_equation_solver,
+                                                            initial_guess = 𝓂.solution.perturbation.qme_solution, 
+                                                            verbose = verbose)
 
     if solved 𝓂.solution.perturbation.qme_solution = qme_sol end