Merge branch 'nonlinear-estimation' of https://github.com/thorek1/Mac…

…roModelling.jl into nonlinear-estimation

src/MacroModelling.jl

@@ -3649,8 +3649,9 @@ function solve_steady_state!(𝓂::ℳ; verbose::Bool = false)
                             end
 
                             # Zero initial value if startin without guess
-                            for i in 1:2:length(closest_solution)
-                                if !isfinite(sum(abs,closest_solution[i+1]))
+                            if !isfinite(sum(abs,closest_solution[2]))
+                                closest_solution = copy(closest_solution)
+                                for i in 1:2:length(closest_solution)
                                     closest_solution[i] = zeros(length(closest_solution[i]))
                                 end
                             end

src/algorithms/quadratic_matrix_equation.jl

@@ -28,7 +28,7 @@ function solve_quadratic_matrix_equation(A::AbstractMatrix{R},
 
         reached_tol = ℒ.norm(AXX) / AXXnorm
 
-        if reached_tol < tol # 1e-12 is too large eps is too small
+        if reached_tol < (tol * length(initial_guess) / 1e6)# 1e-12 is too large eps is too small; if you use the low tol it can be that a small change in the parameters still yields an acceptable solution but as a better tol can be reached it is actually not accurate
             if verbose println("Quadratic matrix equation solver previous solution has tolerance: $reached_tol") end
 
             return initial_guess, true

src/macros.jl

@@ -1544,14 +1544,14 @@ macro parameters(𝓂,ex...)
                 # println("Find SS solver parameters which solve for the NSSS:\t",round(time() - start_time, digits = 3), " seconds")
             end
 
-            if !found_solution
-                @warn "Could not find non-stochastic steady state. Consider setting bounds on variables or calibrated parameters in the `@parameters` section (e.g. `k > 10`)."
-            end
-
             if !$silent 
                 println(round(time() - start_time, digits = 3), " seconds") 
             end
 
+            if !found_solution
+                @warn "Could not find non-stochastic steady state. Consider setting bounds on variables or calibrated parameters in the `@parameters` section (e.g. `k > 10`)."
+            end
+
             mod.$𝓂.solution.non_stochastic_steady_state = SS_and_pars
             mod.$𝓂.solution.outdated_NSSS = false
         end

src/moments.jl

@@ -130,17 +130,31 @@ function calculate_mean(parameters::Vector{T},
     # return mean_of_variables, 𝐒₁, ∇₁, 𝐒₂, ∇₂, true
 end
 
-
-
+function calculate_second_order_moments(parameters::Vector{R}, 
+                                        𝓂::ℳ;
+                                        covariance::Bool = true,
+                                        verbose::Bool = false, 
+                                        sylvester_algorithm::Symbol = :doubling,
+                                        lyapunov_algorithm::Symbol = :doubling,
+                                        tol::AbstractFloat = eps()) where R <: Real
+    calculate_second_order_moments(
+        parameters, 
+        𝓂,
+        Val(covariance);
+        verbose = verbose, 
+        sylvester_algorithm = sylvester_algorithm,
+        lyapunov_algorithm = lyapunov_algorithm,
+        tol = tol)
+end
 
 function calculate_second_order_moments(
     parameters::Vector{R}, 
-    𝓂::ℳ; 
-    covariance::Bool = true,
+    𝓂::ℳ,
+    ::Val{false}; # covariance; 
     verbose::Bool = false, 
     sylvester_algorithm::Symbol = :doubling,
     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
+    tol::AbstractFloat = eps())::Tuple{Vector{R}, Vector{R}, Matrix{R}, Matrix{R}, Vector{R}, Matrix{R}, Matrix{R}, AbstractSparseMatrix{R}, AbstractSparseMatrix{R}, Bool} where R <: Real
 
     Σʸ₁, 𝐒₁, ∇₁, SS_and_pars, solved = calculate_covariance(parameters, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm)
 
@@ -244,10 +258,122 @@ function calculate_second_order_moments(
     Δμˢ₂ = vec((ℒ.I - s_to_s₁) \ (s_s_to_s₂ * vec(Σᶻ₁) / 2 + (v_v_to_s₂ + e_e_to_s₂ * vec(ℒ.I(nᵉ))) / 2))
     μʸ₂  = SS_and_pars[1:𝓂.timings.nVars] + ŝ_to_y₂ * μˢ⁺₂ + yv₂
 
-    if !covariance
-        return μʸ₂, Δμˢ₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂
+    return μʸ₂, Δμˢ₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂, (solved && solved2)
+end
+
+
+
+function calculate_second_order_moments(
+    parameters::Vector{R}, 
+    𝓂::ℳ,
+    ::Val{true}; # covariance
+    verbose::Bool = false, 
+    sylvester_algorithm::Symbol = :doubling,
+    lyapunov_algorithm::Symbol = :doubling,
+    tol::AbstractFloat = eps())::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}, Bool} where R <: Real
+
+    Σʸ₁, 𝐒₁, ∇₁, SS_and_pars, solved = calculate_covariance(parameters, 𝓂, verbose = verbose, lyapunov_algorithm = lyapunov_algorithm)
+
+    nᵉ = 𝓂.timings.nExo
+
+    nˢ = 𝓂.timings.nPast_not_future_and_mixed
+
+    iˢ = 𝓂.timings.past_not_future_and_mixed_idx
+
+    Σᶻ₁ = Σʸ₁[iˢ, iˢ]
+
+    # precalc second order
+    ## mean
+    I_plus_s_s = sparse(reshape(ℒ.kron(vec(ℒ.I(nˢ)), ℒ.I(nˢ)), nˢ^2, nˢ^2) + ℒ.I)
+
+    ## covariance
+    E_e⁴ = zeros(nᵉ * (nᵉ + 1)÷2 * (nᵉ + 2)÷3 * (nᵉ + 3)÷4)
+
+    quadrup = multiplicate(nᵉ, 4)
+
+    comb⁴ = reduce(vcat, generateSumVectors(nᵉ, 4))
+
+    comb⁴ = comb⁴ isa Int64 ? reshape([comb⁴],1,1) : comb⁴
+
+    for j = 1:size(comb⁴,1)
+        E_e⁴[j] = product_moments(ℒ.I(nᵉ), 1:nᵉ, comb⁴[j,:])
     end
 
+    e⁴ = quadrup * E_e⁴
+
+    # second order
+    ∇₂ = calculate_hessian(parameters, SS_and_pars, 𝓂)# * 𝓂.solution.perturbation.second_order_auxilliary_matrices.𝐔∇₂
+
+    𝐒₂, solved2 = calculate_second_order_solution(∇₁, ∇₂, 𝐒₁, 
+                                                𝓂.solution.perturbation.second_order_auxilliary_matrices; 
+                                                T = 𝓂.timings, 
+                                                tol = tol, 
+                                                initial_guess = 𝓂.solution.perturbation.second_order_solution,
+                                                sylvester_algorithm = sylvester_algorithm, 
+                                                verbose = verbose)
+
+    if eltype(𝐒₂) == Float64 && solved2 𝓂.solution.perturbation.second_order_solution = 𝐒₂ end
+
+    𝐒₂ *= 𝓂.solution.perturbation.second_order_auxilliary_matrices.𝐔₂
+
+    𝐒₂ = sparse(𝐒₂)
+
+    s_in_s⁺ = BitVector(vcat(ones(Bool, nˢ), zeros(Bool, nᵉ + 1)))
+    e_in_s⁺ = BitVector(vcat(zeros(Bool, nˢ + 1), ones(Bool, nᵉ)))
+    v_in_s⁺ = BitVector(vcat(zeros(Bool, nˢ), 1, zeros(Bool, nᵉ)))
+
+    kron_s_s = ℒ.kron(s_in_s⁺, s_in_s⁺)
+    kron_e_e = ℒ.kron(e_in_s⁺, e_in_s⁺)
+    kron_v_v = ℒ.kron(v_in_s⁺, v_in_s⁺)
+    kron_s_e = ℒ.kron(s_in_s⁺, e_in_s⁺)
+
+    # first order
+    s_to_y₁ = 𝐒₁[:, 1:nˢ]
+    e_to_y₁ = 𝐒₁[:, (nˢ + 1):end]
+
+    s_to_s₁ = 𝐒₁[iˢ, 1:nˢ]
+    e_to_s₁ = 𝐒₁[iˢ, (nˢ + 1):end]
+
+
+    # second order
+    s_s_to_y₂ = 𝐒₂[:, kron_s_s]
+    e_e_to_y₂ = 𝐒₂[:, kron_e_e]
+    v_v_to_y₂ = 𝐒₂[:, kron_v_v]
+    s_e_to_y₂ = 𝐒₂[:, kron_s_e]
+
+    s_s_to_s₂ = 𝐒₂[iˢ, kron_s_s] |> collect
+    e_e_to_s₂ = 𝐒₂[iˢ, kron_e_e]
+    v_v_to_s₂ = 𝐒₂[iˢ, kron_v_v] |> collect
+    s_e_to_s₂ = 𝐒₂[iˢ, kron_s_e]
+
+    s_to_s₁_by_s_to_s₁ = ℒ.kron(s_to_s₁, s_to_s₁) |> collect
+    e_to_s₁_by_e_to_s₁ = ℒ.kron(e_to_s₁, e_to_s₁)
+    s_to_s₁_by_e_to_s₁ = ℒ.kron(s_to_s₁, e_to_s₁)
+
+    # # Set up in pruned state transition matrices
+    ŝ_to_ŝ₂ = [ s_to_s₁             zeros(nˢ, nˢ + nˢ^2)
+                zeros(nˢ, nˢ)       s_to_s₁             s_s_to_s₂ / 2
+                zeros(nˢ^2, 2*nˢ)   s_to_s₁_by_s_to_s₁                  ]
+
+    ê_to_ŝ₂ = [ e_to_s₁         zeros(nˢ, nᵉ^2 + nᵉ * nˢ)
+                zeros(nˢ,nᵉ)    e_e_to_s₂ / 2       s_e_to_s₂
+                zeros(nˢ^2,nᵉ)  e_to_s₁_by_e_to_s₁  I_plus_s_s * s_to_s₁_by_e_to_s₁]
+
+    ŝ_to_y₂ = [s_to_y₁  s_to_y₁         s_s_to_y₂ / 2]
+
+    ê_to_y₂ = [e_to_y₁  e_e_to_y₂ / 2   s_e_to_y₂]
+
+    ŝv₂ = [ zeros(nˢ) 
+            vec(v_v_to_s₂) / 2 + e_e_to_s₂ / 2 * vec(ℒ.I(nᵉ))
+            e_to_s₁_by_e_to_s₁ * vec(ℒ.I(nᵉ))]
+
+    yv₂ = (vec(v_v_to_y₂) + e_e_to_y₂ * vec(ℒ.I(nᵉ))) / 2
+
+    ## Mean
+    μˢ⁺₂ = (ℒ.I - ŝ_to_ŝ₂) \ ŝv₂
+    Δμˢ₂ = vec((ℒ.I - s_to_s₁) \ (s_s_to_s₂ * vec(Σᶻ₁) / 2 + (v_v_to_s₂ + e_e_to_s₂ * vec(ℒ.I(nᵉ))) / 2))
+    μʸ₂  = SS_and_pars[1:𝓂.timings.nVars] + ŝ_to_y₂ * μˢ⁺₂ + yv₂
+
     # Covariance
     Γ₂ = [ ℒ.I(nᵉ)             zeros(nᵉ, nᵉ^2 + nᵉ * nˢ)
             zeros(nᵉ^2, nᵉ)    reshape(e⁴, nᵉ^2, nᵉ^2) - vec(ℒ.I(nᵉ)) * vec(ℒ.I(nᵉ))'     zeros(nᵉ^2, nᵉ * nˢ)
@@ -265,7 +391,7 @@ function calculate_second_order_moments(
 
     autocorr_tmp = ŝ_to_ŝ₂ * Σᶻ₂ * ŝ_to_y₂' + ê_to_ŝ₂ * Γ₂ * ê_to_y₂'
 
-    return Σʸ₂, Σᶻ₂, μʸ₂, Δμˢ₂, autocorr_tmp, ŝ_to_ŝ₂, ŝ_to_y₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂
+    return Σʸ₂, Σᶻ₂, μʸ₂, Δμˢ₂, autocorr_tmp, ŝ_to_ŝ₂, ŝ_to_y₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂, (solved && solved2 && info)
 end
 
 
@@ -286,15 +412,16 @@ function calculate_third_order_moments(parameters::Vector{T},
                                             tol::AbstractFloat = eps()) where {U, T <: Real}
 
     second_order_moments = calculate_second_order_moments(parameters, 
-                                                            𝓂, 
+                                                            𝓂,
+                                                            Val(true);
                                                             verbose = verbose, 
                                                             sylvester_algorithm = sylvester_algorithm, 
                                                             lyapunov_algorithm = lyapunov_algorithm)
 
-    Σʸ₂, Σᶻ₂, μʸ₂, Δμˢ₂, autocorr_tmp, ŝ_to_ŝ₂, ŝ_to_y₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂ = second_order_moments
+    Σʸ₂, Σᶻ₂, μʸ₂, Δμˢ₂, autocorr_tmp, ŝ_to_ŝ₂, ŝ_to_y₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂, solved = second_order_moments
 
     if !covariance && !autocorrelation
-        return μʸ₂, Δμˢ₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂
+        return μʸ₂, Δμˢ₂, Σʸ₁, Σᶻ₁, SS_and_pars, 𝐒₁, ∇₁, 𝐒₂, ∇₂, solved
     end
 
     ∇₃ = calculate_third_order_derivatives(parameters, SS_and_pars, 𝓂)# * 𝓂.solution.perturbation.third_order_auxilliary_matrices.𝐔∇₃
@@ -354,6 +481,8 @@ function calculate_third_order_moments(parameters::Vector{T},
         autocorr = zeros(T, size(Σʸ₂,1), length(autocorrelation_periods))
     end
 
+    solved_lyapunov = true
+
     # Threads.@threads for ords in orders 
     for ords in orders 
         variance_observable, dependencies_all_vars = ords
@@ -497,6 +626,8 @@ function calculate_third_order_moments(parameters::Vector{T},
 
         Σᶻ₃, info = solve_lyapunov_equation(ŝ_to_ŝ₃, C, lyapunov_algorithm = lyapunov_algorithm, verbose = verbose)
 
+        solved_lyapunov = solved_lyapunov && info
+
         Σʸ₃tmp = ŝ_to_y₃ * Σᶻ₃ * ŝ_to_y₃' + ê_to_y₃ * Γ₃ * ê_to_y₃' + ê_to_y₃ * Eᴸᶻ * ŝ_to_y₃' + ŝ_to_y₃ * Eᴸᶻ' * ê_to_y₃'
 
         for obs in variance_observable
@@ -533,9 +664,9 @@ function calculate_third_order_moments(parameters::Vector{T},
     end
 
     if autocorrelation
-        return Σʸ₃, μʸ₂, autocorr, SS_and_pars
+        return Σʸ₃, μʸ₂, autocorr, SS_and_pars, solved && solved3 && solved_lyapunov
     else
-        return Σʸ₃, μʸ₂, SS_and_pars
+        return Σʸ₃, μʸ₂, SS_and_pars, solved && solved3 && solved_lyapunov
     end
 
 end