Skip to content

Instantly share code, notes, and snippets.

@chriselrod

chriselrod/quant_econ_julia.ipynb

Created February 12, 2018 00:08
Show Gist options
  • Save chriselrod/f30fdfadccd8413000f0703a7cc6f8a5 to your computer and use it in GitHub Desktop.
Save chriselrod/f30fdfadccd8413000f0703a7cc6f8a5 to your computer and use it in GitHub Desktop.
Download ZIP
Quant Econ - Numba and Julia
Display the source blob
Display the rendered blob
Raw
quant_econ_julia.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
quant_econ_numba.ipynb
{
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment