Chapter 31: Implementing DSGE Models

// ============================================================
// Minimal NK Model: 3-equation system
// Variables: pi (inflation), x (output gap), r (interest rate)
// Shocks: eps_u (cost-push), eps_r (demand/natural rate)
// Parameters: beta, kappa, sigma, phi_pi, phi_y, rho_u, rho_r
// ============================================================

var pi x r u_shock r_n;  // endogenous variables
varexo eps_u eps_r;       // structural shocks

parameters beta kappa sigma phi_pi phi_y rho_u rho_r;

// Calibration / prior means
beta    = 0.99;   // household discount factor
kappa   = 0.15;   // NKPC slope (from Calvo pricing)
sigma   = 1.0;    // inverse EIS
phi_pi  = 1.5;    // Taylor rule: inflation response
phi_y   = 0.5;    // Taylor rule: output gap response
rho_u   = 0.5;    // cost-push shock persistence
rho_r   = 0.8;    // demand shock persistence

model(linear);    // declare this is a linearized model

// [1] New Keynesian Phillips Curve
pi = beta * pi(+1) + kappa * x + u_shock;

// [2] Dynamic IS equation
x = x(+1) - sigma * (r - pi(+1) - r_n);

// [3] Taylor rule (with interest rate smoothing = 0 here)
r = phi_pi * pi + phi_y * x;

// [4] Cost-push shock process
u_shock = rho_u * u_shock(-1) + eps_u;

// [5] Natural rate process (demand shock)
r_n = rho_r * r_n(-1) + eps_r;

end;

// Steady state: all log-deviations = 0
initval;
  pi = 0; x = 0; r = 0; u_shock = 0; r_n = 0;
end;

// Shock variances
shocks;
  var eps_u = 0.01^2;  // std dev 1% for cost-push
  var eps_r = 0.01^2;  // std dev 1% for demand
end;

// Solve and simulate
stoch_simul(order=1, irf=20, periods=200);
// This calls gensys internally, checks Blanchard-Kahn,
// computes IRFs, and simulates 200 periods.

// For estimation, replace stoch_simul with:
/*
estimated_params;
  kappa,    beta_pdf,    0.15, 0.10;     // Beta(mean=0.15, std=0.10)
  phi_pi,   normal_pdf,  1.50, 0.25;     // N(1.5, 0.25) truncated [1,inf)
  rho_u,    beta_pdf,    0.50, 0.20;     // Beta(mean=0.5, std=0.2)
  stderr eps_u, inv_gamma_pdf, 0.01, 2;  // IG(0.01, 2)
  stderr eps_r, inv_gamma_pdf, 0.01, 2;
end;

varobs pi x r;   // observed variables

estimation(datafile='us_data', mh_replic=100000, mh_nblocks=2,
           mh_jscale=0.3, bayesian_irf);
*/

import numpy as np
from scipy.linalg import ordqz
from scipy.stats import norm, beta as beta_dist
import matplotlib.pyplot as plt

class MinimalNKModel:
    """Three-equation New Keynesian model — complete implementation."""
    
    def __init__(self, beta=0.99, kappa=0.15, sigma=1.0,
                 phi_pi=1.5, phi_y=0.5, rho_u=0.5, rho_r=0.8,
                 sig_u=0.01, sig_r=0.01):
        self.params = dict(beta=beta, kappa=kappa, sigma=sigma,
                           phi_pi=phi_pi, phi_y=phi_y,
                           rho_u=rho_u, rho_r=rho_r,
                           sig_u=sig_u, sig_r=sig_r)
        self._build_matrices()
    
    def _build_matrices(self):
        p = self.params
        self.G0 = np.array([[1.0, -p['kappa']],
                             [p['sigma']*p['phi_pi'], 1+p['sigma']*p['phi_y']]])
        self.G1 = np.array([[p['beta'], 0.0],
                             [-p['sigma'], 1.0]])
        self.Psi = np.array([[1.0, 0.0],
                              [0.0, p['sigma']]])
        self.Phi = np.diag([p['rho_u'], p['rho_r']])
        self.Q   = np.diag([p['sig_u']**2, p['sig_r']**2])
    
    def solve(self):
        """Compute MSV solution via Sylvester equation."""
        from numpy import kron
        A = np.linalg.inv(self.G0) @ self.G1
        eigs = np.linalg.eigvals(A)
        self.determinate = np.all(np.abs(eigs) > 1.0)
        if not self.determinate:
            return False
        C_mat = np.linalg.inv(self.G0) @ self.Psi
        KPA = kron(self.Phi.T, A)
        I4 = np.eye(4)
        vec_Omega = np.linalg.solve(I4 - KPA, C_mat.flatten(order='F'))
        self.Omega = vec_Omega.reshape(2, 2, order='F')
        return True
    
    def irfs(self, H=20):
        """Impulse response functions for both shocks."""
        if not hasattr(self, 'Omega'):
            self.solve()
        irf = np.zeros((H, 2, 2))  # [horizon, variable, shock]
        for j in range(2):
            shock = np.zeros(2); shock[j] = 1.0
            state = shock.copy()
            for h in range(H):
                irf[h, :, j] = self.Omega @ state
                state = self.Phi @ state
        return irf
    
    def theoretical_moments(self):
        """Theoretical variances via Lyapunov equation."""
        from numpy import kron
        from scipy.linalg import solve_discrete_lyapunov
        Sigma_z = solve_discrete_lyapunov(self.Phi, self.Q)
        Sigma_y = self.Omega @ Sigma_z @ self.Omega.T
        return {'var_pi': Sigma_y[0,0], 'var_x': Sigma_y[1,1],
                'cov_pi_x': Sigma_y[0,1], 'Sigma_y': Sigma_y}
    
    def log_likelihood(self, Y_obs):
        """Evaluate likelihood via Kalman filter (simplified)."""
        if not self.solve():
            return -np.inf
        moments = self.theoretical_moments()
        # Simplified: use theoretical variance for Gaussian lik
        # Full version would run Kalman filter on observed series
        var_y = moments['var_pi']
        T = len(Y_obs)
        return -T/2 * np.log(2*np.pi*var_y) - np.sum(Y_obs**2) / (2*var_y)
    
    def plot_irfs(self, H=20):
        irf = self.irfs(H)
        fig, axes = plt.subplots(2, 2, figsize=(11, 7))
        names = ['Inflation (π̂)', 'Output gap (x̂)']
        shocks = ['Cost-push shock (εᵤ)', 'Demand shock (εᵣ)']
        for j, shock_name in enumerate(shocks):
            for i, var_name in enumerate(names):
                ax = axes[i, j]
                ax.bar(range(H), irf[:, i, j]*100, color='steelblue', alpha=0.8)
                ax.axhline(0, color='k', lw=0.5)
                ax.set_title(f'{var_name}\nResponse to {shock_name}')
                ax.set_xlabel('Quarters')
                ax.set_ylabel('% deviation')
        plt.tight_layout()
        return fig

