Created
February 12, 2018 00:08
-
-
Save chriselrod/f30fdfadccd8413000f0703a7cc6f8a5 to your computer and use it in GitHub Desktop.
Revisions
-
chriselrod created this gist
Feb 12, 2018 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,437 @@ { "cells": [ { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "mainOld (generic function with 2 methods)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Compat, BenchmarkTools\n", "\n", "\n", "function mainNew(verbose = true)\n", " ## 1. Calibration\n", "\n", " α = 1/3 # Elasticity of output w.r.t. capital\n", " β = 0.95 # Discount factor\n", "\n", " # Productivity values\n", " vProductivity = [0.9792 0.9896 1.0000 1.0106 1.0212]\n", "\n", " # Transition matrix\n", " mTransition = [0.9727 0.0273 0.0000 0.0000 0.0000;\n", " 0.0041 0.9806 0.0153 0.0000 0.0000;\n", " 0.0000 0.0082 0.9837 0.0082 0.0000;\n", " 0.0000 0.0000 0.0153 0.9806 0.0041;\n", " 0.0000 0.0000 0.0000 0.0273 0.9727]\n", "\n", " # 2. Steady State\n", "\n", " capitalSteadyState = (α*β)^(1/(1-α))\n", " outputSteadyState = capitalSteadyState^α\n", " consumptionSteadyState = outputSteadyState-capitalSteadyState\n", "\n", " verbose && println(\"Output = \",outputSteadyState,\" Capital = \",capitalSteadyState,\" Consumption = \",consumptionSteadyState)\n", "\n", " # We generate the grid of capital\n", " vGridCapital = collect(0.5*capitalSteadyState:0.00001:1.5*capitalSteadyState)\n", " #vGridCapital = 0.5*capitalSteadyState:0.00001:1.5*capitalSteadyState\n", "\n", " nGridCapital = length(vGridCapital)\n", " nGridProductivity = length(vProductivity)\n", "\n", " # 3. Required matrices and vectors\n", "\n", " mOutput = Matrix{Float64}(uninitialized, nGridCapital, nGridProductivity)\n", " mValueFunction = fill(0.0, nGridCapital, nGridProductivity)\n", " mValueFunctionNew = Matrix{Float64}(uninitialized, nGridCapital,nGridProductivity)\n", " mPolicyFunction = Matrix{Float64}(uninitialized, nGridCapital,nGridProductivity)\n", " expectedValueFunction = Matrix{Float64}(uninitialized, nGridCapital, nGridProductivity)\n", "\n", " # 4. We pre-build output for each point in the grid\n", "\n", " @. mOutput = vGridCapital^ α * vProductivity;\n", "\n", " # 5. Main iteration\n", "\n", " maxDifference = 10.0\n", " tolerance = 0.0000001\n", " iteration = 0\n", "\n", " while(maxDifference > tolerance)\n", " #A_mul_Bt! for Julia < 0.7\n", " A_mul_Bt!(expectedValueFunction, mValueFunction, mTransition)\n", "# mul!(expectedValueFunction, mValueFunction, mTransition')\n", "\n", " @inbounds for nProductivity in 1:nGridProductivity\n", "\n", " # We start from previous choice (monotonicity of policy function)\n", " gridCapitalNextPeriod = 1\n", "\n", " for nCapital in 1:nGridCapital\n", "\n", " valueHighSoFar = -1000.0\n", " capitalChoice = vGridCapital[1]\n", "\n", " for nCapitalNextPeriod in gridCapitalNextPeriod:nGridCapital\n", "\n", " consumption = mOutput[nCapital,nProductivity]-vGridCapital[nCapitalNextPeriod]\n", " valueProvisional = (1-β)*log(consumption)+β*expectedValueFunction[nCapitalNextPeriod,nProductivity]\n", "\n", " if (valueProvisional>valueHighSoFar)\n", " \t valueHighSoFar = valueProvisional\n", " \t capitalChoice = vGridCapital[nCapitalNextPeriod]\n", " \t gridCapitalNextPeriod = nCapitalNextPeriod\n", " else\n", " \t break # We break when we have achieved the max\n", " end\n", "\n", " end\n", "\n", " mValueFunctionNew[nCapital,nProductivity] = valueHighSoFar\n", " mPolicyFunction[nCapital,nProductivity] = capitalChoice\n", "\n", " end\n", "\n", " end\n", " \n", " # I can write into expectedValueFunction, because we are about to\n", " # overwrite it anyway on the start of the next iteration.\n", " @. expectedValueFunction = abs(mValueFunctionNew-mValueFunction)\n", " maxDifference = maximum(expectedValueFunction)\n", " mValueFunction, mValueFunctionNew = mValueFunctionNew, mValueFunction\n", " \n", " iteration = iteration+1\n", " if verbose && (mod(iteration,10)==0 || iteration == 1)\n", " println(\" Iteration = \", iteration, \" Sup Diff = \", maxDifference)\n", " end\n", "\n", " end\n", "\n", " if verbose\n", " println(\" Iteration = \", iteration, \" Sup Diff = \", maxDifference)\n", " println(\" \")\n", " println(\" My check = \", mPolicyFunction[1000,3])\n", " println(\" My check = \", mValueFunction[1000,3])\n", " end\n", " return iteration, maxDifference, mPolicyFunction[1000,3], mValueFunction[1000,3]\n", "\n", "end\n", "\n", "\n", "function mainOld(verbose = true)\n", "\n", " ## 1. Calibration\n", "\n", " α = 1/3 # Elasticity of output w.r.t. capital\n", " β = 0.95 # Discount factor\n", "\n", " # Productivity values\n", " vProductivity = [0.9792 0.9896 1.0000 1.0106 1.0212]\n", "\n", " # Transition matrix\n", " mTransition = [0.9727 0.0273 0.0000 0.0000 0.0000;\n", " 0.0041 0.9806 0.0153 0.0000 0.0000;\n", " 0.0000 0.0082 0.9837 0.0082 0.0000;\n", " 0.0000 0.0000 0.0153 0.9806 0.0041;\n", " 0.0000 0.0000 0.0000 0.0273 0.9727]\n", "\n", " # 2. Steady State\n", "\n", " capitalSteadyState = (α*β)^(1/(1-α))\n", " outputSteadyState = capitalSteadyState^α\n", " consumptionSteadyState = outputSteadyState-capitalSteadyState\n", "\n", " verbose && println(\"Output = \",outputSteadyState,\" Capital = \",capitalSteadyState,\" Consumption = \",consumptionSteadyState)\n", "\n", " # We generate the grid of capital\n", " vGridCapital = collect(0.5*capitalSteadyState:0.00001:1.5*capitalSteadyState)\n", "\n", " nGridCapital = length(vGridCapital)\n", " nGridProductivity = length(vProductivity)\n", "\n", " # 3. Required matrices and vectors\n", "\n", " mOutput = zeros(nGridCapital,nGridProductivity)\n", " mValueFunction = zeros(nGridCapital,nGridProductivity)\n", " mValueFunctionNew = zeros(nGridCapital,nGridProductivity)\n", " mPolicyFunction = zeros(nGridCapital,nGridProductivity)\n", " expectedValueFunction = zeros(nGridCapital,nGridProductivity)\n", "\n", " # 4. We pre-build output for each point in the grid\n", "\n", " mOutput = (vGridCapital.^α)*vProductivity;\n", "\n", " # 5. Main iteration\n", "\n", " maxDifference = 10.0\n", " tolerance = 0.0000001\n", " iteration = 0\n", "\n", " while(maxDifference > tolerance)\n", " expectedValueFunction = mValueFunction*mTransition';\n", "\n", " for nProductivity in 1:nGridProductivity\n", "\n", " # We start from previous choice (monotonicity of policy function)\n", " gridCapitalNextPeriod = 1\n", "\n", " for nCapital in 1:nGridCapital\n", "\n", " valueHighSoFar = -1000.0\n", " capitalChoice = vGridCapital[1]\n", "\n", " for nCapitalNextPeriod in gridCapitalNextPeriod:nGridCapital\n", "\n", " consumption = mOutput[nCapital,nProductivity]-vGridCapital[nCapitalNextPeriod]\n", " valueProvisional = (1-β)*log(consumption)+β*expectedValueFunction[nCapitalNextPeriod,nProductivity]\n", "\n", " if (valueProvisional>valueHighSoFar)\n", " \t valueHighSoFar = valueProvisional\n", " \t capitalChoice = vGridCapital[nCapitalNextPeriod]\n", " \t gridCapitalNextPeriod = nCapitalNextPeriod\n", " else\n", " \t break # We break when we have achieved the max\n", " end\n", "\n", " end\n", "\n", " mValueFunctionNew[nCapital,nProductivity] = valueHighSoFar\n", " mPolicyFunction[nCapital,nProductivity] = capitalChoice\n", "\n", " end\n", "\n", " end\n", " \n", " maxDifference = maximum(abs.(mValueFunctionNew-mValueFunction))\n", " mValueFunction, mValueFunctionNew = mValueFunctionNew, mValueFunction\n", " \n", " iteration = iteration+1\n", " if verbose && (mod(iteration,10)==0 || iteration == 1)\n", " println(\" Iteration = \", iteration, \" Sup Diff = \", maxDifference)\n", " end\n", "\n", " end\n", "\n", "\n", " if verbose\n", " println(\" Iteration = \", iteration, \" Sup Diff = \", maxDifference)\n", " println(\" \")\n", " println(\" My check = \", mPolicyFunction[1000,3])\n", " println(\" My check = \", mValueFunction[1000,3])\n", " end\n", " return iteration, maxDifference, mPolicyFunction[1000,3], mValueFunction[1000,3]\n", "\n", "\n", "end" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output = 0.5627314338711378 Capital = 0.178198287392527 Consumption = 0.3845331464786108\n", " Iteration = 1 Sup Diff = 0.05274159340733661\n", " Iteration = 10 Sup Diff = 0.031346949265852075\n", " Iteration = 20 Sup Diff = 0.01870345989335709\n", " Iteration = 30 Sup Diff = 0.01116551203397076\n", " Iteration = 40 Sup Diff = 0.00666854170813258\n", " Iteration = 50 Sup Diff = 0.003984292748717033\n", " Iteration = 60 Sup Diff = 0.0023813118039327508\n", " Iteration = 70 Sup Diff = 0.0014236586450983024\n", " Iteration = 80 Sup Diff = 0.0008513397747205165\n", " Iteration = 90 Sup Diff = 0.0005092051752288995\n", " Iteration = 100 Sup Diff = 0.00030462324421465237\n", " Iteration = 110 Sup Diff = 0.00018226485632300005\n", " Iteration = 120 Sup Diff = 0.00010906950872624499\n", " Iteration = 130 Sup Diff = 6.527643736320421e-5\n", " Iteration = 140 Sup Diff = 3.907108211997912e-5\n", " Iteration = 150 Sup Diff = 2.3388077119990136e-5\n", " Iteration = 160 Sup Diff = 1.4008644637186762e-5\n", " Iteration = 170 Sup Diff = 8.391317202871562e-6\n", " Iteration = 180 Sup Diff = 5.026474537817016e-6\n", " Iteration = 190 Sup Diff = 3.010899863653549e-6\n", " Iteration = 200 Sup Diff = 1.8035522479920019e-6\n", " Iteration = 210 Sup Diff = 1.080340915837752e-6\n", " Iteration = 220 Sup Diff = 6.471316943423844e-7\n", " Iteration = 230 Sup Diff = 3.876361938104367e-7\n", " Iteration = 240 Sup Diff = 2.3219657907525004e-7\n", " Iteration = 250 Sup Diff = 1.3908720930544405e-7\n", " Iteration = 257 Sup Diff = 9.716035664908418e-8\n", " \n", " My check = 0.1465491436962635\n", " My check = -0.9714880021887344\n" ] }, { "data": { "text/plain": [ "(257, 9.716035664908418e-8, 0.1465491436962635, -0.9714880021887344)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mainNew()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output = 0.5627314338711378 Capital = 0.178198287392527 Consumption = 0.3845331464786108\n", " Iteration = 1 Sup Diff = 0.05274159340733661\n", " Iteration = 10 Sup Diff = 0.031346949265852075\n", " Iteration = 20 Sup Diff = 0.01870345989335709\n", " Iteration = 30 Sup Diff = 0.01116551203397076\n", " Iteration = 40 Sup Diff = 0.00666854170813258\n", " Iteration = 50 Sup Diff = 0.003984292748717033\n", " Iteration = 60 Sup Diff = 0.0023813118039327508\n", " Iteration = 70 Sup Diff = 0.0014236586450983024\n", " Iteration = 80 Sup Diff = 0.0008513397747205165\n", " Iteration = 90 Sup Diff = 0.0005092051752288995\n", " Iteration = 100 Sup Diff = 0.00030462324421465237\n", " Iteration = 110 Sup Diff = 0.00018226485632300005\n", " Iteration = 120 Sup Diff = 0.00010906950872624499\n", " Iteration = 130 Sup Diff = 6.527643736320421e-5\n", " Iteration = 140 Sup Diff = 3.907108211997912e-5\n", " Iteration = 150 Sup Diff = 2.3388077119990136e-5\n", " Iteration = 160 Sup Diff = 1.4008644637186762e-5\n", " Iteration = 170 Sup Diff = 8.391317202871562e-6\n", " Iteration = 180 Sup Diff = 5.026474537817016e-6\n", " Iteration = 190 Sup Diff = 3.010899863653549e-6\n", " Iteration = 200 Sup Diff = 1.8035522479920019e-6\n", " Iteration = 210 Sup Diff = 1.080340915837752e-6\n", " Iteration = 220 Sup Diff = 6.471316943423844e-7\n", " Iteration = 230 Sup Diff = 3.876361938104367e-7\n", " Iteration = 240 Sup Diff = 2.3219657907525004e-7\n", " Iteration = 250 Sup Diff = 1.3908720930544405e-7\n", " Iteration = 257 Sup Diff = 9.716035664908418e-8\n", " \n", " My check = 0.1465491436962635\n", " My check = -0.9714880021887344\n" ] }, { "data": { "text/plain": [ "(257, 9.716035664908418e-8, 0.1465491436962635, -0.9714880021887344)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mainOld()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 3.54 MiB\n", " allocs estimate: 15\n", " --------------\n", " minimum time: 1.333 s (0.00% GC)\n", " median time: 1.335 s (0.00% GC)\n", " mean time: 1.336 s (0.00% GC)\n", " maximum time: 1.342 s (0.00% GC)\n", " --------------\n", " samples: 4\n", " evals/sample: 1" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark mainNew(false)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 528.54 MiB\n", " allocs estimate: 1563\n", " --------------\n", " minimum time: 1.569 s (0.80% GC)\n", " median time: 1.598 s (0.83% GC)\n", " mean time: 1.605 s (1.69% GC)\n", " maximum time: 1.656 s (4.21% GC)\n", " --------------\n", " samples: 4\n", " evals/sample: 1" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark mainOld(false)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" } }, "nbformat": 4, "nbformat_minor": 2 } This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,434 @@ { "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Basic RBC model with full depreciation (Alternate 1)\n", "#\n", "# Jesus Fernandez-Villaverde\n", "# Haverford, July 3, 2013\n", "\n", "import numpy as np\n", "import math\n", "import time\n", "from numba import autojit\n", "\n", "# - Start Inner Loop - #\n", "# - bbeta float\n", "# - nGridCapital: int64\n", "# - gridCapitalNextPeriod: int64\n", "# - mOutput: float (17820 x 5)\n", "# - nProductivity: int64\n", "# - vGridCapital: float (17820, )\n", "# - mValueFunction: float (17820 x 5)\n", "# - mPolicyFunction: float (17820 x 5)\n", "\n", "@autojit\n", "def innerloop(bbeta, nGridCapital, gridCapitalNextPeriod, mOutput, nProductivity, vGridCapital, expectedValueFunction, mValueFunction, mValueFunctionNew, mPolicyFunction):\n", "\n", " for nCapital in range(nGridCapital):\n", " valueHighSoFar = -100000.0\n", " capitalChoice = vGridCapital[0]\n", " \n", " for nCapitalNextPeriod in range(gridCapitalNextPeriod, nGridCapital):\n", " consumption = mOutput[nCapital,nProductivity] - vGridCapital[nCapitalNextPeriod]\n", " valueProvisional = (1-bbeta)*np.log(consumption)+bbeta*expectedValueFunction[nCapitalNextPeriod,nProductivity];\n", "\n", " if valueProvisional > valueHighSoFar:\n", " valueHighSoFar = valueProvisional\n", " capitalChoice = vGridCapital[nCapitalNextPeriod]\n", " gridCapitalNextPeriod = nCapitalNextPeriod\n", " else:\n", " break \n", "\n", " mValueFunctionNew[nCapital,nProductivity] = valueHighSoFar\n", " mPolicyFunction[nCapital,nProductivity] = capitalChoice\n", "\n", " return mValueFunctionNew, mPolicyFunction\n", "\n", "def main_func():\n", "\n", " # 1. Calibration\n", "\n", " aalpha = 1.0/3.0 # Elasticity of output w.r.t. capital\n", " bbeta = 0.95 # Discount factor\n", "\n", " # Productivity values\n", " vProductivity = np.array([0.9792, 0.9896, 1.0000, 1.0106, 1.0212],float)\n", "\n", " # Transition matrix\n", " mTransition = np.array([[0.9727, 0.0273, 0.0000, 0.0000, 0.0000],\n", " [0.0041, 0.9806, 0.0153, 0.0000, 0.0000],\n", " [0.0000, 0.0082, 0.9837, 0.0082, 0.0000],\n", " [0.0000, 0.0000, 0.0153, 0.9806, 0.0041],\n", " [0.0000, 0.0000, 0.0000, 0.0273, 0.9727]],float)\n", "\n", " ## 2. Steady State\n", "\n", " capitalSteadyState = (aalpha*bbeta)**(1/(1-aalpha))\n", " outputSteadyState = capitalSteadyState**aalpha\n", " consumptionSteadyState = outputSteadyState-capitalSteadyState\n", "\n", " print(\"Output = \", outputSteadyState, \" Capital = \", capitalSteadyState, \" Consumption = \", consumptionSteadyState)\n", "\n", " # We generate the grid of capital\n", " vGridCapital = np.arange(0.5*capitalSteadyState,1.5*capitalSteadyState,0.00001)\n", "\n", " nGridCapital = len(vGridCapital)\n", " nGridProductivity = len(vProductivity)\n", "\n", " ## 3. Required matrices and vectors\n", "\n", " mOutput = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " mValueFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " mValueFunctionNew = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " mPolicyFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " expectedValueFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", "\n", " # 4. We pre-build output for each point in the grid\n", "\n", " for nProductivity in range(nGridProductivity):\n", " mOutput[:,nProductivity] = vProductivity[nProductivity]*(vGridCapital**aalpha)\n", "\n", " ## 5. Main iteration\n", "\n", " maxDifference = 10.0\n", " tolerance = 0.0000001\n", " iteration = 0\n", "\n", " log = math.log\n", " zeros = np.zeros\n", " dot = np.dot\n", "\n", " while(maxDifference > tolerance):\n", "\n", " expectedValueFunction = dot(mValueFunction,mTransition.T)\n", "\n", " for nProductivity in range(nGridProductivity):\n", "\n", " # We start from previous choice (monotonicity of policy function)\n", " gridCapitalNextPeriod = 0\n", "\n", " # - Start Inner Loop - #\n", "\n", " mValueFunctionNew, mPolicyFunction = innerloop(bbeta, nGridCapital, gridCapitalNextPeriod, mOutput, nProductivity, vGridCapital, expectedValueFunction, mValueFunction, mValueFunctionNew, mPolicyFunction)\n", "\n", " # - End Inner Loop - #\n", "\n", " maxDifference = (abs(mValueFunctionNew-mValueFunction)).max()\n", "\n", " mValueFunction = mValueFunctionNew\n", " mValueFunctionNew = zeros((nGridCapital,nGridProductivity),dtype=float)\n", "\n", " iteration += 1\n", " if(iteration%10 == 0 or iteration == 1):\n", " print(\" Iteration = \", iteration, \", Sup Diff = \", maxDifference)\n", "\n", " return (maxDifference, iteration, mValueFunction, mPolicyFunction)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output = 0.5627314338711378 Capital = 0.178198287392527 Consumption = 0.3845331464786108\n", " Iteration = 1 , Sup Diff = 0.05274159340733661\n", " Iteration = 10 , Sup Diff = 0.031346949265852186\n", " Iteration = 20 , Sup Diff = 0.0187034598933572\n", " Iteration = 30 , Sup Diff = 0.01116551203397076\n", " Iteration = 40 , Sup Diff = 0.006668541708132469\n", " Iteration = 50 , Sup Diff = 0.003984292748717144\n", " Iteration = 60 , Sup Diff = 0.0023813118039327508\n", " Iteration = 70 , Sup Diff = 0.0014236586450981914\n", " Iteration = 80 , Sup Diff = 0.0008513397747206275\n", " Iteration = 90 , Sup Diff = 0.0005092051752287885\n", " Iteration = 100 , Sup Diff = 0.0003046232442147634\n", " Iteration = 110 , Sup Diff = 0.00018226485632311107\n", " Iteration = 120 , Sup Diff = 0.00010906950872624499\n", " Iteration = 130 , Sup Diff = 6.527643736298216e-05\n", " Iteration = 140 , Sup Diff = 3.907108211997912e-05\n", " Iteration = 150 , Sup Diff = 2.338807712010116e-05\n", " Iteration = 160 , Sup Diff = 1.4008644637408807e-05\n", " Iteration = 170 , Sup Diff = 8.391317202982584e-06\n", " Iteration = 180 , Sup Diff = 5.026474538039061e-06\n", " Iteration = 190 , Sup Diff = 3.0108998638755935e-06\n", " Iteration = 200 , Sup Diff = 1.8035522481030242e-06\n", " Iteration = 210 , Sup Diff = 1.0803409159487742e-06\n", " Iteration = 220 , Sup Diff = 6.471316944534067e-07\n", " Iteration = 230 , Sup Diff = 3.876361940324813e-07\n", " Iteration = 240 , Sup Diff = 2.3219657929729465e-07\n", " Iteration = 250 , Sup Diff = 1.3908720952748865e-07\n", " Iteration = 257 , Sup Duff = 9.716035664908418e-08\n", " \n", " My Check = 0.14654914369624122\n", " \n", "Elapse time = is 1.9711215496063232\n" ] } ], "source": [ "t1=time.time()\n", "# - Call Main Function - #\n", "maxDiff, iterate, mValueF, mPolicyFunction = main_func()\n", "# - End Timer - #\n", "t2 = time.time()\n", "print(\" Iteration = \", iterate, \", Sup Duff = \", maxDiff)\n", "print(\" \")\n", "print(\" My Check = \", mPolicyFunction[1000-1,3-1])\n", "print(\" \")\n", "print(\"Elapse time = is \", t2-t1)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# same as above except without print statements\n", "\n", "def main_func2():\n", "\n", " # 1. Calibration\n", "\n", " aalpha = 1.0/3.0 # Elasticity of output w.r.t. capital\n", " bbeta = 0.95 # Discount factor\n", "\n", " # Productivity values\n", " vProductivity = np.array([0.9792, 0.9896, 1.0000, 1.0106, 1.0212],float)\n", "\n", " # Transition matrix\n", " mTransition = np.array([[0.9727, 0.0273, 0.0000, 0.0000, 0.0000],\n", " [0.0041, 0.9806, 0.0153, 0.0000, 0.0000],\n", " [0.0000, 0.0082, 0.9837, 0.0082, 0.0000],\n", " [0.0000, 0.0000, 0.0153, 0.9806, 0.0041],\n", " [0.0000, 0.0000, 0.0000, 0.0273, 0.9727]],float)\n", "\n", " ## 2. Steady State\n", "\n", " capitalSteadyState = (aalpha*bbeta)**(1/(1-aalpha))\n", " outputSteadyState = capitalSteadyState**aalpha\n", " consumptionSteadyState = outputSteadyState-capitalSteadyState\n", "\n", "# print(\"Output = \", outputSteadyState, \" Capital = \", capitalSteadyState, \" Consumption = \", consumptionSteadyState)\n", "\n", " # We generate the grid of capital\n", " vGridCapital = np.arange(0.5*capitalSteadyState,1.5*capitalSteadyState,0.00001)\n", "\n", " nGridCapital = len(vGridCapital)\n", " nGridProductivity = len(vProductivity)\n", "\n", " ## 3. Required matrices and vectors\n", "\n", " mOutput = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " mValueFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " mValueFunctionNew = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " mPolicyFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", " expectedValueFunction = np.zeros((nGridCapital,nGridProductivity),dtype=float)\n", "\n", " # 4. We pre-build output for each point in the grid\n", "\n", " for nProductivity in range(nGridProductivity):\n", " mOutput[:,nProductivity] = vProductivity[nProductivity]*(vGridCapital**aalpha)\n", "\n", " ## 5. Main iteration\n", "\n", " maxDifference = 10.0\n", " tolerance = 0.0000001\n", " iteration = 0\n", "\n", " log = math.log\n", " zeros = np.zeros\n", " dot = np.dot\n", "\n", " while(maxDifference > tolerance):\n", "\n", " expectedValueFunction = dot(mValueFunction,mTransition.T)\n", "\n", " for nProductivity in range(nGridProductivity):\n", "\n", " # We start from previous choice (monotonicity of policy function)\n", " gridCapitalNextPeriod = 0\n", "\n", " # - Start Inner Loop - #\n", "\n", " mValueFunctionNew, mPolicyFunction = innerloop(bbeta, nGridCapital, gridCapitalNextPeriod, mOutput, nProductivity, vGridCapital, expectedValueFunction, mValueFunction, mValueFunctionNew, mPolicyFunction)\n", "\n", " # - End Inner Loop - #\n", "\n", " maxDifference = (abs(mValueFunctionNew-mValueFunction)).max()\n", "\n", " mValueFunction = mValueFunctionNew\n", " mValueFunctionNew = zeros((nGridCapital,nGridProductivity),dtype=float)\n", "\n", " iteration += 1\n", "# if(iteration%10 == 0 or iteration == 1):\n", "# print(\" Iteration = \", iteration, \", Sup Diff = \", maxDifference)\n", "\n", " return (maxDifference, iteration, mValueFunction, mPolicyFunction)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(9.716035664908418e-08,\n", " 257,\n", " array([[-0.9972862 , -0.98551961, -0.97408008, -0.96027188, -0.94819411],\n", " [-0.99728346, -0.98551687, -0.97407734, -0.96026914, -0.94819137],\n", " [-0.99728072, -0.98551414, -0.9740746 , -0.96026641, -0.94818864],\n", " ...,\n", " [-0.97049336, -0.95872676, -0.947286 , -0.93347903, -0.92140127],\n", " [-0.97049244, -0.95872585, -0.94728509, -0.93347812, -0.92140036],\n", " [-0.97049153, -0.95872494, -0.94728418, -0.93347721, -0.92139945]]),\n", " array([[0.13848914, 0.13996914, 0.14144914, 0.14293914, 0.14443914],\n", " [0.13849914, 0.13996914, 0.14145914, 0.14293914, 0.14443914],\n", " [0.13850914, 0.13997914, 0.14145914, 0.14294914, 0.14444914],\n", " ...,\n", " [0.19973914, 0.20185914, 0.20399914, 0.20613914, 0.20830914],\n", " [0.19973914, 0.20185914, 0.20399914, 0.20614914, 0.20830914],\n", " [0.19973914, 0.20185914, 0.20399914, 0.20614914, 0.20830914]]))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "main_func2()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.64 s ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit main_func2()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.66 s ± 42.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit main_func2()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.63 s ± 5.87 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit main_func2()" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.64 s ± 27.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit main_func2()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.63 s ± 1.58 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit main_func2()" ] }, { "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.3" } }, "nbformat": 4, "nbformat_minor": 2 }