diff --git a/Project.toml b/Project.toml index d52ccf3c..ea6813ee 100644 --- 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 index 0611c3a8..7e27e931 100644 --- 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 index b009237e..d9ae8ba0 100644 --- 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 index 6dcd6931..f7f3643d 100644 --- 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 diff --git a/src/algorithms/sylvester.jl b/src/algorithms/sylvester.jl index 3df4562d..5576e6b0 100644 --- a/src/algorithms/sylvester.jl +++ b/src/algorithms/sylvester.jl @@ -1404,46 +1404,46 @@ function solve_sylvester_equation(A::AbstractMatrix{Float64}, end -function solve_sylvester_equation(A::AbstractMatrix{Float64}, - B::AbstractMatrix{Float64}, - C::AbstractMatrix{Float64}, - ::Val{:speedmapping}; - initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), - timer::TimerOutput = TimerOutput(), - verbose::Bool = false, - tol::AbstractFloat = 1e-14) - # guess_provided = true +# function solve_sylvester_equation(A::AbstractMatrix{Float64}, +# B::AbstractMatrix{Float64}, +# C::AbstractMatrix{Float64}, +# ::Val{:speedmapping}; +# initial_guess::AbstractMatrix{<:AbstractFloat} = zeros(0,0), +# timer::TimerOutput = TimerOutput(), +# verbose::Bool = false, +# tol::AbstractFloat = 1e-14) +# # guess_provided = true - if length(initial_guess) == 0 - # guess_provided = false - initial_guess = zero(C) - end +# if length(initial_guess) == 0 +# # guess_provided = false +# initial_guess = zero(C) +# end - 𝐂 = A * initial_guess * B + C - initial_guess - # 𝐂 = copy(C) +# 𝐂 = A * initial_guess * B + C - initial_guess +# # 𝐂 = copy(C) - if !(C isa DenseMatrix) - C = collect(C) - end +# if !(C isa DenseMatrix) +# C = collect(C) +# end - CB = similar(C) +# CB = similar(C) - soll = speedmapping(-𝐂; - m! = (X, x) -> begin - β„’.mul!(CB, x, B) - β„’.mul!(X, A, CB) - β„’.axpy!(1, 𝐂, X) - end, stabilize = false, maps_limit = 10000, tol = tol) +# soll = speedmapping(-𝐂; +# m! = (X, x) -> begin +# β„’.mul!(CB, x, B) +# β„’.mul!(X, A, CB) +# β„’.axpy!(1, 𝐂, X) +# end, stabilize = false, maps_limit = 10000, tol = tol) - 𝐂 = soll.minimizer +# 𝐂 = soll.minimizer - 𝐂 += initial_guess +# 𝐂 += initial_guess - reached_tol = β„’.norm(A * 𝐂 * B + C - 𝐂) / max(β„’.norm(𝐂), β„’.norm(C)) +# reached_tol = β„’.norm(A * 𝐂 * B + C - 𝐂) / max(β„’.norm(𝐂), β„’.norm(C)) - # if reached_tol > tol - # println("Sylvester: speedmapping $reached_tol") - # end +# # if reached_tol > tol +# # println("Sylvester: speedmapping $reached_tol") +# # end - return 𝐂, soll.maps, reached_tol -end +# return 𝐂, soll.maps, reached_tol +# end diff --git a/src/get_functions.jl b/src/get_functions.jl index de894cea..ca1acbe2 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -718,7 +718,7 @@ function get_conditional_forecast(𝓂::β„³, Y[:,i] = pruning ? sum(initial_state) : initial_state end - elseif algorithm ∈ [:first_order, :riccati, :quadratic_iteration, :linear_time_iteration, :binder_pesaran, :first_order_doubling] + elseif algorithm ∈ [:first_order, :first_order_doubling] C = @views 𝓂.solution.perturbation.first_order.solution_matrix[:,𝓂.timings.nPast_not_future_and_mixed+1:end] CC = C[cond_var_idx,free_shock_idx] @@ -1595,7 +1595,7 @@ function get_solution(𝓂::β„³; solve!(𝓂, parameters = parameters, verbose = verbose, dynamics = true, silent = silent, algorithm = algorithm) - if algorithm ∈ [:riccati, :first_order, :first_order_doubling, :quadratic_iteration, :binder_pesaran, :linear_time_iteration] + if algorithm ∈ [:first_order, :first_order_doubling] solution_matrix = 𝓂.solution.perturbation.first_order.solution_matrix end @@ -2127,7 +2127,7 @@ function get_correlation(𝓂::β„³; lyapunov_algorithm::Symbol = :doubling, verbose::Bool = false) - @assert algorithm ∈ [:first_order,:linear_time_iteration,:quadratic_iteration, :first_order_doubling, :binder_pesaran, :pruned_second_order,:pruned_third_order] "Correlation can only be calculated for first order perturbation or second and third order pruned perturbation solutions." + @assert algorithm ∈ [:first_order, :first_order_doubling, :pruned_second_order,:pruned_third_order] "Correlation can only be calculated for first order perturbation or second and third order pruned perturbation solutions." solve!(𝓂, parameters = parameters, algorithm = algorithm, verbose = verbose) @@ -2221,7 +2221,7 @@ function get_autocorrelation(𝓂::β„³; lyapunov_algorithm::Symbol = :doubling, verbose::Bool = false) - @assert algorithm ∈ [:first_order,:linear_time_iteration,:quadratic_iteration, :first_order_doubling, :binder_pesaran, :pruned_second_order,:pruned_third_order] "Autocorrelation can only be calculated for first order perturbation or second and third order pruned perturbation solutions." + @assert algorithm ∈ [:first_order, :first_order_doubling, :pruned_second_order, :pruned_third_order] "Autocorrelation can only be calculated for first order perturbation or second and third order pruned perturbation solutions." solve!(𝓂, parameters = parameters, algorithm = algorithm, verbose = verbose) @@ -2840,7 +2840,7 @@ function get_statistics(𝓂, lyapunov_algorithm::Symbol = :doubling, verbose::Bool = false) where {U,T} - @assert algorithm ∈ [:first_order,:linear_time_iteration,:quadratic_iteration, :first_order_doubling, :binder_pesaran,:pruned_second_order,:pruned_third_order] "Statistics can only be provided for first order perturbation or second and third order pruned perturbation solutions." + @assert algorithm ∈ [:first_order, :first_order_doubling, :pruned_second_order, :pruned_third_order] "Statistics can only be provided for first order perturbation or second and third order pruned perturbation solutions." @assert !(non_stochastic_steady_state == Symbol[]) || !(standard_deviation == Symbol[]) || !(mean == Symbol[]) || !(variance == Symbol[]) || !(covariance == Symbol[]) || !(autocorrelation == Symbol[]) "Provide variables for at least one output." diff --git a/src/macros.jl b/src/macros.jl index d11c83fe..70ebe277 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -1,7 +1,7 @@ import MacroTools: unblock, postwalk, @capture, flatten -const all_available_algorithms = [:linear_time_iteration, :riccati, :first_order, :first_order_doubling, :quadratic_iteration, :binder_pesaran, :second_order, :pruned_second_order, :third_order, :pruned_third_order] +const all_available_algorithms = [:first_order, :first_order_doubling, :second_order, :pruned_second_order, :third_order, :pruned_third_order] """ diff --git a/src/moments.jl b/src/moments.jl index 0360faf0..c0b2a7b5 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -40,11 +40,11 @@ function calculate_mean(parameters::Vector{T}, # Matrix{T}, Matrix{T}, AbstractSparseMatrix{T}, AbstractSparseMatrix{T}, Bool} where T <: Real # Theoretical mean identical for 2nd and 3rd order pruned solution. - @assert algorithm ∈ [:linear_time_iteration, :riccati, :first_order, :quadratic_iteration, :binder_pesaran, :pruned_second_order, :pruned_third_order] "Theoretical mean available only for first order, pruned second and pruned third order perturbation solutions." + @assert algorithm ∈ [:first_order, :first_order_doubling, :pruned_second_order, :pruned_third_order] "Theoretical mean available only for first order, pruned second and pruned third order perturbation solutions." SS_and_pars, (solution_error, iters) = get_NSSS_and_parameters(𝓂, parameters, verbose = verbose) - if algorithm ∈ [:linear_time_iteration, :riccati, :first_order, :quadratic_iteration, :binder_pesaran] + if algorithm ∈ [:first_order, :first_order_doubling] return SS_and_pars[1:𝓂.timings.nVars], solution_error < tol end diff --git a/test/runtests.jl b/test/runtests.jl index 7dcdd013..c5d77642 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -669,13 +669,13 @@ if test_set == "basic" sol = get_solution(RBC_CME) - sol1 = get_solution(RBC_CME, algorithm = :linear_time_iteration, verbose = true) + # sol1 = get_solution(RBC_CME, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol, sol1, atol = 1e-4) + # @test isapprox(sol, sol1, atol = 1e-4) - sol2 = get_solution(RBC_CME, algorithm = :quadratic_iteration, verbose = true) + # sol2 = get_solution(RBC_CME, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol, sol2, atol = 1e-4) + # @test isapprox(sol, sol2, atol = 1e-4) sol3 = get_solution(RBC_CME, algorithm = :first_order_doubling, verbose = true) @@ -733,16 +733,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -804,16 +804,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -874,16 +874,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -945,16 +945,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1015,16 +1015,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1085,16 +1085,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1155,16 +1155,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1225,16 +1225,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1294,16 +1294,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1365,16 +1365,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1438,16 +1438,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -1511,16 +1511,16 @@ if test_set == "basic" sol_exo = get_solution(m) - sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) + # sol_exo1 = get_solution(m, algorithm = :linear_time_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) - @test isapprox(sol_exo, sol_exo1, atol = 1e-4) + # @test isapprox(sol_exo, sol_exo1, atol = 1e-4) - sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) + # sol_exo2 = get_solution(m, algorithm = :quadratic_iteration, verbose = true) - @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) + # @test isapprox(sol(axiskeys(sol, 1)[1:end-1],:), sol_exo2(axiskeys(sol, 1)[1:end-1], axiskeys(sol, 2)), atol = 1e-4) sol_exo3 = get_solution(m, algorithm = :first_order_doubling, verbose = true) @@ -2224,84 +2224,84 @@ if test_set == "basic" - @testset verbose = true "First order: linear time iteration" begin - # Numerical test with calibration targets - @model RBC_CME begin - y[0]=A[0]*k[-1]^alpha - 1/c[0]=beta*1/c[1]*(alpha*A[1]*k[0]^(alpha-1)+(1-delta)) - 1/c[0]=beta*1/c[1]*(R[0]/Pi[+1]) - R[0] * beta =(Pi[0]/Pibar)^phi_pi - A[0]*k[-1]^alpha=c[0]+k[0]-(1-delta*z_delta[0])*k[-1] - z_delta[0] = 1 - rho_z_delta + rho_z_delta * z_delta[-1] + std_z_delta * delta_eps[x] - A[0] = 1 - rhoz + rhoz * A[-1] + std_eps * eps_z[x] - end + # @testset verbose = true "First order: linear time iteration" begin + # # Numerical test with calibration targets + # @model RBC_CME begin + # y[0]=A[0]*k[-1]^alpha + # 1/c[0]=beta*1/c[1]*(alpha*A[1]*k[0]^(alpha-1)+(1-delta)) + # 1/c[0]=beta*1/c[1]*(R[0]/Pi[+1]) + # R[0] * beta =(Pi[0]/Pibar)^phi_pi + # A[0]*k[-1]^alpha=c[0]+k[0]-(1-delta*z_delta[0])*k[-1] + # z_delta[0] = 1 - rho_z_delta + rho_z_delta * z_delta[-1] + std_z_delta * delta_eps[x] + # A[0] = 1 - rhoz + rhoz * A[-1] + std_eps * eps_z[x] + # end - @parameters RBC_CME verbose = true begin - alpha | k[ss] / (4 * y[ss]) = cap_share - cap_share = 1.66 - # alpha = .157 + # @parameters RBC_CME verbose = true begin + # alpha | k[ss] / (4 * y[ss]) = cap_share + # cap_share = 1.66 + # # alpha = .157 - beta | R[ss] = R_ss # beta needs to enter into function: block in order to solve - R_ss = 1.0035 - # beta = .999 + # beta | R[ss] = R_ss # beta needs to enter into function: block in order to solve + # R_ss = 1.0035 + # # beta = .999 - # delta | c[ss]/y[ss] = 1 - I_K_ratio - delta | delta * k[ss] / y[ss] = I_K_ratio #check why this doesnt solve for y; because delta is not recognised as a free parameter here. - I_K_ratio = .15 - # delta = .0226 + # # delta | c[ss]/y[ss] = 1 - I_K_ratio + # delta | delta * k[ss] / y[ss] = I_K_ratio #check why this doesnt solve for y; because delta is not recognised as a free parameter here. + # I_K_ratio = .15 + # # delta = .0226 - Pibar | Pi[ss] = Pi_ss - Pi_ss = 1.0025 - # Pibar = 1.0008 + # Pibar | Pi[ss] = Pi_ss + # Pi_ss = 1.0025 + # # Pibar = 1.0008 - phi_pi = 1.5 - rhoz = .9 - std_eps = .0068 - rho_z_delta = .9 - std_z_delta = .005 - end + # phi_pi = 1.5 + # rhoz = .9 + # std_eps = .0068 + # rho_z_delta = .9 + # std_z_delta = .005 + # end - get_solution(RBC_CME, algorithm = :linear_time_iteration) + # get_solution(RBC_CME, algorithm = :linear_time_iteration) - @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]], [ 0.0 0.0068 - 6.73489e-6 0.000168887 - 1.01124e-5 0.000253583 - -0.000365783 0.00217203 - -0.00070019 0.00749279 - 0.0 0.00966482 - 0.005 0.0], atol = 1e-6) + # @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]], [ 0.0 0.0068 + # 6.73489e-6 0.000168887 + # 1.01124e-5 0.000253583 + # -0.000365783 0.00217203 + # -0.00070019 0.00749279 + # 0.0 0.00966482 + # 0.005 0.0], atol = 1e-6) - get_solution(RBC_CME, algorithm = :linear_time_iteration, parameters = :I_K_ratio => .1) + # get_solution(RBC_CME, algorithm = :linear_time_iteration, parameters = :I_K_ratio => .1) - @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]],[ 0.0 0.0068 - 3.42408e-6 0.000111417 - 5.14124e-6 0.000167292 - -0.000196196 0.00190741 - -0.000430554 0.0066164 - 0.0 0.00852381 - 0.005 0.0], atol = 1e-6) + # @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]],[ 0.0 0.0068 + # 3.42408e-6 0.000111417 + # 5.14124e-6 0.000167292 + # -0.000196196 0.00190741 + # -0.000430554 0.0066164 + # 0.0 0.00852381 + # 0.005 0.0], atol = 1e-6) - get_solution(RBC_CME, algorithm = :linear_time_iteration, parameters = :cap_share => 1.5) + # get_solution(RBC_CME, algorithm = :linear_time_iteration, parameters = :cap_share => 1.5) - @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]],[ 0.0 0.0068 - 4.00629e-6 0.000118171 - 6.01543e-6 0.000177434 - -0.000207089 0.00201698 - -0.00041124 0.00639229 - 0.0 0.00840927 - 0.005 0.0], atol = 1e-6) + # @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]],[ 0.0 0.0068 + # 4.00629e-6 0.000118171 + # 6.01543e-6 0.000177434 + # -0.000207089 0.00201698 + # -0.00041124 0.00639229 + # 0.0 0.00840927 + # 0.005 0.0], atol = 1e-6) - RBC_CME = nothing - end + # RBC_CME = nothing + # end - @testset verbose = true "First order: quadratic iteration" begin + @testset verbose = true "First order: doubling" begin # Numerical test with calibration targets @model RBC_CME begin y[0]=A[0]*k[-1]^alpha @@ -2339,7 +2339,7 @@ if test_set == "basic" std_z_delta = .005 end - get_solution(RBC_CME, algorithm = :quadratic_iteration) + get_solution(RBC_CME, algorithm = :first_order_doubling) @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]], [ 0.0 0.0068 6.73489e-6 0.000168887 @@ -2350,7 +2350,7 @@ if test_set == "basic" 0.005 0.0], atol = 1e-6) - get_solution(RBC_CME, algorithm = :quadratic_iteration, parameters = :I_K_ratio => .1) + get_solution(RBC_CME, algorithm = :first_order_doubling, parameters = :I_K_ratio => .1) @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]],[ 0.0 0.0068 3.42408e-6 0.000111417 @@ -2361,7 +2361,7 @@ if test_set == "basic" 0.005 0.0], atol = 1e-6) - get_solution(RBC_CME, algorithm = :quadratic_iteration, parameters = :cap_share => 1.5) + get_solution(RBC_CME, algorithm = :first_order_doubling, parameters = :cap_share => 1.5) @test isapprox(RBC_CME.solution.perturbation.first_order.solution_matrix[:,[(end-RBC_CME.timings.nExo+1):end...]],[ 0.0 0.0068 4.00629e-6 0.000118171