# Run the model
model = MinimalNKModel()
model.solve()
moments = model.theoretical_moments()
print(f"Determinacy: {model.determinate}")
print(f"Var(π̂): {moments['var_pi']*100:.4f}%², Var(x̂): {moments['var_x']*100:.4f}%²")
print(f"Welfare loss: {(moments['var_pi'] + 0.1*moments['var_x'])*1e4:.4f} × 10⁻⁴")

irf = model.irfs(H=20)
print(f"\nCost-push IRF: π̂ at h=0: {irf[0,0,0]*100:.2f}%, x̂ at h=0: {irf[0,1,0]*100:.2f}%")
print(f"Demand IRF:    π̂ at h=0: {irf[0,0,1]*100:.2f}%, x̂ at h=0: {irf[0,1,1]*100:.2f}%")

# Sensitivity: vary phi_pi and compute welfare loss
phi_pi_vals = np.linspace(1.01, 4.0, 50)
welfare = []
for phi in phi_pi_vals:
    m = MinimalNKModel(phi_pi=phi)
    if m.solve():
        mom = m.theoretical_moments()
        welfare.append(mom['var_pi'] + 0.1*mom['var_x'])
    else:
        welfare.append(np.nan)

plt.figure(figsize=(8,4))
plt.plot(phi_pi_vals, [w*1e4 if w is not np.nan else np.nan for w in welfare])
plt.xlabel('φ_π'); plt.ylabel('Welfare loss (×10⁻⁴)')
plt.title('NK Welfare Loss vs Taylor Rule Coefficient')
plt.axvline(1.0, color='r', ls='--', label='Taylor principle boundary')
plt.legend(); plt.tight_layout(); plt.show()

using LinearAlgebra, Statistics

struct NKModel
    beta::Float64; kappa::Float64; sigma::Float64
    phi_pi::Float64; phi_y::Float64
    rho_u::Float64; rho_r::Float64
    sig_u::Float64; sig_r::Float64
end

NKModel(; beta=0.99, kappa=0.15, sigma=1.0,
          phi_pi=1.5, phi_y=0.5, rho_u=0.5, rho_r=0.8,
          sig_u=0.01, sig_r=0.01) =
    NKModel(beta, kappa, sigma, phi_pi, phi_y, rho_u, rho_r, sig_u, sig_r)

function build_matrices(m::NKModel)
    G0  = [1.0 -m.kappa; m.sigma*m.phi_pi 1+m.sigma*m.phi_y]
    G1  = [m.beta 0.0; -m.sigma 1.0]
    Psi = [1.0 0.0; 0.0 m.sigma]
    Phi = diagm([m.rho_u, m.rho_r])
    Q   = diagm([m.sig_u^2, m.sig_r^2])
    return G0, G1, Psi, Phi, Q
end

