rm SpeedMapping, quad_mat_iter, lin_time_iter, riccati, binder_pesaran

thorek1 · Dec 18, 2024 · 3e9fcb1 · 3e9fcb1
1 parent afccd8f
commit 3e9fcb1
Show file tree

Hide file tree

Showing 9 changed files with 268 additions and 278 deletions.
diff --git a/Project.toml b/Project.toml
@@ -35,7 +35,6 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
@@ -88,7 +87,6 @@ Requires = "1"
 RuntimeGeneratedFunctions = "0.5"
 SparseArrays = "1"
 SpecialFunctions = "2"
-SpeedMapping = "0.3, 0.4"
 StatsPlots = "0.15"
 Subscripts = "0.1.3, 0.2"
 Suppressor = "0.2"

diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl
@@ -44,7 +44,7 @@ import Subscripts: super, sub
 import Krylov
 import LinearOperators
 import DataStructures: CircularBuffer
-import SpeedMapping: speedmapping
+# import SpeedMapping: speedmapping
 import Suppressor: @suppress
 import REPL
 import Unicode
@@ -148,7 +148,7 @@ export get_equations, get_steady_state_equations, get_dynamic_equations, get_cal
 export irf, girf
 
 # Remove comment for debugging
-# export riccati_forward, block_solver, remove_redundant_SS_vars!, write_parameters_input!, parse_variables_input_to_index, undo_transformer , transformer, calculate_third_order_stochastic_steady_state, calculate_second_order_stochastic_steady_state, filter_and_smooth
+# export block_solver, remove_redundant_SS_vars!, write_parameters_input!, parse_variables_input_to_index, undo_transformer , transformer, calculate_third_order_stochastic_steady_state, calculate_second_order_stochastic_steady_state, filter_and_smooth
 # export create_symbols_eqs!, solve_steady_state!, write_functions_mapping!, solve!, parse_algorithm_to_state_update, block_solver, block_solver_AD, calculate_covariance, calculate_jacobian, calculate_first_order_solution, expand_steady_state, get_symbols, calculate_covariance_AD, parse_shocks_input_to_index
 
 
