diff --git a/NKModel_CLMM_Model.ipynb b/NKModel_CLMM_Model.ipynb new file mode 100644 index 0000000..fa17511 --- /dev/null +++ b/NKModel_CLMM_Model.ipynb @@ -0,0 +1,113 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# A stylized New Keynesian Model\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 this notebook we summarize the key equations for the stylized New Keynesian model we solved in the paper.\n", + "\n", + "For more information on the model itself see section 5 of \n", + "\n", + "> Maliar, L., & Maliar, S. (2015). Merging simulation and projection approaches to solve high-dimensional problems with an application to a new Keynesian model. Quantitative Economics, 15(7), 424. http://doi.org/10.1186/s13059-014-0424-0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The model features Calvo-type price frictions and a\n", + "Taylor (1993) rule. The economy is populated by households, final-good\n", + "firms, intermediate-good firms, a monetary authority and government.\n", + "\n", + "In the equations below, the following variables appear:\n", + "\n", + "- $\\pi$: profit of intermediate good firms\n", + "- $S$, $F$: intermediate variables introduced to help notation below\n", + "- $C$: Consumption\n", + "- $G$: Level of government spending\n", + "- $Y$: Final goods production\n", + "- $L$: Labor supply\n", + "- $\\Delta$: measure of price dispersion across firms (also referred to as efficiency distortion)\n", + "- $Y_N$: Natural level of final goods production -- it is the level of $Y$ in the planner's solution, or the level of $Y$ in the model without distorting taxes.\n", + "- $\\bar{G}$: steady-state share of government spendint in output\n", + "- $R$: gross nominal interest rate\n", + "- $\\eta_u$: exogenous preference shock to utility\n", + "- $\\eta_L$: exogenous preference shock to labor dis-utility\n", + "- $\\eta_B$: exogenous premium on bond returns\n", + "- $\\eta_a$: log of productivity of intermediate goods firms\n", + "- $\\eta_G$: government spending shock\n", + "- $\\eta_R$: monetary shock to interest rate\n", + "\n", + "We also see the following parameters:\n", + "\n", + "- $\\beta$: discount factor\n", + "- $\\varepsilon$: Parameter in Dixit-Stiglitz aggregator over intermediate goods\n", + "- $\\theta$: Calvo-parameter -- a fraction (1-$\\theta$) of firms set prices optimally each period.\n", + "- $\\gamma$, $\\vartheta$: utility function parameters\n", + "- $\\pi^*$: target inflation\n", + "- $\\phi_y$, $\\phi_{\\pi}$, $\\mu$: Taylor rule parameters" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In our computation we will approximate $S_t$, $F_t$ and $C_t$ using a complete monomial of degree N. Given these three variables, we express the equilibrium conditions of the model in the following way:\n", + "\n", + "\\begin{eqnarray}\n", + "\\pi_{t} &=& \\left(\\frac{1-(1-\\theta)}{\\theta}\\left(\\frac{S_t}{F_t}\\right)^{1- \\varepsilon}\\right)^{\\frac{1}{\\varepsilon-1}} \\\\\n", + "\\Delta _{t} &=&\\left[ \\left( 1-\\theta \\right) \\left[ \\frac{1-\\theta \\pi\n", + "_{t}^{\\varepsilon -1}}{1-\\theta }\\right] ^{\\frac{\\varepsilon }{\\varepsilon -1%\n", + "}}+\\theta \\frac{\\pi _{t}^{\\varepsilon }}{\\Delta _{t-1}}\\right] ^{-1} \\\\\n", + "Y_t &=& \\frac{C_t}{1-\\frac{\\bar{G}}{\\exp\\left(\\eta_{G,t}\\right)}} \\\\\n", + "L_t &=& \\frac{Y_t }{\\exp\\left(\\eta_{a,t}\\right) \\Delta_t} \\\\\n", + "Y_{N,t} &=& \\left[ \\frac{\\exp \\left( \\eta _{a,t}\\right) ^{1+\\vartheta }\\left[1- \\frac{\\bar{G}}{\\exp \\left( \\eta _{G,t}\\right)} \\right] ^{-\\gamma}}{\\exp \\left( \\eta_{L,t}\\right) }\\right] ^{\\frac{1}{\\vartheta +\\gamma }}\\\\\n", + "R_{t} &=& \\max \\left\\{1, \\frac{\\pi^*}{\\beta} \\left(R_{t-1} \\frac{\\beta}{\\pi^*} \\right)^{\\mu} \\left(\\left(\\frac{\\pi_t}{\\pi^*}\\right)^{\\phi_{\\pi}} \\left(\\frac{Y_t}{Y_{N,t}} \\right)^{\\phi_y} \\right)^{1-\\mu}\\exp\\left(\\eta_{R,t}\\right) \\right\\}.\n", + "\\end{eqnarray}\n", + "\n", + "The Euler equations are given by\n", + "\n", + "\\begin{eqnarray}\n", + "S_{t} &=&\\frac{\\exp \\left( \\eta _{u,t}+\\eta _{L,t}\\right) }{\\left[ \\exp\n", + "\\left( \\eta _{a,t}\\right) \\right] ^{\\vartheta +1}}\\frac{\\left(\n", + "G_{t}^{-1}C_{t}\\right) ^{1+\\vartheta }}{\\left( \\Delta _{t}\\right)\n", + "^{\\vartheta }}+\\beta \\theta E_{t}\\left \\{ \\pi _{t+1}^{\\varepsilon\n", + "}S_{t+1}\\right \\} \\\\\n", + "F_{t} &=&\\exp \\left( \\eta _{u,t}\\right) C_{t}^{1-\\gamma }G_{t}^{-1}+\\beta\n", + "\\theta E_{t}\\left \\{ \\pi _{t+1}^{\\varepsilon -1}F_{t+1}\\right \\} \\\\\n", + "C_{t}^{-\\gamma } &=&\\beta \\frac{\\exp \\left( \\eta _{B,t}\\right) }{\\exp \\left(\n", + "\\eta _{u,t}\\right) }R_{t}E_{t}\\left[ \\frac{C_{t+1}^{-\\gamma }\\exp \\left(\n", + "\\eta _{u,t+1}\\right) }{\\pi _{t+1}}\\right]\n", + "\\end{eqnarray}" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/NKModel_CLMM_julia.ipynb b/NKModel_CLMM_julia.ipynb new file mode 100644 index 0000000..0278c36 --- /dev/null +++ b/NKModel_CLMM_julia.ipynb @@ -0,0 +1,968 @@ +{ + "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 +} diff --git a/NKModel_CLMM_matlab.ipynb b/NKModel_CLMM_matlab.ipynb new file mode 100644 index 0000000..4db85fb --- /dev/null +++ b/NKModel_CLMM_matlab.ipynb @@ -0,0 +1,1203 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solving a New Keynesian model with Matlab\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 Matlab and the a jupyter Matlab kernel. We assume you already have Matlab installed and show how to install the jupyter matlab kernel in a few steps:\n", + "\n", + "1. Install Python version 3.5. We recommend following the instructions on [quantecon](https://lectures.quantecon.org/py/getting_started.html) and using the Anaconda Python distribution. Make sure you don't get a python version above version 3.5.\n", + "2. Install the [MATLAB Engine API for Python](https://www.mathworks.com/help/matlab/matlab_external/install-the-matlab-engine-for-python.html). You can follow the official instructions from MathWorks at the link in the previous sentence.\n", + "3. Install the `imatlab` kernel by doing the following from the command prompt (or terminal prompt for OSX/Linux userse):\n", + "```shell\n", + "python -m pip install matlab_kernel\n", + "python -m matlab_kernel install\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Numerical Tools\n", + "\n", + "\n", + "In order to implement our routines, we need a few functions defining numerical tools.\n", + "\n", + "Do to limitations in the [Matlab Juptyer kernel](https://github.com/Calysto/matlab_kernel/issues/35), we cannot define functions in the notebook itself. So, we will download the file we need from online, save them to the same directory as this notebook, and then add this directory to the Matlab path.\n", + "\n", + "So that we know what these files are doing, we will print the text of these files here in the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "ans =\n", + "\n", + " 2\n", + "\n" + ] + } + ], + "source": [ + "1+1" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "mon1_txt = fileread(websave('Monomials_1.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/Monomials_1.m'));\n", + "mon2_txt = fileread(websave('Monomials_2.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/Monomials_2.m'));\n", + "poly_txt = fileread(websave('Ord_Polynomial_N.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/Ord_Polynomial_N.m'));" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "mon1_txt =\n", + "\n", + " '% Monomials_1.m is a routine that constructs integration nodes and weights \n", + " % under N-dimensional monomial (non-product) integration rule with 2N nodes;\n", + " % see Judd, Maliar and Maliar, (2010), \"A Cluster-Grid Projection Method:\n", + " % Solving Problems with High Dimensionality\", NBER Working Paper 15965\n", + " % (henceforth, JMM, 2010).\n", + " % -------------------------------------------------------------------------\n", + " % Inputs: \"N\" is the number of random variables; N>=1;\n", + " % \"vcv\" is the variance-covariance matrix; N-by-N\n", + " \n", + " % Outputs: \"n_nodes\" is the total number of integration nodes; 2*N;\n", + " % \"epsi_nodes\" are the integration nodes; n_nodes-by-N;\n", + " % \"weight_nodes\" are the integration weights; n_nodes-by-1\n", + " % -------------------------------------------------------------------------\n", + " % Copyright � 2011 by Lilia Maliar and Serguei Maliar. All rights reserved.\n", + " % The code may be used, modified and redistributed under the terms provided\n", + " % in the file \"License_Agreement.txt\".\n", + " % -------------------------------------------------------------------------\n", + " \n", + " function [n_nodes,epsi_nodes,weight_nodes] = Monomials_1(N,vcv)\n", + " \n", + " n_nodes = 2*N; % Total number of integration nodes\n", + " \n", + " % 1. N-dimensional integration nodes for N uncorrelated random variables with\n", + " % zero mean and unit variance\n", + " % ---------------------------------------------------------------------------\n", + " z1 = zeros(n_nodes,N); % A supplementary matrix for integration nodes;\n", + " % n_nodes-by-N\n", + " \n", + " for i = 1:N % In each node, random variable i takes value either\n", + " % 1 or -1, and all other variables take value 0\n", + " z1(2*(i-1)+1:2*i,i) = [1; -1];\n", + " end % For example, for N = 2, z1 = [1 0; -1 0; 0 1; 0 -1]\n", + " \n", + " % z = z1*sqrt(N); % Integration nodes\n", + " \n", + " % 2. N-dimensional integration nodes and weights for N correlated random\n", + " % variables with zero mean and variance-covaraince matrix vcv\n", + " % ----------------------------------------------------------------------\n", + " sqrt_vcv = chol(vcv); % Cholesky decomposition of the variance-covariance\n", + " % matrix\n", + " \n", + " R = sqrt(N)*sqrt_vcv; % Variable R; see condition (20) in JMM (2010)\n", + " \n", + " epsi_nodes = z1*R; % Integration nodes; see condition (20) in JMM (2010);\n", + " % n_nodes-by-N\n", + " \n", + " % 3. Integration weights\n", + " %-----------------------\n", + " weight_nodes = ones(n_nodes,1)/n_nodes;\n", + " % Integration weights are equal for all integration\n", + " % nodes; n_nodes-by-1; the weights are the same for\n", + " % the cases of correlated and uncorrelated random\n", + " % variables\n", + " '\n", + "\n" + ] + } + ], + "source": [ + "mon1_txt" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "mon2_txt =\n", + "\n", + " '% Monomials_2.m is a routine that constructs integration nodes and weights \n", + " % under N-dimensional monomial (non-product) integration rule with 2N^2+1\n", + " % nodes; see Judd, Maliar and Maliar, (2010), \"A Cluster-Grid Projection\n", + " % Method: Solving Problems with High Dimensionality\", NBER Working Paper\n", + " % 15965 (henceforth, JMM, 2010).\n", + " % -------------------------------------------------------------------------\n", + " % Inputs: \"N\" is the number of random variables; N>=1;\n", + " % \"vcv\" is the variance-covariance matrix; N-by-N;\n", + " \n", + " % Outputs: \"n_nodes\" is the total number of integration nodes; 2*N^2+1;\n", + " % \"epsi_nodes\" are the integration nodes; n_nodes-by-N;\n", + " % \"weight_nodes\" are the integration weights; n_nodes-by-1\n", + " % -------------------------------------------------------------------------\n", + " % Copyright � 2011 by Lilia Maliar and Serguei Maliar. All rights reserved.\n", + " % The code may be used, modified and redistributed under the terms provided\n", + " % in the file \"License_Agreement.txt\".\n", + " % -------------------------------------------------------------------------\n", + " \n", + " \n", + " function [n_nodes,epsi_nodes,weight_nodes] = Monomials_2(N,vcv)\n", + " \n", + " n_nodes = 2*N^2+1; % Total number of integration nodes\n", + " \n", + " % 1. N-dimensional integration nodes for N uncorrelated random variables with\n", + " % zero mean and unit variance\n", + " % ---------------------------------------------------------------------------\n", + " \n", + " % 1.1 The origin point\n", + " % --------------------\n", + " z0 = zeros(1,N); % A supplementary matrix for integration nodes: the\n", + " % origin point\n", + " \n", + " % 1.2 Deviations in one dimension\n", + " % -------------------------------\n", + " z1 = zeros(2*N,N); % A supplementary matrix for integration nodes;\n", + " % n_nodes-by-N\n", + " \n", + " for i = 1:N % In each node, random variable i takes value either\n", + " % 1 or -1, and all other variables take value 0\n", + " z1(2*(i-1)+1:2*i,i) = [1; -1];\n", + " end % For example, for N = 2, z1 = [1 0; -1 0; 0 1; 0 -1]\n", + " \n", + " % 1.3 Deviations in two dimensions\n", + " % --------------------------------\n", + " z2 = zeros(2*N*(N-1),N); % A supplementary matrix for integration nodes;\n", + " % 2N(N-1)-by-N\n", + " \n", + " i=0; % In each node, a pair of random variables (p,q)\n", + " % takes either values (1,1) or (1,-1) or (-1,1) or\n", + " % (-1,-1), and all other variables take value 0\n", + " for p = 1:N-1\n", + " for q = p+1:N\n", + " i=i+1;\n", + " z2(4*(i-1)+1:4*i,p) = [1;-1;1;-1];\n", + " z2(4*(i-1)+1:4*i,q) = [1;1;-1;-1];\n", + " end\n", + " end % For example, for N = 2, z2 = [1 1;1 -1;-1 1;-1 1]\n", + " \n", + " % z = [z0;z1*sqrt(N+2);z2*sqrt((N+2)/2)]; % Integration nodes\n", + " \n", + " % 2. N-dimensional integration nodes and weights for N correlated random\n", + " % variables with zero mean and variance-covaraince matrix vcv\n", + " % ----------------------------------------------------------------------\n", + " sqrt_vcv = chol(vcv); % Cholesky decomposition of the variance-\n", + " % covariance matrix\n", + " \n", + " R = sqrt(N+2)*sqrt_vcv; % Variable R; see condition (21) in JMM\n", + " % (2010)\n", + " \n", + " S = sqrt((N+2)/2)* sqrt_vcv; % Variable S; see condition (21) in JMM\n", + " % (2010)\n", + " \n", + " epsi_nodes = [z0;z1*R;z2*S];\n", + " % Integration nodes; see condition (21)\n", + " % in JMM (2010); n_nodes-by-N\n", + " \n", + " % 3. Integration weights\n", + " %-----------------------\n", + " weight_nodes = [2/(N+2)*ones(size(z0,1),1);(4-N)/2/(N+2)^2*ones(size(z1,1),1);1/(N+2)^2*ones(size(z2,1),1)];\n", + " % See condition (21) in JMM (2010);\n", + " % n_nodes-by-1; the weights are the same\n", + " % for the cases of correlated and\n", + " % uncorrelated random variables\n", + " '\n", + "\n" + ] + } + ], + "source": [ + "mon2_txt" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "poly_txt =\n", + "\n", + " '% Ord_Polynomial_N.m is a routine that constructs the basis functions of \n", + " % complete ordinary polynomial of the degrees from one to five for the\n", + " % multi-dimensional case; see \"Numerically Stable and Accurate Stochastic\n", + " % Simulation Approaches for Solving Dynamic Economic Models\" by Kenneth L.\n", + " % Judd, Lilia Maliar and Serguei Maliar, (2011), Quantitative Economics 2/2,\n", + " % 173�210 (henceforth, JMM, 2011).\n", + " %\n", + " % This version: July 14, 2011. First version: August 27, 2009.\n", + " % -------------------------------------------------------------------------\n", + " % Inputs: \"z\" is the data points on which the polynomial basis functions\n", + " % must be constructed; n_rows-by-dimen;\n", + " % \"D\" is the degree of the polynomial whose basis functions must\n", + " % be constructed; (can be 1,2,3,4 or 5)\n", + " %\n", + " % Output: \"basis_fs\" is the matrix of basis functions of a complete\n", + " % polynomial of the given degree\n", + " % -------------------------------------------------------------------------\n", + " % Copyright � 2011 by Lilia Maliar and Serguei Maliar. All rights reserved.\n", + " % The code may be used, modified and redistributed under the terms provided\n", + " % in the file \"License_Agreement.pdf\".\n", + " % -------------------------------------------------------------------------\n", + " \n", + " function basis_fs = Ord_Polynomial_N(z,D)\n", + " \n", + " \n", + " % A polynomial is given by the sum of polynomial basis functions, phi(i),\n", + " % multiplied by the coefficients; see condition (13) in JMM (2011). By\n", + " % convention, the first basis function is one.\n", + " \n", + " [n_rows,dimen] = size(z); % Compute the number of rows, n_rows, and the\n", + " % number of variables (columns), dimen, in the\n", + " % data z on which the polynomial basis functions\n", + " % must be constructed\n", + " \n", + " % 1. The matrix of the basis functions of the first-degree polynomial\n", + " % (the default option)\n", + " % -------------------------------------------------------------------\n", + " basis_fs = [ones(n_rows,1) z]; % The matrix includes a column of ones\n", + " % (the first basis function is one for\n", + " % n_rows points) and linear basis\n", + " % functions\n", + " i = dimen+1; % Index number of a polynomial basis function; the first\n", + " % basis function (equal to one) and linear basis functions\n", + " % are indexed from 1 to dimen+1, and subsequent polynomial\n", + " % basis functions will be indexed from dimen+2 and on\n", + " \n", + " % 2. The matrix of the basis functions of the second-degree polynomial\n", + " % --------------------------------------------------------------------\n", + " if D == 2\n", + " \n", + " % Version one (not vectorized):\n", + " for j1 = 1:dimen\n", + " for j2 = j1:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2);\n", + " end\n", + " end\n", + " \n", + " % Version 2 (vectorized): Note that this version works only for a second-degree\n", + " % polynomial in which all state variables take non-zero values\n", + " % for r = 1:n_rows\n", + " % basis_fs(r,2+dimen:1+dimen+dimen*(dimen+1)/2) = [nonzeros(tril(z(r,:)'*z(r,:)))'];\n", + " % Compute linear and quadratic polynomial basis functions for\n", + " % each row r; \"tril\" extracts a lower triangular part of z'z\n", + " % (note that matrix z'z is symmetric so that an upper triangular\n", + " % part is redundant); \"nonzeros\" forms a column vector by\n", + " % stacking the columns of the original matrix one after another\n", + " % and by eliminating zero terms\n", + " % end\n", + " \n", + " % 3. The matrix of the basis functions of the third-degree polynomial\n", + " % -------------------------------------------------------------------\n", + " elseif D == 3\n", + " \n", + " for j1 = 1:dimen\n", + " for j2 = j1:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2);\n", + " for j3 = j2:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3);\n", + " end\n", + " end\n", + " end\n", + " \n", + " % 4. The matrix of the basis functions of the fourth-degree polynomial\n", + " % -------------------------------------------------------------------\n", + " elseif D == 4\n", + " \n", + " for j1 = 1:dimen\n", + " for j2 = j1:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2);\n", + " for j3 = j2:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3);\n", + " for j4 = j3:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3).*z(:,j4);\n", + " end\n", + " end\n", + " end\n", + " end\n", + " \n", + " % 5. The matrix of the basis functions of the fifth-degree polynomial\n", + " % -------------------------------------------------------------------\n", + " \n", + " elseif D == 5\n", + " \n", + " for j1 = 1:dimen\n", + " for j2 = j1:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2);\n", + " for j3 = j2:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3);\n", + " for j4 = j3:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3).*z(:,j4);\n", + " for j5 = j4:dimen\n", + " i = i+1;\n", + " basis_fs(:,i) = z(:,j1).*z(:,j2).*z(:,j3).*z(:,j4).*z(:,j5);\n", + " end\n", + " end\n", + " end\n", + " end\n", + " end\n", + " \n", + " end % end of \"if\"/\"elseif\"\n", + " '\n", + "\n" + ] + } + ], + "source": [ + "poly_txt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Accuracy and Simulation\n", + "\n", + "Before proceeding to the main code, we will also need to have routines for simulating and checking the accuracy of our solution to the model.\n", + "\n", + "We also download these from online as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "nkacc_txt = fileread(websave('NK_accuracy.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/NK_accuracy.m'));\n", + "nksim_txt = fileread(websave('NK_simulation.m', 'https://s3.amazonaws.com/clmm-resources/mfiles/NK_simulation.m'));" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "nkacc_txt =\n", + "\n", + " '% NK_accuracy.m is a routine for evaluating accuracy of the solutions \n", + " % to the new Keynesian model considerd in the article \"Merging Simulation\n", + " % and Projection Approaches to Solve High-Dimensional Problems with an\n", + " % application to a New Keynesian Model\" by Lilia Maliar and Serguei Maliar,\n", + " % Quantitative Economics 6, 1-47 (2015) (henceforth, MM, 2015).\n", + " % -------------------------------------------------------------------------\n", + " % Inputs: \"nua\", \"nuL\", \"nuR\", \"nuG\", \"nuB\", \"nuu\" are the time series\n", + " % of shocks generated with the innovations for the test;\n", + " % \"delta\", \"L\", \"Y\", \"Yn\", \"pie\", \"S\", \"F\", \"C\" are the time\n", + " % series solution generated with the innovations for test;\n", + " % \"rho_nua\", \"rho_nuL\", \"rho_nuR\", \"rho_nuu\", \"rho_nuB\", \"rho_nuG\"\n", + " % are the parameters of the laws of motion for shocks;\n", + " % \"mu\", \"gam\", \"epsil\", \"vartheta\", \"beta\", \"A\", \"tau\", \"rho\",\n", + " % \"vcv\", \"beta\", \"phi_y\", \"phi_pie\", \"theta\", \"piestar\" and \"Gbar\"\n", + " % are the parameters of the model;\n", + " % \"vk_2d\" is the matrix of coefficients of the GSSA solution;\n", + " % \"discard\" is the number of data points to discard;\n", + " % \"Degree\" is the degree of polynomial approximation\n", + " %\n", + " % Output: \"Residuals_mean\" and \"Residuals_max\" are, respectively, the mean\n", + " % and maximum absolute residuals across all points and all equi-\n", + " % librium conditions; \"Residuals_max_E\" is the maximum absolute\n", + " % residuals across all points, disaggregated by optimality conditions;\n", + " % \"Residuals\" are absolute residuals disaggregated by the equi-\n", + " % librium conditions\n", + " % -------------------------------------------------------------------------\n", + " % Copyright � 2013-2016 by Lilia Maliar and Serguei Maliar. All rights reserved.\n", + " % The code may be used, modified and redistributed under the terms provided\n", + " % in the file \"License_Agreement.txt\".\n", + " % -------------------------------------------------------------------------\n", + " \n", + " function [Residuals_mean, Residuals_max, Residuals_max_E, Residuals] = NK_accuracy(nua,nuL,nuR,nuG,nuB,nuu,R,delta,L,Y,Yn,pie,S,F,C,rho_nua,rho_nuL,rho_nuR,rho_nuu,rho_nuB,rho_nuG,gam,vartheta,epsil,beta,phi_y,phi_pie,mu,theta,piestar,vcv,discard,vk_2d,Gbar,zlb,Degree)\n", + " \n", + " tic % Start counting time for running the test\n", + " \n", + " [T] = size(nua,1); % Infer the number of points on which accuracy is\n", + " % evaluated\n", + " Residuals = zeros(T,6); % Allocate memory to the matrix of residuals; T-by-6\n", + " \n", + " % Integration method for evaluating accuracy\n", + " % ------------------------------------------\n", + " [n_nodes,epsi_nodes,weight_nodes] = Monomials_2(6,vcv);\n", + " % Monomial integration rule with 2N^2+1 nodes\n", + " \n", + " % Compute variables on the given set of points\n", + " %---------------------------------------------\n", + " \n", + " for t = 1:T; % For each given point,\n", + " t;\n", + " \n", + " % Take the corresponding value for shocks at t\n", + " %---------------------------------------------\n", + " nuR0 = nuR(t,1); % nuR(t)\n", + " nua0 = nua(t,1); % nua(t)\n", + " nuL0 = nuL(t,1); % nuL(t)\n", + " nuu0 = nuu(t,1); % nuu(t)\n", + " nuB0 = nuB(t,1); % nuB(t)\n", + " nuG0 = nuG(t,1); % nuG(t)\n", + " \n", + " % Compute shocks at t+1 in all future nodes using their laws of motion\n", + " %---------------------------------------------------------------------\n", + " % Note that we do not premultiply by standard deviations as epsi_nodes\n", + " % already include them\n", + " nuR1(1:n_nodes,1) = (ones(n_nodes,1)*nuR0)*rho_nuR + epsi_nodes(:,1);\n", + " % nuR(t+1); n_nodes-by-1\n", + " nua1(1:n_nodes,1) = (ones(n_nodes,1)*nua0)*rho_nua + epsi_nodes(:,2);\n", + " % nua(t+1); n_nodes-by-1\n", + " nuL1(1:n_nodes,1) = (ones(n_nodes,1)*nuL0)*rho_nuL + epsi_nodes(:,3);\n", + " % nuL(t+1); n_nodes-by-1\n", + " nuu1(1:n_nodes,1) = (ones(n_nodes,1)*nuu0)*rho_nuu + epsi_nodes(:,4);\n", + " % nuu(t+1); n_nodes-by-1\n", + " nuB1(1:n_nodes,1) = (ones(n_nodes,1)*nuB0)*rho_nuB + epsi_nodes(:,5);\n", + " % nuB(t+1); n_nodes-by-1\n", + " nuG1(1:n_nodes,1) = (ones(n_nodes,1)*nuG0)*rho_nuG + epsi_nodes(:,6);\n", + " % nuG(t+1); n_nodes-by-1\n", + " \n", + " R0 = R(t,1); % R(t-1)\n", + " delta0 = delta(t,1); % delta(t-1)\n", + " R1 = R(t+1,1); % R(t)\n", + " delta1 = delta(t+1,1); % delta(t)\n", + " \n", + " L0 = L(t,1); % L(t)\n", + " Y0 = Y(t,1); % Y(t)\n", + " Yn0 = Yn(t,1); % Yn(t)\n", + " pie0 = pie(t,1); % pie(t)\n", + " S0 = S(t,1); % S(t)\n", + " F0 = F(t,1); % F(t)\n", + " C0 = C(t,1); % C(t)\n", + " \n", + " % Future choices at t+1\n", + " %----------------------\n", + " delta1_dupl = ones(n_nodes,1)*delta1;\n", + " R1_dupl = ones(n_nodes,1)*R1;\n", + " % Duplicate \"delta1\" and \"R1\" n_nodes times to create a matrix with\n", + " % n_nodes identical rows; n_nodes-by-1\n", + " \n", + " X1 = Ord_Polynomial_N([log(R1_dupl) log(delta1_dupl) nuR1 nua1 nuL1 nuu1 nuB1 nuG1],Degree);\n", + " % Form a complete polynomial of degree \"Degree\" (at t+1) on future state\n", + " % variables; n_nodes-by-npol_2d\n", + " \n", + " S1 = X1*vk_2d(:,1); % Compute S(t+1) in all nodes using vk_2d\n", + " F1 = X1*vk_2d(:,2); % Compute F(t+1) in all nodes using vk_2d\n", + " C1 = (X1*vk_2d(:,3)).^(-1/gam); % Compute C(t+1) in all nodes using vk_2d\n", + " pie1 = ((1-(1-theta)*(S1./F1).^(1-epsil))/theta).^(1/(epsil-1));\n", + " % Compute pie(t+1) using condition (35)\n", + " % in MM (2015)\n", + " \n", + " % Compute residuals for each of the 9 equilibrium conditions\n", + " %-----------------------------------------------------------\n", + " Residuals(t,1) = 1-weight_nodes'*(exp(nuu0)*exp(nuL0)*L0^vartheta*Y0/exp(nua0) + beta*theta*pie1.^epsil.*S1)/S0;\n", + " Residuals(t,2) = 1-weight_nodes'*(exp(nuu0)*C0^(-gam)*Y0 + beta*theta*pie1.^(epsil-1).*F1)/F0;\n", + " Residuals(t,3) = 1-weight_nodes'*(beta*exp(nuB0)/exp(nuu0)*R1*exp(nuu1).*C1.^(-gam)./pie1)/C0^(-gam);\n", + " Residuals(t,4) = 1-((1-theta*pie0^(epsil-1))/(1-theta))^(1/(1-epsil))*F0/S0;\n", + " Residuals(t,5) = 1-((1-theta)*((1-theta*pie0^(epsil-1))/(1-theta))^(epsil/(epsil-1)) + theta*pie0^epsil/delta0)^(-1)/delta1;\n", + " Residuals(t,6) = 1-exp(nua0)*L0*delta1/Y0;\n", + " Residuals(t,7) = 1-(1-Gbar/exp(nuG0))*Y0/C0;\n", + " Residuals(t,8) = 1-(exp(nua0)^(1+vartheta)*(1-Gbar/exp(nuG0))^(-gam)/exp(nuL0))^(1/(vartheta+gam))/Yn0;\n", + " Residuals(t,9) = 1-piestar/beta*(R0*beta/piestar)^mu*((pie0/piestar)^phi_pie * (Y0/Yn0)^phi_y)^(1-mu)*exp(nuR0)/R1; % Taylor rule\n", + " if (zlb==1); Residuals(t,9) = Residuals(t,9)*(R1>1);end\n", + " % If the ZLB is imposed and R>1, the residuals in the Taylor rule (the\n", + " % 9th equation) are zero\n", + " \n", + " end\n", + " \n", + " % Residuals across all the equilibrium conditions and all test points\n", + " %--------------------------------------------------------------------\n", + " Residuals_mean = log10(mean(mean(abs(Residuals(1+discard:end,:)))));\n", + " % Mean absolute residuals computed after discarding the first\n", + " % \"discard\" observations\n", + " \n", + " Residuals_max = log10(max(max(abs(Residuals(1+discard:end,:)))));\n", + " % Maximum absolute residuals computed after discarding the first\n", + " % \"discard\" observations\n", + " \n", + " % Residuals disaggregated by the eqiulibrium conditions\n", + " %------------------------------------------------------\n", + " Residuals_max_E = log10(max(abs(Residuals(1+discard:end,:))))';\n", + " % Maximum absolute residuals across all test points for each of the 9\n", + " % equilibrium conditions computed after discarding the first \"discard\"\n", + " % observations;\n", + " \n", + " time_test = toc; % Time needed to run the test\n", + " '\n", + "\n" + ] + } + ], + "source": [ + "nkacc_txt" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "nksim_txt =\n", + "\n", + " '% NK_simulation.m is a routine for simulating the solution to the new \n", + " % Keynesian model considerd in the article \"Merging Simulation and Projection\n", + " % Approaches to Solve High-Dimensional Problems with an application to a New\n", + " % Keynesian Model\" by Lilia Maliar and Serguei Maliar, Quantitative Economics\n", + " % 6, 1-47 (2015) (henceforth, MM, 2015).\n", + " % -------------------------------------------------------------------------\n", + " \n", + " % -------------------------------------------------------------------------\n", + " % Inputs: \"vk\" is the matrix of coefficients of the GSSA solution;\n", + " % \"nua\", \"nuL\", \"nuR\", \"nuG\", \"nuB\", \"nuu\" are the time series\n", + " % of shocks generated with the innovations for the test;\n", + " % \"R_init\" and \"delta_init\" are initial values for R and delta\n", + " % \"gam\", \"vartheta\", \"epsil\", \"betta\", \"phi_y\", \"phi_pie\", \"mu\",\n", + " % \"theta\", \"piestar\", \"Gbar\" are the parameters of the model;\n", + " % \"zlb\" is a dummy parameter which is equal to 0 when ZLB is not\n", + " % imposed and is equal to 1 when it is imposed;\n", + " % \"Degree\" is the degree of polynomial approximation\n", + " %\n", + " % Output: \"S\", \"F\", \"delta\", \"C\", \"Y\", \"Yn\", \"L\", \"R\" and \"pie2\" are the\n", + " % simulated series\n", + " % -------------------------------------------------------------------------\n", + " % Copyright � 2013-2016 by Lilia Maliar and Serguei Maliar. All rights reserved.\n", + " % The code may be used, modified and redistributed under the terms provided\n", + " % in the file \"License_Agreement.txt\".\n", + " % -------------------------------------------------------------------------\n", + " \n", + " \n", + " function [S F delta C Y Yn L R pie w] = NK_simulation(vk,nuR,nua,nuL,nuu,nuB,nuG,R_init,delta_init,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,Gbar,zlb,Degree)\n", + " \n", + " \n", + " [T] = size(nua,1); % Infer the number of points on which the accuracy\n", + " % is evaluated\n", + " \n", + " delta = ones(T+1,1); % Allocate memory for the time series of delta(t)\n", + " R = ones(T+1,1); % Allocate memory for the time series of R(t)\n", + " S = ones(T,1); % Allocate memory for the time series of S(t)\n", + " F = ones(T,1); % Allocate memory for the time series of F(t)\n", + " C = ones(T,1); % Allocate memory for the time series of C(t)\n", + " pie = ones(T,1); % Allocate memory for the time series of pie(t)\n", + " Y = ones(T,1); % Allocate memory for the time series of Y(t)\n", + " L = ones(T,1); % Allocate memory for the time series of L(t)\n", + " Yn = ones(T,1); % Allocate memory for the time series of Yn(t)\n", + " w = ones(T,1);\n", + " state = ones(T,8); % Allocate memory for the time series of the state\n", + " % variables; T-by-8\n", + " \n", + " delta(1,1) = delta_init; % Initial condition for delta(t-1), i.e., delta(-1)\n", + " R(1,1) = R_init; % Initial condition for R(t-1)\n", + " \n", + " for t = 1:T\n", + " pol_bases = Ord_Polynomial_N([log(R(t,1)) log(delta(t,1)) nuR(t,1) nua(t,1) nuL(t,1) nuu(t,1) nuB(t,1) nuG(t,1)],Degree);\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 \"Degree\"\n", + " S(t,1) = pol_bases*vk(:,1); % Compute S(t) using vk\n", + " F(t,1) = pol_bases*vk(:,2); % Compute F(t) using vk\n", + " C(t,1) = (pol_bases*vk(:,3)).^(-1/gam); % Compute C(t) using vk\n", + " pie(t,:) = ((1-(1-theta)*(S(t,:)/F(t,:))^(1-epsil))/theta)^(1/(epsil-1));\n", + " % Compute pie(t) from condition (35) in MM (2015)\n", + " delta(t+1,:) = ((1-theta)*((1-theta*pie(t,1)^(epsil-1))/(1-theta))^(epsil/(epsil-1))+theta*pie(t,1)^epsil/delta(t,1))^-1;\n", + " % Compute delta(t) from condition (36) in MM (2015)\n", + " Y(t,:) = C(t,1)/(1-Gbar/exp(nuG(t,1)));\n", + " % Compute Y(t) from condition (38) in MM (2015)\n", + " L(t,1) = Y(t,1)/delta(t+1,1)/exp(nua(t,1));\n", + " % Compute L(t) from condition (37) in MM (2015)\n", + " Yn(t,:) = (exp(nua(t,1))^(1+vartheta)*(1-Gbar/exp(nuG(t,1)))^(-gam)/exp(nuL(t,1)))^(1/(vartheta+gam));\n", + " % Compute Yn(t) from condition (31) in MM (2015)\n", + " R(t+1,:) = piestar/betta*(R(t,:)*betta/piestar)^mu*((pie(t,:)/piestar)^phi_pie*(Y(t,:)/Yn(t,:))^phi_y)^(1-mu)*exp(nuR(t,1));\n", + " % Compute R(t) from conditions (27), (39) in MM (2015)\n", + " w(t,1) = exp(nuL(t,1))*(L(t,1)^vartheta)*(C(t,1)^gam);\n", + " % Compute real wage\n", + " \n", + " if zlb == 1\n", + " R(t+1,:) = max(R(t+1,:),1);\n", + " % If ZLB is imposed, set R(t)=1 if ZLB binds\n", + " end;\n", + " \n", + " end\n", + " '\n", + "\n" + ] + } + ], + "source": [ + "nksim_txt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Main routine" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% Start counting time\n", + "% -------------------\n", + "tic;" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 1. Parameter values\n", + "% -------------------\n", + "zlb = 0; % Impose ZLB on nominal interest rate\n", + "gam = 1; % Utility-function parameter\n", + "betta = 0.99; % Discount factor\n", + "vartheta = 2.09; % Utility-function parameter\n", + "epsil = 4.45; % Parameter in the Dixit-Stiglitz aggregator\n", + "phi_y = 0.07; % Parameter of the Taylor rule\n", + "phi_pie = 2.21; % Parameter of the Taylor rule\n", + "mu = 0.82; % Parameter of the Taylor rule\n", + "theta = 0.83; % Share of non-reoptimizing firms (Calvo's pricing)\n", + "piestar = 1.0; % Target (gross) inflation rate\n", + "Gbar = 0.23; % Steady-state share of government spending in output\n", + "\n", + "% Autocorrelation coefficients in the processes for shocks\n", + "%----------------------------------------------------------\n", + "rho_nua = 0.95; % See process (22) in MM (2015)\n", + "rho_nuL = 0.25; % See process (16) in MM (2015)\n", + "rho_nuR = 0.0; % See process (28) in MM (2015)\n", + "rho_nuu = 0.92; % See process (15) in MM (2015)\n", + "rho_nuB = 0.0; % See process (17) in MM (2015)\n", + "rho_nuG = 0.95; % See process (26) in MM (2015)\n", + "\n", + "\n", + "% Standard deviations of the innovations in the processes for shocks\n", + "%--------------------------------------------------------------------\n", + "sigma_nua = 0.0045; % See process (22) in MM (2015)\n", + "sigma_nuL = 0.4054; % See process (16) in MM (2015)\n", + "sigma_nuR = 0.0028; % See process (28) in MM (2015)\n", + "sigma_nuu = 0.0054; % See process (15) in MM (2015)\n", + "sigma_nuB = 0.0010; % See process (17) in MM (2015)\n", + "sigma_nuG = 0.0038; % See process (26) in MM (2015)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 2. Steady state of the model (see page 19 in Supplement to MM (2015))\n", + "% -------------------------------------------------------------------\n", + "Yn_ss = exp(Gbar)^(gam/(vartheta+gam));\n", + "Y_ss = Yn_ss;\n", + "pie_ss = 1;\n", + "delta_ss = 1;\n", + "L_ss = Y_ss/delta_ss;\n", + "C_ss = (1-Gbar)*Y_ss;\n", + "F_ss = C_ss^(-gam)*Y_ss/(1-betta*theta*pie_ss^(epsil-1));\n", + "S_ss = L_ss^vartheta*Y_ss/(1-betta*theta*pie_ss^epsil);\n", + "R_ss = pie_ss/betta;\n", + "w_ss = (L_ss^vartheta)*(C_ss^gam);" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 3. Construct a grid for computing a solution\n", + "% --------------------------------------------\n", + "m = 200; % Choose the number of grid points\n", + "grid_type = 1; % Choose a grid type; grid_type = 1 corresponds to a uniformely\n", + " % distribued random grid, and grid_type = 2 corresponds to\n", + " % a quasi Monte Carlo grid\n", + "\n", + "% Random grid\n", + "% -----------\n", + "if grid_type == 1 % Choose the grid type\n", + "nuR0 = (-2*sigma_nuR+4*sigma_nuR*rand(m,1))/sqrt(1-rho_nuR^2);\n", + "nua0 = (-2*sigma_nua+4*sigma_nua*rand(m,1))/sqrt(1-rho_nua^2);\n", + "nuL0 = (-2*sigma_nuL+4*sigma_nuL*rand(m,1))/sqrt(1-rho_nuL^2);\n", + "nuu0 = (-2*sigma_nuu+4*sigma_nuu*rand(m,1))/sqrt(1-rho_nuu^2);\n", + "nuB0 = (-2*sigma_nuB+4*sigma_nuB*rand(m,1))/sqrt(1-rho_nuB^2);\n", + "nuG0 = (-2*sigma_nuG+4*sigma_nuG*rand(m,1))/sqrt(1-rho_nuG^2);\n", + " % Values of exogenous state variables are distributed uniformly\n", + " % in the interval +/- std/sqrt(1-rho_nu^2)\n", + "\n", + "R0 = 1+0.05*rand(m,1);\n", + "delta0 = 0.95+0.05*rand(m,1);\n", + " % Values of endogenous state variables are distributed uniformly\n", + " % in the intervals [1 1.05] and [0.95 1], respectively\n", + "end\n", + "\n", + "% Quasi Monte-Carlo grid\n", + "%-----------------------\n", + "if grid_type == 2 % Choose the grid type\n", + "dimensionality = 8; % The number of state variables (exogenous and endogenous)\n", + "Def_sobol = sobolset(dimensionality);\n", + " % Constructs a Sobol sequence point set in\n", + " % \"dimensionality\" dimensions\n", + "\n", + "Sob = net(Def_sobol,m);\n", + " % Get the first m points\n", + "\n", + "nuR0 = (-2*sigma_nuR+4*(max(Sob(:,1))-Sob(:,1))/(max(Sob(:,1))-min(Sob(:,1)))*sigma_nuR)/sqrt(1-rho_nuR^2);\n", + "nua0 = (-2*sigma_nua+4*(max(Sob(:,2))-Sob(:,2))/(max(Sob(:,2))-min(Sob(:,2)))*sigma_nua)/sqrt(1-rho_nua^2);\n", + "nuL0 = (-2*sigma_nuL+4*(max(Sob(:,3))-Sob(:,3))/(max(Sob(:,3))-min(Sob(:,3)))*sigma_nuL)/sqrt(1-rho_nuL^2);\n", + "nuu0 = (-2*sigma_nuu+4*(max(Sob(:,4))-Sob(:,4))/(max(Sob(:,4))-min(Sob(:,4)))*sigma_nuu)/sqrt(1-rho_nuu^2);\n", + "nuB0 = (-2*sigma_nuB+4*(max(Sob(:,5))-Sob(:,5))/(max(Sob(:,5))-min(Sob(:,5)))*sigma_nuB)/sqrt(1-rho_nuB^2);\n", + "nuG0 = (-2*sigma_nuG+4*(max(Sob(:,6))-Sob(:,6))/(max(Sob(:,6))-min(Sob(:,6)))*sigma_nuG)/sqrt(1-rho_nuG^2);\n", + " % Values of exogenous state variables are in the interval +/- std/sqrt(1-rho^2)\n", + "\n", + "R0 = 1+0.05*(max(Sob(:,7))-Sob(:,7))/(max(Sob(:,7))-min(Sob(:,7)));\n", + "delta0 = 0.95+0.05*(max(Sob(:,8))-Sob(:,8))/(max(Sob(:,8))-min(Sob(:,8)));\n", + " % Values of endogenous state variables are in the intervals [1 1.05] and\n", + " % [0.95 1], respectively\n", + "end\n", + "\n", + "if zlb == 1; R0=max(R0,1); end\n", + " % If ZLB is imposed, set R(t)=1 if ZLB binds\n", + "\n", + "Grid = [log(R0(1:m,1)) log(delta0(1:m,1)) nuR0 nua0 nuL0 nuu0 nuB0 nuG0];\n", + " % Construct the matrix of grid points; m-by-dimensionality" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 4. Constructing polynomial on the grid\n", + "% --------------------------------------\n", + "Degree = 2; % Degree of polynomial approximation\n", + "X0_Gs{1} = Ord_Polynomial_N(Grid, 1);\n", + "X0_Gs{Degree} = Ord_Polynomial_N(Grid, Degree);\n", + " % Construct the matrix of explanatory variables X0_Gs\n", + " % on the grid of state variables; the columns of X0_Gs\n", + " % are given by the basis functions of polynomial of\n", + " % degree \"Degree\"\n", + "npol = size(X0_Gs{Degree},2); % Number of coefficients in polynomial of degree\n", + " % \"Degree\"; it must be smaller than the number of grid\n", + " % points" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 5. Integration in the solution procedure\n", + "% ----------------------------------------\n", + "N = 6; % Total number of exogenous shocks\n", + "vcv = diag([sigma_nuR^2 sigma_nua^2 sigma_nuL^2 sigma_nuu^2 sigma_nuB^2 sigma_nuG^2]);\n", + " % Variance covariance matrix\n", + "\n", + "% Compute the number of integration nodes, their values and weights\n", + "%------------------------------------------------------------------\n", + "[n_nodes,epsi_nodes,weight_nodes] = Monomials_1(N,vcv);\n", + " % Monomial integration rule with 2N nodes\n", + "%[n_nodes,epsi_nodes,weight_nodes] = Monomials_2(N,vcv);\n", + " % Monomial integration rule with 2N^2+1 nodes\n", + "\n", + "nuR1(:,1:n_nodes) = (nuR0*ones(1,n_nodes)).*rho_nuR + ones(m,1)*epsi_nodes(:,1)';\n", + "nua1(:,1:n_nodes) = (nua0*ones(1,n_nodes)).*rho_nua + ones(m,1)*epsi_nodes(:,2)';\n", + "nuL1(:,1:n_nodes) = (nuL0*ones(1,n_nodes)).*rho_nuL + ones(m,1)*epsi_nodes(:,3)';\n", + "nuu1(:,1:n_nodes) = (nuu0*ones(1,n_nodes)).*rho_nuu + ones(m,1)*epsi_nodes(:,4)';\n", + "nuB1(:,1:n_nodes) = (nuB0*ones(1,n_nodes)).*rho_nuB + ones(m,1)*epsi_nodes(:,5)';\n", + "nuG1(:,1:n_nodes) = (nuG0*ones(1,n_nodes)).*rho_nuG + ones(m,1)*epsi_nodes(:,6)';\n", + " % Compute future shocks in all grid points and all integration\n", + " % nodes; the size of each of these matrices is m-by-n_nodes" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 6. Allocate memory\n", + "%--------------------\n", + "% TODO: pick up here!!!!\n", + "e = zeros(m,3);\n", + "% Allocate memory to integrals in the right side of 3 Euler equations\n", + "\n", + "S0_old_G = ones(m,1);\n", + "F0_old_G = ones(m,1);\n", + "C0_old_G = ones(m,1);\n", + "% Allocate memory to S, F and C from the previous iteration (to check\n", + "% convergence)\n", + "\n", + "S0_new_G = ones(m,1);\n", + "F0_new_G = ones(m,1);\n", + "C0_new_G = ones(m,1);\n", + "% Allocate memory to S, F, C from the current iteration (to check\n", + "% convergence)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%--------------------------------------------------------------------------\n", + "%\n", + "% The main iterative cycle\n", + "%\n", + "% -------------------------------------------------------------------------\n", + "\n", + "% 7. Algorithm parameters\n", + "%------------------------\n", + "damp = 0.1; % Damping parameter for (fixed-point) iteration on\n", + " % the coefficients of 3 decision functions (for\n", + " % S, F and C^(-gam))" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 8. The loop over the polynomial coefficients\n", + "% --------------------------------------------\n", + "for deg = [1, Degree]\n", + " diff = 1e+10;\n", + " it = 0;\n", + " X0_G = X0_Gs{deg};\n", + "\n", + " % 9. Initial guess for coefficients of the decision functions for the\n", + " % variables S and F and marginal utility MU\n", + " % -------------------------------------------------------------------\n", + " if deg == 1\n", + " vk = ones(size(Grid, 2)+1, 3)*1e-5; % Initialize first all the coefficients\n", + " % at 1e-5\n", + " vk(1,:) = [S_ss F_ss C_ss.^(-gam)]; % Set the initial values of the constant\n", + " % terms in the decision rules for S,\n", + " % F and MU to values that give the\n", + " % deterministic steady state\n", + " else\n", + " % For degree > 1, initial guess for coefficients is given by\n", + " % regressing final state matrix from degree 1 solution (e) on the\n", + " % complete polynomial basis matrix\n", + " vk = X0_G\\e;\n", + " end\n", + "\n", + " while diff > 1e-7 % The convergence criterion (which is unit free\n", + " % because diff is unit free)\n", + "\n", + " % Current choices (at t)\n", + " % ------------------------------\n", + " S0 = X0_G*vk(:,1); % Compute S(t) using vk\n", + " F0 = X0_G*vk(:,2); % Compute F(t) using vk\n", + " C0 = (X0_G*vk(:,3)).^(-1/gam); % Compute C(t) using vk\n", + "\n", + " pie0 = ((1-(1-theta)*(S0./F0).^(1-epsil))/theta).^(1/(epsil-1));\n", + " % Compute pie(t) from condition (35) in MM (2015)\n", + " delta1 = ((1-theta)*((1-theta*pie0.^(epsil-1))/(1-theta)).^(epsil/(epsil-1))+theta*pie0.^epsil./delta0).^(-1);\n", + " % Compute delta(t) from condition (36) in MM (2015)\n", + " Y0 = C0./(1-Gbar./exp(nuG0));\n", + " % Compute Y(t) from condition (38) in MM (2015)\n", + " L0 = Y0./exp(nua0)./delta1;\n", + " % Compute L(t) from condition (37) in MM (2015)\n", + " Yn0 = (exp(nua0).^(1+vartheta).*(1-Gbar./exp(nuG0)).^(-gam)./exp(nuL0)).^(1/(vartheta+gam));\n", + " % Compute Yn(t) from condition (31) in MM (2015)\n", + " R1 = piestar/betta*(R0*betta./piestar).^mu.*((pie0./piestar).^phi_pie .* (Y0./Yn0).^phi_y).^(1-mu).*exp(nuR0); % Taylor rule\n", + " % Compute R(t) from conditions (27), (39) in MM (2015)\n", + " if zlb == 1; R1=max(R1,1); end\n", + " % If ZLB is imposed, set R(t)=1 if ZLB binds\n", + "\n", + " % Future choices (at t+1)\n", + " %--------------------------------\n", + " for u = 1:n_nodes\n", + "\n", + " X1 = Ord_Polynomial_N([log(R1) log(delta1) nuR1(:,u) nua1(:,u) nuL1(:,u) nuu1(:,u) nuB1(:,u) nuG1(:,u)],deg);\n", + " % Form complete polynomial of degree \"Degree\" (at t+1) on future state\n", + " % variables; n_nodes-by-npol\n", + "\n", + " S1(:,u) = X1*vk(:,1); % Compute S(t+1) in all nodes using vk\n", + " F1(:,u) = X1*vk(:,2); % Compute F(t+1) in all nodes using vk\n", + " C1(:,u) = (X1*vk(:,3)).^(-1/gam); % Compute C(t+1) in all nodes using vk\n", + "\n", + " end\n", + "\n", + " pie1 = ((1-(1-theta)*(S1./F1).^(1-epsil))/theta).^(1/(epsil-1));\n", + " % Compute next-period pie using condition\n", + " % (35) in MM (2015)\n", + "\n", + "\n", + " % Evaluate conditional expectations in the Euler equations\n", + " %---------------------------------------------------------\n", + " e(:,1) = exp(nuu0).*exp(nuL0).*L0.^vartheta.*Y0./exp(nua0) + (betta*theta*pie1.^epsil.*S1)*weight_nodes;\n", + " e(:,2) = exp(nuu0).*C0.^(-gam).*Y0 + (betta*theta*pie1.^(epsil-1).*F1)*weight_nodes;\n", + " e(:,3) = betta*exp(nuB0)./exp(nuu0).*R1.*((exp(nuu1).*C1.^(-gam)./pie1)*weight_nodes);\n", + "\n", + "\n", + " % Variables of the current iteration\n", + " %-----------------------------------\n", + " S0_new_G(:,1) = S0(:,1);\n", + " F0_new_G(:,1) = F0(:,1);\n", + " C0_new_G(:,1) = C0(:,1);\n", + "\n", + " % Compute and update the coefficients of the decision functions\n", + " % -------------------------------------------------------------\n", + " vk_hat_2d = X0_G\\e; % Compute the new coefficients of the decision\n", + " % functions using a backslash operator\n", + "\n", + " vk = damp*vk_hat_2d + (1-damp)*vk;\n", + " % Update the coefficients using damping\n", + "\n", + " % Evaluate the percentage (unit-free) difference between the values\n", + " % on the grid from the previous and current iterations\n", + " % -----------------------------------------------------------------\n", + " diff = mean(mean(abs(1-S0_new_G./S0_old_G)))+mean(mean(abs(1-F0_new_G./F0_old_G)))+mean(mean(abs(1-C0_new_G./C0_old_G)));\n", + " % The convergence criterion is adjusted to the damping\n", + " % parameters\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", + " S0_old_G = S0_new_G;\n", + " F0_old_G = F0_new_G;\n", + " C0_old_G = C0_new_G; \n", + " end\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 10. Finish counting time\n", + "% ------------------------\n", + "running_time = toc;" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "% 11. Simualating a time-series solution\n", + "%---------------------------------------\n", + "T = 10201; % The length of stochastic simulation\n", + "\n", + "\n", + "% Initialize the values of 6 exogenous shocks and draw innovations\n", + "%-----------------------------------------------------------------\n", + "nuR = zeros(T,1); eps_nuR = randn(T,1)*sigma_nuR;\n", + "nua = zeros(T,1); eps_nua = randn(T,1)*sigma_nua;\n", + "nuL = zeros(T,1); eps_nuL = randn(T,1)*sigma_nuL;\n", + "nuu = zeros(T,1); eps_nuu = randn(T,1)*sigma_nuu;\n", + "nuB = zeros(T,1); eps_nuB = randn(T,1)*sigma_nuB;\n", + "nuG = zeros(T,1); eps_nuG = randn(T,1)*sigma_nuG;\n", + "\n", + "% Generate the series for shocks\n", + "%-------------------------------\n", + "for t = 1:T-1\n", + " nuR(t+1,1) = rho_nuR*nuR(t,1) + eps_nuR(t);\n", + " nua(t+1,1) = rho_nua*nua(t,1) + eps_nua(t);\n", + " nuL(t+1,1) = rho_nuL*nuL(t,1) + eps_nuL(t);\n", + " nuu(t+1,1) = rho_nuu*nuu(t,1) + eps_nuu(t);\n", + " nuB(t+1,1) = rho_nuB*nuB(t,1) + eps_nuB(t);\n", + " nuG(t+1,1) = rho_nuG*nuG(t,1) + eps_nuG(t) ;\n", + "end\n", + "\n", + "% Initial values of two endogenous state variables\n", + "%-------------------------------------------------\n", + "R_initial = 1; % Nominal interest rate R\n", + "delta_initial = 1; % Price dispersion \"delta\"\n", + "\n", + "tic;\n", + "% Simulate the model\n", + "%-------------------\n", + "[S, F, delta, C, Y, Yn, L, R, pie, w] = NK_simulation(vk,nuR,nua,nuL,nuu,nuB,nuG,R_initial,delta_initial,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,Gbar,zlb,Degree);\n", + "\n", + "simulation_time = toc;" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Solver time: 3.988508 seconds\n", + "Simulation time: 0.194495 seconds\n", + "Accuracy time: 0.922081 seconds\n", + "Total time: 5.105085 seconds\n", + "pistar: 1.0000 sigma_L: 0.4054\n", + "zlb : 0 grid: 1\n", + "tex line: -1.08 & -1.27 & 5.11\n" + ] + } + ], + "source": [ + "% 12. Compute unit free Euler equation residuals on simulated points\n", + "%-------------------------------------------------------------------\n", + "discard = 200; % The number of observations to discard\n", + "time0 = tic;\n", + "\n", + "[Residuals_mean(1), Residuals_max(1), Residuals_max_E(1:9,1), Residuals] = NK_accuracy(nua,nuL,nuR,nuG,nuB,nuu,R,delta,L,Y,Yn,pie,S,F,C,rho_nua,rho_nuL,rho_nuR,rho_nuu,rho_nuB,rho_nuG,gam,vartheta,epsil,betta,phi_y,phi_pie,mu,theta,piestar,vcv,discard,vk,Gbar,zlb,Degree);\n", + "\n", + "accuracy_time = toc;\n", + "% Compute the residuals; Residuals_mean, Residuals_max and Residuals_max_E(1:9,5)\n", + "% are filled into the 5h columns of the corresponding vectors/matrix, while\n", + "% the first 4 columns correspond to PER1, PER2 with and without ZLB\n", + "\n", + "% ------------------------------------------------------------------------\n", + "% disp(' '); disp(' OUTPUT:'); disp(' ');\n", + "% disp('RUNNING TIME (in seconds):'); disp('');\n", + "% display(running_time);\n", + "% disp('SIMULATION TIME (in seconds):'); disp('');\n", + "% display(simulation_time);\n", + "% disp('ACCURACY TIME (in seconds):'); disp('');\n", + "% display(accuracy_time);\n", + "% disp('APPROXIMATION ERRORS (log10):'); disp('');\n", + "% disp('a) mean error in the model equations');\n", + "% disp(Residuals_mean)\n", + "% disp('b) max error in the model equations');\n", + "% disp(Residuals_max)\n", + "% disp('b) max error in by equation');\n", + "% disp(Residuals_max_E(1:9,1))\n", + "\n", + "l1 = log10(sum(max(abs(Residuals), [], 1)));\n", + "tot_time = running_time + simulation_time + accuracy_time;\n", + "fprintf('Solver time: %4.6f seconds\\n', running_time);\n", + "fprintf('Simulation time: %4.6f seconds\\n', simulation_time);\n", + "fprintf('Accuracy time: %4.6f seconds\\n', accuracy_time);\n", + "fprintf('Total time: %4.6f seconds\\n', running_time + simulation_time + accuracy_time)\n", + "fprintf('pistar: %3.4f sigma_L: %3.4f\\n', piestar, sigma_nuL);\n", + "fprintf('zlb : %d grid: %d\\n', zlb, grid_type);\n", + "fprintf('tex line: %0.2f & %0.2f & %0.2f\\n', l1, Residuals_max, tot_time);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Matlab", + "language": "matlab", + "name": "matlab" + }, + "language_info": { + "codemirror_mode": "octave", + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-matlab", + "name": "matlab", + "version": "0.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/NKModel_CLMM_python.ipynb b/NKModel_CLMM_python.ipynb new file mode 100644 index 0000000..de6c845 --- /dev/null +++ b/NKModel_CLMM_python.ipynb @@ -0,0 +1,828 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Solving a New Keynesian model with Python\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 Python packages. We recommend following the instructions on [quantecon.org](https://lectures.quantecon.org/jl/getting_started.html) for getting a base python installation set up. Then to acquire additional packages used in this notebook, uncomment the lines in the cell below (delete the `#` and space at the beginning of the line) and then run the cell:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# !pip install git+https://github.com/EconForge/interpolation.py.git\n", + "# !pip install git+https://github.com/naught101/sobol_seq.git" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "outputExpanded": false + }, + "source": [ + "## Python Code\n", + "\n", + "The Python version of our algorithm is implemented as a few methods defined on\n", + "a core class named `Model`. This class is itself composed of instances of three\n", + "different classes that hold the model parameters, steady state, and grids\n", + "needed to describe the numerical model. Before we get to the classes, we need\n", + "to bring in some dependencies:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true, + "outputExpanded": false + }, + "outputs": [], + "source": [ + "import math\n", + "from math import sqrt\n", + "import time as time\n", + "from collections import namedtuple\n", + "\n", + "import numpy as np\n", + "from numpy import exp\n", + "from interpolation.complete_poly import (_complete_poly_impl_vec,\n", + " _complete_poly_impl,\n", + " complete_polynomial)\n", + "\n", + "import sobol_seq\n", + "\n", + "# set seed on random number generator to make results reproducible\n", + "np.random.seed(42)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "outputExpanded": false + }, + "source": [ + "We will also need the following two functions, which use monomial rules to\n", + "compute quadrature nodes and weights:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true, + "outputExpanded": false + }, + "outputs": [], + "source": [ + "def qnwmonomial1(vcv):\n", + " n = vcv.shape[0]\n", + " n_nodes = 2*n\n", + "\n", + " z1 = np.zeros((n_nodes, n))\n", + "\n", + " # In each node, random variable i takes value either 1 or -1, and\n", + " # all other variables take value 0. For example, for N = 2,\n", + " # z1 = [1 0; -1 0; 0 1; 0 -1]\n", + " for i in range(n):\n", + " z1[2*i:2*(i+1), i] = [1, -1]\n", + "\n", + "\n", + " sqrt_vcv = np.linalg.cholesky(vcv)\n", + " R = np.sqrt(n)*sqrt_vcv\n", + " ϵj = z1 @ R\n", + " ωj = np.ones(n_nodes) / n_nodes\n", + "\n", + " return ϵj, ωj\n", + "\n", + "\n", + "def qnwmonomial2(vcv):\n", + " n = vcv.shape[0]\n", + " assert n == vcv.shape[1], \"Variance covariance matrix must be square\"\n", + " z0 = np.zeros((1, n))\n", + "\n", + " z1 = np.zeros((2*n, n))\n", + " # In each node, random variable i takes value either 1 or -1, and\n", + " # all other variables take value 0. For example, for N = 2,\n", + " # z1 = [1 0; -1 0; 0 1; 0 -1]\n", + " for i in range(n):\n", + " z1[2*i:2*(i+1), i] = [1, -1]\n", + "\n", + " z2 = np.zeros((2*n*(n-1), n))\n", + " i = 0\n", + "\n", + " # In each node, a pair of random variables (p,q) takes either values\n", + " # (1,1) or (1,-1) or (-1,1) or (-1,-1), and all other variables take\n", + " # value 0. For example, for N = 2, `z2 = [1 1; 1 -1; -1 1; -1 1]`\n", + " for p in range(n-1):\n", + " for q in range(p+1, n):\n", + " z2[4*i:4*(i+1), p] = [1, -1, 1, -1]\n", + " z2[4*i:4*(i+1), q] = [1, 1, -1, -1]\n", + " i += 1\n", + "\n", + " sqrt_vcv = np.linalg.cholesky(vcv)\n", + " R = np.sqrt(n+2)*sqrt_vcv\n", + " S = np.sqrt((n+2)/2)*sqrt_vcv\n", + " ϵj = np.row_stack([z0, z1 @ R, z2 @ S])\n", + "\n", + " ωj = np.concatenate([2/(n+2) * np.ones(z0.shape[0]),\n", + " (4-n)/(2*(n+2)**2) * np.ones(z1.shape[0]),\n", + " 1/(n+2)**2 * np.ones(z2.shape[0])])\n", + " return ϵj, ωj" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "outputExpanded": false + }, + "source": [ + "## Classes\n", + "\n", + "First we have the `Params` class, which holds all the model parameters as well\n", + "as the paramters that drive the algorithm." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true, + "outputExpanded": false + }, + "outputs": [], + "source": [ + "SteadyState = namedtuple(\"SteadyState\",\n", + " [\"Yn\", \"Y\", \"π\", \"δ\", \"L\", \"C\", \"F\", \"S\", \"R\", \"w\"])\n", + "\n", + "class Params(object):\n", + " def __init__(self, zlb=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, σηu=0.0054, σηB=0.0010,\n", + " σηG=0.0038, degree=2):\n", + "\n", + " self.zlb = zlb # whether or not the zlb should be imposed\n", + " self.γ = γ # Utility-function parameter\n", + " self.β = β # Discount factor\n", + " self.ϑ = ϑ # Utility-function parameter\n", + " self.ϵ = ϵ # Parameter in the Dixit-Stiglitz aggregator\n", + " self.ϕ_y = ϕ_y # Parameter of the Taylor rule\n", + " self.ϕ_π = ϕ_π # Parameter of the Taylor rule\n", + " self.μ = μ # Parameter of the Taylor rule\n", + " self.Θ = Θ # Share of non-reoptimizing firms (Calvo's pricing)\n", + " self.πstar = πstar # Target (gross) inflation rate\n", + " self.gbar = gbar # Steady-state share of government spending in output\n", + "\n", + " # autocorrelation coefficients\n", + " self.ρηR = ρηR # See process (28) in MM (2015)\n", + " self.ρηa = ρηa # See process (22) in MM (2015)\n", + " self.ρηL = ρηL # See process (16) in MM (2015)\n", + " self.ρηu = ρηu # See process (15) in MM (2015)\n", + " self.ρηB = ρηB # See process (17) in MM (2015)\n", + " self.ρηG = ρηG # See process (26) in MM (2015)\n", + "\n", + " # standard deviations\n", + " self.σηR = σηR # See process (28) in MM (2015)\n", + " self.σηa = σηa # See process (22) in MM (2015)\n", + " self.σηL = σηL # See process (16) in MM (2015)\n", + " self.σηu = σηu # See process (15) in MM (2015)\n", + " self.σηB = σηB # See process (17) in MM (2015)\n", + " self.σηG = σηG # See process (26) in MM (2015)\n", + "\n", + " self.degree = degree\n", + "\n", + " @property\n", + " def vcov(self):\n", + " return np.diag([self.σηR**2, self.σηa**2, self.σηL**2,\n", + " self.σηu**2, self.σηB**2, self.σηG**2])\n", + "\n", + " @property\n", + " def steady_state(self):\n", + " Yn_ss = exp(self.gbar)**(self.γ/(self.ϑ+self.γ))\n", + " Y_ss = Yn_ss\n", + " π_ss = 1.0\n", + " δ_ss = 1.0\n", + " L_ss = Y_ss/δ_ss\n", + " C_ss = (1-self.gbar)*Y_ss\n", + " F_ss = C_ss**(-self.γ)*Y_ss/(1-self.β*self.Θ*π_ss**(self.ϵ-1))\n", + " S_ss = L_ss**self.ϑ*Y_ss/(1-self.β*self.Θ*π_ss**self.ϵ)\n", + " R_ss = π_ss/self.β\n", + " w_ss = (L_ss**self.ϑ)*(C_ss**self.γ)\n", + "\n", + " return SteadyState(Yn_ss, Y_ss, π_ss, δ_ss, L_ss, C_ss, F_ss, S_ss, R_ss, w_ss)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "outputExpanded": false + }, + "source": [ + "Notice that we have a namedtuple to hold the steady state of the model. Using\n", + "the namedtuple infrastructure allows us to have convenient \"dot-style\" access\n", + "to the steady state, without defining a full class.\n", + "\n", + "Given an instance of `Params` class, we can construct the grid on which we will\n", + "solve the model.\n", + "\n", + "The `Grids` class holds this grid as well as matrices used to compute\n", + "expectations." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": true, + "outputExpanded": false + }, + "outputs": [], + "source": [ + "class Grids(object):\n", + "\n", + " def __init__(self, p, m=200, kind=\"rands\"):\n", + "\n", + " if kind == \"sobol\":\n", + " ub = np.array([\n", + " 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[[6, 7]] = [1.0, 0.95] # adjust lower bound for R and δ\n", + " s = sobol_seq.i4_sobol_generate(8, m)\n", + " s *= (ub - lb)\n", + " s += lb\n", + "\n", + " ηR = s[:, 0]\n", + " ηa = s[:, 1]\n", + " ηL = s[:, 2]\n", + " ηu = s[:, 3]\n", + " ηB = s[:, 4]\n", + " ηG = s[:, 5]\n", + " R = s[:, 6]\n", + " δ = s[:, 7]\n", + " else:\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*np.random.rand(m)) / sqrt(1-p.ρηR**2)\n", + " ηa = (-2*p.σηa + 4*p.σηa*np.random.rand(m)) / sqrt(1-p.ρηa**2)\n", + " ηL = (-2*p.σηL + 4*p.σηL*np.random.rand(m)) / sqrt(1-p.ρηL**2)\n", + " ηu = (-2*p.σηu + 4*p.σηu*np.random.rand(m)) / sqrt(1-p.ρηu**2)\n", + " ηB = (-2*p.σηB + 4*p.σηB*np.random.rand(m)) / sqrt(1-p.ρηB**2)\n", + " ηG = (-2*p.σηG + 4*p.σηG*np.random.rand(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*np.random.rand(m)\n", + " δ = 0.95 + 0.05*np.random.rand(m)\n", + "\n", + " self.ηR = ηR\n", + " self.ηa = ηa\n", + " self.ηL = ηL\n", + " self.ηu = ηu\n", + " self.ηB = ηB\n", + " self.ηG = ηG\n", + " self.R = R\n", + " self.δ = δ\n", + "\n", + " # shape (8, m)\n", + " self.X = np.vstack([np.log(R), np.log(δ), ηR, ηa, ηL, ηu, ηB, ηG])\n", + "\n", + " # shape (n_complete(8, 2), m)\n", + " self.X0_G = {\n", + " 1: complete_polynomial(self.X, 1),\n", + " p.degree: complete_polynomial(self.X, p.degree)\n", + " }\n", + "\n", + " # shape (2*n=12, n=6)\n", + " self.ϵ_nodes, self.ω_nodes = qnwmonomial1(p.vcov)\n", + "\n", + " # all shape (len(ϵ_nodes), m)\n", + " self.ηR1 = p.ρηR * ηR[None, :] + self.ϵ_nodes[:, None, 0]\n", + " self.ηa1 = p.ρηa * ηa[None, :] + self.ϵ_nodes[:, None, 1]\n", + " self.ηL1 = p.ρηL * ηL[None, :] + self.ϵ_nodes[:, None, 2]\n", + " self.ηu1 = p.ρηu * ηu[None, :] + self.ϵ_nodes[:, None, 3]\n", + " self.ηB1 = p.ρηB * ηB[None, :] + self.ϵ_nodes[:, None, 4]\n", + " self.ηG1 = p.ρηG * ηG[None, :] + self.ϵ_nodes[:, None, 5]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "outputExpanded": false + }, + "source": [ + "Finally, we construct the Model class, which has an instance of Params,\n", + "SteadyState and Grids as its three attributes.\n", + "\n", + "This block of code will be longer than the others because we also include\n", + "routines to solve and simulate the model as methods on the Model class. These\n", + "methods will be clearly marked and commented." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": true, + "outputExpanded": false + }, + "outputs": [], + "source": [ + "class Model(object):\n", + " def __init__(self, p=Params(), g=None):\n", + " if g is None:\n", + " g = Grids(p)\n", + "\n", + " self.p = p\n", + " self.g = g\n", + " self.s = self.p.steady_state\n", + "\n", + " def init_coefs(self):\n", + " \"Iniital guess for coefs. We evaluate interpoland as coefs @ basis_mat\"\n", + " npol = self.g.X0_G[1].shape[0]\n", + " coefs = np.full((3, npol), 1e-5)\n", + " coefs[:, 0] = [self.s.S, self.s.F, self.s.C**(-self.p.γ)]\n", + " return coefs\n", + "\n", + " def step(self, S, F, C, δ0, R0, ηG, ηa, ηL, ηR):\n", + " # simplify notation\n", + " Θ, ϵ, gbar, ϑ, γ = self.p.Θ, self.p.ϵ, self.p.gbar, self.p.ϑ, self.p.γ\n", + " β, μ, ϕ_π, ϕ_y, πs = self.p.β, self.p.μ, self.p.ϕ_π, self.p.ϕ_y, self.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)\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", + " return π0, δ1, Y0, L0, Yn0, R1\n", + "\n", + " def solve(self, damp=0.1, tol=1e-7):\n", + " # rename self to m to make code below readable\n", + " m = self\n", + "\n", + " n = len(m.g.ηR)\n", + " n_nodes = len(m.g.ω_nodes)\n", + "\n", + " ## allocate memory\n", + " # euler equations\n", + " e = np.zeros((3, n))\n", + "\n", + " # previous iteration S, F, C\n", + " S0_old_G = np.ones(n)\n", + " F0_old_G = np.ones(n)\n", + " C0_old_G = np.ones(n)\n", + "\n", + " # current iteration S, F, C\n", + " S0_new_G = np.ones(n)\n", + " F0_new_G = np.ones(n)\n", + " C0_new_G = np.ones(n)\n", + "\n", + " # future S, F, C\n", + " S1 = np.zeros((n_nodes, n))\n", + " F1 = np.zeros((n_nodes, n))\n", + " C1 = np.zeros((n_nodes, n))\n", + "\n", + " for deg in [1, self.p.degree]:\n", + " # housekeeping\n", + " err = 1.0\n", + " X0_G = m.g.X0_G[deg]\n", + "\n", + " if deg > 1:\n", + " # compute degree coefs using degree 1 coefs as guess\n", + " coefs = np.linalg.lstsq(X0_G.T, e.T)[0].T\n", + " else:\n", + " coefs = self.init_coefs()\n", + "\n", + " while err > tol:\n", + " # Current choices (at t)\n", + " # ------------------------------\n", + " SFC0 = coefs @ X0_G\n", + " S0 = SFC0[0, :] # Compute S(t) using coefs\n", + " F0 = SFC0[1, :] # Compute F(t) using coefs\n", + " C0 = (SFC0[2, :])**(-1/m.p.γ) # Compute C(t) using coefs\n", + " π0, δ1, Y0, L0, Yn0, R1 = self.step(S0, F0, C0, m.g.δ, m.g.R,\n", + " m.g.ηG, m.g.ηa, m.g.ηL, m.g.ηR)\n", + "\n", + " if self.p.zlb:\n", + " R1 = np.maximum(R1, 1.0)\n", + "\n", + " for u in range(n_nodes):\n", + " # Form complete polynomial of degree \"Degree\" (at t+1) on future state\n", + " grid1 = [np.log(R1), np.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", + "\n", + " X1 = complete_polynomial(grid1, deg)\n", + "\n", + " S1[u, :] = coefs[0, :] @ X1 # Compute S(t+1)\n", + " F1[u, :] = coefs[1, :] @ X1 # Compute F(t+1)\n", + " C1[u, :] = (coefs[2, :] @ X1)**(-1/m.p.γ) # Compute C(t+1)\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[0, :] = exp(m.g.ηu)*exp(m.g.ηL)*L0**m.p.ϑ*Y0/exp(m.g.ηa) + m.g.ω_nodes @ (m.p.β*m.p.Θ*π1**m.p.ϵ*S1)\n", + " e[1, :] = exp(m.g.ηu)*C0**(-m.p.γ)*Y0 + m.g.ω_nodes @ (m.p.β*m.p.Θ*π1**(m.p.ϵ-1)*F1)\n", + " e[2, :] = m.p.β*exp(m.g.ηB)/exp(m.g.ηu)*R1 * (m.g.ω_nodes @ ((exp(m.g.ηu1)*C1**(-m.p.γ)/π1)))\n", + "\n", + " # Variables of the current iteration\n", + " #-----------------------------------\n", + " np.copyto(S0_new_G, S0)\n", + " np.copyto(F0_new_G, F0)\n", + " np.copyto(C0_new_G, C0)\n", + "\n", + " # Compute and update the coefficients of the decision functions\n", + " # -------------------------------------------------------------\n", + " coefs_hat = np.linalg.lstsq(X0_G.T, e.T)[0].T\n", + "\n", + " # Update the coefficients using damping\n", + " coefs = damp*coefs_hat + (1-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 = (np.mean(abs(1-S0_new_G/S0_old_G)) +\n", + " np.mean(abs(1-F0_new_G/F0_old_G)) +\n", + " np.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", + " np.copyto(S0_old_G, S0_new_G)\n", + " np.copyto(F0_old_G, F0_new_G)\n", + " np.copyto(C0_old_G, C0_new_G)\n", + "\n", + " return coefs\n", + "\n", + " def simulate(self, coefs=None, capT=10201):\n", + " if coefs is None:\n", + " coefs = self.solve()\n", + "\n", + " # rename self to m to make code below readable\n", + " m = self\n", + "\n", + " # create namedtuple to hold simulation results in an organized container\n", + " Simulation = namedtuple(\"Simulation\",\n", + " [\"nuR\", \"nua\", \"nuL\", \"nuu\", \"nuB\", \"nuG\",\n", + " \"δ\", \"R\", \"S\", \"F\", \"C\", \"π\", \"Y\", \"L\", \"Yn\", \"w\"])\n", + "\n", + " # 11. Simualating a time-series solution\n", + " #---------------------------------------\n", + "\n", + " # Initialize the values of 6 exogenous shocks and draw innovations\n", + " #-----------------------------------------------------------------\n", + " nuR = np.zeros(capT)\n", + " nua = np.zeros(capT)\n", + " nuL = np.zeros(capT)\n", + " nuu = np.zeros(capT)\n", + " nuB = np.zeros(capT)\n", + " nuG = np.zeros(capT)\n", + "\n", + " # Generate the series for shocks\n", + " #-------------------------------\n", + " rands = np.random.randn(capT-1, 6)\n", + "\n", + " for t in range(capT-1):\n", + " nuR[t+1] = self.p.ρηR*nuR[t] + self.p.σηR*rands[t, 0]\n", + " nua[t+1] = self.p.ρηa*nua[t] + self.p.σηa*rands[t, 1]\n", + " nuL[t+1] = self.p.ρηL*nuL[t] + self.p.σηL*rands[t, 2]\n", + " nuu[t+1] = self.p.ρηu*nuu[t] + self.p.σηu*rands[t, 3]\n", + " nuB[t+1] = self.p.ρηB*nuB[t] + self.p.σηB*rands[t, 4]\n", + " nuG[t+1] = self.p.ρηG*nuG[t] + self.p.σηG*rands[t, 5]\n", + "\n", + " δ = np.ones(capT+1) # Allocate memory for the time series of delta(t)\n", + " R = np.ones(capT+1) # Allocate memory for the time series of R(t)\n", + " S = np.ones(capT) # Allocate memory for the time series of S(t)\n", + " F = np.ones(capT) # Allocate memory for the time series of F(t)\n", + " C = np.ones(capT) # Allocate memory for the time series of C(t)\n", + " π = np.ones(capT) # Allocate memory for the time series of π(t)\n", + " Y = np.ones(capT) # Allocate memory for the time series of Y(t)\n", + " L = np.ones(capT) # Allocate memory for the time series of L(t)\n", + " Yn = np.ones(capT) # Allocate memory for the time series of Yn(t)\n", + " w = np.ones(capT)\n", + "\n", + " pol_bases = np.empty(coefs.shape[1])\n", + " states = np.empty(8)\n", + "\n", + " for t in range(capT):\n", + " states[0] = math.log(R[t]); states[1] = math.log(δ[t])\n", + " states[2] = nuR[t]; states[3] = nua[t]\n", + " states[4] = nuL[t]; states[5] = nuu[t]\n", + " states[6] = nuB[t]; states[7] = nuG[t]\n", + "\n", + " _complete_poly_impl_vec(states, self.p.degree, pol_bases)\n", + "\n", + " vals = coefs @ pol_bases\n", + " S[t] = vals[0]\n", + " F[t] = vals[1]\n", + " C[t] = (vals[2])**(-1/m.p.γ)\n", + "\n", + " π[t], δ[t+1], Y[t], L[t], Yn[t], R[t+1] = self.step(S[t], F[t], C[t],\n", + " δ[t], R[t], nuG[t],\n", + " nua[t], nuL[t],\n", + " nuR[t])\n", + "\n", + " # Compute real wage\n", + " w[t] = exp(nuL[t])*(L[t]**m.p.ϑ)*(C[t]**m.p.γ)\n", + "\n", + " # If ZLB is imposed, set R(t)=1 if ZLB binds\n", + " if self.p.zlb:\n", + " R[t+1] = max(R[t+1], 1.0)\n", + "\n", + " return Simulation(nuR, nua, nuL, nuu, nuB, nuG, δ, R, S, F, C, π, Y, L, Yn, w)\n", + "\n", + " def residuals(self, coefs, sim, burn=200):\n", + " m = self # rename self to m so the rest of this code is more readable\n", + " capT = len(sim.w)\n", + " resids = np.zeros((capT, 9))\n", + "\n", + " # Integration method for evaluating accuracy\n", + " # ------------------------------------------\n", + " # Monomial integration rule with 2N**2+1 nodes\n", + " ϵ_nodes, ω_nodes = qnwmonomial2(m.p.vcov)\n", + " n_nodes = len(ω_nodes)\n", + "\n", + " # Allocate for arrays needed in the loop\n", + " basis_mat = np.empty((8, n_nodes))\n", + " X1 = np.empty((coefs.shape[1], n_nodes))\n", + "\n", + " nuR1 = np.empty(n_nodes)\n", + " nua1 = np.empty(n_nodes)\n", + " nuL1 = np.empty(n_nodes)\n", + " nuu1 = np.empty(n_nodes)\n", + " nuB1 = np.empty(n_nodes)\n", + " nuG1 = np.empty(n_nodes)\n", + "\n", + " for t in range(capT): # For each given point,\n", + " # Take the corresponding value for shocks at t\n", + " #---------------------------------------------\n", + " nuR0 = sim.nuR[t] # nuR(t)\n", + " nua0 = sim.nua[t] # nua(t)\n", + " nuL0 = sim.nuL[t] # nuL(t)\n", + " nuu0 = sim.nuu[t] # nuu(t)\n", + " nuB0 = sim.nuB[t] # nuB(t)\n", + " nuG0 = sim.nuG[t] # nuG(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", + " nuR1[:] = nuR0*m.p.ρηR + ϵ_nodes[:, 0]\n", + " nua1[:] = nua0*m.p.ρηa + ϵ_nodes[:, 1]\n", + " nuL1[:] = nuL0*m.p.ρηL + ϵ_nodes[:, 2]\n", + " nuu1[:] = nuu0*m.p.ρηu + ϵ_nodes[:, 3]\n", + " nuB1[:] = nuB0*m.p.ρηB + ϵ_nodes[:, 4]\n", + " nuG1[:] = nuG0*m.p.ρηG + ϵ_nodes[:, 5]\n", + "\n", + " basis_mat[0, :] = np.log(R1)\n", + " basis_mat[1, :] = np.log(δ1)\n", + " basis_mat[2, :] = nuR1\n", + " basis_mat[3, :] = nua1\n", + " basis_mat[4, :] = nuL1\n", + " basis_mat[5, :] = nuu1\n", + " basis_mat[6, :] = nuB1\n", + " basis_mat[7, :] = nuG1\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_poly_impl(basis_mat, self.p.degree, X1)\n", + "\n", + " # Compute S(t+1), F(t+1) and C(t+1) in all nodes using coefs\n", + " S1 = coefs[0, :] @ X1\n", + " F1 = coefs[1, :] @ X1\n", + " C1 = (coefs[2, :] @ X1)**(-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[t, 0] = 1-(ω_nodes @\n", + " (exp(nuu0)*exp(nuL0)*L0**m.p.ϑ*Y0/exp(nua0) +\n", + " m.p.β*m.p.Θ*π1**m.p.ϵ*S1)/S0\n", + " )\n", + " resids[t, 1] = 1 - (ω_nodes @\n", + " (exp(nuu0)*C0**(-m.p.γ)*Y0 + m.p.β*m.p.Θ*π1**(m.p.ϵ-1)*F1)/F0\n", + " )\n", + " resids[t, 2] = 1.0 -(ω_nodes @\n", + " (m.p.β*exp(nuB0)/exp(nuu0)*R1*exp(nuu1)*C1**(-m.p.γ)/π1)/C0**(-m.p.γ)\n", + " )\n", + "\n", + " resids[t, 3] = 1-((1-m.p.Θ*π0**(m.p.ϵ-1))/(1-m.p.Θ))**(1/(1-m.p.ϵ))*F0/S0\n", + " resids[t, 4] = 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[t, 5] = 1-exp(nua0)*L0*δ1/Y0\n", + " resids[t, 6] = 1-(1-m.p.gbar/exp(nuG0))*Y0/C0\n", + " resids[t, 7] = 1-(exp(nua0)**(1+m.p.ϑ)*(1-m.p.gbar/exp(nuG0))**(-m.p.γ)/exp(nuL0))**(1/(m.p.ϑ+m.p.γ))/Yn0\n", + " resids[t, 8] = 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(nuR0)/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 and R1 <= 1:\n", + " resids[t, 8] = 0.0\n", + "\n", + " return resids[burn:, :]" + ] + }, + { + "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": 6, + "metadata": { + "collapsed": true, + "outputExpanded": false + }, + "outputs": [], + "source": [ + "def main(m=Model(), file=None):\n", + " if file is None:\n", + " mprint = print\n", + " else:\n", + " def mprint(*x):\n", + " print(*x, file=file)\n", + "\n", + " # solve the model\n", + " t1 = time.time()\n", + " coefs = m.solve()\n", + " solve_time = time.time() - t1\n", + "\n", + " # simulate the model\n", + " t1 = time.time()\n", + " sim = m.simulate(coefs)\n", + " sim_time = time.time() - t1\n", + "\n", + " # check accuracy\n", + " t1 = time.time()\n", + " resids = m.residuals(coefs, sim)\n", + " resids_time = time.time() - t1\n", + "\n", + " tot_time = solve_time + sim_time + resids_time\n", + " mean_err = np.log10(abs(resids).mean())\n", + " max_err = np.log10(abs(resids).max())\n", + " max_err_eqn = np.log10(abs(resids).max(1) + 1e-16)\n", + " l1 = np.log10(abs(resids).max(0).sum())\n", + "\n", + " mprint(\"Solver time (in seconds): \", solve_time)\n", + " mprint(\"Simulation time (in seconds): \", sim_time)\n", + " mprint(\"Residuals time (in seconds): \", resids_time)\n", + " mprint(\"Total time (in seconds): \", tot_time)\n", + " mprint(\"\\nAPPROXIMATION ERRORS (log10):\")\n", + " mprint(\"\\ta) mean error in the model equations: {:0.3f}\".format(mean_err))\n", + " mprint(\"\\tb) sum of max error per equation: {:0.3f}\".format(l1));\n", + " mprint(\"\\tc) max error in the model equations: {:0.3f}\".format(max_err))\n", + " mprint(\"\\td) max error by equation: \", max_err_eqn)\n", + " mprint(\"tex row:\", \"{:.2f} & {:.2f} & {:.2f}\".format(l1, max_err, tot_time))\n", + "\n", + " return solve_time, sim_time, resids_time, coefs, sim, resids\n", + "\n", + "\n", + "def build_paper_table():\n", + " with open(\"output.log\", \"w\") as 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", + "\n", + "\n", + " for grid_kind in [\"sobol\", \"random\"]:\n", + " p = Params(**params)\n", + " g = Grids(p, kind=grid_kind)\n", + " m = Model(p, g)\n", + "\n", + " print(\"working with params:\", params, file=f)\n", + " print(\"And grid type:\", grid_kind, file=f)\n", + " main(m, f)\n", + " print(\"\\n\"*5, file=f)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "build_paper_table()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Solver time (in seconds): 1.620750904083252\n", + "Simulation time (in seconds): 0.671532154083252\n", + "Residuals time (in seconds): 0.826570987701416\n", + "Total time (in seconds): 3.11885404586792\n", + "\n", + "APPROXIMATION ERRORS (log10):\n", + "\ta) mean error in the model equations: -4.342\n", + "\tb) sum of max error per equation: -1.375\n", + "\tc) max error in the model equations: -1.704\n", + "\td) max error by equation: [-3.73802289 -3.31801009 -3.40755216 ..., -4.56959248 -3.7830201 -3.6766427 ]\n", + "tex row: -1.37 & -1.70 & 3.12\n" + ] + } + ], + "source": [ + "main();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.0" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}