function solve_msv(m::NKModel)
    G0, G1, Psi, Phi, Q = build_matrices(m)
    A = inv(G0) * G1
    eigs = eigvals(A)
    all(abs.(eigs) .> 1.0) || return nothing, false
    
    C_mat = inv(G0) * Psi
    KPA   = kron(Phi', A)
    I4    = I(4)
    vec_O = (I4 - KPA) \ vec(C_mat)
    Omega = reshape(vec_O, 2, 2)
    return Omega, true
end

function compute_irfs(m::NKModel; H=20)
    Omega, ok = solve_msv(m)
    ok || return nothing
    _, _, _, Phi, _ = build_matrices(m)
    irf = zeros(H, 2, 2)
    for j in 1:2
        shock = zeros(2); shock[j] = 1.0
        state = copy(shock)
        for h in 1:H
            irf[h,:,j] = Omega * state
            state = Phi * state
        end
    end
    return irf
end

function welfare_loss(m::NKModel; lambda_x=0.1)
    Omega, ok = solve_msv(m)
    ok || return Inf
    _, _, _, Phi, Q = build_matrices(m)
    # Lyapunov: Sigma_z = Phi * Sigma_z * Phi' + Q
    from_vec = kron(Phi, Phi)
    vec_Sz = (I(4) - from_vec) \ vec(Q)
    Sigma_z = reshape(vec_Sz, 2, 2)
    Sigma_y = Omega * Sigma_z * Omega'
    return Sigma_y[1,1] + lambda_x * Sigma_y[2,2]
end

# Run
m = NKModel()
Omega, ok = solve_msv(m)
println("Determinate: $ok")
println("Omega (MSV solution):\n", round.(Omega, digits=4))

irf = compute_irfs(m, H=20)
println("\nCost-push shock IRF (h=1..5):")
println("  π̂: ", round.(irf[1:5,1,1].*100, digits=3))
println("  x̂: ", round.(irf[1:5,2,1].*100, digits=3))

wl = welfare_loss(m)
println("\nWelfare loss (φ_π=1.5): $(round(wl*1e4, digits=4)) × 10⁻⁴")

# Optimal phi_pi
opt_phi = 1.5
min_wl = wl
for phi in 1.01:0.05:5.0
    m_test = NKModel(phi_pi=phi)
    wl_test = welfare_loss(m_test)
    if wl_test < min_wl; min_wl = wl_test; opt_phi = phi; end
end
println("Optimal φ_π ≈ $opt_phi, min welfare loss = $(round(min_wl*1e4, digits=4)) × 10⁻⁴")

# R — NK model: full pipeline
library(Matrix)

nk_model <- function(beta=0.99, kappa=0.15, sigma=1.0,
                     phi_pi=1.5, phi_y=0.5, rho_u=0.5, rho_r=0.8,
                     sig_u=0.01, sig_r=0.01) {
  G0 <- matrix(c(1, sigma*phi_pi, -kappa, 1+sigma*phi_y), 2, 2)
  G1 <- matrix(c(beta, -sigma, 0, 1), 2, 2)
  Psi<- matrix(c(1,0,0,sigma),2,2)
  Phi<- diag(c(rho_u, rho_r))
  Q  <- diag(c(sig_u^2, sig_r^2))
  
  A    <- solve(G0) %*% G1
  eigs <- eigen(A)$values
  if(!all(Mod(eigs) > 1)) return(list(ok=FALSE))
  
  # MSV solution
  C_mat <- solve(G0) %*% Psi
  KPA   <- kronecker(t(Phi), A)
  I4    <- diag(4)
  vec_O <- solve(I4 - KPA, as.vector(C_mat))
  Omega <- matrix(vec_O, 2, 2)
  
  # IRFs
  H <- 20; irf <- array(0, c(H, 2, 2))
  for(j in 1:2) {
    shock <- c(0,0); shock[j] <- 1; state <- shock
    for(h in 1:H) {
      irf[h,,j] <- Omega %*% state; state <- Phi %*% state
    }
  }
  
  # Lyapunov for variances
  KPP <- kronecker(Phi, Phi)
  vec_Sz <- solve(diag(4)-KPP, as.vector(Q))
  Sigma_z <- matrix(vec_Sz,2,2)
  Sigma_y <- Omega %*% Sigma_z %*% t(Omega)
  
  list(ok=TRUE, Omega=Omega, irf=irf, Sigma_y=Sigma_y,
       welfare=Sigma_y[1,1]+0.1*Sigma_y[2,2])
}

m <- nk_model()
cat(sprintf("Determinate: %s\n", m$ok))
cat(sprintf("Var(π̂)=%.4f%%, Var(x̂)=%.4f%%\n",
            m$Sigma_y[1,1]*100, m$Sigma_y[2,2]*100))
cat(sprintf("Welfare loss: %.4f × 10⁻⁴\n", m$welfare*1e4))
cat("\nCost-push IRF (π̂, h=1..5):",
    round(m$irf[1:5,1,1]*100, 3), "\n")

{
  "model": "minimal_NK",
  "version": "1.0",
  "calibrated": {"beta": 0.99, "delta": 0.025},
  "estimated_posterior_means": {"kappa": 0.148, "phi_pi": 1.523, "rho_u": 0.487},
  "estimated_posterior_stds": {"kappa": 0.032, "phi_pi": 0.121, "rho_u": 0.093}
}