diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl index 81f155d5..89197a75 100644 --- a/src/MacroModelling.jl +++ b/src/MacroModelling.jl @@ -4579,10 +4579,8 @@ function solve!(𝓂::β„³; end # timeit_debug - if solution_error > tol - @warn "Could not find non stochastic steady steady." - end - + @assert solution_error < tol "Could not find non stochastic steady steady." + @timeit_debug timer "Calculate Jacobian" begin βˆ‡β‚ = calculate_jacobian(𝓂.parameter_values, SS_and_pars, 𝓂)# |> Matrix diff --git a/src/get_functions.jl b/src/get_functions.jl index 805a60e1..b77bac5c 100644 --- a/src/get_functions.jl +++ b/src/get_functions.jl @@ -2135,7 +2135,9 @@ function get_correlation(𝓂::β„³; elseif algorithm == :pruned_second_order covar_dcmp, Ξ£αΆ»β‚‚, state_ΞΌ, Δμ˒₂, autocorr_tmp, sΜ‚_to_sΜ‚β‚‚, sΜ‚_to_yβ‚‚, Σʸ₁, Σᢻ₁, SS_and_pars, 𝐒₁, βˆ‡β‚, 𝐒₂, βˆ‡β‚‚ = calculate_second_order_moments(𝓂.parameter_values, 𝓂, sylvester_algorithm = sylvester_algorithm, lyapunov_algorithm = lyapunov_algorithm, verbose = verbose) else - covar_dcmp, sol, _, SS_and_pars = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, sol, _, SS_and_pars, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." end std = sqrt.(β„’.diag(covar_dcmp)) @@ -2236,7 +2238,9 @@ function get_autocorrelation(𝓂::β„³; sΜ‚_to_ŝ₂ⁱ *= sΜ‚_to_sΜ‚β‚‚ end else - covar_dcmp, sol, _, SS_and_pars = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, sol, _, SS_and_pars, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." A = @views sol[:,1:𝓂.timings.nPast_not_future_and_mixed] * β„’.diagm(ones(𝓂.timings.nVars))[𝓂.timings.past_not_future_and_mixed_idx,:] @@ -2380,6 +2384,8 @@ function get_moments(𝓂::β„³; NSSS, (solution_error, iters) = 𝓂.solution.outdated_NSSS ? get_NSSS_and_parameters(𝓂, 𝓂.parameter_values, verbose = verbose) : (copy(𝓂.solution.non_stochastic_steady_state), (eps(), 0)) + @assert solution_error < 1e-12 "Could not find non-stochastic steady state." + if length_par * length(NSSS) > 200 && derivatives @info "Most of the time is spent calculating derivatives wrt parameters. If they are not needed, add `derivatives = false` as an argument to the function call." maxlog = 3 end @@ -2471,7 +2477,9 @@ function get_moments(𝓂::β„³; var_means = KeyedArray(state_ΞΌ[var_idx]; Variables = axis1) end else - covar_dcmp, ___, __, _ = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, ___, __, _, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." # dvariance = π’œ.jacobian(𝒷(), x -> covariance_parameter_derivatives(x, param_idx, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm), 𝓂.parameter_values[param_idx])[1] dvariance = π’Ÿ.jacobian(x -> max.(β„’.diag(calculate_covariance(x, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm)[1]),eps(Float64)), backend, 𝓂.parameter_values)[:,param_idx] @@ -2536,8 +2544,10 @@ function get_moments(𝓂::β„³; var_means = KeyedArray(state_ΞΌ[var_idx]; Variables = axis1) end else - covar_dcmp, ___, __, _ = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, ___, __, _, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + @assert solved "Could not find covariance matrix." + # dst_dev = π’œ.jacobian(𝒷(), x -> sqrt.(covariance_parameter_derivatives(x, param_idx, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm)), 𝓂.parameter_values[param_idx])[1] dst_dev = π’Ÿ.jacobian(x -> sqrt.(max.(β„’.diag(calculate_covariance(x, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm)[1]),eps(Float64))), backend, 𝓂.parameter_values)[:,param_idx] end @@ -2612,7 +2622,9 @@ function get_moments(𝓂::β„³; var_means = KeyedArray(state_ΞΌ[var_idx]; Variables = axis1) end else - covar_dcmp, ___, __, _ = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, ___, __, _, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." end varr = convert(Vector{Real},max.(β„’.diag(covar_dcmp),eps(Float64))) @@ -2636,7 +2648,9 @@ function get_moments(𝓂::β„³; var_means = KeyedArray(state_ΞΌ[var_idx]; Variables = axis1) end else - covar_dcmp, ___, __, _ = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, ___, __, _, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." end st_dev = KeyedArray(sqrt.(convert(Vector{Real},max.(β„’.diag(covar_dcmp),eps(Float64))))[var_idx]; Variables = axis1) end @@ -2653,7 +2667,9 @@ function get_moments(𝓂::β„³; var_means = KeyedArray(state_ΞΌ[var_idx]; Variables = axis1) end else - covar_dcmp, ___, __, _ = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, ___, __, _, solved = calculate_covariance(𝓂.parameter_values, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." end end end @@ -2874,7 +2890,9 @@ function get_statistics(𝓂, end else - covar_dcmp, sol, _, SS_and_pars = calculate_covariance(all_parameters, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + covar_dcmp, sol, _, SS_and_pars, solved = calculate_covariance(all_parameters, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + + @assert solved "Could not find covariance matrix." end SS = SS_and_pars[1:end - length(𝓂.calibration_equations)] diff --git a/src/moments.jl b/src/moments.jl index 32579da7..489c32c5 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -2,9 +2,13 @@ function calculate_covariance(parameters::Vector{R}, 𝓂::β„³; lyapunov_algorithm::Symbol = :doubling, - verbose::Bool = false)::Tuple{Matrix{R}, Matrix{R}, Matrix{R}, Vector{R}} where R <: Real + verbose::Bool = false)::Tuple{Matrix{R}, Matrix{R}, Matrix{R}, Vector{R}, Bool} where R <: Real SS_and_pars, (solution_error, iters) = get_NSSS_and_parameters(𝓂, parameters, verbose = verbose) + if solution_error > 1e-12 + return zeros(0,0), zeros(0,0), zeros(0,0), SS_and_pars, solution_error < 1e-12 + end + βˆ‡β‚ = calculate_jacobian(parameters, SS_and_pars, 𝓂) sol, qme_sol, solved = calculate_first_order_solution(βˆ‡β‚; T = 𝓂.timings, initial_guess = 𝓂.solution.perturbation.qme_solution, verbose = verbose) @@ -17,9 +21,13 @@ function calculate_covariance(parameters::Vector{R}, CC = C * C' - covar_raw, _ = solve_lyapunov_equation(A, CC, lyapunov_algorithm = lyapunov_algorithm, verbose = verbose) + if !solved + return CC, sol, βˆ‡β‚, SS_and_pars, solved + end + + covar_raw, solved = solve_lyapunov_equation(A, CC, lyapunov_algorithm = lyapunov_algorithm, verbose = verbose) - return covar_raw, sol , βˆ‡β‚, SS_and_pars + return covar_raw, sol , βˆ‡β‚, SS_and_pars, solved end @@ -134,7 +142,7 @@ function calculate_second_order_moments( lyapunov_algorithm::Symbol = :doubling, tol::AbstractFloat = eps())::Union{Tuple{Matrix{R}, Matrix{R}, Vector{R}, Vector{R}, Matrix{R}, Matrix{R}, Matrix{R}, Matrix{R}, Matrix{R}, Vector{R}, Matrix{R}, Matrix{R}, AbstractSparseMatrix{R}, AbstractSparseMatrix{R}}, Tuple{Vector{R}, Vector{R}, Matrix{R}, Matrix{R}, Vector{R}, Matrix{R}, Matrix{R}, AbstractSparseMatrix{R}, AbstractSparseMatrix{R}}} where R <: Real - Σʸ₁, 𝐒₁, βˆ‡β‚, SS_and_pars = calculate_covariance(parameters, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) + Σʸ₁, 𝐒₁, βˆ‡β‚, SS_and_pars, solved = calculate_covariance(parameters, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm) nᡉ = 𝓂.timings.nExo