From 57738b31fc94946c16882b511fc08805e9446464 Mon Sep 17 00:00:00 2001 From: JohnRGreen Date: Wed, 10 May 2023 14:46:37 -0400 Subject: [PATCH] Semester project --- REMARKs/HumanCapital_EGM.ipynb | 1274 ++++++++++++++++++++++++++++++++ 1 file changed, 1274 insertions(+) create mode 100644 REMARKs/HumanCapital_EGM.ipynb diff --git a/REMARKs/HumanCapital_EGM.ipynb b/REMARKs/HumanCapital_EGM.ipynb new file mode 100644 index 0000000..88a8d2c --- /dev/null +++ b/REMARKs/HumanCapital_EGM.ipynb @@ -0,0 +1,1274 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Implementing the EGM in a Model of Human Capital\n", + "\n", + " John Green
May 2023 \n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this notebook I will demonstrate the solution method for my second year paper, \"Skills-Biased Technological Change and the Human Capital Production Function.\" In the paper I estimate a model of life-cycle labor supply and consumption with human capital similar to that of Imai and Keane (2004). I estimate this on data from 1979 and 1997, and compare the parameters for the human capital production function (HCPF) in order to test the skills-biased technological change (SBTC) hypothesis. I solve the model using the Endogenous Grid Method of Carroll (2006), which provides a significant speedup (by about a factor of 50) over using backwards-recursion of the value function. The gain in computation time is in addition to a gain in accuracy, since the value function iteration seems prone to getting stuck at local minima. However, because of the ``twisted grid\" problem which occurs in the endogenous grid step, implementation is not straightforward.\n", + "\n", + "I will first present the model, then show my first order conditions and give some intuition for the solution, and then demonstrate how to use the EGM to solve the agent's problem. Full details as well as the results of the estimation are provided in the full paper, which is available upon request." + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### The model\n", + "In each annual period agents choose their consumption and their labor supply. Agents derive current-period utility from consumption and disutility from labor, following the classic additive CRRA form, as used in MaCurdy (1981) and Altonji (1986). Utility from consumption is:\n", + "\n", + "\\begin{equation}\n", + "\\tag{1}\n", + "u(c) = \\dfrac{c^{1-\\gamma}}{1 - \\gamma}\n", + "\\end{equation}\n", + "\n", + "Disutility from labor is:\n", + "\\begin{equation}\n", + "\\tag{2}\n", + "v(h) = b \\dfrac{h^{1+ \\eta}}{1 + \\eta}\n", + "\\end{equation}\n", + "\n", + "Labor income is determined by the number of hours worked as well as the level of human capital, which determines workers' wages:\n", + "\\begin{equation}\n", + "\\tag{3}\n", + "Y_t = K_t h_t\n", + "\\end{equation}\n", + "\n", + "Assets are determined by the intertemporal budget constraint:\n", + "\\begin{equation}\n", + "\\tag{4}\n", + "A_{t+1} = \\left( A_t + Y_t - c_t \\right) R\n", + "\\end{equation}\n", + "Where $R$ is the risk-free interest rate and is set to 1.05.\n", + "\n", + "Human capital is accumulated by learning-by-doing. That is, agents increase their next period human capital by working in the current period. At the same time, the existing human capital stock depreciates. There is also a stochastic shock. The deterministic portion of human capital is:\n", + "\\begin{equation}\n", + "\\tag{5}\n", + "K_{\\bar{t}} = \\delta K_t + G(h_t,s) K_t\n", + "\\end{equation}\n", + "\n", + "where $G(\\cdot,\\cdot)$ is a function increasing in hours worked and decreasing in age ($h$ and $s$ respectively). The notation $K_{\\bar{t}}$ indicates the level of human capital at the *end* of period $t$, after the control variables have been realized but before the stochastic shock. Then human capital at the beginning of the next period is:\n", + "\\begin{equation}\n", + "\\tag{6}\n", + "K_{\\underline{t+1}} = \\epsilon_{t+1} K_{\\bar{t}}\n", + "\\end{equation}\n", + "\n", + "This distinction between end-of-period deterministic human capital and the beginning-period human capital will be useful when we implement the EGM. \n", + "\n", + "The HCPF $G(\\cdot,\\cdot)$ is increasing in hours worked and decreasing in age. It is specified as:\n", + "\\begin{equation}\n", + "\\tag{7}\n", + "G(h,t) = \\omega_0(1 - \\omega_1(s-19)) \\left[ (h+d)^{\\alpha} - \\omega_2(h+d)\\right]\n", + "\\end{equation}\n", + "\n", + "This form is similar to that of Imai and Keane (2004), and is meant to capture several features of the data:\n", + "* Wages grow at a slower rate as workers get older, captured by $-\\omega_1(s-19)$\n", + "* Returns to hours worked are not linear, which is why we include $\\alpha<1$\n", + "* At high levels of hours, returns may be negative, which is reflected by $-\\omega_2(h+d)$\n", + "* $d$ captures a lower-bound to wages\n", + "\n", + "The lower bound is likely an artifact of the minimum wage, and in future work I would like to model this explicitly.\n", + "\n", + "Putting all of this together, agents' problem can be written in the Bellman equation form:\n", + "\\begin{equation}\n", + "\\tag{8}\n", + " \\begin{split}\n", + " & V_t(A_t,K_t) = \\text{max}_{c_t,h_t} \\left\\{ u(c_t) - v(h_t) + \\beta \\mathbb{E}_t V_{t+1} \\left(A_{t+1}, K_{t+1} \\right) \\right\\} \\\\\n", + " & = \\text{max}_{c_t,h_t} \\left\\{ u(c_t) - v(h_t) + \n", + " \\beta \\mathbb{E}_t V_{t+1} \\left( R(A_t + K_t h_t - c_t), \\epsilon_{t+1}( \\delta K_{\\bar{t}} + G(h_t,s) K_{\\bar{t}}) \\right) \\right\\}\n", + " \\end{split}\n", + "\\end{equation}\n", + "\n", + "Where $\\beta$ gives the agent's discount factor." + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### First order conditions\n", + "\n", + "The solution method will take advantage of the problem's first order conditions. The FOC with respect to consumption is:\n", + "\\begin{equation}\n", + "\\tag{9}\n", + "u'(c_t) = \\beta R \\mathbb{E} \\left[ u'(c_{t+1}) \\right]\n", + "\\end{equation}\n", + "which is the classic consumption Euler equation.\n", + "\n", + "The FOC with respect to labor supply is slightly more complicated:\n", + "\\begin{equation}\n", + "\\tag{10}\n", + "v'(h_t) = K_t u'(c_t) + \\beta \\kappa(h_t, K_t, s) \\mathbb{E} \\left[ \\epsilon_t \\dfrac{\\partial V_{t+1}}{\\partial K_{t+1}} \\right]\n", + "\\end{equation}\n", + "\n", + "$\\frac{\\partial V_{t}}{\\partial K_{t}}$ takes the form:\n", + "\\begin{equation}\n", + "\\tag{11}\n", + "\\dfrac{\\partial V_{t}}{\\partial K_{t}} = h_t u'(c_{t}) + \\beta \\lambda (h_t,s) \\mathbb{E} \\left[ \\epsilon_{t+1} \\dfrac{\\partial V_{t+1}}{\\partial K_{t+1}} \\right]\n", + "\\end{equation}\n", + "\n", + "Note that in our terminal period $T$, the RHS term $\\mathbb{E} \\left[ \\epsilon_{t+1} \\frac{\\partial V_{t+1}}{\\partial K_{t+1}} \\right] = 0$, which will be useful.\n", + "\n", + "Finally, we will need to invert the law of motion for human capital. Define:\n", + "\\begin{equation}\n", + " \\tag{12}\n", + " \\begin{split}\n", + " K_{\\underline{t}} & = K_{\\bar{t}}( \\delta + \\omega_0(1 - \\omega_1(s-19))[(h_t + d)^{\\alpha} - \\omega_2(h_t+d)])^{-1} \\\\\n", + " & \\equiv \\Psi(K_{\\bar{t}}, h_t,s)\n", + " \\end{split}\n", + "\\end{equation}\n", + "\n", + "So that eqn. (10) can be rewritten as:\n", + "\\begin{equation}\n", + "\\tag{13}\n", + "v'(h_t) = \\Psi(K_{\\bar{t}}, h_t,s) u'(c_t) + \\beta \\kappa(h_t, \\Psi(K_{\\bar{t}}, h_t, s), s) \\mathbb{E} \\left[ \\epsilon_t \\dfrac{\\partial V_{t+1}}{\\partial K_{t+1}} \\right]\n", + "\\end{equation}" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Implementing the EGM\n", + "\n", + "I will now solve the agent's problem using the EGM, in order to get an approximation of the policy function for consumption and labor supply. First, some basic setup." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# Import necessary packages here\n", + "import math as mat\n", + "import numpy as np\n", + "import scipy\n", + "import pandas as pd\n", + "import xarray as xr\n", + "from scipy.optimize import root_scalar\n", + "from scipy import interpolate\n", + "from numpy import vectorize\n", + "from scipy.optimize import minimize_scalar\n", + "from scipy.interpolate import griddata\n", + "from scipy.interpolate import LinearNDInterpolator\n", + "import statsmodels.api as sm\n", + "from scipy.interpolate import interp2d\n", + "from scipy.interpolate import RegularGridInterpolator\n", + "\n", + "# Now declare the global variables\n", + "# risk-free rate\n", + "r = 0.05\n", + "R = 1 + r\n", + "# Number of grid points for \n", + "nA = 20 # assets\n", + "nK = 20 # human capital\n", + "nq = 6 # quadrature integration\n", + "# minimum consumption\n", + "cmin = 1e-5\n", + "# and maximum work\n", + "hmax = 5096\n", + "# Agents live from age 20 to 65\n", + "T = 65\n", + "t0 = 20\n", + "lifespan = T- t0\n", + "\n", + "# Preset variables to be estimated\n", + "# note: in the paper I estimate for 4 separate education groups. Here I will use the values for college graduates.\n", + "\n", + "# utility function parameters\n", + "gamma = 0.68 # curvature of utility from consumption\n", + "eta = 0.87 # curvature of disutility from labor\n", + "b = 5.1e-5 # scale parameter on labor disutility\n", + "beta = .986\n", + "# HCPF parameters\n", + "delta = .303\n", + "omega0 = .13\n", + "omega1 = 1e-3\n", + "alpha = .24\n", + "omega2 = 4.1e-4\n", + "d=512\n", + "# sd of wage shocks\n", + "sigma = 8.06e-2" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now I will define some of the functions from above, both from the structural model and those that are used in calculating the FOCs and in other intermediate steps." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Utility functions:\n", + "def utility(c,gamma):\n", + " if c < cmin:\n", + " # illegal value of c; give a large negative utility\n", + " print(\"invalid value of c\")\n", + " return -10000\n", + " elif c >= cmin and gamma != 1:\n", + " return (c**(1 - gamma))/(1-gamma)\n", + " else:\n", + " # case of log utility\n", + " return mat.log(c)\n", + " \n", + "def getMargUtilityC(c,gamma):\n", + " if c < cmin:\n", + " # illegal value of c; give an error\n", + " # print(\"invalid value of c\")\n", + " return 999\n", + " elif c >= cmin and gamma != 1:\n", + " # print(\"normal value of c\")\n", + " return c**(-gamma)\n", + " else:\n", + " # case of log utility\n", + " # print(\"log utility\")\n", + " return 1/c\n", + "\n", + "def disutility(h,b,eta):\n", + " return b*(h**(1+eta))/(1+eta)\n", + "\n", + "def getMargDisutilityH(h,b,eta):\n", + " return b*(h**eta)\n", + "\n", + "def gFunc(h,K,age,delta,omega0,omega1,alpha,omega2,d):\n", + " G = omega0*(1+omega1*(age-19))*((h+d)**alpha - omega2*(h+d))\n", + " return delta*K + G*K\n", + "\n", + "# Now define functions relating to the FOCs:\n", + "# Tau\n", + "def taufunc(age,omega0,omega1):\n", + " return omega0*(1+omega1*(age-19))\n", + "\n", + "# Psi (inverse law of motion for human capital)\n", + "def psifunc(Kend, h, this_tau, d, alpha, omega2, delta):\n", + " this_chi = ((h+d)**alpha - omega2*(h+d))\n", + " return Kend/(delta + this_tau*this_chi)\n", + "\n", + "# lambda\n", + "def lambdafunc(h, this_tau, omega2, d, alpha, delta):\n", + " return delta + this_tau*((h + d)**alpha - omega2*(h+d))\n", + "\n", + "# kappa\n", + "def kappafunc(Kstart, h, this_tau, d, omega2, alpha):\n", + " return this_tau*(alpha*((h+d))**(alpha-1)-omega2)*Kstart" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our final piece of setup is to establish the exogenous grid of state variables, assets and human capital, over which we would like to get our policy function approximation. First, I define a function to calculate the grids, given minimum and maximum values, number of points, and a method for the steps." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Write a function to get the grids\n", + "def GetGrid(min, max, points, method):\n", + "\n", + " # MMake sure min and max are arrays of the same length\n", + " if len(min) != len(max):\n", + " print(\"Error: min and max must be arrays of the same length.\")\n", + " exit()\n", + "\n", + " # get span in each period\n", + " span = max - min\n", + "\n", + " # Get the grid using log steps\n", + " if method == \"log\":\n", + " loggrid = np.full((points, len(span)), np.nan)\n", + " for i in range(len(min)):\n", + " vec = np.linspace(np.log(1), np.log(1 + span[i]), points)\n", + " loggrid[:, i] = vec\n", + "\n", + " # then exponentiate\n", + " grid = np.exp(loggrid) - 1\n", + "\n", + " # add back in the min\n", + " min_mat = np.full((points, len(span)), min)\n", + " grid = grid + min_mat\n", + "\n", + " # grid[1,] == asset_lb\n", + " # round(grid[points,], 2) == round(asset_ub, 2)\n", + " # all good\n", + " elif method == \"linear\": # if we want to do linear interpolation\n", + " grid = np.full((points, len(span)), np.nan)\n", + " for i in range(len(min)):\n", + " vec = np.linspace(min[i], max[i], points)\n", + " grid[:, i] = vec\n", + " else:\n", + " print(\"error: method invalid. Choose 'log' or 'linear'\")\n", + "\n", + " return grid" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I set \\$4 as the minimum value of the human capital grid. This is slightly below the minimum wage in the time period we are studying, allowing for some below-minimum earnings (due, for example, to black market work). I set the minimum as \\$30 in the first period, which increases to \\$100 by period $T$. This is arbitrary, but should capture well the range of earnings we see in the data." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "kmin = np.repeat([4],lifespan)\n", + "kmax = np.linspace(30, 100, 45)\n", + "k_grid = GetGrid(kmin,kmax,nK,\"linear\")\n", + "# k_grid" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I impose some mild liquidity constraints. Individuals are allowed to have negative assets up to twice their full-time earnings in all periods except for $T-1$, when they are allowed to have one year of full time earnings. Thus the lower bound for assets is defined by the borrowing constraint for those with the highest level of human capital on our grid. The upper bound of assets is set to \\$250 thousand; this is arbitrary, but captures most of the observations in the data." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# minimum level of assets on our grid \n", + "amin = np.full(lifespan, np.nan)\n", + "amin[len(amin)-1] = -k_grid[nK-1,len(amin)-1]*2040\n", + "for i in (range(len(amin)-1)):\n", + " amin[i] = -k_grid[nK-1,i]*2040*2\n", + "amax = np.repeat([250000],lifespan)\n", + "a_grid = np.ascontiguousarray(GetGrid(amin, amax, nA, \"linear\"))\n", + "# a_grid[0,:]" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are two final pieces of setup. First, we need to discretize the distribution for our wage shocks, and then we need establish the arrays and vectors we will use in our solution." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "def discretize_log_distribution(points, mean, sd):\n", + " var = sd ** 2\n", + " Z = np.full(points + 1, np.nan)\n", + " EVbetweenZ = np.full(points, np.nan)\n", + "\n", + " # manually set upper bounds at +/- 3\n", + " Z[0] = -3\n", + " Z[points] = 3\n", + "\n", + " # Get the cutoff points\n", + " prob = 1 / points\n", + " for i in range(points - 1, 0, -1):\n", + " Z[i] = scipy.stats.norm.ppf(prob * i)\n", + "\n", + " # get mean value in each interval (see Adda & Cooper page 58)\n", + " for i in range(points):\n", + " ev = (scipy.stats.norm.pdf(Z[i]) - scipy.stats.norm.pdf(Z[i + 1])) / (scipy.stats.norm.cdf(Z[i + 1]) - scipy.stats.norm.cdf(Z[i]))\n", + " pt = sd * ev + mean\n", + " EVbetweenZ[i] = pt\n", + "\n", + " quad2 = np.exp(EVbetweenZ)\n", + " return quad2\n", + "\n", + "quad = discretize_log_distribution(nq, -(1/2)*sigma**2, sigma)\n", + "\n", + "# Now set up empty matrices and vectors that we will use\n", + " # set up empty matrices \n", + "policyC = np.full((nA,nK,lifespan), np.nan)\n", + "policyH = np.full((nA,nK,lifespan), np.nan)\n", + "dVdK = np.full((nA,nK,lifespan), np.nan)\n", + "E_dVdK = np.full((nA,nK,lifespan), np.nan)\n", + "E_eps_dVdK= np.full((nA,nK,lifespan), np.nan)\n", + "uPrime = np.full((nA,nK,lifespan), np.nan)\n", + "E_uPrime = np.full((nA,nK,lifespan), np.nan)\n", + "lc = np.full((nA,nK,lifespan), np.nan)\n", + "extrap = np.full((nA,nK,lifespan), np.nan)\n", + "replace = np.full((nA,nK,lifespan-1), np.nan)\n", + " \n", + "# vectors to store input for each E_{} value\n", + "Uprime_nextvec = np.full(nq, np.nan) \n", + "dVdK_nextvec = np.full(nq, np.nan) \n", + "eps_dVdK_nextvec = np.full(nq, np.nan) \n", + "E_dVdK_nextvec = np.full(nq, np.nan) \n", + " \n", + "# In periods T-1 and prior, we will not have a grid, but a data frame of \n", + "# nA*nK points\n", + "# rows for R^{-1}A_endofperiod, K_endofperiod, policyC, policyH,\n", + "# uPrime, dVdK, E_dVdK, E_eps_dVdK, E_uPrime\n", + "policy_df = np.full((nA*nK,13,lifespan-1), np.nan)\n", + "\n", + "# Get minimum utility and marginal value at the constraint\n", + "umin = utility(cmin,gamma) - disutility(hmax,b,eta);\n", + "uprimemin = getMargUtilityC(cmin,gamma); " + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now have our grids and all of our variables set, and are ready to solve our problem. First we will solve the problem in period $T$. \n", + "\n", + "#### Period $T$\n", + "\n", + "The agent's problem ends in period $T$, and for now we use no terminal value function. Thus the agent will consume all of their total resources (assets and labor income), and in practice only needs to choose the optimal labor supply (in the future I plan to implement a terminal period value function, representing retirement, bequest motive, etc.).\n", + "\n", + "We will first set up a function from eqn. (10) which we can use to solve for $h_T^*$. Recalling that the second term on the RHS of (10) is equal to 0 since there is no value from human capital past period $T$, this is:\n", + "\n", + "\\begin{equation}\n", + "\\tag{14}\n", + "v'(h_T) - K_T u'(A_T + h_T K_T) = h_T^{\\eta} - K_T (A_T + h_T K_T)^{-\\gamma} = 0\n", + "\\end{equation}\n", + "\n", + "Unfortunately, eqn. (14) has no analytical solution, so we will use a root-finder. We loop through every point on our grid, plug $A$ and $K$ into eqn. (14), and then get the optimal labor supply. We will first check whether the agent is liquidity constrained." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# Define function for eqn. (14)\n", + "def hFuncTerminal(h,A,K,b,gamma,eta):\n", + " c = A + h*K\n", + " dis = getMargDisutilityH(h,b,eta)\n", + " util = getMargUtilityC(c,gamma)\n", + " return dis - util*K\n", + "\n", + "for ixa in range(nA): # looping through asset points\n", + " # First, find the optimal labor supply for each K grid point at this A\n", + " for ixk in range(nK): # looping through human capital points\n", + " ub = hmax*k_grid[ixk,lifespan-1]+a_grid[ixa,lifespan-1]\n", + " # print(ub)\n", + " if ub < cmin:\n", + " # liquidity constrained\n", + " policyC[ixa,ixk,lifespan-1] = cmin\n", + " policyH[ixa,ixk,lifespan-1] = hmax\n", + " lc[ixa,ixk,lifespan-1] = 1\n", + " \n", + " # fill in uPrime and dVdK\n", + " uPrime[ixa,ixk,lifespan-1] = uprimemin\n", + " dVdK[ixa,ixk,lifespan-1] = hmax*uprimemin\n", + " else:\n", + " # not liquidity constrained\n", + " # get minimum hours\n", + " hmin = (2*cmin-a_grid[ixa,lifespan-1])/k_grid[ixk,lifespan-1]+.1\n", + " hmin = hmin if hmin > 0 else 0.1\n", + " \n", + " p1 = hFuncTerminal(hmin,a_grid[ixa,lifespan-1],k_grid[ixk,lifespan-1],b,gamma,eta)\n", + " p2 = hFuncTerminal(hmax,a_grid[ixa,lifespan-1],k_grid[ixk,lifespan-1],b,gamma,eta)\n", + "\n", + " if p1*p2>0: # either both positive or both negative\n", + " h = hmin if abs(p1)= k_grid[nK-1,lifespan-1], k_grid[nK-1,lifespan-1]-0.01, nextK)\n", + "\n", + " f = interpolate.interp1d(k_grid[:,lifespan-1], uPrime[ixa,:,lifespan-1])\n", + " Uprime_nextvec = f(nextK)\n", + " f = interpolate.interp1d(k_grid[:,lifespan-1], dVdK[ixa,:,lifespan-1])\n", + " dVdK_nextvec = f(nextK)\n", + " for i in range(len(dVdK_nextvec)):\n", + " eps_dVdK_nextvec[i] = quad[i] * dVdK_nextvec[i]\n", + " \n", + " E_uPrime[ixa,ixk,lifespan-1] = np.mean(Uprime_nextvec)\n", + " E_dVdK[ixa,ixk,lifespan-1] = np.mean(dVdK_nextvec)\n", + " E_eps_dVdK[ixa,ixk,lifespan-1] = np.mean(eps_dVdK_nextvec)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Write out this last step in markdown." + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now have the problem in period $T$ solved, and can move on to period $T-1$.\n", + "\n", + "#### Period $T-1$\n", + "\n", + "As discussed, our use of the EGM will result in the ``twisted grid\" problem. In order to fill in the exogenous grid from our endogenous grid, we will use a triangulation based on our endogenous values. One issue that we will confront is the necessity of filling in values of our exogenous grid not covered by our endogenous values. This is required, because if we simply drop grid points outside of the endogenous values, our grid will continuously narrow. In order to do this, I use a simple polynomial approximation, where the policy function is a function of $A$, $K$, $A\\times K$, $A^2$, and $K^2$. This is not ideal, as there are well known shortcomings to using a polynomial approximation. The approximation may be very accurate on the domain of the endogenous values ($R^2 > 0.99$), outside of the domain it may not be accurate. A better solution is needed in the future.\n", + "\n", + "First, I set up a data frame for the polynomial function." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# set up the polynomial df\n", + "polynomial_df = np.empty((nA*nK, 5, lifespan-1))\n", + "for ixt in range(0, lifespan-1):\n", + " polynomial_df[:, 0, ixt] = np.tile(a_grid[:, ixt], len(k_grid[:, ixt]))\n", + " polynomial_df[:, 1, ixt] = np.kron(k_grid[:, ixt], np.repeat(1, len(a_grid[:, ixt]))) # x2\n", + " polynomial_df[:, 2, ixt] = polynomial_df[:, 0, ixt] * polynomial_df[:, 1, ixt] # x3\n", + " polynomial_df[:, 3, ixt] = (polynomial_df[:, 0, ixt]**2) / 1e5 # x4\n", + " polynomial_df[:, 4, ixt] = polynomial_df[:, 1, ixt]**2 # x5" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we enter period $T-1$ , first melting our values from the beginning of period $T$ to be our values at the end of period $T-1$" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# set our time indicator\n", + "ixt = lifespan-1-1\n", + "\n", + "# end of period assets\n", + "policy_df[:,0,ixt] = np.tile(a_grid[:,ixt+1]*(R**(-1)),len(k_grid[:,ixt+1])) # unwind the interest payment to get assets at end of T-1\n", + "\n", + "# end of period human capital\n", + "policy_df[:,1,ixt] = np.kron(k_grid[:,ixt+1], np.repeat(1,len(a_grid[:,ixt+1])))\n", + "# here we are treating realized HC in T as the deterministic HC at the end of T-1\n", + "\n", + "# now melt the expected values we calculated\n", + "df = pd.DataFrame(E_uPrime[:, :, ixt+1])\n", + "policy_df[:,2, ixt] = pd.melt(df)[\"value\"]\n", + "df = pd.DataFrame(E_dVdK[:, :, ixt+1])\n", + "policy_df[:,3, ixt] = pd.melt(df)[\"value\"]\n", + "df = pd.DataFrame(E_eps_dVdK[:, :, ixt+1])\n", + "policy_df[:,4, ixt] = pd.melt(df)[\"value\"]\n", + "\n", + "# a few other values\n", + "df = pd.DataFrame(extrap[:, :, ixt+1])\n", + "policy_df[:,5, ixt] = pd.melt(df)[\"value\"]\n", + "df = pd.DataFrame(lc[:, :, ixt+1])\n", + "policy_df[:,6, ixt] = pd.melt(df)[\"value\"]" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now have a series of points $\\left\\{A_{\\overline{T-1}},K_{\\overline{T-1}} \\right\\}$, with the associated values of $\\mathbb{E} \\left( u'(c_T) \\right)$, $\\mathbb{E} \\left( \\frac{\\partial V}{\\partial K} \\right)$, and $\\mathbb{E} \\left( \\epsilon \\frac{\\partial V}{\\partial K} \\right)$. Essentially we will now work backwards to determine what the optimal controls $c_{T-1}^*$ and $h_{T-1}^*$, as well as the beginning of period $T-1$ state variables $\\left\\{A_{\\underline{T-1}},K_{\\underline{T-1}} \\right\\}$ *must have been* in order to arrive at the exogenous grid points we set for period $T$.\n", + "\n", + "Note that we can solve eqn. (9) for $c_t$:\n", + "\n", + "\\begin{equation}\n", + "\\tag{15}\n", + "c_t = (\\beta R \\mathbb{E} u'(c_{t+1}))^{-1/\\gamma}\n", + "\\end{equation}\n", + "\n", + "Since we have $ \\mathbb{E} u(c_{t+1})$ for each endogenous point, we can calculate $c_{T-1}^*$ very quickly for each of our points, and then fill in the marginal utility from that consumption (which we will use when we solve period $T-2$)." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Now get that period's optimal consumption using the inverted FOC\n", + "policy_df[:, 7, ixt] = (beta*R*policy_df[:,2,ixt])**(1/(-gamma));\n", + "# policy_df[:, 7, lifespan-1-1]\n", + "# and get marginal utility:\n", + "for i in range(policy_df.shape[0]):\n", + " # print(i)\n", + " policy_df[i,8, ixt] = getMargUtilityC(policy_df[i, 7, ixt], gamma)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have the optimal level of consumption, we can use our first order condition to solve for labor supply. Recall eqn. (13):\n", + "\n", + "\\begin{equation}\n", + "\\tag{13}\n", + "v'(h_t) = \\Psi(K_{\\bar{t}}, h_t,t) u'(c_t) + \\beta \\kappa(h_t, \\Psi(K_{\\bar{t}}, h_t,t), t) \\mathbb{E} \\left[ \\epsilon_t \\dfrac{\\partial V_{t+1}}{\\partial K_{t+1}} \\right]\n", + "\\end{equation}\n", + "\n", + "For each of our endogenous points, we have all of the elements except for $h$. Thus we can take the difference and set up a function we can minimize to solve for optimal labor supply." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Define a function from eqn. 13 to use as the root finder.\n", + "def hFunc(h, Ktilde, uprime, age, E_eps_dvdk, b, eta, beta, d, alpha, omega2, omega0, omega1, delta):\n", + " # print(h)\n", + " tau = taufunc(age, omega0, omega1)\n", + " kstart = psifunc(Ktilde, h, tau, d, alpha, omega2, delta)\n", + " kappa = kappafunc(kstart, h, tau, d, omega2, alpha)\n", + " dis = getMargDisutilityH(h, b, eta)\n", + " ret = dis - uprime * kstart - beta * kappa * E_eps_dvdk\n", + " # ret should be zero; we want to minimize, so take absolute value, and then multiply\n", + " # since this number may be quite small\n", + " ret = abs(ret) * 1e8\n", + " return ret\n", + "\n", + "# Now solve for optimal labor supply\n", + "for ixn in range(policy_df[:, :, ixt].shape[0]):\n", + " if policy_df[ixn, 6, ixt] == 1:\n", + " policy_df[ixn, 7, ixt] = np.nan\n", + " policy_df[ixn, 9, ixt] = np.nan \n", + " else:\n", + " #print(ixn)\n", + " Ktilde = policy_df[ixn, 1, ixt]\n", + " this_E_eps_dVdK = policy_df[ixn, 4, ixt]\n", + " thisuprime = policy_df[ixn, 8, ixt]\n", + " out = scipy.optimize.minimize_scalar(hFunc, bounds=(1, hmax), method=\"bounded\", args=(Ktilde, thisuprime, ixt+1+20,\n", + " this_E_eps_dVdK, b, eta, beta,\n", + " d, alpha, omega2, omega0, omega1, delta))\n", + " policy_df[ixn, 9, ixt] = out.x" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "# Now get the subset that is not liquidity constrained as a data frame\n", + "col_names = [\"A_end\", \"K_end\",\"E_uprime\",\"E_dVdK\",\"E_eps_dVdK\",\"extrap\",\"lc\",\"C\",\"uprime\",\"H\",\"K_start\",\"A_start\",\"dVdK\"]\n", + "sub_policy_df = pd.DataFrame(policy_df[:, :, ixt], columns=col_names)\n", + "sub_policy_df = sub_policy_df[sub_policy_df['lc'] != 1]" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the last step, I drop all of the rows where the agent was liquidity constrained in period $T$. If this was the case, this means they were not following their FOCs, and thus the endogenous step we take is not valid. For purposes of filling in our exogenous grid, then, these observations are not helpful.\n", + "\n", + "Here I need to comment on another shortcoming of my implementation of the EGM. If there are too few observations that are not liquidity constrained, either due to a poorly-set exogenous grid or to a strange parameterization of the problem, we cannot use the polynomial approximation to help fill in our exogenous grid. In my full implementation in R, all of the code I am presenting here is inside of a function called \"SolveProblem()\". If it is the case that during one of the endogenous steps there are fewer than 1/10 of the grid points that are done by interpolation, I stop the process and return a large number for the likelihood function. This is problematic inasmuch as it renders the likelihood function discontinuous, preventing me from using gradient-based optimization methods. However, as I currently implement the EGM, there may be some parameterizations where it is not possible to solve the agent's problem. My argument is that because I am setting a \"reasonable\" exogenous grid (matching what we find in the data), any parameterizaiton at which no solution can be reached cannot be a good one, and so my short-circuiting of the likelihood function only causes a problem for optimization, but will not prevent me from finding the true parameters.\n", + "\n", + "With that out of the way, let's fill in the rest of the necessary information for our endogenous grid step. First, we can use our inverse laws of motion to get $\\left\\{A_{\\underline{T-1}},K_{\\underline{T-1}} \\right\\}$, and then will calculate $\\mathbb{E} \\epsilon \\frac{\\partial V}{\\partial K}$." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# Use the inverted law of motion for K to get start-of-period HC\n", + "thistau = taufunc(ixt+1+20, omega0, omega1)\n", + "for i in range(sub_policy_df.shape[0]):\n", + " # print(i)\n", + " sub_policy_df.iloc[i, 10] = psifunc(sub_policy_df.iloc[i, 1], sub_policy_df.iloc[i, 9],thistau,d,alpha,omega2,delta)\n", + "\n", + "# Now use the inverted IBC to get start-of-period A\n", + "for i in range(sub_policy_df.shape[0]):\n", + " # print(i)\n", + " sub_policy_df.iloc[i, 11] = sub_policy_df.iloc[i, 0] - sub_policy_df.iloc[i, 9]*sub_policy_df.iloc[i, 10] + sub_policy_df.iloc[i, 7]\n", + "\n", + "# Now get dVdK\n", + "for i in range(sub_policy_df.shape[0]):\n", + " # print(i)\n", + " thislambda = lambdafunc(sub_policy_df.iloc[i, 9], thistau, omega2, d, alpha, delta)\n", + " sub_policy_df.iloc[i, 12] = sub_policy_df.iloc[i, 9]*sub_policy_df.iloc[i, 8] + beta*thislambda*sub_policy_df.iloc[i, 4]\n", + "\n", + "# policy_df index:\n", + "# 0 = end-of-period assets\n", + "# 1 = end-of-period capital\n", + "# 2 = E_uPrime\n", + "# 3 = E_dVdK\n", + "# 4 = E_eps_dVdK\n", + "# 5 = extrap\n", + "# 6 = lc\n", + "# 7 = C\n", + "# 8 = u'(C)\n", + "# 9 = H\n", + "# 10 = start-of-period capital\n", + "# 11 = start-of-period assets\n", + "# 12 = E_eps_dVdK" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our endogenous grid step has left us with a cross-section of data for the policy function. We now want to construct an exogenous, rectangular consumption grid from this data. The method currently used is essentially a Delauney triangulation." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "# Now create a grid from our endogenous data points\n", + "y = a_grid[:,ixt]\n", + "x = k_grid[:,ixt]\n", + "xx, yy = np.meshgrid(x, y)\n", + "cgrid = griddata((sub_policy_df['K_start'], sub_policy_df['A_start']), sub_policy_df['C'], (xx, yy), method='linear')\n", + "policyC[:,:,ixt] = cgrid" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's create the polynomial approximation we'll use to fill in the values of our grid that are not covered by the endogenous grid points and are not liquidity constrained." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# create a data frame xx\n", + "xx = pd.DataFrame({'y1': sub_policy_df[\"C\"],\n", + " 'x1': sub_policy_df[\"A_start\"],\n", + " 'x2': sub_policy_df[\"K_start\"]})\n", + "# Set up the polynomials of A and K\n", + "xx['x3'] = xx['x1'] * xx['x2']\n", + "xx['x4'] = (xx['x1'] ** 2) / 1e5 # assets\n", + "xx['x5'] = xx['x2'] ** 2\n", + "\n", + "# get the model\n", + "out_C2 = sm.OLS(xx['y1'], xx[['x1', 'x2', 'x3', 'x4', 'x5']]).fit()\n", + "\n", + "# get the relevant portion of the polynomial df\n", + "# set up the data frame with all the polynomials ahead of time\n", + "col_names = [\"x1\",\"x2\",\"x3\",\"x4\",\"x5\"]\n", + "polynomialdf = pd.DataFrame(polynomial_df[:,:,ixt], columns=col_names)\n", + "# polynomialdf" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we will fill in the rest of our exogenous grid: those points where the agent is liquidity constrained, and those points which the endogenous grid did not cover. In order to get optimal H from optimal C, we will define a function using the FOC for labor supply which we can plug into a root finder." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "# First, cast our slices for interpolation as contiguous arrays\n", + "K_slice = np.ascontiguousarray(k_grid[:,ixt+1])\n", + "A_slice = np.ascontiguousarray(a_grid[:,ixt+1])\n", + "E_eps_dVdK_slice = np.ascontiguousarray(E_eps_dVdK[:,:,ixt+1])\n", + "\n", + "# and define our interpolator function\n", + "interp_func_eps_dVdK = RegularGridInterpolator((K_slice, A_slice), E_eps_dVdK_slice, method='linear')\n", + "\n", + "def hFuncZero(K,c,h,age,E_eps_dvdk,b,eta,gamma,beta,d, omega2,alpha,omega0,omega1):\n", + " thistau = taufunc(age, omega0, omega1)\n", + " kappa = kappafunc(K, h, thistau, d, omega2, alpha)\n", + " dis = getMargDisutilityH(h, b, eta)\n", + " util = getMargUtilityC(c,gamma)\n", + " ret = dis - util * K - beta * kappa * E_eps_dvdk\n", + " return ret\n", + "\n", + "def FOCzero(h,c,age,A,K,gamma,b,eta,omega0,omega1,d,alpha,omega2,delta,beta,E_eps_dVdK,sub_agrid,sub_kgrid):\n", + " hclb =np.min(sub_kgrid)\n", + " hcub =np.max(sub_kgrid)\n", + " alb =np.min(sub_agrid)\n", + " aub =np.max(sub_agrid)\n", + "\n", + " # next period assets given this choice\n", + " Aend = R*(A+h*K - c)\n", + " # Check it doesn't violate our bounds so that we can interpolate\n", + " # Aend = np.where(Aend > alb, Aend, alb+0.1)\n", + " # Aend = np.where(Aend < aub, Aend, aub-0.1)\n", + "\n", + " # deterministic part of human capital:\n", + " Kend = gFunc(h,K,age,delta,omega0,omega1,alpha,omega2,d).item()\n", + " # Check it doesn't violate our bounds so that we can interpolate\n", + " Kend = np.where(Kend > hclb, Kend, hclb+0.1).item()\n", + " Kend = np.where(Kend < hcub, Kend, hcub-0.1).item()\n", + "\n", + " if (Aend < aub) & (Aend > alb):\n", + " # get the value of E_eps_dVdK for next period\n", + " # this_epsdVdKnext = interp2d(sub_kgrid, sub_agrid, E_eps_dVdK, kind='linear')(Kend, Aend)\n", + " thispoint = (Kend, Aend)\n", + " this_epsdVdKnext = interp_func_eps_dVdK(thispoint)\n", + "\n", + " # solve problem\n", + " val = hFuncZero(K, c, h, age, this_epsdVdKnext, b, eta, gamma, beta, d, omega2, alpha, omega0, omega1)\n", + " else:\n", + " val = np.where(Aend <= alb, -(Aend*1.2)*1e2, (Aend*1.2)*1e2)\n", + " \n", + " return val*10\n", + "\n", + "# Now let's go through our grid points\n", + "for ixa in range(nA): # looping through asset points\n", + " for ixk in range(nK): # looping through human capital points\n", + " if ixt == 43:\n", + " cub = hmax*k_grid[ixk,ixt] + a_grid[ixa,ixt] + (2040*k_grid[ixk,ixt])*(1/R)\n", + " elif ixt < 43:\n", + " cub = hmax*k_grid[ixk,ixt] + a_grid[ixa,ixt] + (2*2040*k_grid[ixk,ixt])*(1/R)\n", + "\n", + " if cub < cmin:\n", + " # liquidity constrained\n", + " policyC[ixa,ixk,ixt] = cmin\n", + " policyH[ixa,ixk,ixt] = hmax\n", + " lc[ixa,ixk,ixt] = 1\n", + " \n", + " # fill in uPrime and dVdK\n", + " E_eps_dVdK_tp1 = np.max(E_eps_dVdK[:,:,ixt+1])\n", + " uPrime[ixa,ixk,ixt] = uprimemin\n", + " dVdK[ixa,ixk,ixt] = hmax*uprimemin + beta*lambdafunc(hmax,taufunc(ixt+1+20,omega0,omega1),omega2,d,alpha,delta)*E_eps_dVdK_tp1\n", + " else:\n", + " if np.isnan(policyC[ixa,ixk,ixt]) or (ixa > 1 and (not np.isnan(policyC[ixa,ixk,ixt]) and policyC[ixa-1,ixk,ixt] > policyC[ixa,ixk,ixt])):\n", + " extrap[ixa,ixk,ixt] = 1\n", + " nxx = polynomialdf[(polynomialdf['x1'] == a_grid[ixa,ixt]) & (polynomialdf['x2'] == k_grid[ixk,ixt])]\n", + " policyC[ixa,ixk,ixt] = out_C2.predict(nxx) \n", + "\n", + " # make sure c is not over our upper bound or below our lower bound\n", + " c = np.where(policyC[ixa,ixk,ixt] > cmin, policyC[ixa,ixk,ixt], cmin+.1).item()\n", + " c = np.where(c < cub, c, cub-.1).item()\n", + " policyC[ixa,ixk,ixt] = c\n", + " \n", + " # catch occasions where the polynomial predicts lower consumption with assets\n", + " if ixa > 0:\n", + " if policyC[ixa-1,ixk,ixt] > policyC[ixa,ixk,ixt]:\n", + " c = policyC[ixa-1,ixk,ixt]*1.01\n", + " policyC[ixa,ixk,ixt] = c\n", + " replace[ixa,ixk,ixt] = 1\n", + "\n", + " # Now we have optimal C and move on to optimal H\n", + " # get minimum hours\n", + " # now use our consumption choice to get minimum hours\n", + " if ixt == 43:\n", + " hmin = (-(2040 * k_grid[ixk, ixt]) * (1 / R) - a_grid[ixa, ixt] + c) / k_grid[ixk, ixt] + 0.1\n", + " elif ixt < 43:\n", + " hmin = (-(2 * 2040 * k_grid[ixk, ixt]) * (1 / R) - a_grid[ixa, ixt] + c) / k_grid[ixk, ixt] + 0.1\n", + "\n", + " # Make sure hmin is positive but less than hmax\n", + " hmin = np.where(hmin > 0, hmin, 0.1).item()\n", + " hmin = np.where(hmin >= hmax, hmax-.1, hmin).item()\n", + " \n", + " p1 = FOCzero(hmin,c,ixt+1+20,a_grid[ixa,ixt],k_grid[ixk,ixt],gamma,b,eta,omega0,omega1,d,alpha,omega2,delta,beta,\n", + " E_eps_dVdK_slice, A_slice, K_slice).item()\n", + " p2 = FOCzero(hmax,c,ixt+1+20,a_grid[ixa,ixt],k_grid[ixk,ixt],gamma,b,eta,omega0,omega1,d,alpha,omega2,delta,beta,\n", + " E_eps_dVdK_slice, A_slice, K_slice).item()\n", + " if p1*p2>0: # either both positive or both negative\n", + " h = hmin if abs(p1) np.min(k_grid[:,ixt+1]), Knext, np.min(k_grid[:,ixt+1])+0.1).item()\n", + " Knext = np.where(Knext < np.max(k_grid[:,ixt+1]), Knext, np.max(k_grid[:,ixt+1])-0.1).item()\n", + "\n", + " # next period A\n", + " Anext = a_grid[ixa,ixt] + h*thisK - c\n", + " Anext = np.where(Anext > np.min(a_grid[:,ixt+1]), Anext, np.min(a_grid[:,ixt+1])+0.1).item()\n", + " Anext = np.where(Anext < np.max(a_grid[:,ixt+1]), Anext, np.max(a_grid[:,ixt+1])-0.1).item()\n", + " \n", + " # now interpolate E[eps_dVdK]_{t+1}\n", + " point = ( Knext, Anext)\n", + " E_eps_dVdK_tp1 = interp_func_eps_dVdK(point).item()\n", + " uPrime[ixa,ixk,ixt] = getMargUtilityC(c,gamma)\n", + " dVdK[ixa,ixk,ixt] = h*uPrime[ixa,ixk,ixt] + beta*lambdafunc(h,taufunc(ixt+1+20,omega0,omega1), omega2, d, alpha,delta)*E_eps_dVdK_tp1\n", + " # Once we have looped through human capital once and filled in the grid for a given level of assets, we can loop through again and fill in our expectations grid.\n", + " \n", + "# set up our interpolator functions for the expectations grid\n", + " interp_func_uprime = interpolate.interp1d(k_grid[:,ixt], uPrime[ixa,:,ixt])\n", + " interp_func_dVdK = interpolate.interp1d(k_grid[:,ixt], dVdK[ixa,:,ixt])\n", + "\n", + " for ixk in range(nK):\n", + " nextK = quad*k_grid[ixk,ixt]\n", + " # If HC is off the grid, force it to stay on.\n", + " nextK = np.where(nextK <= k_grid[0,ixt], k_grid[0,ixt]+0.01, nextK)\n", + " nextK = np.where(nextK >= k_grid[nK-1,ixt], k_grid[nK-1,ixt]-0.01, nextK)\n", + "\n", + " Uprime_nextvec = interp_func_uprime(nextK)\n", + " dVdK_nextvec = interp_func_dVdK(nextK)\n", + " for i in range(len(dVdK_nextvec)):\n", + " eps_dVdK_nextvec[i] = quad[i] * dVdK_nextvec[i]\n", + " \n", + " E_uPrime[ixa,ixk,ixt] = np.mean(Uprime_nextvec)\n", + " E_dVdK[ixa,ixk,ixt] = np.mean(dVdK_nextvec)\n", + " E_eps_dVdK[ixa,ixk,ixt] = np.mean(eps_dVdK_nextvec)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are now done solving the problem in period $T-1$. We have an approximation of our policy function on an exogenous grid over which we can interpolate when we want to simulate our data, and we have all of the other objects we will need to perform the endogenous grid step in period $T-2$, and iterate this process backwards to the initial period. In the next cell, I will loop over all the previous time periods and do just this." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "42\n", + "41\n", + "40\n", + "39\n", + "38\n", + "37\n", + "36\n", + "35\n", + "34\n", + "33\n", + "32\n", + "31\n", + "30\n", + "29\n", + "28\n", + "27\n", + "26\n", + "25\n", + "24\n", + "23\n", + "22\n", + "21\n", + "20\n", + "19\n", + "18\n", + "17\n", + "16\n", + "15\n", + "14\n", + "13\n", + "12\n", + "11\n", + "10\n", + "9\n", + "8\n", + "7\n", + "6\n", + "5\n", + "4\n", + "3\n", + "2\n", + "1\n", + "0\n" + ] + } + ], + "source": [ + "for ixt in range(lifespan-3, -1, -1):\n", + " # Melt our previous values from nexy period to use in this one\n", + " # end of period assets\n", + " policy_df[:,0,ixt] = np.tile(a_grid[:,ixt+1]*(R**(-1)),len(k_grid[:,ixt+1])) # unwind the interest payment to get assets at end of T-1\n", + "\n", + " # end of period human capital\n", + " policy_df[:,1,ixt] = np.kron(k_grid[:,ixt+1], np.repeat(1,len(a_grid[:,ixt+1])))\n", + " # here we are treating realized HC in T as the deterministic HC at the end of T-1\n", + "\n", + " # now melt the expected values we calculated\n", + " df = pd.DataFrame(E_uPrime[:, :, ixt+1])\n", + " policy_df[:,2, ixt] = pd.melt(df)[\"value\"]\n", + " df = pd.DataFrame(E_dVdK[:, :, ixt+1])\n", + " policy_df[:,3, ixt] = pd.melt(df)[\"value\"]\n", + " df = pd.DataFrame(E_eps_dVdK[:, :, ixt+1])\n", + " policy_df[:,4, ixt] = pd.melt(df)[\"value\"]\n", + "\n", + " # a few other values\n", + " df = pd.DataFrame(extrap[:, :, ixt+1])\n", + " policy_df[:,5, ixt] = pd.melt(df)[\"value\"]\n", + " df = pd.DataFrame(lc[:, :, ixt+1])\n", + " policy_df[:,6, ixt] = pd.melt(df)[\"value\"]\n", + "\n", + " # Now get that period's optimal consumption using the inverted FOC\n", + " policy_df[:, 7, ixt] = (beta*R*policy_df[:,2,ixt])**(1/(-gamma));\n", + " # and get marginal utility:\n", + " for i in range(policy_df.shape[0]):\n", + " policy_df[i,8, ixt] = getMargUtilityC(policy_df[i, 7, ixt], gamma)\n", + "\n", + " # Now solve for optimal labor supply\n", + " for ixn in range(policy_df[:, :, ixt].shape[0]):\n", + " if policy_df[ixn, 6, ixt] == 1:\n", + " policy_df[ixn, 7, ixt] = np.nan\n", + " policy_df[ixn, 9, ixt] = np.nan \n", + " else:\n", + " Ktilde = policy_df[ixn, 1, ixt]\n", + " this_E_eps_dVdK = policy_df[ixn, 4, ixt]\n", + " thisuprime = policy_df[ixn, 8, ixt]\n", + " out = scipy.optimize.minimize_scalar(hFunc, bounds=(1, hmax), method=\"bounded\", args=(Ktilde, thisuprime, ixt+1+20,\n", + " this_E_eps_dVdK, b, eta, beta,\n", + " d, alpha, omega2, omega0, omega1, delta))\n", + " policy_df[ixn, 9, ixt] = out.x\n", + "\n", + " # Turn into a data frame\n", + " col_names = [\"A_end\", \"K_end\",\"E_uprime\",\"E_dVdK\",\"E_eps_dVdK\",\"extrap\",\"lc\",\"C\",\"uprime\",\"H\",\"K_start\",\"A_start\",\"dVdK\"]\n", + " # sub_policy_df = pd.DataFrame(policy_df[policy_df[:, 6, lifespan-1-1] != 1, :, lifespan-1-1], columns=col_names)\n", + " sub_policy_df = pd.DataFrame(policy_df[:, :, ixt], columns=col_names)\n", + " sub_policy_df = sub_policy_df[sub_policy_df['lc'] != 1]\n", + "\n", + " # Work backwards to get start-of-period K\n", + " thistau = taufunc(ixt+1+20, omega0, omega1)\n", + " for i in range(sub_policy_df.shape[0]):\n", + " # print(i)\n", + " sub_policy_df.iloc[i, 10] = psifunc(sub_policy_df.iloc[i, 1], sub_policy_df.iloc[i, 9],thistau,d,alpha,omega2,delta)\n", + "\n", + " # Now use the inverted IBC to get start-of-period A\n", + " for i in range(sub_policy_df.shape[0]):\n", + " sub_policy_df.iloc[i, 11] = sub_policy_df.iloc[i, 0] - sub_policy_df.iloc[i, 9]*sub_policy_df.iloc[i, 10] + sub_policy_df.iloc[i, 7]\n", + "\n", + " # Now get dVdK\n", + " for i in range(sub_policy_df.shape[0]):\n", + " thislambda = lambdafunc(sub_policy_df.iloc[i, 9], thistau, omega2, d, alpha, delta)\n", + " sub_policy_df.iloc[i, 12] = sub_policy_df.iloc[i, 9]*sub_policy_df.iloc[i, 8] + beta*thislambda*sub_policy_df.iloc[i, 4]\n", + "\n", + " # Now create a grid from our endogenous data points\n", + " y = a_grid[:,ixt]\n", + " x = k_grid[:,ixt]\n", + " xx, yy = np.meshgrid(x, y)\n", + "\n", + " cgrid = griddata((sub_policy_df['K_start'], sub_policy_df['A_start']), sub_policy_df['C'], (xx, yy), method='linear')\n", + " policyC[:,:,ixt] = cgrid\n", + "\n", + " # Now setup the polynomial approximation for the policy function\n", + " xx = pd.DataFrame({'y1': sub_policy_df[\"C\"],\n", + " 'x1': sub_policy_df[\"A_start\"],\n", + " 'x2': sub_policy_df[\"K_start\"]})\n", + " # Set up the polynomials of A and K\n", + " xx['x3'] = xx['x1'] * xx['x2']\n", + " xx['x4'] = (xx['x1'] ** 2) / 1e5 # assets\n", + " xx['x5'] = xx['x2'] ** 2\n", + "\n", + " # get the model\n", + " out_C2 = sm.OLS(xx['y1'], xx[['x1', 'x2', 'x3', 'x4', 'x5']]).fit()\n", + "\n", + " # get the relevant portion of the polynomial df\n", + " # set up the data frame with all the polynomials ahead of time\n", + " col_names = [\"x1\",\"x2\",\"x3\",\"x4\",\"x5\"]\n", + " polynomialdf = pd.DataFrame(polynomial_df[:,:,ixt], columns=col_names)\n", + "\n", + " # Now go through all of our grid points\n", + " # First, cast our slices for interpolation as contiguous arrays\n", + " K_slice = np.ascontiguousarray(k_grid[:,ixt+1])\n", + " A_slice = np.ascontiguousarray(a_grid[:,ixt+1])\n", + " E_eps_dVdK_slice = np.ascontiguousarray(E_eps_dVdK[:,:,ixt+1])\n", + " # and define our interpolator function\n", + " interp_func_eps_dVdK = RegularGridInterpolator((K_slice, A_slice), E_eps_dVdK_slice, method='linear')\n", + "\n", + " for ixa in range(nA): # looping through asset points\n", + " for ixk in range(nK): # looping through human capital points\n", + " if ixt == 43:\n", + " cub = hmax*k_grid[ixk,ixt] + a_grid[ixa,ixt] + (2040*k_grid[ixk,ixt])*(1/R)\n", + " elif ixt < 43:\n", + " cub = hmax*k_grid[ixk,ixt] + a_grid[ixa,ixt] + (2*2040*k_grid[ixk,ixt])*(1/R)\n", + "\n", + " if cub < cmin:\n", + " # liquidity constrained\n", + " policyC[ixa,ixk,ixt] = cmin\n", + " policyH[ixa,ixk,ixt] = hmax\n", + " lc[ixa,ixk,ixt] = 1\n", + " \n", + " # fill in uPrime and dVdK\n", + " E_eps_dVdK_tp1 = np.max(E_eps_dVdK[:,:,ixt+1])\n", + " uPrime[ixa,ixk,ixt] = uprimemin\n", + " dVdK[ixa,ixk,ixt] = hmax*uprimemin + beta*lambdafunc(hmax,taufunc(ixt+1+20,omega0,omega1),omega2,d,alpha,delta)*E_eps_dVdK_tp1\n", + " else:\n", + " if np.isnan(policyC[ixa,ixk,ixt]) or (ixa > 1 and (not np.isnan(policyC[ixa,ixk,ixt]) and policyC[ixa-1,ixk,ixt] > policyC[ixa,ixk,ixt])):\n", + " extrap[ixa,ixk,ixt] = 1\n", + " nxx = polynomialdf[(polynomialdf['x1'] == a_grid[ixa,ixt]) & (polynomialdf['x2'] == k_grid[ixk,ixt])]\n", + " policyC[ixa,ixk,ixt] = out_C2.predict(nxx) \n", + "\n", + " # make sure c is not over our upper bound or below our lower bound\n", + " c = np.where(policyC[ixa,ixk,ixt] > cmin, policyC[ixa,ixk,ixt], cmin+.1).item()\n", + " c = np.where(c < cub, c, cub-.1).item()\n", + " policyC[ixa,ixk,ixt] = c\n", + " \n", + " # catch occasions where the polynomial predicts lower consumption with assets\n", + " if ixa > 0:\n", + " if policyC[ixa-1,ixk,ixt] > policyC[ixa,ixk,ixt]:\n", + " c = policyC[ixa-1,ixk,ixt]*1.01\n", + " policyC[ixa,ixk,ixt] = c\n", + " replace[ixa,ixk,ixt] = 1\n", + "\n", + " # Now we have optimal C and move on to optimal H\n", + " # get minimum hours\n", + " # now use our consumption choice to get minimum hours\n", + " if ixt == 43:\n", + " hmin = (-(2040 * k_grid[ixk, ixt]) * (1 / R) - a_grid[ixa, ixt] + c) / k_grid[ixk, ixt] + 0.1\n", + " elif ixt < 43:\n", + " hmin = (-(2 * 2040 * k_grid[ixk, ixt]) * (1 / R) - a_grid[ixa, ixt] + c) / k_grid[ixk, ixt] + 0.1\n", + "\n", + " # Make sure hmin is positive but less than hmax\n", + " hmin = np.where(hmin > 0, hmin, 0.1).item()\n", + " hmin = np.where(hmin >= hmax, hmax-.1, hmin).item()\n", + " \n", + " p1 = FOCzero(hmin,c,ixt+1+20,a_grid[ixa,ixt],k_grid[ixk,ixt],gamma,b,eta,omega0,omega1,d,alpha,omega2,delta,beta,\n", + " E_eps_dVdK_slice, A_slice, K_slice).item()\n", + " p2 = FOCzero(hmax,c,ixt+1+20,a_grid[ixa,ixt],k_grid[ixk,ixt],gamma,b,eta,omega0,omega1,d,alpha,omega2,delta,beta,\n", + " E_eps_dVdK_slice, A_slice, K_slice).item()\n", + " if p1*p2>0: # either both positive or both negative\n", + " h = hmin if abs(p1) np.min(k_grid[:,ixt+1]), Knext, np.min(k_grid[:,ixt+1])+0.1).item()\n", + " Knext = np.where(Knext < np.max(k_grid[:,ixt+1]), Knext, np.max(k_grid[:,ixt+1])-0.1).item()\n", + "\n", + " # next period A\n", + " Anext = a_grid[ixa,ixt] + h*thisK - c\n", + " Anext = np.where(Anext > np.min(a_grid[:,ixt+1]), Anext, np.min(a_grid[:,ixt+1])+0.1).item()\n", + " Anext = np.where(Anext < np.max(a_grid[:,ixt+1]), Anext, np.max(a_grid[:,ixt+1])-0.1).item()\n", + " \n", + " # now interpolate E[eps_dVdK]_{t+1}\n", + " point = ( Knext, Anext) # (Knext,Anext)\n", + " E_eps_dVdK_tp1 = interp_func_eps_dVdK(point).item()\n", + " uPrime[ixa,ixk,ixt] = getMargUtilityC(c,gamma)\n", + " dVdK[ixa,ixk,ixt] = h*uPrime[ixa,ixk,ixt] + beta*lambdafunc(h,taufunc(ixt+1+20,omega0,omega1), omega2, d, alpha,delta)*E_eps_dVdK_tp1\n", + " # Once we have looped through human capital once and filled in the grid for a given level of assets, we can loop through again and fill in our expectations grid.\n", + " \n", + " # set up our interpolator functions for the expectations grid\n", + " interp_func_uprime = interpolate.interp1d(k_grid[:,ixt], uPrime[ixa,:,ixt])\n", + " interp_func_dVdK = interpolate.interp1d(k_grid[:,ixt], dVdK[ixa,:,ixt])\n", + "\n", + " for ixk in range(nK):\n", + " nextK = quad*k_grid[ixk,ixt]\n", + " # If HC is off the grid, force it to stay on.\n", + " nextK = np.where(nextK <= k_grid[0,ixt], k_grid[0,ixt]+0.01, nextK)\n", + " nextK = np.where(nextK >= k_grid[nK-1,ixt], k_grid[nK-1,ixt]-0.01, nextK)\n", + "\n", + " Uprime_nextvec = interp_func_uprime(nextK)\n", + " dVdK_nextvec = interp_func_dVdK(nextK)\n", + " for i in range(len(dVdK_nextvec)):\n", + " eps_dVdK_nextvec[i] = quad[i] * dVdK_nextvec[i]\n", + " \n", + " E_uPrime[ixa,ixk,ixt] = np.mean(Uprime_nextvec)\n", + " E_dVdK[ixa,ixk,ixt] = np.mean(dVdK_nextvec)\n", + " E_eps_dVdK[ixa,ixk,ixt] = np.mean(eps_dVdK_nextvec)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "econ-ark", + "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.9" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "d7363176c81bf58bfee85a7a82bd18d1cb91d56304d8eb5248e7f7cc27dde8b4" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}