Fix bug in NK with Capital Taylor rule and real excess returns calcul…

…ation.

examples/nk_with_capital/example_nk_with_capital.jl

@@ -1,20 +1,25 @@
-# This script actually solves the WachterDisasterRisk model with a risk-adjusted linearization
-# and times the methods, if desired
+# This script solves the NKCapital model with a risk-adjusted linearization.
+# Note that, under default parameters, forward difference equations can be
+# approximated up to four periods ahead. A five-period ahead
+# approximation does not yield a saddle-path stable
+# deterministic or stochastic steady state.
 using RiskAdjustedLinearizations, JLD2, LinearAlgebra, Test
 
-# Number of quadrature points
-n_GH = 5
-include("nk_with_capital.jl")
-
 # Settings
+define_functions      = true
 testing               = true         # check model's solution under default parameters against saved output
 autodiff              = false
 algorithm             = :relaxation
 euler_equation_errors = false
 test_price_dispersion = false        # check if price dispersion in steady state is always bounded below by 1
 plot_irfs             = false
 horizon               = 40    # horizon for IRFs
-N_approx              = 5     # Number of periods ahead used for forward-difference equations
+N_approx              = 4     # Number of periods ahead used for forward-difference equations
+n_GH                  = 5     # number of nodes for Gauss-Hermite quadrature
+
+if define_functions
+    include("nk_with_capital.jl")
+end
 
 # Set up
 m_nk = NKCapital(; N_approx = N_approx) # create parameters
@@ -31,6 +36,9 @@ if testing
     test_m = nk_capital(test_m_nk) # instantiate risk-adjusted linearization
 
     solve!(test_m; algorithm = :deterministic, verbose = :none)
+    zdet = copy(test_m.z)
+    ydet = copy(test_m.y)
+    Psidet = copy(test_m.Ψ)
     @test test_m.z ≈ out["z_det"]
     @test test_m.y ≈ out["y_det"]
     @test test_m.Ψ ≈ out["Psi_det"]
@@ -40,12 +48,12 @@ if testing
     @test test_m.y ≈ out["y"]
     @test test_m.Ψ ≈ out["Psi"]
 
-    test_5_m_nk = NKCapital(; N_approx = 5) # create parameters
-    test_5_m = nk_capital(test_5_m_nk) # instantiate risk-adjusted linearization
-    solve!(test_5_m; algorithm = :relaxation, verbose = :none)
-    @test test_5_m.z ≈ out["z_5"]
-    @test test_5_m.y ≈ out["y_5"]
-    @test test_5_m.Ψ ≈ out["Psi_5"]
+    test_4_m_nk = NKCapital(; N_approx = 4) # create parameters
+    test_4_m = nk_capital(test_4_m_nk) # instantiate risk-adjusted linearization
+    solve!(test_4_m; algorithm = :relaxation, verbose = :none)
+    @test test_4_m.z ≈ out["z_4"]
+    @test test_4_m.y ≈ out["y_4"]
+    @test test_4_m.Ψ ≈ out["Psi_4"]
 end
 
 if test_price_dispersion
@@ -177,9 +185,11 @@ if plot_irfs
         EₜRₖₜ₊₁ = exp.(y_irfs[k][m_nk.J[:rk], 2:end] .+ m.y[m_nk.J[:rk]])
         EₜQₜ₊₁ = exp.(y_irfs[k][m_nk.J[:q], 2:end] .+ m.y[m_nk.J[:q]])
         EₜΩₜ₊₁ = exp.(y_irfs[k][m_nk.J[:ω], 2:end] .+ m.y[m_nk.J[:ω]])
