1st.ass.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIGMA = 1.5 # utility parameter\n",
    "BETA = 0.95 # \n",
    "DELTA = 0.1 # \n",
    "ALPHA = 0.3 # \n",
    "DEV = 0.9 # maximal deviation from steady state\n",
    "EPSILON = 1e-6 # convergence parameter\n",
    "N = 1000 # number of data points in the grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "analytical soulition of stationary kapital and consumtion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.6257456456982027"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# steady state value of k\n",
    "kss = ((DELTA-1+(1/BETA))/ALPHA)**(1/(ALPHA-1))\n",
    "kss\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.073331114820493"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# steady consumption\n",
    "css = kss*(kss**(ALPHA-1)-DELTA)\n",
    "css"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "compute the grid search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2625745645698202, 4.988916726826585)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kmin = (1 - DEV)*kss # grid's lower bound\n",
    "kmax = (1 + DEV)*kss # grid's upper bound\n",
    "kmin, kmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dk = (kmax - kmin)/(N-1) # compute distance between two consecutive data points on the grid (step size)\n",
    "kgrid = np.arange(kmin, kmax+dk, dk) # build the grid as a COLUMN vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "vectors for value function, decision rule and Bellman operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Define initial condition for value function. As we know from Contraction Mapping Theorem, this could be anything.\n",
    "To start with, you could leave the vector of zeros as initial condition and later on think about something better\n",
    "'''\n",
    "v = np.zeros(N) # value function\n",
    "dr = np.zeros(N) # decision rule\n",
    "Bv = np.zeros(N) # Bellman operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "functions for utility, production and distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u(c) : return (c**(1-SIGMA))/(1-SIGMA) # \n",
    "def F(k) : return k**ALPHA + (1-DELTA)*k #\n",
    "def dist(a,b) : return max(abs(a-b)) #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "crit = 1 # could be anything bigger than stopping criterion - the algorithm will keep iterating until crit < epsi\n",
    "iter = 0 # iteration counter (not really needed, just FYI)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "non-optimized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\martin.ranieri\\AppData\\Local\\Temp\\ipykernel_6828\\3902570932.py:1: RuntimeWarning: invalid value encountered in power\n",
      "  def u(c) : return (c**(1-SIGMA))/(1-SIGMA) #\n"
     ]
    }
   ],
   "source": [
    "while crit > EPSILON:\n",
    "    v = Bv.copy() # update the value function\n",
    "\n",
    "    #tmp = [u(F(k) - kgrid) + v*BETA for k in kgrid]\n",
    "    tmp = np.zeros(N) \n",
    "    for i, k in enumerate(kgrid):\n",
    "        c = F(k) - kgrid\n",
    "        tmp[i] = u(c) + v*BETA\n",
    "\n",
    "    rd = np.argmax(tmp, axis=-1)\n",
    "    Bv = np.max(tmp, axis=-1)\n",
    "    # Bv = np.take_along_axis(tmp,np.expand_dims(rd, axis=-1),-1).squeeze(axis=-1)\n",
    "    \n",
    "    crit = dist(Bv,v) # compute convergence criterion\n",
    "\n",
    "    print(crit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "c =  [ F(kgrid[i]) - kgrid[rd[i]] for i in range(N)]\n",
    "kk = [ F(kgrid[i]) - c[i]  for i in range(N)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\martin.ranieri\\OneDrive - Università degli Studi di Milano\\Desktop\\MASSARO CM\\venv\\lib\\site-packages\\seaborn\\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n",
      "c:\\Users\\martin.ranieri\\OneDrive - Università degli Studi di Milano\\Desktop\\MASSARO CM\\venv\\lib\\site-packages\\seaborn\\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjMUlEQVR4nO3de3CkV5nf8e/pi9S6jzXSaK6aC3OzMdgY4QvjxcTGZgzGyx9bgYQslQ2pYb1LCipLtpKqJNRu5Q9S2Wwtf+yGOIZkHchCChZCNkCgiMEGfJsZMPgCxnPTaEa3GY3U3ZJafXlP/jj9St0tzUhjv939dvfvUzWl7ul295EsPTrzvM/zHGOtRUREwitS7wWIiMi1KVCLiIScArWISMgpUIuIhJwCtYhIyMWq8aIDAwN2z5491XhpEZGmdOLEiUvW2sG1HqtKoN6zZw/Hjx+vxkuLiDQlY8y5qz22odSHMWaTMeZrxphfGWNeMcbcFdzyRETkWja6o/4c8F1r7e8YY9qAziquSURESqwbqI0xfcC7gH8MYK3NAtnqLktERHwbSX3sBaaB/2qM+Zkx5jFjTFflk4wxx4wxx40xx6enpwNfqIhIq9pIoI4BtwH/yVr7NmAe+JeVT7LWPmqtHbHWjgwOrnnhUkREXoeNBOoxYMxa+2zx/tdwgVtERIB83uOF81f47ovjvHB+lnzeC/T1181RW2snjDHnjTGHrLW/Bu4DXg50FSIiDSif9zg3k+TkaIp/879eJJPzSMQj/LsP3swHb9lBLBZMT+FGqz7+GfDlYsXHaeD3Anl3EZEGlM97nL+S5OxMhvZobDlIA2RyHv/6my9yYLCbW4ZvCOT9NhSorbU/B0YCeUcRkQbl76DHk0tEiBI1EaZSS8tB2pfJeVyYW+QWahioRURamedZxpPzvDqZpj0aI2qiTCaXAOjvipOIR8qCdSIeYaCrPbD311AmEZGr8DzL+Stpnjk9zblLmeUd9HRqif6uOP1dcUYvz/OZD7yZRNyF00Q8wmc+8GaG+oIL1NpRi4hU8DzLxeQ8pybniUQMsUikbAcdj0Y4NZVi66YOejriLOUK/Nnv3MJ8Nk9XW4yeRJRdN6xqN3ndFKhFRIo8z3Jhbp7zlxcoWDDGMLuQx1pLf1ccgNHL8+za3ElvR5yZ9BI7buigPRZlPlsgHo2wb6CLN23pJhIxga1LgVpEBMhk8jw3ehlj3A56dj5LwbP0JmJ0tEXLdtCTcxl2D3QSM1HS2Tzz6Rx7BjrZv6Un0ADtU6AWkZaWzRY4P5tiIpkDVnbQ3e0xYlFDJu+RXFxacwedyWarsoOupEAtIi0pmy1wYS7F2csZOtpiJBfzLOW95R10ailPeiFLWyzKYs4Si1j2DfYwn82TTufYO9DB/i29VQ3QPgVqEWkp/g56fG6JWDRKNBJhYm6JGzrjdLVHl3fQYChYQzqTY9+WHjI5j8lkhh2bOjiyry+wrsONUKAWkZaw3E1Y3EFHI1Em5korOQzTqUzJDtpj7+ZuFnIFxmcz7NiU4PabttY0QPsUqEWkqfndhKNXMnTE48s76IhhuZLjS8+c4cO378ZSvoO+MLtY1wDtU8OLiDSlfN7j1NQsT52aYjKZd80qyQwz8zn6u+IM9rTzpWfOkPc87rtxG1957hy7+jvY2d/NxSsZ2qKG+w8Pcetwf12DNChQi0iT8TzLhdk0T52aYiqVX+4mnJnPMdjTTn9XvBigC7ynGKAPDHXxT+7ez8UrS8Qi8MCNQ4zs2UxbW7Tenw6g1IeINInSZhWIEDURLqVdLbSf4rhwZZ7+7nbuu3Eb//3ps/zhvQfYvXk/F64ssbW3nfcc2kIiEb6wGL4ViYhcp8XFHM+Nzqxq997S00bes4xenmfrpg7isSjpTI6DQ13s3ryf0csZhnrauT+kAdoX3pWJiFyD51nGZufJ5fNcnMut2e6dLXjkPUtPRbNKMlMgk/U4PNRT9WaVIChQi0jD8du9856lqy3O7EJuzXbvfMHDszDQ3bbcTeiaVarX7l0NCtQi0jAq271jEcNkMnPNdu+uthgLOY/J5FJdmlWCoEAtIqGXz3uMzSY5c6m83RtgqLf9mu3ek8ksewc6uHPv5obZQVdSoBaR0Co9+ioWWd3uDfCNk6Pcd+PW0LR7V4MCtYiETmk3oX/01Vrt3t2JOG/fM8APXpngw3fsIZPz6t7uXQ0K1CISGqUBOhGLLTerXK3deya9xOFtveze3MXFKxm29rVz/+Gh0DSqBKU5ft2ISENzpXbpsnbvS+nshtq9z19eDGU3YZC0oxaRuvG7CU9PzdMWiy7voKG8WWXX5s7ldu/ff/d+dm927d5b+9p54Mbm20FXUqAWkbooPfrKGLOq3bu0WcU/+upjd7t276Gedu47OEhHR7zOn0VtKFCLSE1lMnnG5tJMplaOvip4dlW791rNKplsjkND3Q3VrBKEDQVqY8xZIAUUgLy1dqSaixKRkCnkYeIFmLsIfTtg61shen37vGsdfRWLmjXbvZuhWSUI1/OV/nvW2ktVW4mIhE8hDzNn4Pyz8J1PQ24R4h3wvv8Ib/3QhoK13004OpOhsy2+5tFXfrNK5Q7aNat0NnSzShCU+hCR1fwAPXPaBWY/SIP7+O0/gsHDsPPtV32J1UdfRZhMZvDs6qOvVjertPYOutJGA7UFvmeMscB/ttY+WvkEY8wx4BjA8PBwcCsUkdopDdBtnW7HnJ5cCdK+3CLMja0ZqD3PMp6c59XJdNnRV+DavXMF76pHXzVjs0oQNhqo77bWXjDGbAG+b4z5lbX2ydInFIP3owAjIyM24HWKSLVlF+HskxDrcAF6fgq8AnRudrvq0mAd74CugbL/PJ/3OHM5yXQqiz+4f6pkBw2u3fvBt2xfroX+/XfvJ5unqZtVgrChX1nW2gvFj1PAN4Dbq7koEamhfBamfwMXjkOkuINeuAydA+7PzFk4+lkXnMF9PPpZ6N0OlB99NZHMYq569JXH2/cM8J1fXuTwtp5QH30VNuvuqI0xXUDEWpsq3n4A+NOqr0xEqiufhSvnimmOrvIdNAYKOfDykOiFfAZ++68gm4a2bkj04fUNc+FKuuzoq8rB/aVHX600q3Rx/vJiyzSrBGEjqY8h4BvGGP/5/8Na+92qrkpEqqcyQEdjkBqH7i0uUM+ccSV4hSxY6/4+moBcGqJx2HyATO8+njt1CWPKj76qHNxfefSV36wS1rMJw2rdr5S19jRwSw3WIiLVVMi7AH35tZUAnZ4EW9xFl+6g56dh0zDEuyE/D+mL0LeL7M67OZ9cZGJsDr9ZpXQHvdbg/kY8+ips9CtNpNmVVXKU7KCx0LPNBejnvwC3fGj1Djp9ETYfgN13M7+U47mzV1Y1q5TuoD3PWzW43x191cH+Lb0K0K+TArVIs6qsha7cQWNg9ry7fegovPBVuOsT4OVcgO7dBXvvgWiMTCbPL8dTazarVO6gt/b63YTNM7i/3hSoRZqN50HyAky+XF4LXbmDvu2jLue8lHTNKzfshfQ4dA/BwQchnlh+yV+Oz605uL9Zj74KGwVqkWaSy8DYcbDe6lroyh30ycfhyKfcc5eSrk566M0u1REp3wFPJJfWHNy/+uirBEf2bdIOOmAK1CLNIJdxFwpTF1ZqobFud1xayVG5g06ed8/ZfaRsB11pa287/+Wp1/jw7bvVrFIHCtQijSyfdbtkv5Jj8Yq7INjZz6paaL+SI5pY2UFvuQkGDq7aQVd6y7Y+3nPj9pYd3F9vCtQijaiQdwH60qvQ3r1SyZHohXgnTL969VropSmX3thAgPYlEjEefss29g50cvHKEkO9qoWuJX2VRRqJX8kxN+bSGNEYpCZWKjnyS25XfY1aaL+S43olEjHesXdzFT4pWY8CtUgjKG1WiXdANHr1WujcostTDxx0Ld/pcejfD7vv3vAOWsJF/9dEwszzYHYUzjzlSu78Urv0lMtDdw64+udCbqUWum8YugYhdRHa++DgURg6rCDdwLSjFgmrXAbOPQ2RaHEHXazk6NwM0TZ45vOuFnr4Dnjlf8PIx1ZqoXu2w6H3Q6yt3p+FBECBWiRs/KFJqYuAB+lLYMxKJcfMGejfW14LfcNemD3jAvThhxSgm4wCtUhYVJbaZWbdxUF/B+1XciR6XRqk/01w1x+s1EIffBDaOur9WUgVKFCL1Ju/gy6r5BiHjhugrWdlB11ZC51LX1cttDQuBWqReimbydFVXsnh76JTE+U76EjUVXIsTcPm/QrQLUKBWqTWPM+lKy6dcsHY30GX5qH9oUngmlY274fsAuQXoP8ADCpAtxL9nxapFc+D2XNw+kdw+RRgS84n7Hd55ue/uFJqd/Jx2LQHNu12gbytCw68V6V2LUg7apFayC7CaLHUDg8WZ92wpNJmlXf809WVHMnzquQQBWqRqvIrOebOAXalkiPRB5F4+djR5x+D3/o09O8rVnJshUPvu+ZUO2kNCtQi1VDa8t3e7eZveAXo6HeVHPkMZJKrx47OnnEpEAVoKaFALRKkNc8nnHDBOhJ3gTk1AbF2N7Dfy7sddCbpLhSq1E7WoEAtEoS1jr8qHZq0lHLBuL3XPX8p5QJybhHSE+58wu2vb6qdND99V4i8EX6p3ZVRN/d5reOvXvgqHHzA7aAXLrldc3Ye5s5Dr2ZyyPr07yuR1yu7CKd/uEap3YD745faDd8Br34Pbtjn5kHPnoVYAm58GIbvVJCWdW14R22MiQLHgQvW2oeqtySRkKus5FicXfv4q6uV2mkHLdfpelIfnwReAXqrtBaRcPNncsycLq/kSPStffzVpmG485Fiqd0WDU2S121DqQ9jzE7g/cBj1V2OSAgV8q7d+9QTLr3hV3K0dbsUh1coP/4q0Qe9O11L+EJxiNKuOxWk5XXb6I76L4A/BnqqtxSREFpagLNPlZTarVHJEWvX8VdSVesGamPMQ8CUtfaEMebd13jeMeAYwPDwcFDrE6k9z4O5UcjnXGD2j7+yV6nkKC21S72xA2RF1rKRX/VHgIeNMWeBrwD3GmO+VPkka+2j1toRa+3I4OBgwMsUqRG/kmPq167Mzq/k6N5y7UqO5JjLUx/+AOwcUZCWQBlr7caf7HbUn16v6mNkZMQeP378ja1MpJZKKzn8n4nFWejYBBiXk+7c7IL3r7/rKjkKS+5+z3bYcZsqOeQNMcacsNaOrPWYfu1La/M8SF6EyZfKKznA5aH9So7KmRz+8VeaySE1cF2B2lr7Q+CHVVmJSC0V8q5RZX4KMKtncoDLQ++7p/z4q6Wkjr+SmtOOWlpL6UyOWNxVaqQmWHMmx/Adbsj/rR+B3II7/qr/TTB4SAFaakqBWlpHLgNjx12lRjTm0hylHYWVMzmGbnZpjvS4m8mx97eUh5a60LZAml8+6xpWRn8KeCuVHG3d5cdfaSaHhJR21NK8/JbvuTF3MbByB53PuL/TTA4JOQVqaT6llRxtXRCNusaVRG/5TA4v7y4MxjvhrkcgOeqOvzpwFNo76/1ZiCxToJbmUTobermSY9zN3Ojsd2cVls7k2DQMPTsgP1+8vxu23apmFQkdfUdK4/MD9KXTYGx5JUfnZoi2wTOfh1s+pJkc0pD0nSmNLZdxJXQzZwHP7ZjTU24HvdzyXZwN/cJXoW8YugbdTI72Pjh4FIYOK0hLqOm7UxpTIQ+Xz6xUcmRm167kOHQUnn8MttwMd/2B20FHYm429C7N5JDGoO9SaSz+Kd+llRxeATr6oa3n6pUcs2dcJcfhh1TJIQ1HgVoag+e5CXWTr0C8o7ySIxJ3QTs1sTIbWpUc0kQUqCXcli8UnnLB2J8NXVrJ4bd8Q/ls6Plp6NsN229VikMamr57JbxyGTj3tKvkwLoBSl5hdSWH3/K95SbIzsPcedfyrYYVaRIK1BI+hbzrKEyO4So5Zl2A7h5yH2fOQP/elUqOuz4BXs61fPfuci3f2kFLE9F3s4SHH6Avv+Y6CjOzLrWR6HN56ELOdRMmet0EvP43rVRydA+5PLQOkJUmpEAt9edfKLx0CiLRlY7CjhtWKjn8Q2StdcdiRROQS7vxo5oNLU1OgVrqy89DR6LFSo5JyjoK/UqOykNk0zpEVlqHvsOlPvzJdqmLuNGjl1YqOTDw/Bfgto+655YG6OSYy0PvfqcCtLQMfadLbRXy7hDZS6+W56H9HbQ/2a60YaWw5Co5NHpUWpQCtdRGZUdhZR7ar+QonWx35yM6RFYEBWqphewinH0SYiUdhZV56NJKjkjUTbbLpXShUAQFaqkmPw+dnnCDkEo7Civz0IUsbN4P2QXIL0D/ARhUgBYBTc+TavAn2516wgXn+Sk32a6zf/Vku5OPw6Y9bmh/atzlrQ+8V6NHRUpoRy3BuVoeunvL6o7Ctc4o1GQ7kTUpUEswrpWHvlpHoX+h8OCD6igUuYZ1A7UxJgE8CbQXn/81a+1nqr0waRAbyUPf8iF1FIq8ARv5CVkC7rXW3gLcChw1xtxZ1VVJ+F1PHvqFr0LXFsC4jsL2PreL3qI8tMhGrLujttZaIF28Gy/+sdVclITY681Dp8fVsCLyOm0oR22MiQIngP3AX1prn13jOceAYwDDw8NBrlHCQnlokbrY0L87rbUFa+2twE7gdmPMzWs851Fr7Yi1dmRwcDDgZUpd5bMw/Ru4cHwlD1150nduwbWCd2+BzQchvwiLMy4PPXyXgrTIG3BdCUJr7SzwBHC0KquRcFEeWiQUNlL1MQjkrLWzxpgO4H7g31d9ZVI/ykOLhMpGctTbgL8u5qkjwP+01v5ddZcldaM8tEjobKTq4xfA22qwFqkn1UOLhJZ+slqd8tAioacW8lZVek6h8tAioaZA3YqyizD6tEtxKA8tEnoK1K0kn3XHYM2dA6zy0CINQj91rcDzYHasmIeegMUrykOLNBDtqJuZ57mUxZVRwBTz0BPQ3g2RuPLQIg1CgbpZ5TJw7mkw1uWiUxOAhZ5tsJSCTFJ5aJEGoUDdbAp5Vw/t56EXZ13O2c9Dv/BVOPgAWM9dNCw9p1B5aJFQ0k9ks/A8d6HwtR+4Kg4/D93WXZ6HHr4DXv0e3LAP+na5nHVCeWiRMNOOuhnkMjB23O2So8WuQj8Pnc+4oH21cwqVhxYJPQXqRlaa5vDbvivz0LF2yC1CvBPuegSSo9C9FQ4chfbOen8GIrIBCtSNyJ9uN3Ma2rrcjrk0Dz173u2orecC9sBBF6znp2HTbth2q9t5i0hD0E9ro1lagLNPQVvnStt3otftmKdfhb4driU8vwSbhl1eem7MNbAozSHSkBSoG0UhX+wqPO8C9PyUm8nRudkF5cUrLmDPT7sAvdxVuAhbblQ1h0gDU6AOO3940uQrLp3h56G7h1yg9tu+c4suTz1wELJpWJqG/v0wqAAt0ugUqMPK7yq8dBoi0ZWuQj8P7Q9P8tu+7/oEeDlIXYTeXbD3HuWhRZqEfpLDyO8qjEQAD9KXWK7m8PPQhay7gLhpGO58xLV9d29x1RzqKhRpKgrUYeJ5rqX78m8A68rr8ksuD+13Fe67ZyUPHe+G/DwsFO+rmkOkKemnOiz8XXQ07i4MegXo6Ie2HleK17fDdRWe/hHc+hE3enRp3M3p2H238tAiTUyBut78ao7ZMyzPiPa7CqNxl5curebwp9t1D8GeIxBP1PszEJEqU6Cul8qjsPxddGlXYXuve24hWzI8KQmDN6maQ6SFKFDXw/LFwujKUViJXreLLu0qXLjkJtpl591zerfD3t9S04pIi1GgriV/NkdyjOVqDv8orPyS20VXdhVeOeuqOQ4/pAAt0qIUqGvBr+aYfNnN5sjMrlRzRNvgmc8Xzyoszoju31cyI1pdhSKtTj/91ZbLwLmfusoNfzZHvAs6B2DmbPlZhX3DmhEtIqusu6M2xuwCHgeGAAs8aq39XLUX1vA8D5IX4fKr5SNI/V20X81RehSWX82hphURKbGR1Ece+CNr7UljTA9wwhjzfWvty1VeW+PyB/nD6hGkz38Bbvuoe6ysmmNB1RwisqZ1I4K1dtxae7J4OwW8Auyo9sIaUj4L07+B0Z8Cnktz+EdhzZx1gdk/aWXTHjcbOjXu8tYH3gtDSnOIyGrXdTHRGLMHeBvw7BqPHQOOAQwPDwextsZRebGwdBe91gjSOx9ZOfH70PvUtCIi17Th7Zsxphv4OvApa22y8nFr7aPW2hFr7cjg4GCQawy3XMa1dV85t3KxsPRAWS/vRpCaiKve8PKu6mPwJhi+S0FaRNa1oR21MSaOC9Jfttb+bXWX1CBKLxb6rd+2YpC/RpCKSAA2UvVhgC8Ar1hr/7z6S2oAlRcL/dbvQq58kL9/oGz6ojtQdv97daCsiFy3jaQ+jgC/C9xrjPl58c/7qryucPI8mB1bfbGwc8C1flfWRHcNurx0327YdYeCtIi8LuvuqK21PwZMDdYSblcruStt/V5KwuDh8gl3Bx9UHlpE3hDVgq2nkHcT7tYqufMvFlrPld717XRBO3PFDVPSxUIRCYCuaF1N6aGyVyu5q7xYmJ5wFwu362KhiARHO+q1bLTkrvRiYSTqLhbuGlGQFpFAKVBXWlpYSXPMT8HCZbeL9grlu2hdLBSRGtHWz5fPuh10enKl5K57yH1UyZ2I1JECtd+4MvmSy0WXnllYyLk0R2Uu2t9Fb79VaQ4RqbrWjjL5LJx/zt32c9E9W1fOLIy1uwuI/nyO9Lg7bUVjSEWkhlo3UC8twNiz7iJgqmRWdOmZhUspN58jtwgLxYFK227VLlpEaqr1Ik4h74Lx3HlYnFk9K/rNHyw/s3BuzO2iD71fZxaKSF20TtWH58HsKLz2AzccKT1ZXnLnt3+/9E33d9l5WJqDwRtd44qCtIjUSWvsqPNZOPsTl+aIxtwxWJWNKycfhyOfcu3f81PQsx32vksBWkTqrvkDdXYRxp4DPEhfAqybdDf9KnQNQD6zUnKXHHUldweOquROREKjeQO1X3Y3c9oN6s8vuYuFGFdqt+8eN0SpZ6s7bDY95VIeO9+hXbSIhEpz5qhzGTj3U7hy1uWi411uFKl/buHwHa5FvH8/EHHVHYM3wvA7FaRFJHSaa0ftn114+TcQibmyu85+t2NOTZSfW1g6inTPEU25E5HQap5A7V8wjMbLJ93NnIW+4qHphSxs3g/ZBcgn3bmFgwd18reIhFrjB2p/HOnMGZbnRSd63QXCZz7v6qL9XbSJuse1ixaRBtLYW8lCHs78yO2aM7NrT7p76Zuwaberk85cWamLVpAWkQbRuIE6n4XRp12rd+kFw7XmRS9MQ7wD9t0HQzcq1SEiDaUxI1Z2Ec79xLWA+7voWLu7rXnRItJkGitHXVobvTgD7X3Q1uPy0/4FQxOBOz/u2sS7hzQvWkQaXuME6tKRpOlJd8Ew1r667M5Eiy3gW2HH25WLFpGG1xiBOrsI559ZGUnqz+lITbjHVXYnIk0s/IE6Mw8XT5SPJJ056+Z0WM9dONy0y6VEuodg9zs11F9Emkp4t5yFPFw+A+M/Wz2SNNHr5nR09kPfTpi/BD1DsOt2BWkRaTrrBmpjzBeNMVPGmBdrsSDApTr8udGp8avXRmNcwNacDhFpYhvZUf834GiV17Eiuwhjz5fPjZ456xpWVBstIi1o3ehmrX0SmKnBWlyQvnDCpToWLrsjsCpTHV0D5SNJVXonIk0usG2oMeaYMea4Meb49PT09b/A0oLbSaeLVR2dA65xRakOEWlxgQVqa+2j1toRa+3I4ODg9f3H+SxcPOmCdPeW8rnRo8/CHR93ddLxDth7r1IdItJSwlGed/HkykXD2fOr50argUVEWlg4tqXJ8ZWLhpGo+7sbihPvFmfczI5ddyhIi0hL2kh53t8ATwOHjDFjxpiPBb6K3u0rFw0LWWjrhNyCy0d3bHJjSaPh2PyLiNTautHPWvsPqr6K7W+DQw+6i4Z3fNwF6/yiWsFFRAhLjjrWBm/9+zCw3+Wqe7bB3nepqkNEhLAEanBBefjOeq9CRCR0lFMQEQk5BWoRkZBToBYRCTkFahGRkFOgFhEJOQVqEZGQU6AWEQk5BWoRkZBToBYRCTkFahGRkFOgFhEJOQVqEZGQU6AWEQk5BWoRkZALz5hTEQm3XAYu/Rps8X4kCl5h9e2gHqvFewS1NoBsCtJT7sSq7bcFOk9fgVpE1pfLwOkfrdyPtrmTmCpvB/VYLd4jqLWBO4D7O/8CcosQ74AH/8wdhhJQsFagFpH1XTxZfm5pNMpy+Ci9HdRjtXiPoNYGK0Ea3MfvfNqdWBXQYSgK1CKyvtQEePmV+yYC1lt9O6jHavEeQa0NVoJ06f3UOEFRoBaR9fVsKw9G0TgUcqtvB/VYLd4jqLWBS3eUfn3iHe5rFhAFahFZ3/a3leeoMSWBypQHrUAeq8V7BLU24MH/sDpHvf02gqJALSLriydg3z2q+rja/Z6t8JGvqepDROosnoBtt9R7FS1JDS8iIiG3oR21MeYo8DkgCjxmrf1s0AvJZgu8dimF50EkAp7n/oUVLbnNOvfD9Fi9378Z113v92/UtQW17vlsgenUEtv6Erx1ex9tbVGkNtYN1MaYKPCXwP3AGPC8MeZb1tqXg1pENlvgJ2emsRbiUUO+YLGU32ad+2F6rN7v34zrrvf7N+raglr3dCrLZ771EpmcRyIe4U8fvpkP3rJdwbpGNrKjvh14zVp7GsAY8xXgt4HAAvUvLs4RNREwEItEsMX6xNLb690P02P1fv9mXHe9379R1xbUuv0gDZDJefzbb73IvsEuRvb0I9W3kUC9Azhfcn8MuKPyScaYY8AxgOHh4etaxEQyQ77gfo9HDHjFX+mlt9e7H6bH6v3+zbjuer9/o64tqHX7QdqXyXlMJjNIbQRW9WGtfRR4FGBkZMSu8/Qy23oTzGdd11M8GiFX8FbdXu9+mB6r9/s347rr/f6Nurag1p2IR8qCdSIeYag3gdTGRgL1BWBXyf2dxb8LzFu29y3nqI2xFDyXFyu9zTr3w/RYvd+/Gddd7/dv1LUFte4/efjNq3LUb93eh9TGRgL188ABY8xeXID+MPAPg1xEW1uUI3sHl6s+ohEoeO5Kc6zktqH8scr713purR+r9/s347rr/f6NurYg1j3YneCvf+92plNLbFXVR82tG6ittXljzCeA/4srz/uitfaloBfS1hblpu2bgn5ZEZGGt6EctbX228C3q7wWERFZQ6TeCxARkWtToBYRCTkFahGRkFOgFhEJOWOtXf9Z1/uixkwD5wJ/4fAbAC7VexF11upfA33+rf35w+v/Guy21g6u9UBVAnWrMsYct9aO1Hsd9dTqXwN9/q39+UN1vgZKfYiIhJwCtYhIyClQB+vRei8gBFr9a6DPXwL/GihHLSISctpRi4iEnAK1iEjIKVAHwBjzRWPMlDHmxXqvpR6MMbuMMU8YY142xrxkjPlkvddUa8aYhDHmOWPMC8WvwZ/Ue031YIyJGmN+Zoz5u3qvpdaMMWeNMb80xvzcGHM80NdWjvqNM8a8C0gDj1trb673emrNGLMN2GatPWmM6QFOAB8M8gDksDPGGKDLWps2xsSBHwOftNY+U+el1ZQx5p8DI0Cvtfaheq+nlowxZ4ERa23gDT/aUQfAWvskMFPvddSLtXbcWnuyeDsFvII7a7NlWCddvBsv/mmpXZAxZifwfuCxeq+l2ShQS6CMMXuAtwHP1nkpNVf8Z//PgSng+9baVvsa/AXwx4C3zvOalQW+Z4w5UTzsOzAK1BIYY0w38HXgU9baZL3XU2vW2oK19lbcuaK3G2NaJg1mjHkImLLWnqj3WurobmvtbcCDwB8WU6KBUKCWQBTzsl8Hvmyt/dt6r6eerLWzwBPA0TovpZaOAA8X87RfAe41xnypvkuqLWvtheLHKeAbwO1BvbYCtbxhxQtpXwBesdb+eb3XUw/GmEFjzKbi7Q7gfuBXdV1UDVlr/5W1dqe1dg/uAOz/Z639R3VeVs0YY7qKF9IxxnQBDwCBVYEpUAfAGPM3wNPAIWPMmDHmY/VeU40dAX4Xt4v6efHP++q9qBrbBjxhjPkF8DwuR91yJWotbAj4sTHmBeA54P9Ya78b1IurPE9EJOS0oxYRCTkFahGRkFOgFhEJOQVqEZGQU6AWEQk5BWoRkZBToBYRCbn/D3DfZN9HNTqxAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(kgrid,c)\n",
    "sns.scatterplot(kgrid,kk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "optimized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "while crit > EPSILON: # Matlab will keep doing whatever you write between \"while\" and \"end\" as long as crit > epsi (i.e. until convergence is reached)\n",
    "    \n",
    "    v = Bv # update the value function\n",
    "\n",
    "    for i, k in enumerate(kgrid):\n",
    "\n",
    "        ceil = math.floor((F(k) - kmin) / dk)\n",
    "        \n",
    "        c = F(k) - kgrid\n",
    "        Bv[i] = max(u(c) + BETA*v)\n",
    "    \n",
    "    crit = dist(Bv,v) # compute convergence criterion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi(k) : return k/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "k1 = 100 # compute optimal next period capital choice at each grid point (recall that you have save the position of optimal k' at each grid point in vector dr)\n",
    "c = phi(k1) # compute optimal consumption at each grid point"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "eaa82c54e884903e943417a673a7150460a0dc8edf0085058b9432fcba54d5ca"
  },
  "kernelspec": {
   "display_name": "Python 3.10.2 ('venv': venv)",
   "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.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}