{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Projection implementation of stochastic growth model adapted from Fabrice Collard's Matlab code, http://fabcol.free.fr/\n", "# tested in Julia 0.6.1.1\n", "# this code is part of chapter 5, \"Advanced Numerical Techniques\" from the book: \"Introduction to Quantitative Macroeconomics using Julia\"\n", "# Academic Press - Elsevier\n", "# for comments, email at: petre(dot)caraiani(at)gmail(dot)com\n", "\n", "# The Optimal Growth Model\n", "# Collocation method (AR(1) case) \n", "global nbk, kmin, ksup, XK, kt;\n", "global nba, amin, asup, XA, at;\n", "global nstate, nodea, nodek, wmat, wgrid;\n", "\n", "# Parameters \n", "nbk = 4; # Degree of polynomials (capital)\n", "nba = 2; # Degree of polynomials (technology shock)\n", "ncoef = (nbk+1)*(nba+1); # of coefficients\n", "nodek = nbk+1; # of Nodes\n", "nodea = nba+1; # of Nodes\n", "nstate= 12;\n", " \n", "# Structural Parameters \n", "delta = 0.1;\n", "beta = 0.95;\n", "alpha = 0.3;\n", "sigma = 1.5;\n", "ysk =(1-beta*(1-delta))/(alpha*beta);\n", "ksy = 1/ysk;\n", "ys = ksy^(alpha/(1-alpha));\n", "ks = ys^(1/alpha);\n", "is = delta*ks;\n", "cs = ys-is;\n", "ab = 0;\n", "\n", "# Technology shock process\n", "rho = 0.8;\n", "se = 0.2;\n", "ma = 0;" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "function hernodes(nstate)\n", "\n", "TOL = sqrt(2.2204e-16);\n", "MAXIT = 30;\n", "PIM4 = pi^(-1/4);\n", "\n", "n = nstate;\n", "m = (n+1)/2;\n", "m = convert(Int64, floor(m));\n", "x = zeros(n,1);\n", "w = zeros(n,1);\n", "\n", "z = 0;\n", "z1 = 0;\n", "pp = 0;\n", "for i= 1:m\n", "\n", "# Initialize the first four roots\n", "if i == 1; \n", " z= sqrt(2*n+1) - 1.85575 * (2*n+1)^(-1/6);\n", "elseif i== 2;\n", " z= z - 1.14 * (n^.426)/z;\n", "elseif i== 3;\n", " z= 1.86 * z - .86 * x[1];\n", "elseif i== 4;\n", " z= 1.91 * z - .91 * x[2];\n", "else;\n", " z= 2 * z - x[i-2];\n", "end;\n", "\n", "for its= 1:MAXIT\n", " p1 = PIM4;\n", " p2 = 0;\n", " for j = 1:n\n", " p3 = p2;\n", " p2 = p1;\n", " p1 = z*sqrt(2/j)*p2 - sqrt((j-1)/j)*p3;\n", " end;\n", "\n", " pp = p2 * sqrt(2*n);\n", " z1 = z;\n", " z = z1 - p1/pp;\n", "if abs(z-z1) < TOL; break; end;\n", "end;\n", "\n", "x[i] = z;\n", "x[n+1-i]= -z;\n", "w[i]= 2/(pp*pp);\n", "w[n+1-i]= w[i];\n", "end; \n", " \n", "return x,w\n", " end;" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "wgrid,wmat=hernodes(nstate);\n", "wgrid=wgrid*sqrt(2)*se;" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "#define functions rcheb transfo itransfo\n", "function rcheb(nn);\n", "mod=nn-floor(nn/2)*2;\n", "n1=floor(nn/2);\n", "k=collect(1:n1).';\n", "r1=cos.((2*k-1)*pi/(2*nn));\n", "r1=[r1 -r1];\n", "if mod==1;\n", " r1=[r1 0];\n", "end;\n", "r1=vec(r1);\n", "rr=real(sort!(r1));\n", "return rr\n", "end;\n", "\n", "function transfo(x,xmin,xmax);\n", " z=(2*(x-xmin)/(xmax-xmin))-1;\n", "return z\n", "end;\n", "\n", "function itransfo(x,xmin,xmax);\n", "z=0.5*(x+1)*(xmax-xmin)+xmin;\n", "return z\n", "end;\n", "\n", "function cheb(xx,nn);\n", "cc=real(cos.(kron(acos.(complex(xx)),nn))) \n", "return cc\n", " end;" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# grid for the income\n", "amin=(ma+wgrid[nstate]);\n", "asup=(ma+wgrid[1]);\n", "ra=rcheb(nodea); #roots\n", "at=itransfo(ra,amin,asup); #grid\n", "vnba = collect(0:nba).';\n", "XA=cheb(ra,vnba); # Polynomials" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# grid for the capital stock\n", "kmin = log(1.2);\n", "ksup = log(6);\n", "rk = rcheb(nodek); #roots\n", "kt = exp.(itransfo(rk,kmin,ksup)) #grid\n", "vrbk = collect(0:nbk).';\n", "XK = cheb(rk,vrbk);" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "#initial conditions\n", "a0=[ -0.23759592487257;\n", " 0.60814488103911;\n", " 0.03677400318790;\n", " 0.69025680170443;\n", " -0.21654209984197;\n", " 0.00551243342828;\n", " 0.03499834613714;\n", " -0.00341171507904;\n", " -0.00449139656933;\n", " 0.00085302605779;\n", " 0.00285737302122;\n", " -0.00002348542016;\n", " -0.00011606672164;\n", " -0.00003323351559;\n", " 0.00018045618825];\n", "param=[alpha beta delta sigma ma rho]';" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# function makepoly\n", "function makepoly(XA,XW);\n", "nba = size(XA,2);\n", "nba1 = size(XA,1);\n", "nbw = size(XW,2);\n", "nmax = max(nba,nbw);\n", "XX = Float64[] ; \n", " \n", "for i=1:nbw\n", " for j=1:nba \n", " XX=[XX; kron(XW[:,i],XA[:,j])];\n", " end\n", "end\n", "return XX\n", " end;" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "f =function(theta);\n", "RHS=Float64[];\n", "LHS=Float64[];\n", "XX =Float64[];\n", "ct = 0.0; \n", "c1 = 0.0; \n", "for i = 1:nodek\n", " for j=1:nodea\n", " XX0 = makepoly(XA[j,:],XK[i,:]);\n", " ct = exp.(dot(XX0,theta));\n", " XX = [XX;XX0];\n", " k1 = exp.(at[j])*kt[i].^alpha+(1-delta)*kt[i]-ct;\n", " rk1 = transfo(log.(k1),kmin,ksup);\n", " #if abs(rk1)>1;disp('problème k(t+1)');end\n", " vrbk = collect(0:nbk).'; \n", " xk1 = cheb(complex(rk1),vrbk); \n", " a1 = rho*at[j]+(1-rho)*ma+wgrid;\n", " ra1 = transfo(a1,amin,asup); # grid\n", " #if abs(ra1)>1;disp('problème a(t+1)');end\n", " XA1 = cheb(complex(ra1),vnba); # Polynomials\n", " XX1 = makepoly(XA1,xk1);\n", " XX1 = reshape(XX1,size(XA1,1), size(XA1,2)*size(xk1,2)) \n", " c1 = exp.(XX1*theta);\n", " aux = wmat'*(beta*(alpha*exp.(a1)*k1.^(alpha-1)+1-delta).*c1.^(-sigma))/sqrt(pi);\n", " RHS = [RHS;log.(aux)];\n", " LHS = [LHS;log.(ct.^(-sigma))];\n", " end; \n", "end;\n", "resid = LHS-RHS;\n", "return vec(resid);\n", "end;\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Display final rule:[0.117961, 0.493804, 0.0406254, 0.361884, -0.104603, -0.00199674, 0.00831576, 0.00119796, -0.000900948, 3.45746e-5, 0.00026351, 5.36035e-5, -4.07343e-6, -1.03949e-5, 3.20045e-6]\n" ] } ], "source": [ "#Pkg.add(\"NLsolve\");\n", "using NLsolve;\n", "sol = nlsolve(not_in_place(f),a0);\n", "println(\"Display final rule:\",sol.zero)\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "#final\n", "lt = length(sol.zero);\n", "nk = 20;\n", "na = 10;\n", "kt = collect(kmin:(ksup-kmin)/(nk-1):ksup);\n", "at = collect(amin:(asup-amin)/(na-1):asup);\n", "rk = transfo(kt,kmin,ksup);\n", "rk = vec(rk);\n", "ra = transfo(at,kmin,asup);\n", "ra = vec(ra);\n", "vrbk= collect(0:nbk).';\n", "XT = cheb(rk,vrbk);\n", "vnba= collect(0:nba).';\n", "XA = cheb(ra,vnba);\n", "kt = exp.(kt);\n", "ct = zeros(nk,na);\n", "k1 = zeros(nk,na);\n", "ii=[];\n", "for i=1:nk\n", " for j=1:na;\n", " XX0 = makepoly(XA[j,:],XT[i,:]);\n", " ct[i,j] = exp(dot(XX0,sol.zero));\n", " k1[i,j] = exp(at[j])*kt[i].^alpha+(1-delta)*kt[i]-ct[i,j];\n", "end;\n", "end;\n" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Julia 0.6.1", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }