diff --git a/examples/ConsIndShockModel/IndShockConsumerType_IntegrationMethod_Example.ipynb b/examples/ConsIndShockModel/IndShockConsumerType_IntegrationMethod_Example.ipynb new file mode 100644 index 000000000..7dff7f367 --- /dev/null +++ b/examples/ConsIndShockModel/IndShockConsumerType_IntegrationMethod_Example.ipynb @@ -0,0 +1,521 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "cfe8f06b", + "metadata": {}, + "source": [ + "# Integration Method by Matt (Example using `IndShockConsumer`)\n", + "\n", + "Constructing End of Period (Marginal) Value Function using Transitionmatrices instead of Expectations.\n", + "This notebook shows the necessary steps to construct the function. For simplicity, we will only focus on transitory shock without unemployment probability nor permanent shocks.\n", + "\n", + "Idea by Matt White\n", + "Code by Adrian Monninger" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "915c0bb6", + "metadata": {}, + "outputs": [], + "source": [ + "from HARK.interpolation import LinearInterp, MargValueFuncCRRA\n", + "from HARK.ConsumptionSaving.ConsIndShockModel import IndShockConsumerType\n", + "import numpy as np\n", + "\n", + "from HARK.distribution import (\n", + " expected,\n", + ")\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from HARK.utilities import make_grid_exp_mult\n", + "import scipy as sp" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c3cc5e1", + "metadata": {}, + "outputs": [], + "source": [ + "T_cycle = 5\n", + "PermShk = 0.0\n", + "TranShk = 0.2\n", + "ShockCount = 21\n", + "UnempPrb = 1e-10\n", + "CRRA = 2\n", + "Rfree = 1.04**0.25\n", + "DiscFac = 0.975\n", + "PermGroFac = 1.0\n", + "LivPrb = 0.99375\n", + "DiscFacEff = DiscFac * LivPrb\n", + "\n", + "Dict = {\n", + " \"CRRA\": CRRA,\n", + " \"Rfree\": [Rfree] * T_cycle, # Interest factor on assets\n", + " \"DiscFac\": DiscFac, # Intertemporal discount factor\n", + " \"LivPrb\": [LivPrb]\n", + " * T_cycle, # [0.999], #[.99375], # Survival probability\n", + " \"PermGroFac\": [PermGroFac] * T_cycle, # Permanent income growth factor\n", + " # Parameters that specify the income distribution over the lifecycle\n", + " \"PermShkStd\": [PermShk]\n", + " * T_cycle, # [.06], # Standard deviation of log permanent shocks to income\n", + " \"PermShkCount\": ShockCount, # Number of points in discrete approximation to permanent income shocks\n", + " \"TranShkStd\": [TranShk]\n", + " * T_cycle, # [.2], # Standard deviation of log transitory shocks to income\n", + " \"TranShkCount\": ShockCount, # Number of points in discrete approximation to transitory income shocks\n", + " \"UnempPrb\": [UnempPrb] * T_cycle, # Probability of unemployment while working\n", + " \"IncUnemp\": [0.0] * T_cycle, # Unemployment benefits replacement rate\n", + " \"UnempPrbRet\": [0.0] * T_cycle, # Probability of \"unemployment\" while retired\n", + " \"IncUnempRet\": [0.0] * T_cycle, # \"Unemployment\" benefits when retired\n", + " \"T_cycle\": T_cycle, # Number of periods in the cycle for this agent type\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e7b6ec8", + "metadata": {}, + "outputs": [], + "source": [ + "IndShock = IndShockConsumerType(**Dict)\n", + "IndShock.update()\n", + "IndShock.solve()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4e55f42", + "metadata": {}, + "outputs": [], + "source": [ + "t = 0\n", + "tNext = 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b893f23", + "metadata": {}, + "outputs": [], + "source": [ + "### Get Last Period Marginal Value Function and grids\n", + "# cFuncNext = IndShock.solution[1].cFunc\n", + "# vPfuncNext = MargValueFuncCRRA(cFuncNext, CRRA)\n", + "vPfuncNext = IndShock.solution[tNext].vPfunc\n", + "aXtraGrid = IndShock.aXtraGrid\n", + "mNrmMinNext = IndShock.solution[tNext].mNrmMin\n", + "\n", + "### Shocks\n", + "IncShkDstn = IndShock.IncShkDstn[t]\n", + "ShkPrbsNext = IncShkDstn.pmv\n", + "PermShkValsNext = IncShkDstn.atoms[0]\n", + "TranShkValsNext = IncShkDstn.atoms[1]\n", + "PermShkMinNext = np.min(PermShkValsNext)\n", + "TranShkMinNext = np.min(TranShkValsNext)\n", + "\n", + "# Borrowing Constraint\n", + "BoroCnstNat = (\n", + " (IndShock.solution[tNext].mNrmMin - TranShkMinNext)\n", + " * (IndShock.PermGroFac[t] * PermShkMinNext)\n", + " / IndShock.Rfree[t]\n", + ")\n", + "\n", + "aNrmNow = np.asarray(aXtraGrid) + BoroCnstNat" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "09419231", + "metadata": {}, + "outputs": [], + "source": [ + "def m_nrm_next(PermGroFac, shocks, a_nrm, Rfree):\n", + " \"\"\"\n", + " Computes normalized market resources of the next period\n", + " from income shocks and current normalized market resources.\n", + "\n", + " Parameters\n", + " ----------\n", + " shocks: [float]\n", + " Permanent and transitory income shock levels.\n", + " a_nrm: float\n", + " Normalized market assets this period\n", + "\n", + " Returns\n", + " -------\n", + " float\n", + " normalized market resources in the next period\n", + " \"\"\"\n", + " return Rfree / (PermGroFac * shocks[\"PermShk\"]) * a_nrm + shocks[\"TranShk\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47da57d1", + "metadata": {}, + "outputs": [], + "source": [ + "### End-of-period Marginal Value Function (Traditional Way)\n", + "\n", + "\n", + "def vp_next(shocks, a_nrm, Rfree):\n", + " return shocks[\"PermShk\"] ** (-CRRA) * vPfuncNext(\n", + " m_nrm_next(PermGroFac, shocks, a_nrm, Rfree)\n", + " )\n", + "\n", + "\n", + "EndOfPrdvP = (\n", + " DiscFacEff\n", + " * Rfree\n", + " * PermGroFac ** (-CRRA)\n", + " * expected(vp_next, IncShkDstn, args=(aNrmNow, Rfree))\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a92737a4", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(aNrmNow, EndOfPrdvP)\n", + "plt.title(\"End Of Period Marginal Value Function\")\n", + "plt.ylim([0.0, 2.0])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c7c25550", + "metadata": {}, + "source": [ + "### Alternative Start" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36bdd9ec", + "metadata": {}, + "outputs": [], + "source": [ + "### Step 1: In a pre-solution step, specify fairly dense grids of m and b, denser than the grid of a.\n", + "mGrid = make_grid_exp_mult(0.0001, 52, 310, timestonest=1)\n", + "bGrid = make_grid_exp_mult(0.0, 48, 300, timestonest=1)\n", + "\n", + "### Let ist start from minimum of aNrm\n", + "# NO, WRONG MINIMUM!\n", + "mGrid = mGrid + mNrmMinNext\n", + "bGrid = bGrid + aNrmNow[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "475cc722", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(mGrid, vPfuncNext(mGrid))\n", + "plt.ylim([0.0, 2.0])\n", + "plt.title(\"Next period marginal value function\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "282a0419", + "metadata": {}, + "outputs": [], + "source": [ + "### Step 2: In that pre-solution step, for each $F_{\\theta, t}$ (i.e. just once if this is an ifinite horizon problem),\n", + "### calculate $f_{\\theta, t}(m - b) on the cross product of the b-grid and m-grid.\n", + "### Put them in a matrix that has b values by row and m values by column.\n", + "### Still in the pre-solution step, for each b-row in the matrix, take the row-wise sum and\n", + "### divide the row by it. This represents a Markov transition matrix from the exogenous\n", + "### b-grid to the exogenous m-grid.\n", + "\n", + "\n", + "def updateBM_TranMatrix(TranShkStd, bGrid, mGrid):\n", + " \"\"\"\n", + " Calculates the Probabillity of a transioty shock for the b times m matrix\n", + " \"\"\"\n", + " # probGrid = np.zeros((len(self.bNrmGrid_income), len(self.mNrmGrid_income))) # b x m\n", + " # ### getting the probability for each transitory shock with the size b - m (remember m = b * transitory shock)\n", + "\n", + " ### Integration 1: No unemployment probability\n", + " s = TranShkStd\n", + " mu = -0.5 * s**2\n", + " lognorm_dist = sp.stats.lognorm(s, scale=np.exp(mu))\n", + "\n", + " ### Create matrix\n", + " # Construct meshgrid of bNrmGrid_income and mNrmGrid_income\n", + " b, m = np.meshgrid(bGrid, mGrid, indexing=\"ij\")\n", + "\n", + " # Calculate differences between corresponding elements\n", + " probGrid = lognorm_dist.pdf(m - b)\n", + "\n", + " for i_b in range(len(bGrid)):\n", + " probGrid[i_b] = probGrid[i_b] / (np.max([np.sum(probGrid[i_b]), 0.000001]))\n", + "\n", + " return probGrid" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93f258bd", + "metadata": {}, + "outputs": [], + "source": [ + "probGrid = updateBM_TranMatrix(TranShk, bGrid, mGrid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7aa79e03-4449-43f1-bfc0-7216143b0d37", + "metadata": {}, + "outputs": [], + "source": [ + "probGrid.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f53a5d71", + "metadata": {}, + "outputs": [], + "source": [ + "### Step 3: To compute expectations during backwards solution, evaluate vt(mt) on the dense m-grid. Put\n", + "### these (marginal) values into a matrix that has one column.\n", + "vPnext_array = vPfuncNext(mGrid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9642f501", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(mGrid, vPnext_array)\n", + "plt.title(\"Next period marginal value function\")\n", + "plt.xlabel(\"m Next grid\")\n", + "plt.ylim([0.0, 2.0])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69bb8a13", + "metadata": {}, + "outputs": [], + "source": [ + "### Step 4: Pre-multiply the Markov matrix by the matrix you just constructed. The resulting\n", + "### matrix of intermediate expected values will have b row-wise.\n", + "\n", + "vPnext_array_reshaped = np.reshape(vPnext_array, (vPnext_array.size, 1))\n", + "Interm_vP = np.dot(probGrid, vPnext_array_reshaped)\n", + "Interm_vPnvrs = Interm_vP.flatten() ** (-1 / CRRA)\n", + "Interm_vPnvrsFunc = LinearInterp(\n", + " np.insert(bGrid, 0, 0.0), np.insert(Interm_vPnvrs, 0, 0.0)\n", + ") # IMPORTANT: GENERALIZE THIS ZERO\n", + "Interm_vPfunc = MargValueFuncCRRA(Interm_vPnvrsFunc, CRRA)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "913ba280", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(bGrid, Interm_vPfunc(bGrid))\n", + "plt.title(\"Intermediate Marginal Value Function\")\n", + "plt.xlabel(\"b Next Grid\")\n", + "plt.ylim([0.0, 2.0])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56c88261", + "metadata": {}, + "outputs": [], + "source": [ + "### Step 5: Using the pre-specified exogenous grid of a, compute expected end-of-period\n", + "### (marginal) value using the typical discretized approximation to the permanent shock\n", + "### distribution. For any (a) end-of-period state, there is a finite set of future (b)\n", + "### points that can be reached under the discretization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b7a487c", + "metadata": {}, + "outputs": [], + "source": [ + "def b_nrm_next(PermGroFac, shocks, a_nrm, Rfree):\n", + " \"\"\"\n", + " Computes normalized bank balances of the next period\n", + " from income shocks and current normalized market resources.\n", + "\n", + " Parameters\n", + " ----------\n", + " shocks: [float]\n", + " Permanent and transitory income shock levels.\n", + " a_nrm: float\n", + " Normalized market assets this period\n", + "\n", + " Returns\n", + " -------\n", + " float\n", + " normalized market resources in the next period\n", + " \"\"\"\n", + " return Rfree / (PermGroFac * shocks[\"PermShk\"]) * a_nrm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "651ab6ae", + "metadata": {}, + "outputs": [], + "source": [ + "def vp_next_Integration(shocks, a_nrm, Rfree):\n", + " return shocks[\"PermShk\"] ** (-CRRA) * Interm_vPfunc(\n", + " b_nrm_next(PermGroFac, shocks, a_nrm, Rfree)\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "036c3a61", + "metadata": {}, + "outputs": [], + "source": [ + "EndOfPrdvP_Integration = (\n", + " DiscFacEff\n", + " * Rfree\n", + " * PermGroFac ** (-CRRA)\n", + " * expected(vp_next_Integration, IncShkDstn, args=(aNrmNow, Rfree))\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "325313b3", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(aNrmNow, EndOfPrdvP)\n", + "plt.plot(aNrmNow, EndOfPrdvP_Integration, \"--\")\n", + "plt.title(\"End Of Period Marginal Value Functions: Comparison\")\n", + "plt.xlabel(\"a Grid\")\n", + "plt.ylim([0.0, 2.0])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a70cd24a-6184-4a29-8040-93132a43892c", + "metadata": {}, + "outputs": [], + "source": [ + "# Construct the consumption function in two ways: standard approach and new approach\n", + "cNrm_old = EndOfPrdvP ** (-1 / CRRA)\n", + "mNrm_old = aNrmNow + cNrm_old\n", + "cFunc_old = LinearInterp(np.insert(mNrm_old, 0, 0.0), np.insert(cNrm_old, 0, 0.0))\n", + "\n", + "cNrm_new = EndOfPrdvP_Integration ** (-1 / CRRA)\n", + "mNrm_new = aNrmNow + cNrm_new\n", + "cFunc_new = LinearInterp(np.insert(mNrm_new, 0, 0.0), np.insert(cNrm_new, 0, 0.0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e67ea4ac-1fff-4edb-af84-36bf50eb7e9e", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot a comparison of the two consumption functions using different integration methods\n", + "M = np.linspace(0.0, 10.0, 201)\n", + "C0 = cFunc_old(M)\n", + "C1 = cFunc_new(M)\n", + "plt.plot(M, C0)\n", + "plt.plot(M, C1)\n", + "plt.xlabel(r\"Market resources $m_t$\")\n", + "plt.ylabel(r\"Consumption $c_t$\")\n", + "plt.xlim(0.0, 10.0)\n", + "plt.ylim(0.0, None)\n", + "plt.title(\"Consumption function comparison\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05ad7a58-0edf-4d04-837d-c27953508e58", + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(M, C1 - C0)\n", + "plt.xlabel(r\"Market resources $m_t$\")\n", + "plt.ylabel(r\"Difference between methods\")\n", + "plt.title(\"Consumption function comparison\")\n", + "plt.xlim(0.0, 10.0)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fc316d2-5b6b-43c4-ad53-3dfd47c92ba1", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}