NKModel_CLMM_julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving a New Keynesian model with Julia\n",
    "\n",
    "This notebook is part of a computational appendix that accompanies the paper.\n",
    "\n",
    "> MATLAB, Python, Julia: What to Choose in Economics? \n",
    ">> Coleman, Lyon, Maliar, and Maliar (2017)\n",
    "\n",
    "In order to run the codes in this notebook you will need to install and configure a few Julia packages. We recommend downloading [JuliaPro](https://juliacomputing.com/products/juliapro.html) and/or following the instructions on [quantecon.org](https://lectures.quantecon.org/py/getting_started.html). Once your Julia installation is up and running, there are a few additional packages you will need in order to run the code here. To do this uncomment the lines in the cell below (by deleting the `#` and space at the beginning of each line) and run the cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Pkg.add(\"Sobol\")\n",
    "# Pkg.add(\"QuantEcon\")\n",
    "# Pkg.add(\"BasisMatrices\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "## Julia Code\n",
    "\n",
    "The Julia version of our algorithm is implemented as a few functions defined on\n",
    "a core type named `Model`. This type is itself composed of three different\n",
    "types that hold the model parameters, steady state, and grids needed to\n",
    "describe the numerical model. Before we get to the types, we need to bring in\n",
    "some dependencies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [],
   "source": [
    "# for constructing Sobol sequences\n",
    "using Sobol\n",
    "# for basis matrix of complete monomials and monomial quadrature rules\n",
    "using BasisMatrices: complete_polynomial, complete_polynomial!, n_complete\n",
    "using QuantEcon: qnwmonomial1, qnwmonomial2\n",
    "    \n",
    "# Set seed so we can replicate results\n",
    "srand(42);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "## Types\n",
    "\n",
    "First we have the `Params` type, which holds all the model parameters as well\n",
    "as the paramters that drive the algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vcov (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "immutable Params\n",
    "    zlb::Bool\n",
    "    γ::Float64      # Utility-function parameter\n",
    "    β::Float64      # Discount factor\n",
    "    ϑ::Float64      # Utility-function parameter\n",
    "    ϵ::Float64      # Parameter in the Dixit-Stiglitz aggregator\n",
    "    ϕ_y::Float64    # Parameter of the Taylor rule\n",
    "    ϕ_π::Float64    # Parameter of the Taylor rule\n",
    "    μ::Float64      # Parameter of the Taylor rule\n",
    "    Θ::Float64      # Share of non-reoptimizing firms (Calvo's pricing)\n",
    "    πstar::Float64  # Target (gross) inflation rate\n",
    "    gbar::Float64   # Steady-state share of government spending in output\n",
    "\n",
    "    # autocorrelation coefficients\n",
    "    ρηR::Float64  # See process (28) in MM (2015)\n",
    "    ρηa::Float64  # See process (22) in MM (2015)\n",
    "    ρηL::Float64  # See process (16) in MM (2015)\n",
    "    ρηu::Float64  # See process (15) in MM (2015)\n",
    "    ρηB::Float64  # See process (17) in MM (2015)\n",
    "    ρηG::Float64  # See process (26) in MM (2015)\n",
    "\n",
    "    # standard deviations\n",
    "    σηR::Float64  # See process (28) in MM (2015)\n",
    "    σηa::Float64  # See process (22) in MM (2015)\n",
    "    σηL::Float64  # See process (16) in MM (2015)\n",
    "    σηu::Float64  # See process (15) in MM (2015)\n",
    "    σηB::Float64  # See process (17) in MM (2015)\n",
    "    σηG::Float64  # See process (26) in MM (2015)\n",
    "\n",
    "    # algorithm parameters\n",
    "    deg::Int       # max degree of complete monomial\n",
    "    damp::Float64  # dampening parameter for coefficient update\n",
    "    tol::Float64   # Tolerance for stopping iterations\n",
    "    m::Int         # number of points in the grid\n",
    "    grid_kind::Symbol  # type of grid (:sobol or :random)\n",
    "end\n",
    "\n",
    "# constructor that takes only keyword arguments and specifies\n",
    "# defualt values\n",
    "function Params(;zlb::Bool=true, γ= 1, β=0.99, ϑ=2.09, ϵ=4.45, ϕ_y=0.07,\n",
    "                 ϕ_π=2.21, μ=0.82, Θ=0.83, πstar=1, gbar=0.23,\n",
    "                 ρηR=0.0, ρηa=0.95, ρηL=0.25, ρηu=0.92, ρηB=0.0, ρηG=0.95,\n",
    "                 σηR=0.0028, σηa=0.0045, σηL=0.0500,\n",
    "                 σηu=0.0054, σηB=0.0010, σηG=0.0038,\n",
    "                 deg=2, damp=0.1, tol=1e-7, m=200, grid_kind=:sobol)\n",
    "    Params(zlb, γ, β, ϑ, ϵ, ϕ_y, ϕ_π, μ, Θ, πstar, gbar,\n",
    "           ρηR, ρηa, ρηL, ρηu, ρηB, ρηG, σηR, σηa, σηL, σηu, σηB, σηG,\n",
    "           deg, damp, tol, m, grid_kind)\n",
    "end\n",
    "\n",
    "# returns the covariance matrix for the 6 shocks in the model\n",
    "vcov(p::Params) = diagm([p.σηR^2, p.σηa^2, p.σηL^2, p.σηu^2, p.σηB^2, p.σηG^2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "Next we have a type called `SteadyState` that is intended to hold the\n",
    "deterministic steady state realization for each variable in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SteadyState"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "immutable SteadyState\n",
    "    Yn::Float64\n",
    "    Y::Float64\n",
    "    π::Float64\n",
    "    δ::Float64\n",
    "    L::Float64\n",
    "    C::Float64\n",
    "    F::Float64\n",
    "    S::Float64\n",
    "    R::Float64\n",
    "    w::Float64\n",
    "end\n",
    "\n",
    "function SteadyState(p::Params)\n",
    "    Yn_ss = exp(p.gbar)^(p.γ/(p.ϑ+p.γ))\n",
    "    Y_ss  = Yn_ss\n",
    "    π_ss  = 1.0\n",
    "    δ_ss  = 1.0\n",
    "    L_ss  = Y_ss/δ_ss\n",
    "    C_ss  = (1-p.gbar)*Y_ss\n",
    "    F_ss  = C_ss^(-p.γ)*Y_ss/(1-p.β*p.Θ*π_ss^(p.ϵ-1))\n",
    "    S_ss  = L_ss^p.ϑ*Y_ss/(1-p.β*p.Θ*π_ss^p.ϵ)\n",
    "    R_ss  = π_ss/p.β\n",
    "    w_ss  = (L_ss^p.ϑ)*(C_ss^p.γ)\n",
    "\n",
    "    SteadyState(Yn_ss, Y_ss, π_ss, δ_ss, L_ss, C_ss, F_ss, S_ss, R_ss, w_ss)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "Given an instance of `Params` and `SteadyState`, we can construct the grid on\n",
    "which we will solve the model.\n",
    "\n",
    "The `Grids` type holds this grid as well as matrices used to compute\n",
    "expectations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Grids"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "immutable Grids\n",
    "    # period t grids\n",
    "    ηR::Vector{Float64}\n",
    "    ηa::Vector{Float64}\n",
    "    ηL::Vector{Float64}\n",
    "    ηu::Vector{Float64}\n",
    "    ηB::Vector{Float64}\n",
    "    ηG::Vector{Float64}\n",
    "    R::Vector{Float64}\n",
    "    δ::Vector{Float64}\n",
    "\n",
    "    # combined matrix and complete polynomial version of it\n",
    "    X::Matrix{Float64}\n",
    "    X0_G::Dict{Int,Matrix{Float64}}\n",
    "\n",
    "    # quadrature weights and nodes\n",
    "    ϵ_nodes::Matrix{Float64}\n",
    "    ω_nodes::Vector{Float64}\n",
    "\n",
    "    # period t+1 grids at all shocks\n",
    "    ηR1::Matrix{Float64}\n",
    "    ηa1::Matrix{Float64}\n",
    "    ηL1::Matrix{Float64}\n",
    "    ηu1::Matrix{Float64}\n",
    "    ηB1::Matrix{Float64}\n",
    "    ηG1::Matrix{Float64}\n",
    "end\n",
    "\n",
    "function Grids(p::Params, ss::SteadyState)\n",
    "    if p.grid_kind == :sobol\n",
    "        # Values of exogenous state variables are in the interval +/- σ/sqrt(1-ρ^2)\n",
    "        ub = [2 * p.σηR / sqrt(1 - p.ρηR^2),\n",
    "              2 * p.σηa / sqrt(1 - p.ρηa^2),\n",
    "              2 * p.σηL / sqrt(1 - p.ρηL^2),\n",
    "              2 * p.σηu / sqrt(1 - p.ρηu^2),\n",
    "              2 * p.σηB / sqrt(1 - p.ρηB^2),\n",
    "              2 * p.σηG / sqrt(1 - p.ρηG^2),\n",
    "              1.05,  # R\n",
    "              1.0    # δ\n",
    "             ]\n",
    "        lb = -ub\n",
    "        lb[[7, 8]] = [1.0, 0.95]  # adjust lower bound for R and δ\n",
    "\n",
    "        # construct SobolSeq\n",
    "        s = SobolSeq(length(ub), ub, lb)\n",
    "        # skip(s, m)  # See note in README of Sobol.jl\n",
    "\n",
    "        # gather points\n",
    "        seq = Array{Float64}(8, p.m)\n",
    "        for i in 1:p.m\n",
    "            seq[:, i] = next(s)\n",
    "        end\n",
    "        seq = seq'  # transpose so variables are in columns\n",
    "\n",
    "        ηR = seq[:, 1]\n",
    "        ηa = seq[:, 2]\n",
    "        ηL = seq[:, 3]\n",
    "        ηu = seq[:, 4]\n",
    "        ηB = seq[:, 5]\n",
    "        ηG = seq[:, 6]\n",
    "        R  = seq[:, 7]\n",
    "        δ  = seq[:, 8]\n",
    "    else  # assume random\n",
    "        # Values of exogenous state variables are distributed uniformly\n",
    "        # in the interval +/- std/sqrt(1-rho_nu^2)\n",
    "        ηR = (-2*p.σηR + 4*p.σηR*rand(p.m)) / sqrt(1 - p.ρηR^2)\n",
    "        ηa = (-2*p.σηa + 4*p.σηa*rand(p.m)) / sqrt(1 - p.ρηa^2)\n",
    "        ηL = (-2*p.σηL + 4*p.σηL*rand(p.m)) / sqrt(1 - p.ρηL^2)\n",
    "        ηu = (-2*p.σηu + 4*p.σηu*rand(p.m)) / sqrt(1 - p.ρηu^2)\n",
    "        ηB = (-2*p.σηB + 4*p.σηB*rand(p.m)) / sqrt(1 - p.ρηB^2)\n",
    "        ηG = (-2*p.σηG + 4*p.σηG*rand(p.m)) / sqrt(1 - p.ρηG^2)\n",
    "\n",
    "        # Values of endogenous state variables are distributed uniformly\n",
    "        # in the intervals [1 1.05] and [0.95 1], respectively\n",
    "        R = 1 + 0.05*rand(p.m)\n",
    "        δ = 0.95 + 0.05*rand(p.m)\n",
    "    end\n",
    "\n",
    "    X = [log.(R) log.(δ) ηR ηa ηL ηu ηB ηG]\n",
    "    X0_G = Dict(\n",
    "        1 => complete_polynomial(X, 1),\n",
    "        p.deg => complete_polynomial(X, p.deg)\n",
    "    )\n",
    "\n",
    "    ϵ_nodes, ω_nodes = qnwmonomial1(vcov(p))\n",
    "\n",
    "    ηR1 = p.ρηR.*ηR .+ ϵ_nodes[:, 1]'\n",
    "    ηa1 = p.ρηa.*ηa .+ ϵ_nodes[:, 2]'\n",
    "    ηL1 = p.ρηL.*ηL .+ ϵ_nodes[:, 3]'\n",
    "    ηu1 = p.ρηu.*ηu .+ ϵ_nodes[:, 4]'\n",
    "    ηB1 = p.ρηB.*ηB .+ ϵ_nodes[:, 5]'\n",
    "    ηG1 = p.ρηG.*ηG .+ ϵ_nodes[:, 6]'\n",
    "\n",
    "    Grids(ηR, ηa, ηL, ηu, ηB, ηG, R, δ, X, X0_G, ϵ_nodes, ω_nodes,\n",
    "          ηR1, ηa1, ηL1, ηu1, ηB1, ηG1)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "Finally, we construct the `Model` type, which has an instance of `Params`,\n",
    "`SteadyState` and `Grids` as its three fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "outputExpanded": false
   },
   "outputs": [],
   "source": [
    "type Model\n",
    "    p::Params\n",
    "    s::SteadyState\n",
    "    g::Grids\n",
    "end\n",
    "\n",
    "function Model(;kwargs...)\n",
    "    p = Params(;kwargs...)\n",
    "    s = SteadyState(p)\n",
    "    g = Grids(p, s)\n",
    "    Model(p, s, g)\n",
    "end\n",
    "\n",
    "Base.show(io::IO, m::Model) = println(io, \"New Keynesian model\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "Now that we have a model, we will construct one more helper function that takes\n",
    "the control variables $(S, F, C)$ and shocks $(\\delta, R, \\eta_G, \\eta_a,\n",
    "\\eta_L, \\eta_R)$ and applies equilibrium conditions to produce $(\\pi, \\delta',\n",
    "Y, L, Y_n, R')$. We will use this function in both the solution and simulation\n",
    "routines below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "outputExpanded": false
   },
   "outputs": [],
   "source": [
    "function Base.step(m::Model, S, F, C, δ0, R0, ηG, ηa, ηL, ηR)\n",
    "    Θ, ϵ, gbar, ϑ, γ = m.p.Θ, m.p.ϵ, m.p.gbar, m.p.ϑ, m.p.γ\n",
    "    β, μ, ϕ_π, ϕ_y, πs = m.p.β, m.p.μ, m.p.ϕ_π, m.p.ϕ_y, m.s.π\n",
    "\n",
    "    # Compute pie(t) from condition (35) in MM (2015)\n",
    "    π0 = ((1-(1-Θ)*(S./F).^(1-ϵ))/Θ).^(1/(ϵ-1))\n",
    "\n",
    "    # Compute delta(t) from condition (36) in MM (2015)\n",
    "    δ1 = ((1-Θ)*((1-Θ*π0.^(ϵ-1))/(1-Θ)).^(ϵ/(ϵ-1))+Θ*π0.^ϵ./δ0).^(-1.0)\n",
    "\n",
    "    # Compute Y(t) from condition (38) in MM (2015)\n",
    "    Y0 = C./(1-gbar./exp.(ηG))\n",
    "\n",
    "    # Compute L(t) from condition (37) in MM (2015)\n",
    "    L0 = Y0./exp.(ηa)./δ1\n",
    "\n",
    "    #  Compute Yn(t) from condition (31) in MM (2015)\n",
    "    Yn0 = (exp.(ηa).^(1+ϑ).*(1-gbar./exp.(ηG)).^(-γ)./exp.(ηL)).^(1/(ϑ+γ))\n",
    "\n",
    "    # Compute R(t) from conditions (27), (39) in MM (2015) -- Taylor rule\n",
    "    R1 = πs ./ β .* (R0*β./πs).^μ.*((π0./πs).^ϕ_π .* (Y0./Yn0).^ϕ_y).^(1-μ).*exp.(ηR)\n",
    "\n",
    "    π0, δ1, Y0, L0, Yn0, R1\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "### Solution routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "solve (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# construct an initial guess for the solution\n",
    "function initial_coefs(m::Model)\n",
    "    npol = size(m.g.X0_G[1], 2)\n",
    "    coefs = fill(1e-5, npol, 3)\n",
    "    coefs[1, :] = [m.s.S, m.s.F, m.s.C^(-m.p.γ)]\n",
    "    coefs\n",
    "end\n",
    "\n",
    "function solve(m::Model)\n",
    "    # simplify notation\n",
    "    n, n_nodes = size(m.g.ηR, 1), length(m.g.ω_nodes)\n",
    "\n",
    "    ## allocate memory\n",
    "    # euler equations\n",
    "    e = zeros(n, 3)\n",
    "\n",
    "    # previous iteration S, F, C\n",
    "    S0_old_G = ones(n)\n",
    "    F0_old_G = ones(n)\n",
    "    C0_old_G = ones(n)\n",
    "\n",
    "    # current iteration S, F, C\n",
    "    S0_new_G = ones(n)\n",
    "    F0_new_G = ones(n)\n",
    "    C0_new_G = ones(n)\n",
    "\n",
    "    # future S, F, C\n",
    "    S1 = zeros(n, n_nodes)\n",
    "    F1 = zeros(n, n_nodes)\n",
    "    C1 = zeros(n, n_nodes)\n",
    "    π1 = zeros(n, n_nodes)\n",
    "\n",
    "    local coefs\n",
    "\n",
    "    for deg in [1, m.p.deg]\n",
    "        # set up matrices for this degree\n",
    "        X0_G = m.g.X0_G[deg]\n",
    "\n",
    "        # future basis matrix for S, F, C\n",
    "        X1 = Array{Float64}(n, n_complete(8, deg))\n",
    "\n",
    "        # initialize coefs\n",
    "        if deg == 1\n",
    "            coefs = initial_coefs(m)\n",
    "        else\n",
    "            # for higher order, fit the polynomial on the states e we left\n",
    "            # off with in the first order solution\n",
    "            coefs = X0_G \\ e\n",
    "        end\n",
    "\n",
    "        err = 1.0\n",
    "\n",
    "        # solve at this degree of complete polynomial\n",
    "        while err > m.p.tol\n",
    "            # Current choices (at t)\n",
    "            # ------------------------------\n",
    "            S0 = X0_G*coefs[:, 1]                # Compute S(t) using coefs\n",
    "            F0 = X0_G*coefs[:, 2]                # Compute F(t) using coefs\n",
    "            C0 = (X0_G*coefs[:, 3]).^(-1/m.p.γ)  # Compute C(t) using coefs\n",
    "\n",
    "            π0, δ1, Y0, L0, Yn0, R1 = step(m, S0, F0, C0, m.g.δ, m.g.R, m.g.ηG,\n",
    "                                           m.g.ηa, m.g.ηL, m.g.ηR)\n",
    "\n",
    "            if m.p.zlb R1 .= max.(R1, 1.0) end\n",
    "\n",
    "            for u in 1:n_nodes\n",
    "                # Form complete polynomial of degree \"Degree\" (at t+1) on future state\n",
    "                complete_polynomial!(\n",
    "                    X1,\n",
    "                    hcat(log.(R1), log.(δ1), m.g.ηR1[:, u], m.g.ηa1[:, u],\n",
    "                         m.g.ηL1[:, u], m.g.ηu1[:, u], m.g.ηB1[:, u], m.g.ηG1[:, u]),\n",
    "                    deg\n",
    "                )\n",
    "\n",
    "                S1[:, u] = X1*coefs[:, 1]                # Compute S(t+1)\n",
    "                F1[:, u] = X1*coefs[:, 2]                # Compute F(t+1)\n",
    "                C1[:, u] = (X1*coefs[:, 3]).^(-1/m.p.γ)  # Compute C(t+1)\n",
    "            end\n",
    "\n",
    "            # Compute next-period π using condition\n",
    "            # (35) in MM (2015)\n",
    "            π1 .= ((1-(1-m.p.Θ).*(S1./F1).^(1-m.p.ϵ))/m.p.Θ).^(1/(m.p.ϵ-1))\n",
    "\n",
    "            # Evaluate conditional expectations in the Euler equations\n",
    "            #---------------------------------------------------------\n",
    "            e[:, 1] = exp.(m.g.ηu).*exp.(m.g.ηL).*L0.^m.p.ϑ.*Y0./exp.(m.g.ηa) + (m.p.β*m.p.Θ*π1.^m.p.ϵ.*S1)*m.g.ω_nodes\n",
    "            e[:, 2] = exp.(m.g.ηu).*C0.^(-m.p.γ).*Y0 + (m.p.β*m.p.Θ*π1.^(m.p.ϵ-1).*F1)*m.g.ω_nodes\n",
    "            e[:, 3] = m.p.β*exp.(m.g.ηB)./exp.(m.g.ηu).*R1.*((exp.(m.g.ηu1).*C1.^(-m.p.γ)./π1)*m.g.ω_nodes)\n",
    "\n",
    "            # Variables of the current iteration\n",
    "            #-----------------------------------\n",
    "            copy!(S0_new_G, S0)\n",
    "            copy!(F0_new_G, F0)\n",
    "            copy!(C0_new_G, C0)\n",
    "\n",
    "            # Compute and update the coefficients of the decision functions\n",
    "            # -------------------------------------------------------------\n",
    "            coefs_hat = X0_G\\e   # Compute the new coefficients of the decision\n",
    "                                 # functions using a backslash operator\n",
    "\n",
    "            # Update the coefficients using damping\n",
    "            coefs = m.p.damp*coefs_hat + (1-m.p.damp)*coefs\n",
    "\n",
    "            # Evaluate the percentage (unit-free) difference between the values\n",
    "            # on the grid from the previous and current iterations\n",
    "            # -----------------------------------------------------------------\n",
    "            # The convergence criterion is adjusted to the damping parameters\n",
    "            err = mean(abs, 1-S0_new_G./S0_old_G) +\n",
    "                  mean(abs, 1-F0_new_G./F0_old_G) +\n",
    "                  mean(abs, 1-C0_new_G./C0_old_G)\n",
    "\n",
    "            # Store the obtained values for S(t), F(t), C(t) on the grid to\n",
    "            # be used on the subsequent iteration in Section 10.2.6\n",
    "            #-----------------------------------------------------------------------\n",
    "            copy!(S0_old_G, S0_new_G)\n",
    "            copy!(F0_old_G, F0_new_G)\n",
    "            copy!(C0_old_G, C0_new_G)\n",
    "        end\n",
    "    end\n",
    "\n",
    "    coefs\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "### Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Simulation"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type Simulation\n",
    "    # shocks\n",
    "    ηR::Vector{Float64}\n",
    "    ηa::Vector{Float64}\n",
    "    ηL::Vector{Float64}\n",
    "    ηu::Vector{Float64}\n",
    "    ηB::Vector{Float64}\n",
    "    ηG::Vector{Float64}\n",
    "\n",
    "    # variables\n",
    "    δ::Vector{Float64}\n",
    "    R::Vector{Float64}\n",
    "    S::Vector{Float64}\n",
    "    F::Vector{Float64}\n",
    "    C::Vector{Float64}\n",
    "    π::Vector{Float64}\n",
    "    Y::Vector{Float64}\n",
    "    L::Vector{Float64}\n",
    "    Yn::Vector{Float64}\n",
    "    w::Vector{Float64}\n",
    "end\n",
    "\n",
    "function Simulation(m::Model, coefs::Matrix, capT::Int=10201)\n",
    "\n",
    "    # 11. Simualating a time-series solution\n",
    "    #---------------------------------------\n",
    "\n",
    "    # Initialize the values of 6 exogenous shocks and draw innovations\n",
    "    #-----------------------------------------------------------------\n",
    "    ηR = zeros(capT)\n",
    "    ηa = zeros(capT)\n",
    "    ηL = zeros(capT)\n",
    "    ηu = zeros(capT)\n",
    "    ηB = zeros(capT)\n",
    "    ηG = zeros(capT)\n",
    "\n",
    "    # Generate the series for shocks\n",
    "    #-------------------------------\n",
    "    @inbounds for t in 1:capT-1\n",
    "        ηR[t+1] = m.p.ρηR*ηR[t] + m.p.σηR*randn()\n",
    "        ηa[t+1] = m.p.ρηa*ηa[t] + m.p.σηa*randn()\n",
    "        ηL[t+1] = m.p.ρηL*ηL[t] + m.p.σηL*randn()\n",
    "        ηu[t+1] = m.p.ρηu*ηu[t] + m.p.σηu*randn()\n",
    "        ηB[t+1] = m.p.ρηB*ηB[t] + m.p.σηB*randn()\n",
    "        ηG[t+1] = m.p.ρηG*ηG[t] + m.p.σηG*randn()\n",
    "    end\n",
    "\n",
    "    δ  = ones(capT+1) # Allocate memory for the time series of delta(t)\n",
    "    R  = ones(capT+1) # Allocate memory for the time series of R(t)\n",
    "    S  = Array{Float64}(capT)   # Allocate memory for the time series of S(t)\n",
    "    F  = Array{Float64}(capT)   # Allocate memory for the time series of F(t)\n",
    "    C  = Array{Float64}(capT)   # Allocate memory for the time series of C(t)\n",
    "    π  = Array{Float64}(capT)   # Allocate memory for the time series of π(t)\n",
    "    Y  = Array{Float64}(capT)   # Allocate memory for the time series of Y(t)\n",
    "    L  = Array{Float64}(capT)   # Allocate memory for the time series of L(t)\n",
    "    Yn = Array{Float64}(capT)   # Allocate memory for the time series of Yn(t)\n",
    "    w  = Array{Float64}(capT)\n",
    "\n",
    "    pol_bases = Array{Float64}(1, size(coefs, 1))\n",
    "    @inbounds for t in 1:capT\n",
    "        # Construct the matrix of explanatory variables \"pol_bases\" on the series\n",
    "        # of state variables; columns of \"pol_bases\" are given by the basis\n",
    "        # functions of the polynomial of degree 2\n",
    "        complete_polynomial!(\n",
    "            pol_bases,\n",
    "            hcat(log(R[t]), log(δ[t]), ηR[t], ηa[t], ηL[t], ηu[t], ηB[t], ηG[t]),\n",
    "            m.p.deg\n",
    "        )\n",
    "        S[t], F[t], MU = pol_bases*coefs\n",
    "        C[t] = (MU).^(-1/m.p.γ)\n",
    "\n",
    "        π[t], δ[t+1], Y[t], L[t], Yn[t], R[t+1] = step(m, S[t], F[t], C[t], δ[t],\n",
    "                                                       R[t], ηG[t], ηa[t], ηL[t], ηR[t])\n",
    "\n",
    "        # Compute real wage\n",
    "        w[t] = exp(ηL[t])*(L[t]^m.p.ϑ)*(C[t]^m.p.γ)\n",
    "\n",
    "        # If ZLB is imposed, set R(t)=1 if ZLB binds\n",
    "        if m.p.zlb; R[t+1] = max(R[t+1],1.0); end\n",
    "    end\n",
    "\n",
    "    Simulation(ηR, ηa, ηL, ηu, ηB, ηG, δ, R, S, F, C, π, Y, L, Yn, w)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "### Accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "max_E (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type Residuals\n",
    "    resids::Matrix{Float64}\n",
    "end\n",
    "\n",
    "function Residuals(m::Model, coefs::Matrix, sim::Simulation; burn::Int=200)\n",
    "    capT = length(sim.w)\n",
    "    resids = zeros(9, capT)\n",
    "\n",
    "    # Integration method for evaluating accuracy\n",
    "    # ------------------------------------------\n",
    "    # Monomial integration rule with 2N^2+1 nodes\n",
    "    ϵ_nodes, ω_nodes = qnwmonomial2(vcov(m.p))\n",
    "    n_nodes = length(ω_nodes)\n",
    "\n",
    "    # Allocate for arrays needed in the loop\n",
    "    basis_mat = Array{Float64}(n_nodes, 8)\n",
    "    X1 = Array{Float64}(n_nodes, size(coefs, 1))\n",
    "\n",
    "    ηR1 = Array{Float64}(n_nodes)\n",
    "    ηa1 = Array{Float64}(n_nodes)\n",
    "    ηL1 = Array{Float64}(n_nodes)\n",
    "    ηu1 = Array{Float64}(n_nodes)\n",
    "    ηB1 = Array{Float64}(n_nodes)\n",
    "    ηG1 = Array{Float64}(n_nodes)\n",
    "\n",
    "    for t = 1:capT                 # For each given point,\n",
    "        # Take the corresponding value for shocks at t\n",
    "        #---------------------------------------------\n",
    "        ηR0 = sim.ηR[t]  # ηR(t)\n",
    "        ηa0 = sim.ηa[t]  # ηa(t)\n",
    "        ηL0 = sim.ηL[t]  # ηL(t)\n",
    "        ηu0 = sim.ηu[t]  # ηu(t)\n",
    "        ηB0 = sim.ηB[t]  # ηB(t)\n",
    "        ηG0 = sim.ηG[t]  # ηG(t)\n",
    "\n",
    "        # Exctract time t values for all other variables (and t+1 for R, δ)\n",
    "        #------------------------------------------------------------------\n",
    "        R0  = sim.R[t]   # R(t-1)\n",
    "        δ0  = sim.δ[t]   # δ(t-1)\n",
    "        R1  = sim.R[t+1] # R(t)\n",
    "        δ1  = sim.δ[t+1] # δ(t)\n",
    "\n",
    "        L0  = sim.L[t]   # L(t)\n",
    "        Y0  = sim.Y[t]   # Y(t)\n",
    "        Yn0 = sim.Yn[t]  # Yn(t)\n",
    "        π0  = sim.π[t]   # π(t)\n",
    "        S0  = sim.S[t]   # S(t)\n",
    "        F0  = sim.F[t]   # F(t)\n",
    "        C0  = sim.C[t]   # C(t)\n",
    "\n",
    "        # Fill basis matrix with R1, δ1 and shocks\n",
    "        #-----------------------------------------\n",
    "        # Note that we do not premultiply by standard deviations as ϵ_nodes\n",
    "        # already include them. All these variables are vectors of length n_nodes\n",
    "        copy!(ηR1, ηR0*m.p.ρηR + ϵ_nodes[:, 1])\n",
    "        copy!(ηa1, ηa0*m.p.ρηa + ϵ_nodes[:, 2])\n",
    "        copy!(ηL1, ηL0*m.p.ρηL + ϵ_nodes[:, 3])\n",
    "        copy!(ηu1, ηu0*m.p.ρηu + ϵ_nodes[:, 4])\n",
    "        copy!(ηB1, ηB0*m.p.ρηB + ϵ_nodes[:, 5])\n",
    "        copy!(ηG1, ηG0*m.p.ρηG + ϵ_nodes[:, 6])\n",
    "\n",
    "        basis_mat[:, 1] = log(R1)\n",
    "        basis_mat[:, 2] = log(δ1)\n",
    "        basis_mat[:, 3] = ηR1\n",
    "        basis_mat[:, 4] = ηa1\n",
    "        basis_mat[:, 5] = ηL1\n",
    "        basis_mat[:, 6] = ηu1\n",
    "        basis_mat[:, 7] = ηB1\n",
    "        basis_mat[:, 8] = ηG1\n",
    "\n",
    "        # Future choices at t+1\n",
    "        #----------------------\n",
    "        # Form a complete polynomial of degree \"Degree\" (at t+1) on future state\n",
    "        # variables; n_nodes-by-npol\n",
    "        complete_polynomial!(X1, basis_mat, m.p.deg)\n",
    "\n",
    "        # Compute S(t+1), F(t+1) and C(t+1) in all nodes using coefs\n",
    "        S1 = X1*coefs[:, 1]\n",
    "        F1 = X1*coefs[:, 2]\n",
    "        C1 = (X1*coefs[:, 3]).^(-1/m.p.γ)\n",
    "\n",
    "        # Compute π(t+1) using condition (35) in MM (2015)\n",
    "        π1 = ((1-(1-m.p.Θ)*(S1./F1).^(1-m.p.ϵ))/m.p.Θ).^(1/(m.p.ϵ-1))\n",
    "\n",
    "        # Compute residuals for each of the 9 equilibrium conditions\n",
    "        #-----------------------------------------------------------\n",
    "        resids[1, t] = 1-dot(ω_nodes,\n",
    "            (exp(ηu0)*exp(ηL0)*L0^m.p.ϑ*Y0/exp(ηa0) + m.p.β*m.p.Θ*π1.^m.p.ϵ.*S1)/S0\n",
    "        )\n",
    "        resids[2, t] = 1 - dot(ω_nodes,\n",
    "            (exp(ηu0)*C0^(-m.p.γ)*Y0 + m.p.β*m.p.Θ*π1.^(m.p.ϵ-1).*F1)/F0\n",
    "        )\n",
    "        resids[3, t] = 1.0 -dot(ω_nodes,\n",
    "            (m.p.β*exp(ηB0)/exp(ηu0)*R1*exp.(ηu1).*C1.^(-m.p.γ)./π1)/C0^(-m.p.γ)\n",
    "        )\n",
    "        resids[4, t] = 1-((1-m.p.Θ*π0^(m.p.ϵ-1))/(1-m.p.Θ))^(1/(1-m.p.ϵ))*F0/S0\n",
    "        resids[5, t] = 1-((1-m.p.Θ)*((1-m.p.Θ*π0^(m.p.ϵ-1))/(1-m.p.Θ))^(m.p.ϵ/(m.p.ϵ-1)) + m.p.Θ*π0^m.p.ϵ/δ0)^(-1)/δ1\n",
    "        resids[6, t] = 1-exp(ηa0)*L0*δ1/Y0\n",
    "        resids[7, t] = 1-(1-m.p.gbar/exp(ηG0))*Y0/C0\n",
    "        resids[8, t] = 1-(exp(ηa0)^(1+m.p.ϑ)*(1-m.p.gbar/exp(ηG0))^(-m.p.γ)/exp(ηL0))^(1/(m.p.ϑ+m.p.γ))/Yn0\n",
    "        resids[9, t] = 1-m.s.π/m.p.β*(R0*m.p.β/m.s.π)^m.p.μ*((π0/m.s.π)^m.p.ϕ_π * (Y0/Yn0)^m.p.ϕ_y)^(1-m.p.μ)*exp(ηR0)/R1   # Taylor rule\n",
    "\n",
    "        # If the ZLB is imposed and R>1, the residuals in the Taylor rule (the\n",
    "        # 9th equation) are zero\n",
    "        if m.p.zlb && R1 <= 1; resids[9, t] = 0.0; end\n",
    "\n",
    "    end\n",
    "    # discard the first burn observations\n",
    "    Residuals(resids[:, burn+1:end])\n",
    "end\n",
    "\n",
    "Base.mean(r::Residuals) = log10(mean(abs, r.resids))\n",
    "Base.max(r::Residuals) = log10(maximum(abs, r.resids))\n",
    "max_E(r::Residuals) = log10.(maximum(abs, r.resids, 2))[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "outputExpanded": false
   },
   "source": [
    "### Running the code\n",
    "\n",
    "Now that we've done all the hard work to define the model, its solution and\n",
    "simulation, and accuracy checks, let's put things together and run the code!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "outputExpanded": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "build_paper_table (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function main(m=Model(); io::IO=STDOUT)\n",
    "    tic(); coefs = solve(m); solve_time = toq()\n",
    "\n",
    "    # simulate the model\n",
    "    tic(); sim = Simulation(m, coefs); simulation_time = toq()\n",
    "\n",
    "    # check accuracy\n",
    "    tic(); resids = Residuals(m, coefs, sim); resids_time = toq()\n",
    "\n",
    "    err_by_eq = max_E(resids)\n",
    "    l1 = log10(sum(maximum(abs, resids.resids, 2)))\n",
    "    l∞ = max(resids)\n",
    "    tot_time = solve_time + simulation_time + resids_time\n",
    "    round3(x) = round(x, 3)\n",
    "    round2(x) = round(x, 2)\n",
    "\n",
    "    println(io, \"Solver time (in seconds): $(solve_time)\")\n",
    "    println(io, \"Simulation time (in seconds): $(simulation_time)\")\n",
    "    println(io, \"Residuals time (in seconds): $(resids_time)\")\n",
    "    println(io, \"total time (in seconds): $(tot_time)\")\n",
    "    println(io, \"\\nAPPROXIMATION ERRORS (log10):\");\n",
    "    println(io, \"\\ta) mean error in the model equations: $(round3(mean(resids)))\");\n",
    "    println(io, \"\\tb) sum of max error per equation: $(round3(l1))\");\n",
    "    println(io, \"\\tc) max error in the model equations: $(round3(l∞))\");\n",
    "    println(io, \"\\td) max error by equation:$(round3.(err_by_eq))\");\n",
    "    println(io, \"tex row\\n\", join(round2.([l1, l∞, tot_time]), \" & \"))\n",
    "\n",
    "    solve_time, simulation_time, resids_time, coefs, sim, resids\n",
    "end\n",
    "\n",
    "\n",
    "function build_paper_table()\n",
    "    # call main once to precompile all routines\n",
    "    main()\n",
    "\n",
    "    open(\"output.log\", \"w\") do f\n",
    "        for params in [Dict(:πstar=>1.0, :σηL=>0.1821, :zlb=>false),\n",
    "                       Dict(:πstar=>1.0, :σηL=>0.4054, :zlb=>false),\n",
    "                       Dict(:πstar=>1.0, :σηL=>0.1821, :zlb=>true)]\n",
    "            for grid_kind in [:sobol, :random]\n",
    "                m = Model(;grid_kind=grid_kind, params...)\n",
    "\n",
    "                println(f, \"Working with params:\")\n",
    "                show(f, MIME\"text/plain\"(), params)\n",
    "                println(f, \"\\nand grid type: $(grid_kind)\")\n",
    "                main(m, io=f)\n",
    "                println(f, \"\\n\"^5)\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solver time (in seconds): 2.461723847\n",
      "Simulation time (in seconds): 0.092354514\n",
      "Residuals time (in seconds): 0.440504231\n",
      "total time (in seconds): 2.9945825920000004\n",
      "\n",
      "APPROXIMATION ERRORS (log10):\n",
      "\ta) mean error in the model equations: -4.486\n",
      "\tb) sum of max error per equation: -1.803\n",
      "\tc) max error in the model equations: -2.139\n",
      "\td) max error by equation:[-2.139, -2.421, -2.329, -15.0, -Inf, -15.654, -15.654, -Inf, -Inf]\n",
      "tex row\n",
      "-1.8 & -2.14 & 2.99\n"
     ]
    }
   ],
   "source": [
    "main();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernel_info": {
   "name": "julia-0.5"
  },
  "kernelspec": {
   "display_name": "Julia 0.6.0-rc1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}