-        exc_ret = (EₜRₖₜ₊₁ + EₜQₜ₊₁ .* EₜΩₜ₊₁) ./ exp.(y_irfs[k][m_nk.J[:q], 1:end - 1] .+ m.y[m_nk.J[:q]]) -
-            exp.(y_irfs[k][m_nk.J[:r], 1:end - 1] - y_irfs[k][m_nk.J[:π], 1:end - 1] .+ (m.y[m_nk.J[:r]] - m.y[m_nk.J[:π]])) .-
-            (exp.(m.y[m_nk.J[:rk]]) + exp.(m.y[m_nk.J[:q]] + m.y[m_nk.J[:ω]])) / exp.(m.y[m_nk.J[:q]])
+        exc_ret = ((EₜRₖₜ₊₁ + EₜQₜ₊₁ .* EₜΩₜ₊₁) ./ exp.(y_irfs[k][m_nk.J[:q], 1:end - 1] .+ m.y[m_nk.J[:q]])) -
+            (exp.(y_irfs[k][m_nk.J[:r], 1:end - 1] - y_irfs[k][m_nk.J[:π], 1:end - 1] .+
+                  (m.y[m_nk.J[:r]] - m.y[m_nk.J[:π]]))) .-
+            (((exp(m.y[m_nk.J[:rk]]) + exp.(m.y[m_nk.J[:q]] + m.y[m_nk.J[:ω]])) / exp.(m.y[m_nk.J[:q]])) -
+             exp(m.y[m_nk.J[:r]] - m.y[m_nk.J[:π]]))
         plot_dicts[k][:real_excess_ret] = plot(1:(horizon - 1), exc_ret, label = "Real Excess Returns",
                                                linewidth = 3, color = :black)
     end

examples/nk_with_capital/nk_with_capital.jl

@@ -45,8 +45,8 @@ function NKCapital(; β::T = .99, γ::T = 3.8, φ::T = 1., ν::T = 1., χ::T = 4
     # state variables, instead of e.g. η_L and η_A, I use η_l and η_a
     # since the uppercase variable will not appear in the jumps/states.
     S_init  = [:k₋₁, :v₋₁, :r₋₁, :output₋₁, :η_β, :η_l, :η_a, :η_r] # State Variables
-    J_init  = [:output, :c, :l, :w, :r, :π, :q, :x, :rk, :ω, :mc,
-               :s₁, :s₂, :v] # Jump variables
+    J_init  = [:output, :c, :l, :w, :r, :π, :q, :x, :rk, :mc,
+               :s₁, :s₂, :v, :ω] # Jump variables
     E_init  = [:wage, :euler, :tobin, :cap_ret,
                :eq_mc, :kl_ratio, :eq_s₁, :eq_s₂,
                :phillips_curve, :price_dispersion,
@@ -139,13 +139,13 @@ function nk_capital(m::NKCapital{T}) where {T <: Real}
         F[kl_ratio]            = z[k₋₁] - y[l] - log(α / (1. - α)) - (y[w] - y[rk])
         F[phillips_curve]      = (1. - ϵ) * y[π] - log((1. - θ) * exp((1. - ϵ) * (pstarv + y[π])) + θ)
         F[price_dispersion]    = y[v] - ϵ * y[π] - log((1. - θ) * exp(-ϵ * (pstarv + y[π])) + θ * exp(z[v₋₁]))
-        F[mp]                  = ϕ_r * z[r₋₁] + (1. - ϕ_r) .* (y[r] + ϕ_π * (y[π] - π_ss) +
-                                                               ϕ_y * (y[output] - z[output₋₁])) + z[η_r] - y[r]
+        F[mp]                  = (1. - ϕ_r) * r_ss + ϕ_r * z[r₋₁] + (1. - ϕ_r) .*
+            (ϕ_π * (y[π] - π_ss) + ϕ_y * (y[output] - z[output₋₁])) + z[η_r] - y[r]
         F[output_market_clear] = y[output] - log(exp(y[c]) + exp(y[x]))
         F[production]          = z[η_a] + α * z[k₋₁] + (1. - α) * y[l] - y[v] - y[output]
 
         ## Forward-difference equations separately handled b/c recursions
-        F[eq_omega] = 1. - δ + Φv - Φ′v * exp(y[x] - z[k₋₁]) - exp(y[ω])
+        F[eq_omega] = log(1. - δ + Φv - Φ′v * exp(y[x] - z[k₋₁])) - y[ω]
         F[cap_ret]  = y[q] - log(sum([exp(y[J[Symbol("dq$(i)")]]) for i in 1:N_approx]) +
                                 exp(y[J[Symbol("pq$(N_approx)")]]))
         F[eq_s₁]    = y[s₁] - log(sum([exp(y[J[Symbol("ds₁$(i)")]]) for i in 0:(N_approx - 1)]) +
@@ -258,7 +258,7 @@ function nk_capital(m::NKCapital{T}) where {T <: Real}
     RK0 = 1. / β + X̄ - 1.
 
     # Guesses
-    L0 = .5548
+    L0 = .5548429
     V0 = 1. # true if π_ss = 0, otherwise this is only a reasonable guess
 
     # Implied values given guesses