@@ -4748,7 +4748,7 @@ end
 function solve!(𝓂::ℳ; 
     parameters::ParameterType = nothing, 
     dynamics::Bool = false, 
-    algorithm::Symbol = :riccati, 
+    algorithm::Symbol = :first_order, 
     obc::Bool = false,
     verbose::Bool = false,
     silent::Bool = false,
@@ -4778,12 +4778,8 @@ function solve!(𝓂::ℳ;
 
     if dynamics
         obc_not_solved = isnothing(𝓂.solution.perturbation.first_order.state_update_obc)
-        if  ((:riccati             == algorithm) && ((:riccati             ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
-            ((:first_order         == algorithm) && ((:first_order         ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
+        if  ((:first_order         == algorithm) && ((:first_order         ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
             ((:first_order_doubling         == algorithm) && ((:first_order_doubling         ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
-            ((:binder_pesaran  == algorithm) && ((:binder_pesaran   ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
-            ((:quadratic_iteration  == algorithm) && ((:quadratic_iteration   ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
-            ((:linear_time_iteration  == algorithm) && ((:linear_time_iteration   ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
             ((:second_order        == algorithm) && ((:second_order        ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
             ((:pruned_second_order == algorithm) && ((:pruned_second_order ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
             ((:third_order         == algorithm) && ((:third_order         ∈ 𝓂.solution.outdated_algorithms) || (obc && obc_not_solved))) ||
@@ -4807,10 +4803,6 @@ function solve!(𝓂::ℳ;
 
             if algorithm == :first_order_doubling 
                 qme_solver = :doubling
-            elseif algorithm ∈ [:quadratic_iteration, :binder_pesaran]
-                qme_solver = :quadratic_iteration
-            elseif algorithm == :linear_time_iteration 
-                qme_solver = :linear_time_iteration
             else # default option
                 qme_solver = :schur
             end
@@ -4850,7 +4842,7 @@ function solve!(𝓂::ℳ;
             end
 
             𝓂.solution.perturbation.first_order = perturbation_solution(S₁, state_update₁, state_update₁̂)
-            𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:riccati, :first_order])
+            𝓂.solution.outdated_algorithms = setdiff(𝓂.solution.outdated_algorithms,[:first_order, :first_order_doubling])
 
             𝓂.solution.non_stochastic_steady_state = SS_and_pars
             𝓂.solution.outdated_NSSS = solution_error > tol
@@ -6925,7 +6917,7 @@ end
 
 function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasionally_binding_constraints::Bool)::Tuple{Function,Bool}
     if occasionally_binding_constraints
-        if algorithm ∈ [:riccati, :first_order, :first_order_doubling, :linear_time_iteration, :quadratic_iteration, :binder_pesaran]
+        if algorithm ∈ [:first_order, :first_order_doubling]
             state_update = 𝓂.solution.perturbation.first_order.state_update_obc
             pruning = false
         elseif :second_order == algorithm
@@ -6942,7 +6934,7 @@ function parse_algorithm_to_state_update(algorithm::Symbol, 𝓂::ℳ, occasiona
             pruning = true
         end
     else
-        if algorithm ∈ [:riccati, :first_order, :linear_time_iteration, :quadratic_iteration, :binder_pesaran]
+        if algorithm ∈ [:first_order, :first_order_doubling]
             state_update = 𝓂.solution.perturbation.first_order.state_update
             pruning = false
         elseif :second_order == algorithm

diff --git a/src/algorithms/lyapunov.jl b/src/algorithms/lyapunov.jl
@@ -519,27 +519,27 @@ function solve_lyapunov_equation(A::AbstractMatrix{Float64},
 end
 
 
-function solve_lyapunov_equation(A::AbstractMatrix{Float64},
-                                    C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
-                                    ::Val{:speedmapping};
-                                    tol::AbstractFloat = 1e-14,
-                                    timer::TimerOutput = TimerOutput())
-    𝐂A = similar(C)
-
-    soll = speedmapping(C; 
-            m! = (X, x) -> begin
-                ℒ.mul!(𝐂A, x, A')
-                ℒ.mul!(X, A, 𝐂A)
-                ℒ.axpy!(1, C, X)
-            end, stabilize = false, maps_limit = 1000, tol = tol)
+# function solve_lyapunov_equation(A::AbstractMatrix{Float64},
+#                                     C::Union{ℒ.Adjoint{Float64,Matrix{Float64}},DenseMatrix{Float64}},
+#                                     ::Val{:speedmapping};
+#                                     tol::AbstractFloat = 1e-14,
+#                                     timer::TimerOutput = TimerOutput())
+#     𝐂A = similar(C)
+
+#     soll = speedmapping(C; 
+#             m! = (X, x) -> begin
+#                 ℒ.mul!(𝐂A, x, A')
+#                 ℒ.mul!(X, A, 𝐂A)
+#                 ℒ.axpy!(1, C, X)
+#             end, stabilize = false, maps_limit = 1000, tol = tol)
 
-    𝐂 = soll.minimizer
+#     𝐂 = soll.minimizer
 
-    reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂)
+#     reached_tol = ℒ.norm(A * 𝐂 * A' + C - 𝐂) / ℒ.norm(𝐂)
 
-    # if reached_tol > tol
-    #     println("Lyapunov: speedmapping $reached_tol")
-    # end
+#     # if reached_tol > tol
+#     #     println("Lyapunov: speedmapping $reached_tol")
+#     # end
 
-    return 𝐂, soll.maps, reached_tol
-end
+#     return 𝐂, soll.maps, reached_tol
+# end
diff --git a/src/algorithms/quadratic_matrix_equation.jl b/src/algorithms/quadratic_matrix_equation.jl
@@ -394,112 +394,112 @@ end
 
 
 
-function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, 
-                                        B::AbstractMatrix{R}, 
-                                        C::AbstractMatrix{R}, 
-                                        ::Val{:linear_time_iteration}, 
-                                        T::timings; 
-                                        initial_guess::AbstractMatrix{R} = zeros(0,0),
-                                        tol::AbstractFloat = 1e-14, # lower tol not possible for NAWM (and probably other models this size)
-                                        timer::TimerOutput = TimerOutput(),
-                                        verbose::Bool = false,
-                                        max_iter::Int = 5000) where R <: Real
-    # Iterate over: A * X * X + C + B * X = (A * X + B) * X + C = X + (A * X + B) \ C
+# function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, 
+#                                         B::AbstractMatrix{R}, 
+#                                         C::AbstractMatrix{R}, 
+#                                         ::Val{:linear_time_iteration}, 
+#                                         T::timings; 
+#                                         initial_guess::AbstractMatrix{R} = zeros(0,0),
+#                                         tol::AbstractFloat = 1e-14, # lower tol not possible for NAWM (and probably other models this size)
+#                                         timer::TimerOutput = TimerOutput(),
+#                                         verbose::Bool = false,
+#                                         max_iter::Int = 5000) where R <: Real
+#     # Iterate over: A * X * X + C + B * X = (A * X + B) * X + C = X + (A * X + B) \ C
 
-    @timeit_debug timer "Preallocate" begin
+#     @timeit_debug timer "Preallocate" begin
 
-    F = similar(C)
-    t = similar(C)
+#     F = similar(C)
+#     t = similar(C)
 
-    initial_guess = length(initial_guess) > 0 ? initial_guess : zero(A)
+#     initial_guess = length(initial_guess) > 0 ? initial_guess : zero(A)
 
-    end # timeit_debug
-    @timeit_debug timer "Loop" begin
+#     end # timeit_debug
+#     @timeit_debug timer "Loop" begin
 
-    sol = @suppress begin 
-        speedmapping(initial_guess; m! = (F̄, F) -> begin 
-            ℒ.mul!(t, A, F)
-            ℒ.axpby!(-1, B, -1, t)
-            t̂ = ℒ.lu!(t, check = false)
-            ℒ.ldiv!(F̄, t̂, C)
-        end,
-        tol = tol, maps_limit = max_iter)#, σ_min = 0.0, stabilize = false, orders = [3,3,2])
-    end
+#     sol = @suppress begin 
+#         speedmapping(initial_guess; m! = (F̄, F) -> begin 
+#             ℒ.mul!(t, A, F)
+#             ℒ.axpby!(-1, B, -1, t)
+#             t̂ = ℒ.lu!(t, check = false)
+#             ℒ.ldiv!(F̄, t̂, C)
+#         end,
+#         tol = tol, maps_limit = max_iter)#, σ_min = 0.0, stabilize = false, orders = [3,3,2])
+#     end
 
-    X = sol.minimizer
+#     X = sol.minimizer
 
-    AXX = A * X^2
+#     AXX = A * X^2
 
-    AXXnorm = ℒ.norm(AXX)
+#     AXXnorm = ℒ.norm(AXX)
 
-    ℒ.mul!(AXX, B, X, 1, 1)
+#     ℒ.mul!(AXX, B, X, 1, 1)
 
-    ℒ.axpy!(1, C, AXX)
+#     ℒ.axpy!(1, C, AXX)
 
-    reached_tol = ℒ.norm(AXX) / AXXnorm
+#     reached_tol = ℒ.norm(AXX) / AXXnorm
 
-    end # timeit_debug
+#     end # timeit_debug
 
-    # if reached_tol > tol
-    #     println("QME: linear time iteration $reached_tol")
-    # end
+#     # if reached_tol > tol
+#     #     println("QME: linear time iteration $reached_tol")
+#     # end
 
-    return X, sol.maps, reached_tol
-end
+#     return X, sol.maps, reached_tol
+# end
 
 
 
-function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, 
-                                        B::AbstractMatrix{R}, 
-                                        C::AbstractMatrix{R}, 
-                                        ::Val{:quadratic_iteration}, 
-                                        T::timings; 
-                                        initial_guess::AbstractMatrix{R} = zeros(0,0),
-                                        tol::AbstractFloat = 1e-14, # lower tol not possible for NAWM (and probably other models this size)
-                                        timer::TimerOutput = TimerOutput(),
-                                        verbose::Bool = false,
-                                        max_iter::Int = 50000) where R <: Real
-    # Iterate over: Ā * X * X + C̄ + X
+# function solve_quadratic_matrix_equation(A::AbstractMatrix{R}, 
+#                                         B::AbstractMatrix{R}, 
+#                                         C::AbstractMatrix{R}, 
+#                                         ::Val{:quadratic_iteration}, 
+#                                         T::timings; 
+#                                         initial_guess::AbstractMatrix{R} = zeros(0,0),
+#                                         tol::AbstractFloat = 1e-14, # lower tol not possible for NAWM (and probably other models this size)
+#                                         timer::TimerOutput = TimerOutput(),
+#                                         verbose::Bool = false,
+#                                         max_iter::Int = 50000) where R <: Real
+#     # Iterate over: Ā * X * X + C̄ + X
 
-    B̂ =  ℒ.lu(B)
+#     B̂ =  ℒ.lu(B)
 
-    C̄ = B̂ \ C
-    Ā = B̂ \ A
+#     C̄ = B̂ \ C
+#     Ā = B̂ \ A
 
-    X = similar(Ā)
-    X̄ = similar(Ā)
+#     X = similar(Ā)
+#     X̄ = similar(Ā)
 
-    X² = similar(X)
+#     X² = similar(X)
 
-    initial_guess = length(initial_guess) > 0 ? initial_guess : zero(A)
+#     initial_guess = length(initial_guess) > 0 ? initial_guess : zero(A)
 
-    sol = @suppress begin
-        speedmapping(initial_guess; m! = (X̄, X) -> begin 
-                                                ℒ.mul!(X², X, X)
-                                                ℒ.mul!(X̄, Ā, X²)
-                                                ℒ.axpy!(1, C̄, X̄)
-                                            end,
-        tol = tol, maps_limit = max_iter)#, σ_min = 0.0, stabilize = false, orders = [3,3,2])
-    end
+#     sol = @suppress begin
+#         speedmapping(initial_guess; m! = (X̄, X) -> begin 
+#                                                 ℒ.mul!(X², X, X)
+#                                                 ℒ.mul!(X̄, Ā, X²)
+#                                                 ℒ.axpy!(1, C̄, X̄)
+#                                             end,
+#         tol = tol, maps_limit = max_iter)#, σ_min = 0.0, stabilize = false, orders = [3,3,2])
+#     end
 
-    X = -sol.minimizer
+#     X = -sol.minimizer
 
-    AXX = A * X^2
+#     AXX = A * X^2
 
-    AXXnorm = ℒ.norm(AXX)
+#     AXXnorm = ℒ.norm(AXX)
 
-    ℒ.mul!(AXX, B, X, 1, 1)
+#     ℒ.mul!(AXX, B, X, 1, 1)
 
-    ℒ.axpy!(1, C, AXX)
+#     ℒ.axpy!(1, C, AXX)
 
-    reached_tol = ℒ.norm(AXX) / AXXnorm
+#     reached_tol = ℒ.norm(AXX) / AXXnorm
 
-    # if reached_tol > tol
-    #     println("QME: quadratic iteration $reached_tol")
-    # end
+#     # if reached_tol > tol
+#     #     println("QME: quadratic iteration $reached_tol")
+#     # end
 
-    return X, sol.maps, reached_tol
-end
+#     return X, sol.maps, reached_tol
+# end