consumer theory - 4. optimal consumption, algebraic solution (CES).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import all functions from python labrary: sympy\n",
    "from sympy import *\n",
    "# import python library: matplotlib (for nice plots)\n",
    "import matplotlib\n",
    "# set sympy output to print pretty math expressions\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declare variables we are going to use as 'symbols'. We will add more along the way, if necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1, q2, p1, p2, y = symbols('q1 q2 p1 p2 y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Solve for optimal consumption analytically"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the optimization program"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our optimization problem is\n",
    "\\begin{align}\n",
    "    & \\max_{q_1,q_2}\\quad u(q_1,q_2) = (q_1^{\\rho} + q_2^{\\rho})^{\\frac{1}{\\rho}} \\\\\n",
    "    \\text{s.t.}\\quad & p_1 q_1 + p_2 q_2 = y\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the utility function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFkAAAAdCAYAAADBwsDVAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEdklEQVRoBe2a8VUUMRCH73gUgFiB2AFqB9gB2IHaATwq8GEHaAnaAXYg0AF2oNIBfl9I1t3bbLjd24X74+a92WQnk8nkl8kke+/md3d3sw2Nh8B8Pt/D2jl8BLa3Wt72saFxEADgfSztwgJd0QbkCorVK0TutVYAu2Fsq/G2eZkEgQ3Ik8DaNLoBuYnHJG+j5mRy0Q5ensK/4Tfwp5SnqK8ljelzPPgOmKiH3xnvF8z/+4xHg2ncW5Qt+07fK3hffUv477J9n0oPH4s+0z4YjzSnRroA+eMIDkU/in0vMRxOWHp7R9xBbnSvJS3p817U65wD7XclrkBG6RArzwHpe6e1csNHmr/VVPasYy9cyGvydao+6DP+/9DhiE/Wd3TmJQ6pgp5G2xWKjdSx7Hvsb47YSX2o+9XzLb2vW9nXZ/RNK9X8+swngSwgH/p0rOvS112A6H6RqBvFN/Awp2YzD4/B/iQ/SiX2e/msP/B5yWZXm2HuVvCAeoaRQUR/F8kT9SdsengFn2BzUKqIW9MD5zN2JqEhPosTzrwozQsdz7WEgYsSziijZnCqwIgguJUOrY/B2oKPx7DVZWOIz7FP5w6j/Qw+cEzI25VXuJkH3xEckjtlllgdT9hz2Luf7K3hGNaQZFm0EbQe8YFvU/jsHN8WpmGgJRxeo1ftZKOwtDpGutukui9S9xaBvWrFVtoJ2qkztleKZPpP4jN2zcvZucYxQ+RGXC6QhW+GbSo78B+4RUYsQgH1y+3XgkLINcgtzcFrQRP7LE5ikiN38y3jGyAe/GcRm6BrlIY84grUGbkHGqKW3D4aaej3fY/2vYUssvblRbnvxWsh7ZP5jG2BbOHhvCEjt9rtdSxsLIFsW7UFojFXSavZhakbH1rH9uB0Qd/JfHbO2s/NC/lNTq7MdGEaaG2BuO2Uu0J1ciCXMyV4r4DqvYM/In+y1LGsz1Hv1HlABo30Ht9v76udT+eZTa3IO+ctyHZKA1Ft0WIu9jZiHg6Ew26hADylTqwDFX3GQVOdn9SBmIMpxgvAy3tJ53OXlkXbQbm0QFtoCJg/SzYodjJaqwXAGU9XAa2iGL1r2I+G7ODIH42W9RmHPjCXFBj65/3WK58BUyKjdXFnl/RDm5FsJ1cyR0btVwb3K0ZKQPYe6L77ozyX8dkovhzgjQvjgvQjVt8OntrhTud7F6NjJNPc1kHuYdWZ/HN9umTR1ihffNjq9DmNj47AFX2n3R1d1En2FkvTheQgVY4KkvzDr50qH+dVRpF6AD10CC07UNHnmCIMkM6DKw50Qtk/iu2cUKdqCsje82o6Xo+y92Pko0VyGm+MEr9KPhudzrv4ayHtQW+oP3WQvRk07sTJKHJX0NPX3GJqaQGNbK1A1ke402faBK76sKHuoZcNMuQPLkTCKldWINvoIHDn7xg5A0lGv7UCOfmVK+M8BU5gE3v4tyIamTk9C37Odk4Wfk/GyGAip7kw5nNPXh3+AntATPZbMPZXInw2hbTu9Pjc/OvPSqP87/wPurW1pmZgyL8AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left(q_{1}^{\\rho} + q_{2}^{\\rho}\\right)^{\\frac{1}{\\rho}}$"
      ],
      "text/plain": [
       "    ___________\n",
       "rho╱   ρ     ρ \n",
       " ╲╱  q₁  + q₂  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define alpha as a rational number of 1/2\n",
    "rho = symbols('rho', positive=True)\n",
    "\n",
    "\n",
    "# Note: in sympy syntax: use double star ** for power\n",
    "u = (q1**(rho) + q2**(rho))**(1/rho)\n",
    "# show the defined function\n",
    "u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the lagrangian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "lam = Symbol('lambda')\n",
    "lagr = u + lam *(y - p1*q1 - p2*q2)\n",
    "lagr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the first-order conditions. \n",
    "\n",
    "Note the syntax of the solve function, especially the use of the brackets:\n",
    "solve( (eq1, eq2, eq3), var1, var2, var3)\n",
    "\n",
    "The equations are by default treated as equal to 0 by the solver.\n",
    "\n",
    "The solver does not seem to be able to solve the CES case..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "foc_sol = solve((diff(lagr,q1), diff(lagr,q2), diff(lagr,lam) ), q1,q2,lam)\n",
    "foc_sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Eq(diff(lagr,q1),0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Eq(diff(lagr,q2),0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to directly use the tangent condition and the budget constraint to solve..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve((diff(u,q1)/diff(u,q2) - p1/p2, diff(lagr,lam) ), q1,q2,lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is a list with one element, because our solution is unique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "foc_sol[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To access the first element in our solution, which is a list by itself, we use the list index again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1star = foc_sol[0][0]\n",
    "q1star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "q2star = foc_sol[0][1]\n",
    "q2star"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substitute optimal consumption into the direct utility function, to calculate the optimal utility level attainable at the given exogenous variables. \n",
    "\n",
    "The resulting function is our **indirect utility function**, $v(p_1,p_2,y)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = u.subs([(q1,q1star), (q2,q2star)])\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate MRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFEAAAAaCAYAAADPELCZAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC/0lEQVRoBd2YjXHUMBCF48wVcKGE0MFl6ODSwZEOIB0kQwWZowOgBOiAUEF+OoAOGNLB8T2j9diOZelsKxa3M3v625OenlYrWcVutzvKQYqiWILjAr0E09kYTFP2FYNjEWOU2oZJrxhj7cYRmYNlyr5iQWRBIp73COBHCNjEAvfZTdmXb4x2/SAS3Xb5QGe/0TfojQPf7v+gyr55DyIRZn6g70Sc2z4qnxwUY92T6Zz3cbetvxbSrmi9r3neE+WlWyX/H//zlr55D/HES/iQmpwqA6ki88h55lb5gOgU/hWwmbx5BD7vvPci0XmbSLuvze4t+W9Wdh56buXc0iH4QvPei0QIWYuUmteJUNWNutepz8yld94FhETjZ0U+YfwKvUO1fUXetZFKfpDQrxZD20VgdWf8jP6k34+ke8mUfdnAwXmLxFil0wd0E2t/KHaheTc80a3iNX8qDwnSG/Q9egshus7IbU/IWzvFl5UQxhRogvM2b2FwbaU/6Gmt7itliuWW1zZ7sLY5UsbvxZgCE2MG511uZQyXqAi8qgOhLBJnJc7wgCNbjEaiDgzwNuMjdSJ2266fowyObDEeA06iJ6jbMud+XOzR6n+v18+YzxbjArJEVBdZdjeqyHW2k7z57bMYsRidnR5GJLo2SfSN//Qvm+Z3Ueu2/QmmLxE9UZUCQAXYklhSkT6H9GIEkEJP9UkKZoUAXcteJwWr+IZoy1aHCnlda9TwLB5Spzc/XYTLePpSKWMGMWIjzGvDRF7eqLqV1aVIzRPldV9YOb3QSGzFc4mHwhSDUV5Y/67X/9JL18owaumJnrZZPLGNpQ+j2WKj16Tku8ZOZ8ZqiF5hqnjYaMmn0IvRxXAtePLHER+JOkCqUzkf3hpIvBghULFQXniGVyY9mYWoQSKDb1GdZjp9NyrLKCcJYaS9JBDyzkUg5ZWrSzcNix+xKUiyiIldeMEmAnUY6jpmqmvOsst+qrrGKw6DecWt5iRvft5BRjaAUZ+pz+6wkFWM7Lr3738Bn5XdoREGxhsAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle q_{1}^{\\rho - 1} q_{2}^{1 - \\rho}$"
      ],
      "text/plain": [
       "  ρ - 1   1 - ρ\n",
       "q₁     ⋅q₂     "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MRS = simplify(u.diff(q1)/u.diff(q2))\n",
    "MRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "could not solve -p1 + p2*q2**(1 - rho)*((-p2*q2 + y)/p1)**(rho - 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-8088f2015384>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mMRS\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mp1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mp2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mp1\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mq1\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mp2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mq2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.7/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(f, *symbols, **flags)\u001b[0m\n\u001b[1;32m   1171\u001b[0m         \u001b[0msolution\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0msymbols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1172\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1173\u001b[0;31m         \u001b[0msolution\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_solve_system\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msymbols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1175\u001b[0m     \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.7/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36m_solve_system\u001b[0;34m(exprs, symbols, **flags)\u001b[0m\n\u001b[1;32m   1935\u001b[0m                     \u001b[0mgot_s\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1936\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhit\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1937\u001b[0;31m                     \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'could not solve %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0meq2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1938\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1939\u001b[0m                 \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnewresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNotImplementedError\u001b[0m: could not solve -p1 + p2*q2**(1 - rho)*((-p2*q2 + y)/p1)**(rho - 1)"
     ]
    }
   ],
   "source": [
    "# try to solve for the optimal consumption choices using tangent condition with budget constraint\n",
    "solve( (MRS - p1/p2, y - p1*q1 - p2*q2), q1,q2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}