consumer theory - 5. optimal consumption, algebraic solution - CD utility 3 goods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import necessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1, q2, q3, p1, p2, p3, y = symbols('q1 q2 q3 p1 p2 p3 y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Solve for optimal consumption analytically (CD utility, 3 goods)"
   ]
  },
  {
   "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,q_3}\\quad u(q_1,q_2,q_3) = q_1^{\\alpha} q_2^{\\beta} q_3^{1-\\alpha-\\beta} \\\\\n",
    "    \\text{s.t.}\\quad & p_1 q_1 + p_2 q_2 + p_3 q_3= y\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algebraic solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Lagrangian is\n",
    "$$ \\mathcal{L} = q_1^{\\alpha} q_2^{\\beta} q_3^{1-\\alpha-\\beta} + \\lambda (y - p_1 q_1 - p_2 q_2 - p_3 q_3)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which implies the first order conditions (w.r.t. each control variable and $\\lambda$):\n",
    "\n",
    "\\begin{align}\n",
    "    \\frac{\\partial \\mathcal{L}}{\\partial q_1} &= \\alpha q_1^{\\alpha - 1} q_2^{\\beta} q_3^{1-\\alpha-\\beta} - \\lambda p_1 = 0 \\\\\n",
    "    \\frac{\\partial \\mathcal{L}}{\\partial q_2} &= \\beta q_1^{\\alpha} q_2^{\\beta - 1} q_3^{1-\\alpha-\\beta} - \\lambda p_2 = 0 \\\\\n",
    "    \\frac{\\partial \\mathcal{L}}{\\partial q_3} &= (1-\\alpha-\\beta) q_1^{\\alpha} q_2^{\\beta} q_3^{1-\\alpha-\\beta-1} - \\lambda p_3 = 0 \\\\\n",
    "    \\frac{\\partial \\mathcal{L}}{\\partial \\lambda} &= y - p_1 q_1 - p_2 q_2 - p_3 q_3 = 0  \n",
    "\\end{align}\n",
    "\n",
    "And we jointly solve for $q_1, q_2, q_3, \\lambda$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rewrite the FOCs as\n",
    "\n",
    "\\begin{align}\n",
    "    \\alpha q_1^{\\alpha - 1} q_2^{\\beta} q_3^{1-\\alpha-\\beta} &= \\lambda p_1 \\tag{1} \\\\\n",
    "    \\beta q_1^{\\alpha} q_2^{\\beta - 1} q_3^{1-\\alpha-\\beta} &= \\lambda p_2 \\tag{2} \\\\\n",
    "    (1-\\alpha-\\beta) q_1^{\\alpha} q_2^{\\beta} q_3^{1-\\alpha-\\beta-1} &= \\lambda p_3 \\tag{3} \\\\\n",
    "   y - p_1 q_1 - p_2 q_2 - p_3 q_3 &= 0  \\tag{4} \n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve the 4 equations with 4 unknowns, we want to write two of the quantities in terms of a third one.\n",
    "\n",
    "As an example, we will try to write both $q_1$ and $q_2$ in terms of $q_3$ so we can substitute them into the budget constraint (4) and solve for a single variable function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To do that, we divide (1) by (3) and get\n",
    "\n",
    "\\begin{align}\n",
    "    && \\frac{\\alpha q_1^{\\alpha - 1} q_2^{\\beta} q_3^{1-\\alpha-\\beta}}{(1-\\alpha-\\beta) q_1^{\\alpha} q_2^{\\beta} q_3^{1-\\alpha-\\beta-1}} &= \\frac{p_1}{p_3} \\\\\n",
    "        \\qquad\\Rightarrow\\qquad &&  \\frac{\\alpha q_1^{-1} }{(1-\\alpha-\\beta) q_3^{-1}} &= \\frac{p_1}{p_3} \\\\\n",
    "        \\qquad\\Rightarrow\\qquad &&  \\frac{\\alpha q_3 }{(1-\\alpha-\\beta) q_1} &= \\frac{p_1}{p_3} \\\\\n",
    "        \\qquad\\Rightarrow\\qquad &&  q_1 &= \\frac{p_3}{p_1} \\frac{\\alpha q_3 }{(1-\\alpha-\\beta)} \\tag{5}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, dividing (2) by (3), we have\n",
    "\n",
    "\\begin{align}\n",
    "    &&  \\frac{\\beta q_1^{\\alpha} q_2^{\\beta - 1} q_3^{1-\\alpha-\\beta}}{(1-\\alpha-\\beta) q_1^{\\alpha} q_2^{\\beta} q_3^{1-\\alpha-\\beta-1}} &= \\frac{p_2}{p_3} \\\\\n",
    "        \\qquad\\Rightarrow\\qquad &&  \\frac{\\beta q_2^{- 1}}{(1-\\alpha-\\beta) q_3^{-1}} &= \\frac{p_2}{p_3} \\\\\n",
    "        \\qquad\\Rightarrow\\qquad &&  \\frac{\\beta q_3}{(1-\\alpha-\\beta) q_2} &= \\frac{p_2}{p_3} \\\\\n",
    "        \\qquad\\Rightarrow\\qquad &&  q_2 &= \\frac{p_3}{p_2} \\frac{\\beta q_3}{(1-\\alpha-\\beta) }   \\tag{6}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to substitute (5) and (6) into (4) and solve for $q_3^*$. \n",
    "\n",
    "After that, we can substitute $q_3^*$ back into (5) and (6) to find $q_1^*$ and $q_2^*$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We turn to the sympy setup to solve through the entire problem. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the utility function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHgAAAAaCAYAAAB8WJiDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE3klEQVRoBd2ajXHUMBBGOYYCwlABoQMIFUA6CKQCSAfJUAEDHQQqYKADoAIgHYQOgHQQ3jOSR/adfZbu7HDeGZ2tlfRJu6td/fgW19fXt/4nWiwW+4znaTKmPcb4NsnP4nUMOcHcQznPSSfo7JGKWvxPBmaADxnTAWN65+AkeH943Id3VTFm8FMiJ23O0cFJl/gBMzqGBn5g3dtdDW6If9wyrgP+PifjBr2WyGlk6yR0dEEy0v1MK91JMzf5zgxUgG88T3neIx2Rrhh0FWp4n4wYwxs6uyQZ8u4xhjM7h+9ysVEkmVrOIg9WAaSXpFPfFV7iXYWUkt76hfSJpHKrcARm78yl3laJ/uz7A4Z8R9Ij9uE56aSX/x4b/U4qZ7aBR1TAA72D9DMoV2N/IHWuOxupeUVjZDuHbai7SIo/835MmYZJ+UmVrNdJ5cwK0T0KOKFM4behgFRbhuqiqMB43LDV0SUFbb27IYnrljvQF61yy4wih9SrQrXlhfg2XUUNOYOe25HrAL6TLaWPOkPKWHqnwq2hicbuaI/S+uSd2fLfpPycd9pqxNN2G3g/SE/b/LHy9OWZcT/FJ+9Eke/6O1hXq+qKQSqSk3afV2G2edRz73IZ+bkh2gG2vfQ3PPmvSaV0QMNqWx8BmK2ue+6gDdVTkbLVnsMYlOvYzhnHVcibLaXJ5cwK0UgVFVCFtFUKUBEF0usl5+DFTYyKdS2e+oLjGf2eMQ6N7MTVm83vkZxwyu0msJRGkzOM2f2KEdWNofuJy6yLjgDiOmTojAp4G8HgaZRsBdDesDW1MRnutLSJnLQ1RB/mjrjhwYA4czVg9ELDrl71BXB3l87gpV0t/CUe9Wpah1tXzHgZAzOj+5uoWkXN7I4xTrVxoGHcLNWbDHgfSVQp31ysw6XcsJW1kVqHucl4x2pbIuc2xlJ5MN7gmqcxXwPanimuu0U0BJf+svCHYBYNduRGuXJubTjOEsgFmdemp8Lb9PizddyxxtqWfS75eEzygN84joQ1Ts+uD9d6D8krSjdZQ2gQ7hCgpM4YmAn8vF7vaDREahgyiOiarFtXhqdetVaGMuv3UgauWK8CWDyDvqDfq3YHGZjiVB8reIr/mOTyk7Uc0GbnKd1Ft9dez4S1QoJyLlCyisuhXlyAvAGrd+HgG9aNEI2Lj1aHazGpX5+jwfQM+5V0t4Uz++zt4Cl6afQe71k9GunBjbCdo40MXEN+FS0CvvfHHtSNGA3KwPR454eKSN71tidFLJv1M3qw3vo+zHQFjsqo199CLQzB1Xu/Z+CvxWQi1JEn4Bp1hnx4yBjGjlRFGdU5OH0ydD0Y1soylVVfZq+q08Xrw41tqKMhBuP3YVJmaDbcL13wx/7m/oy7aHTQIK/E2l7QqFCY6cUNYdkJlPMvjk5MjOf15xPSYRKdCoe+m826DLzR+tujik5cDOAeQO99hGGWdtAlmLYJWF6/+i+UpXW9B3cWRQ0DowCVYEjzaHFkfhtSrsOlvDIuxvCjup/lHgZeZ/ddmPA9q/8hpcaME+agE3CuBSh0aQ1exUP+6GFOAK+/PM5svLYFXDdzGiQmsYs/sNO2+sNclIO8ewpv5YoxI9auPbM+F6KgrZPeBqgRo0EoctFgZGSC9/qh/ldo9pjnGZjxdJCBtttV/wK4QgiaM9w3/gAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle q_{1}^{\\alpha} q_{2}^{\\beta} q_{3}^{- \\alpha - \\beta + 1}$"
      ],
      "text/plain": [
       "  α   β   -α - β + 1\n",
       "q₁ ⋅q₂ ⋅q₃          "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define alpha as a rational number of 1/2\n",
    "#alpha = Rational(1,2)\n",
    "alpha, beta = symbols('alpha beta')\n",
    "\n",
    "\n",
    "# Note: in sympy syntax: use double star ** for power\n",
    "u = q1**(alpha) * q2**(beta) * q3**(1-alpha-beta)\n",
    "# show the defined function\n",
    "u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the lagrangian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAAaCAYAAACpbxvnAAAACXBIWXMAAA7EAAAOxAGVKw4bAAALpUlEQVR4Ae2ci5EUNxCGfVcEQJkIjDPgEQGQgYEIDBlAEcEVZABEgHEG4AgwzgBngLkM8P/NSYOk0ei1M7N3u+oq7ejZ6m5J3VKPdk6+f//+U4cugS6BLoFDksDJyclN8XPf4em6dN0rJ30Q0TX4FM7rEs4jhaeS2W0EddINBWLo0CXQJXAoEpCiuyVe7ki3vbE8Ke+b4r8o79zmXfVnC59q81oyeDrHu8FpDSyG4lfqns416PldAl0CXQJXVAKPAyOB4vv7kIyEGZcWPjlpzYJk9I8CJ69/3UrX3ESPdwl0CXQJXGUJaEeMIvyk5zM9byj8pnAu5Te4UBTfDETDS3X2RQFXzg3R8JzOlY8bbKeTzdZ8rnKiMEwgkw5dAl0CG0jgKq450fxS4YnCM+JWTIqjWFuB08NHhT8VUNKDm0U4kztp1VsU1B99v5NBeKPADv2m8jBewJOLx06/m/K5uKEwwsBH2KFLoEtgOwm4imi7Xht7WlGR/spuXeFfo6QxGu8UZv3yjSzMNhNvr1WIC+cfp9IHxR+rDAXv5jtVqqKb8jm6nsQAVvy9AkJ+WEWyqSwcHPM4Yh3c7YIWefQ2XQJbSUBr7qPW3y3WoOLspi8tiMY5RfpUZSjRJRSpyz8uqKZTCjJV2/G04yIN4rz4tX59bgz9HpRTxqnmgeoNLijKG/HTNAYen0bO4UnqjvIxWi68F03ji3+3wMZHQ6GK50KAgfim50uXGVs59VQbBuKF2m3uC0zR1cu6BI5FAlp7r7QOPyt8ZD1fYr6LFWkND0YHfY20YRc/KudI+WyW5IjRejBbIV6ALgyN3X/KI//MbdKCv5RP4Z6cotT2g/Jr+fFvPQkBk4vdSIsPDavLTqFDl0CXgJGAFuZ9hZb11CpD1mDJDjiKfyN6ixVplMj5zDsqGq5z2irih/cC3HjCBbUVYCTGnbxogN/HdI6ONWmSrbA5n+OJwqGYifabmKk9wj6SECYWzMHbo10CxygBlARhE9AafKO1i1fgOUqpodMt6LWKdHDVGMXpKdJG2nETvRY+a5jhhXcVW7vC8cw8Fx0YC04SN0UD6esKGC743sU9uBqfhmb0OKcw3nthD75MDIUYwtfJBKNyETOqD1Lrn1O0Q5dAl8AeJcBaxL2T9Dvvkb7VFKn0F0aIsDcQDch/smlW/iSvlci1+DS0T9x0E0NhCGeCcWWt9L4vA7/K0c4YIXxqWFEEzS4BwwTcVfhdzGHYqkB4sfYIxLY9U5ydCP7dvU400TDAMfOOADr/g8uCOYk75Yvm5bgzlmxYB2+Vx9oLgbXImtmLocitLdHcpEhzeEMhlKTXwFnS7x7rtG3oNWj4zbwgJlCifATqWVgWS6veZ4UnsbJd84SXF+t6DDeyuJuMS2yglzIFXs549OfSaoOh4S/9HAktLm58KVmHa836oudoeUeuh8C/eOAmYNE6CueSwz84vrnlSmNAlDWdr6bsc6wsl6e2zfSCW7DK2srhVTkbyfs5/tzyHE637mWJt/C5BO1zJwr8agA7+HEXM+TEf9jd2DbxGg25xtp/Mk0xXtxNdt1hGA58fgOoPnRw5B4/ZnVR8uPX1MEonAlXaF3Hk0QJrh9Yl4/V8G5ofWGoQE7A5KRVwnsprosu1vut5B+eUXCcDpkDnDQZ33E8lR7+EatHcuxNv1lc4FsTAv45HYRzlby5UzxrETlsCmbuJOXbQlAJ3nCsc/2U4Mzh2Ed5LZ+L0aiOvd24EDPBuGeLAmaLcCusE6ZVh9151JorH3ycOFDqpWHok7a2L8WhZTxNkC/gREGEPthRQDOBY7rHl02rjJczSvrlyoMHu4MvwhXB0cxrDJfNg16FFO986GvkV3WHF1BuHnGTr+iPuibf5b0IVwTHYrwbmmrGHuU0jJ1pyxzwduAmv2Tsi3BF+B9krn7DOY5sCWE+ae6vj+PmxlXm8k9771Si9CTPtlcZ81fJOG7yBYvS6+Cc9Ku+oHUcnxRdsTJD66J418AZo/1Q8rxJKuGxM0OpD5NUT2aUpzhijKsOEyFqKGL1a/PArQAt4+IBhwCD5ikEpdkNpgwFtHruKqXhG4QeD0oncdXy0VIfmgxts7yb8pF2pS0/npFXfpb3UlwtvLS0ET0l/KMYR14VZwMxcb0or4T/IlylvKhP5pCn5EvbUk8w4V95gyHgGcNl2njrIlYvlqe2zfSqbVa+hic2FbjOJmM0Q1MR3ljbubxSWufaH1v+NQlsAB3FmHzsph5KCOcXucPLsEeK597Wcyxm8NeC4ejt0GXdCCyi4hd25rgJneE/E8GD0Zw7ylO8LyjhnfH5O0VgBe9ZXKl+VijL8q9x81xMogFl5/2XoJT/Elwr8JhCyctqXK52TVKX+XoeoZUygDmO+2kzKJWv6qFnhvWmZ1ZnVOAF1wvDMBslYOJ+JbMCJ3iYS8ge/FGXpvIPHk7h0AzeX4ri23cXHcdT7v4irBQwKe3gpOq1ljGxQiXOjoQBnFzlKugk9PcOi7Gg3T6qZHnXmPHhMWRhAWXP/XF3LG1ZkvdKXBbnms8s/7ZzzVNu6nEi5hQ8t4FI8l+Jy1Zf88m6CmlOvZ+Alp8jbcjfAkJavbXFnFTgvWdYL0dbWN/Dq8a4tvjvCIEydBJzIQVZnGrMJ4nsh/0+KY2ePDo41cJiISJQrK+njJVG0SDM3ImCencVFgfRhyVnFzIaIuWRZvdwTzS6ClJZ82DqwqOLC4MTU0bziDYqaeHdyAbD7n1KpYX3OVwbsW93fsVjLx5RQPcUHoh23lOMUMt/CteIdJuIp8zEF3OVEJ6KXWoY+1S5W3eReK18SzutwMtXaJGLBU6U/GGM+eNBBc4zNXznNOZbSt54OGUHHb0m7tiRY4Xd20Qu05S/l8C5SjonJCYlp481wA4+nyS2i5975bcT9KToYLfx1sFledp0YaUIdMqqeGeM1JYFgmxiBrSY9wJcDpmrRav4hwr4Fu3MWfvNI/dUVcx/AS6qbAHwwnxlXL8q2E9UeJu6gBDkRv2toUq+FcSV4K11mWZxai65cwdyJy7NCh6udlUJw3uh3ZqWFLjFEX251oqTdgIM0OzL6RC36jKYNfVT99GrcIW07Jqu4V11MRLjLRrFh514igbVifLegivVT2uZ6MiOvepw4uRl5zj3FEcWTJ7kf3sop56lT/FmXBZH+BRO5lDzy+wIPgzA7EtglcF78fyP4F+MXtHhydftS2XNayuF1/ahOsipWA4pnCpjg4rXZbFxtHReleepmF8KGJici6qlL3ZHqd1TC063Df7ecOfglu8zXsS72f2jVM8U51PTGAnG4r8M8RPed8CV6aqpOMu/FhonJ/i0J0M6oh35f5BIgMf/jrjmuoEOQjVoLPhu0ehnV/y6kKB8cYnMASeQXU4TzfRGCPLkGylvzUrilZyY/xgiz/2a6WwWp+bFrEszg/Nwipe0aJIK7pvx38674hYuFgY7Q+8/BDG89Ksw7LZMGxRndgegOuxGvTveSjfhitHVmicaaniHB+TkhVzfqh/jvQlXrq/a8kr+UQyMPTs/Arf3svNQdWL8N+Gq5a+kvujjlD7OYcW9/3iEOFTOvPWufod1tkyLlol8bf8q2+VEkcI7yED4vevktt+5Z4pW20Z1mBussfH0assO/bmI28kKSQJEuS0yUYWHicRuioHhmXQjWBpKn8LnGhUWpGcsSvGsUU+0HC3vyPPY+bdzysjBGj82PuN/ZWwd96lyNmpVCtJtv1RcNGTXlupUG4ocXpVjJGrdr1FahQddhkEajYLi4GeCLqqLlpL7mnhOQL4k6NiHMJnQc9cTl+yu4+oS6BKQBLTucEnxQUvXBXfpZGP0A25R3IPs0NETvEvAvdMMBi/GFNebBfpp/dw6MmUDOV4MMTLGsPwies9tJ8fwXNxQHIPQOo9dAl0Cl0sCUuLs/jkFeCCFfuJlVCSEE0P2WOGraXZXTwzPpTbGhtZFH/8DgzrVSOMVzkcAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\lambda \\left(- p_{1} q_{1} - p_{2} q_{2} - p_{3} q_{3} + y\\right) + q_{1}^{\\alpha} q_{2}^{\\beta} q_{3}^{- \\alpha - \\beta + 1}$"
      ],
      "text/plain": [
       "                                   α   β   -α - β + 1\n",
       "λ⋅(-p₁⋅q₁ - p₂⋅q₂ - p₃⋅q₃ + y) + q₁ ⋅q₂ ⋅q₃          "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lam = symbols('lambda')\n",
    "lagr = u + lam *(y - p1*q1 - p2*q2 - p3*q3)\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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAABZCAYAAABoisMRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dWxIUt5KGDwRP88TBKxi8A2xWYHsH2KwA2AEOVkDADjDvE8GBFYzx27zZZgETY88KDocVjOf/qpWNurqqW7e6daciqqtKJaUyf6WkVEpVfeOvv/76mwdHwBFwBBwBR6CPwI0bN+4q7tso/rbGjJfRvV86AqtDYAq9Fc3bEvQHHU/UBr5aWugbYuCejt9HGHknJr8feZYULYGh/aPovE/K4IkcAUfAEXAEFkdAfTdjw9fqu38yZhT3L13/u+I+WZyfHYE1IVCit8rzSjr9ZEyOQNMmMRhvXw6lVTomOz/r+Kq2jYjWH6IDvaMg2jduRbHMpkgchz/jm9xrAFGeP1WQG2654Hl6R8ARcASWReCh+u4fjQX15wxev9UOSkbPz47ARAiU6O2gkWT8Sec/6PqD2sADixs6K92fSvNOz17rqHJ8Kf+LgTK+U1zHQ2y8YXlWGWtxQUFIXIz/Hsf7tSPgCDgCa0dA/RcdJ5NZlkq+UN/YGTGKZ9nw4r1OkpPB7FednyK/DgaMT5J98eUi8eFhwwhM2bbWoLf0FeLjXzoe63rvtc6tsqG8ogmZzni7mUswJb0KoMPD8mS59OI7uhRMPI0j4AhsAwH1Xxhtb+g8dbAicVdxGDGEx7vTxf/iZWPFBC8CeHRLSsLhpIdC6Tw4AqMIzNC21qK3eN1eTdleJjHexDSGG8ulxVbnaO37A0fAEdgsAurM6NBmMYBUzm0dv3NOBUxp2erxQX0XyyQW2MPyUM8YGOJ4e551pgwdk2OgMrLljwT5Uhjgaev6cZ0x5N7oGN0XFOUdvZxL9lEG/MERApV6ckRvLIK617NJ25boT6K3YzKNxYf2wsQHmScJ8bJpkwJUQWxyxa3H2qwHR8ARcAQ6BELnfUcd2yyTOpXzSWU+UuG/6PwN9wlVwVYP8sSB7SR4nL4Tjf0esDhB6vWcGBTKf0oUlk+TDeE+oTll75ft9+MIlOiJ6pJxfmhPVr8gNvfbdqzktlVIv1+23R/obdDDvgf5a8UzSYvDW/Fe01fRV/whug9EB0OuaWhuvIk7vG5Y1/6SQtOqcmKOwHYRUAeGp4k3F2fdM6Xy2GTM7PetjpQJJcZJ37v2UXHEP9dRHJbAoED+v4lPZP3ngKB4HouM1yVkH+Dfo0YQyNUT0otUSnuKS0xuWyX0U/VWtI+8x8r7s+Jz5YllO7oWPV5ewA7CyG1uvDVdNhWjeNywyIsauPJ5cAQcgQtDQP2CzdK/WUI0daLd7Fl82L61U2wwKO1n5crDgPOQDKKDJ4/77LAkBpnyI9vXOg4+hRCw403T7En5krKPVZR4eqBjX89j6bYeH+otSYwCPUmiGyWapG1F9JvqbUS35hJbKN4zW0PrIG/8nTfWis29eZAo9UaKwjfdeBvroOGn5vd0joAjcHkIhH7hvfqFxSZ1YaBm4/3Jfi6kg0/6Mjxud8X3S8XjvSM/e8CyZ9FLY5Aqv+TD84aRi5HGYEjAYC2Sm8xLyw4PcRA/eIGRJ9sQjemUXoe6wItpgTGz+MPHokf9sCQ5+PFYPT/5DTNjgnPg7Ww7ifOkXgfazduWlS/6xXqrvCc9b4F3PHbUG5NRJoR/pNQbtJWWVYe/61wVRAvdpT4PvvNWSxSB3OtWhaJndgQuCwF1NuaNX8TrZmiqs2MJgw4XI2x0eYR0en60rKL4ozijfe5cgoHy4BVqNsCnym+yKD1eEo6qUCJ7VYFnMosfDB32LtZ+g+ugJNFNMpCUjjHy6MPHiufN5k8HRBNuAj3TE2QbCm+V7oXon508KU1SOxkq5FwctJXmqB0p/ijuHK2x56LVRG/79APvZ/Hr5wv39Dnfqg6a7n27NVJYSfSzkKlmg19JuZ7HEXAE1osA+z1eqvPLHpgmEAle2ECMN61qlSGTtywMwoDcbICPeF1C/izZI16nuoQfBtPWIXUJtuQDsqO8mrGCYTCWSGne6/kLHXj4UtrhEnoyxv5c8ZP1B8L8nbBHDmykbK/9GAA3xx7kxKMUSo/ysDSSohw55D2tI+AIbBAB9Qt4GRjUnq+BffVNdNAslTWb6Z+TqxADBvj9JFg08KxU/7NBivyB32qPG7gUyn4O0sHnYKTjqQ5eTNkH3fffIPxWOIwulyo9Rs5jHdDCiOmCrhnjqoJo0Ba6Dx9DWwfLk3jERj3BVQUeZkZm89AdPundpehJL8vit8KSvqZYbyXz1H0CRtu9wGcTvJoYb+KE9XbCFDOaHWX/dQQcga0hwEwza0IXBrWTg6fSpA7UQ3gxuE/+jbWo4CwMJNvZAX5K+TWItfxSQJbsEWYllyyFsm+MFxE6Q0vnbonSiAVsR50Leo4xNeXHmTGeMKIYyCmrMxgCX7qdNPwq6jlG4tztpEr4xnpbxctIZrONHo48z45uZbzZ/oHRGU02Z57BEXAEto4A3vgDT8gpgTIGz7MD9Yly6KP4eO3oMtOJvCWPsjBQASkD/Fbkz5W9BF/bZP9cdQp2vIhgBhr3v0VEMYw/Rvf7S+VlcL3kD8iCCfKnhrnbSSpfm0wnnTTbqNnE8VYtElJ6Zjn9RlNL1vM7Ao7AhhFQv4DXg2Cd1u5u5PfE4PlEz/abkHXNAJQyUA+WpE6UDdksn+KFaLb/ZKgwlZOFQaDRfSFe1wy23dKp6LzRNV4a/jNxE/IXyi4R8wN1Si6VCUZxnd7XfbxsylgFrkOB1aNHvQfQBW+M5f1mdZWDodc3hEo/8rr/gGzAbL9U2+Mlvo0/fBvHn7oeNFrHMszZTsZ4uMB4+kJWDe4J3+IlXsOl2ngTIVsyjRuN0fezI+AIXCcCXyO2DawJECQNnkYvYaA+VSSdaMffqUQNnuViMFbkfoDfkPytZB/DZCgeY9mWp3iOU+E5FyFgjGHADQXi+wMqBg/xMQ10GiPxIEgfz31qAjonP3wsupSfs7R5wMOZG4zNzsg9ky5+PFc7icu85GsmEugkS6d9XcuW+2Z2juMMpmzxDOc4lcc4Ao7ANSFAv5DTQSUPngFEBuq4z6FTTPLyKR37jcwrpsvJQhYGMgBODfD95ee1y58le8MawODCC4c+8HblXgfDdd9jZkWTbv8s1EW3P0n5ij/ObMR1xpg9+P6pyniquKIPH0d0Uy+RDb3PCXO1kxyetpzW+id0szrcqqawsySZjRhjDUg6CUfAEdg4AjnGFKLa4GlLYBgyB4Mng2gPk9GBupeuf2tlNFm+6BOP7nMxyB3g1yx/ruwRbMWXLG2yzM4HltmHjU71A8vmQ5+KIb0tS4MrabhnfyRGFjpTs7rUeQVFy/Y8od/szyv+OK/y234/PIHgzZf88TwOfTwWY/rIY6i4U2GudnKKh4t5prrmr/qQp8nE8VYNMmIEBeyUsIaO53UEHIGLQyC3X8gdPFMG6jFQu0FJD/FGDA3wY/ly43MxyBng1y5/ruy52B6kD2MRLxx0Bpbu0ad4CdXSs6cMIwb89kH50Ikj40bxR3H7TJkXDN7K0lTfAt8HsvTZCtjkbGEwEnO1EyvvGs7UP58MOfnJmhQgqow3FYC1T7BK3t35ryMwgIAUdmjGO5DSo2oQWBHOycs0OYNnGIxSBuoxGK2/ujOWoGF8MgaUmTLAb0j+LNkrMX+t/HwOg78yY1y6Iyz338oz2orjg7Xf62jdF5lOWVFrOj8TMxizucFkmqOd5PK21fS/iXEmaRxVq5U3KxHAFUtoOpvYkfTfS0JAnSVLDyish+kRmOSPkFPZZmAMaT+m5slMx0D9gDwqa3SgHqOpAfxTeIZ3aJIwMQarln9i2cfqizdy2ZtGP8PboV+NJdQzvGksrzar/0BzsEiVQ7+3yBipslmm5WUK0/lBHociozzNcBoq58riWNIn3N+dyn9rPW/dG0UqnhnPVQY1DgYqOgPeJLINqT9Gin+VuMRCCyMG2i+ESdX+jpjm3Ndrq2fxw4DAJvav+rqme7wLuOab/pdeBuZTd/YM1HckHwM1ejU6UJ/hmbc4pwpZGIT6TB3g1y5/luwtKkA6kLUfTem7/Wwtyj5HQ2Wl1us5UiXP/6HyP5VkjPJM2U6iYq7iEs8bwVYtd3cFv7XGm82wzb1awMJ2s6jDRX7Wrvd7DhTHTIeZcYmbertgjHAuPOjInwmj0gF2hPJ80Wup54AluoVHi4mTtb8jMIQ3y0e/68j6h4MjQmURd0K2pEFDPP5VVswul/JjxO2DZO92Be8jxi+mNDKyMBCL3YxcsoxzO/KkLz/JEjGYSv5k2cU7PPyiI4eX7yXf3iASjWL9KcEbfLcUhmRM1A8TM6duLI+fhxEwW6ka01rjzRgwhobZvcDY0Ong2eh7kxhY8Yp42CHABuGhjcOT4qP6YWbDvpajfS85Ba+pniXLJ/HeTQrEV8oyNLiDf7NN14nY0QYI1j/s7kZ+MweSESqri75mDJJlDzpdNbG7UP1ZnUI7Q/UIoO/quztCOlftu7xZyg4FW97QAO32Ws6PJXffcEP2hzo+XAsICXL+IJyqDKiEMoaSYDgkGQ9DmaO4zdZzwP0HtdUWOESQbPsywmPODfWrAe3a5V9NRaycEdeTySrInF17G6qkpGLjTYWZd8kYKSn/ovJI2fH2gMujixKsUJiAx8Xpx8bqGfx/KKzC0myfQkZbPiulM1U+48v4nKIco21lTVFGKU3jyXgspTOWz+haOWPpPH7dCFj9WX2um9vtcGd4VhlvNcumVrAxUgSdBkKWdQi24f+VPAaDniulva107CnjxYCDDxGGZ6+Vd/K9ZmHwZlM4GLA0BQaGx8EG8rXwLP6KQiX/1MXZ16FVRrIOFAlRmEl8YYwn1XNhEXNkA3/eCp/T+5m8bDYHAANl0I8QjM/dXdtfo21ltaVeR814Mh7rqB3nNrpWznEKj9kCAlZ/Vp9b4HkLPDKhvqfDXnAs4vlmUa5dJnsDpahiNTDyOQOWLXiNmbczWYJk4/8vPKOIkKb7JAD3Cmx8J93PXBMRBbwLcdroUfNLXkXnO1NUAm/8wQ/eNpTdjDhddmEtPBs/ueca/tlUP7o0VagDufzXpM+p55pypswL/n2dnLI8Nst/UgEcVZ3ThExO7lFYOQaTyr9y2SdUq4sjPameXBxa6QKZzVTVLy/peePtqp/U0PeeGRq9BnRepcWIY5M1xljnMQgGnX2SBE9CfzmOuD0tXc8SxPPeoyEe8SDxNmC3AXetPKcC04B/jFlT1KFis3RgiMBccafqeS4eCssBf+ph7kD7zOqcpG941Wk7GJwYfwTa9XPh/6G7a/NjfNlr+22oHlNZKwZzyJ8t+zF864+ZUWeXAGMOPVlCrqXLtL6tql+uMd6s4FOD8yBIUniWGsn/fCABjR6PDYFvOJmgH3Vt3/LBy9bPyxJXPw4aTYN4h2+WeIcCvOMOtbAKno2ZgnMt/8zcrP4Oii/UgQMa3AQ61H0/dLNGPWcS0A94TU8ur2fWc59+dx9o/KIbdCY1HHwGITXTiXRLDaIYRiV77cyTzfYD/gsQOsS19OJhIPK/koO6qWetwloxmEP+UtlbYT8nnTl0dk55rKw59MTKusazeTaLZK8x3moKplMf+/4UHSpLqnix9saYdbSKZ6BmMIw9XvdC3ByeNwzLsXLgYz8grIhnsZUfJuY/WweGJBCPQ8YZRh1eW17FHnojeIhUPy65nvsZ7T7gRwe4ZKCd7nVyRkbwqj5WPdw2PTpVttIxy/8HeXS8Ux7ztHX9QZw3pGUCiKFbIht9iNGPSbe+ngQD020x+4UOcGOf8FifNCTTHPJnyT7E5NrjUnVW6RgXzCExhSd5Kqjm0JOpeF8zXXP+3K5h8mZN5pA3q/OUIsMwx6nOk+e/jnTMeEzwnMTlomSfFHeKppI0Cd0+qBFKY0u3S/M8wm5ydCn/eH2oy4PQQAcO6E10U1LPE7FSRRb8P1ZRKMuMIUagbZ4NarvmCaMN2fYI8h1MiMJAiMHCkR2C7pGXfapTh0kwENO8mPVSB9tLHulIlmVG+bNkn7oipqCfqrMqG0fEb0qPw+GNDrbWrDrMqCerxmEi5mLbpbiIGuPtaFBO4UIKbIyPbmQXHQyxdyP06HgxCuIwZjTFaVpdDw5GYVBhdkWH2g9L89znJ/e+lH+MBvIehAY6cEBvopuSep6IlSqyeN767aWKYErmUMd4g2ibOaHv8XyozGYIdH/cLto5XqZ+2VavNTT6NAfvp8JAhcXeXNqX9amDfPQiZ5G/QvYeu5u4PamzkqDbChEMItrDfkVpxdLNoicrln8O1uibi8PN4pyfM5oL8HPM+StmIEedupSbJRNCZxgGg2gX8/n3YCBSGpSMI3n2+ZlU3lVofH9GfHYEAg/se6CRHvAXSjiIm5PnPAlHU5fyjyf0/gjVGh0YIdkmuqKe2zDQlgoD/eRtY4Rl2kTyvjfhjiFC2++MkXCP521oQqTookC/g+f+QKeLKKVlao5Bj3e8OnjGU8Oc8mfJnirAmtKl6CyGbEhnbWEu3auBak49qeFzi3ltJaSzc0oFuFWaUfnMauw62hw6UuYnUuYXHMqH8dfRUHy3j03xt6NnGABxoCN/HT23jcyTz6RVLkYisyb25D3VmcC+k3+Kd+Oji+z9LMlzj5Wi21L+MRrYm3QUKnXgiF7jiNJ6bszGMTnpHQMijd720PD/pbQRthm8PM7R6SztbImAx+yV+OP/f1PaJ7gjywPlgV/a1MF3E4msDAygjypp5GSfDANhRNviM0v9PvIUf3PKnyv7Kb7X+uyszqqe7G+QftI1bfdfOv9d9ZY9ds4Iwpx6MqNYl1PUrQpRUMLiIMUdnU3r2eAmdAoLCr+faaoRMDDNNZO+r/JZzs3pLJfmGdiqQinmyscHbnl54J6ujzBTXJEOVAmTlrmontNI16USZnvdP0dJuHdL1sqzyExf5eJxwKCkjcZLfWOsk4YXmcjTPIiXbsIl+mNbMpqXORUGkgVMMdzA+IXOo23JhJpb/gLZjdUtnU/qrDBnnGSL0O4PLT87PVYr49x6slogpmOsidHeYtn043QyHlKWUjGL5y2mLoSGwVLravcQbJHnAK9hXIs5g8yoMR6X1fiaBtKkkTTma05yDOjgv1gIRgWeapY/zwW8GL+eSgQdHRhheP7x3j84lb737Jnu5/S6dcVPgAF9IBjgwflL51QMZpc/U/YOr439nNRZDFjJw0oTb15TZ7RHXoZac980u55srM5r2d3bTNKJ26XE8LxB6D90/F8mkf9U+n/T8T+Z+WqS01BYKrDAWzt86HfymXQAuWR/32I8G0iV5yr+VTcsFfysw5YOKtlJy65yU5bpjohV1PMRrSUjwFvlg/kShnNfdIxI2uqo9038MqjB83e6HvuMUPfCgtLgxc3yzokmAydvtE7eV6icodASg1EchwombmH5z8o+xvea4zN0ttsOtGZZjLeF9cTYuPQzNtd/IaT6o1wj/r+VDXutc+UyI2Ym96UIJS+vqJKZ8RHYpD9Lh6gymWF2S0E6sx/mrcouGqSVNyuobIyYjyqPgSM5LMlzMpMnErbgXzSYXVBXbIJddSit57UJJTnYb0jbzO0cJhFF/NDHvBE/WUZXC2ZUNn0G5bN/LrmPa1F2TGMpDNYg/1Kyx/j79WkE1qAnpzm8jKfC2WwuBMra+6i8rDTyXccbt7YEhxiexUgcwkRlFxmJS/I8JEduXAv+RYN9Od3Sga5XPQsVf0X1nIvrlOlDA38C7lOWk0n7G6Xn5Qo+wDubAaXymDhgyD6as9wRbGbHYEXyzy77SB149AACK9KTAe48agiBm0ORHucItEaAgVPHqg231jIvRQ+cV2CoHIgvfjAk8by+DQPFwfMJb3g7lw39i038TLaFMFiF/AvJbtD7+TwCq9CT82x6CkPAjTdDws+OgCMwKQIawPG48basbX2YtLxgJGK4rWbSMCcGa5N/TtknVawLI742PbkweCcTZ1PLppOh4IQdAUdgFgTCAD5XWXj7VrcMPhcGKmd18s8l+ywKdiGFrFFPLgTaScVw421SeJchHmZSv6h09vukBja3Z72MkUrY0zkCjoAj4Ag4Ao5AOwRmN95kWNhbqu2k6FGSEWIfROw9Kb+dmu+WPIeZVPbnBProTC1zvzy/nweBlrpmHLuuGBJ+dgQcgWtEYIp+9RSOsxtvcwt4SvicZ1vlO0fGftprlLmPgd+nIeC6koaTp3IEHAFHoAUCV/HCgrwCOcuHLXCtprFFnquFdgKLIOC6tgjsXqgj4Ag4AsUIXLzxFgYm/kaGj3RuImyR500A60weIeC6dgSJRzgCjoAjsHoELt5403IOb1zxiYLbYaBafaVskefVg+oMDiLgujYIi0c6Ao6AI7BqBGbf87YEGhqgvpThlvrnzUuweFTmFnk+EsIjNoGA69omqsmZdAQcAUdgj0ALzxuerS2E+8HLsAVejcct8my8+3lbCLiubau+nFtHwBG4YgRaGG+rfxkgLJf+c0v1vEWet4Sv8/oZAde1z1j4lSPgCDgCcyFQ41CqMd624nGjHn4QSC/nqpBG5WyR50aiO5mZEXBdmxlwL84RcASuFoE7LSSvMd4+BgaaMNJCmDEaMtxW89+GYzz247fIc18Gv98GAq5r26gn59IRcAQcAUOgxngzGn52BBwBR8ARcAQcAUfAEZgJgZq3TZssm2q/zbeS9Tsd93Q80cEeOuII93U8kmegSVkdxQY/W+S5gdibJeH1tdmqc8YdgSMEwh7Nx3rwpY4/ND7st8SEZ68V9/1RRo9wBNaBADYOocquqfG82bKpMbJjJ//3OzW0H4MgP+t8l8YYGiTfZ3ubT3LyHFvkeXJQVlyA19eKK8dZcwQyEXgWxgfGi2e9vD/oflOfherx77eXj0CTrWY1njeD+Au7yD1rlnRXeX4N+bj+oEb5LqLzh66f2n2YVdE4nyhd9R+vG92ccw7PgV/rXJCPsDpP4o6ty/zNrC/qiI6fGRGTEjy/z6VrH3T24Ag4Agsj0GvPrNgwwY8Dce/jCL92BFaGgDm8zAFWxN7Noly7TAxwBGNkd5f3+zEy1lg2fdPLjlv8bxhBOniOq5zyaspU9qqQzLNKeSH5fgwHbnwqazN/01WF0noyZ9WX2P5C9fWTDpZimFj8sh5RnBNH4OoRiNszE/n+mMGWGzxyHhyBtSJgDi+zoYr4rDHezGosdgFqgOyYl2Fme9z6MyYMtk+k0/EhDKj9mVaR4KWZcnhWGY8j2SjyhY67wRDl3sPECGTW13OxEw8GNLJF9W1ieJy8I7ApBHrtmUn8/ksCoV8lrj+ObEpGZ/biEUBHCWZD7e4yf2uMN5Y0CcbI7q7st3N/W8OEhBoidDHq/sH9CkMKz7yA8dsKeb9Gls7Wl/SPCcKHCByWUF9F937pCDgC60CAlQza66eIHcaLbqIfxfmlI7A2BMzhVeUYqDHerGBjpAYgGl1/tsQSKQ2TlxnWGM7yrI6F5be4c8GY+7NnIKxRtkvk6Wx9mdCaODzVwfL2K+rQ4v3sCDgCq0HgrjixMciY8v1uhoSf14wAjimCOcB2d5m/LYw3GlFxCB42lkf3dBTHPRv9v+kZP8XltMxYwnOQCU/OIi9atJR/a7Ry60s6x363b3R8p7z7F2a2Jrfz6whcMAIHhpvaKZMzDt/vdsGVfiGimcPrQIdzZbuVmyFKbwWbFRk9yrqkwRFeRAMlLyp8pUHUyugSrOgni2fJhWHKfjdkij1xKxLpolnJqi+QoJ5Ub3h9f9f5ve7j5dSLBsuFcwQ2gABt87XaJv0q/1vdvdymc38FR1EeHIFVIWCOqir7pth4C4NbhwjGie5LGbG9SDS6rTS8ZJ7BRnLx1il52MuHV5F9GaV4QcZDHgJn60v1wiTkf3Xg7TVDzQztrxVvcXkle2pHwBFojoDaKG1z/yFetV+MOPbAeb/aHG0n2AqBMM4w1hCqdPXmjkbxrxVulmQJIbwiWzHaTL4knlVR4MKG9+e6vseha/a9fdThYT4EztZXGAyoF9NpuCMfg8RaX5qBRw+OwFUhoH70lY79J5d0zWDIHmneFvfgCKwZAbOVuq9o1DB6qyaz8uKNgBljKItcaHTkPbtPQWlJh+HDgMrnNjCKDv4aRfeThxyexQwdDB3LvqOBQRkKyOFhBgQy64uZ/DPlYRmGcF+HL3V3UPiPI7AaBBgD4rfAX+uel8PiD7yvhllnxBGIEGAVhxA7CXYxmb+1xhsfMS3ahK8Bkny8lEBgwLyjxjf6Zp+eIeyib54W8Pz3Tjr/WQSBgvpiMsLhwRFwBNaLAOMAE3heJmKvG2+Fb231Zr3oOmdTItBsb2at8WYNxqzJZKHDLGlTM6Ut8pxcIReY0OvrAivVRbp6BEK7vnocHIBNIsDWKYL9LejuruD3ZkGefRY1IvNSGEP7Z37hCDgCjoAj4Ag4Ao6AI7BHwLaYmeNr/yD3osp4C4V1a7dyYRtTuTx4ekfAEXAEHAFHwBFwBC4WAdlI7H/HTqp+WQGQWhhvtvTJJlIPjoAj4Ag4Ao6AI+AIOAKHCJiNVO11g2wL483eFO2+Y3bIq985Ao6AI+AIOAKOgCNw9Qg8DAi8aYFEtfGmfW9YkXwLy6zKFnw5DUfAEXAEHAFHwBFwBC4Fgc5GavXCTbXxFlDlI6a3tabrLy5cipq5HI6AI+AIOAKOgCNQjUCwjdjzZtvMqmm2Mt7eBk7MLVjNmBNwBBwBR8ARcAQcAUfgAhCwlckmS6bg0cR4i5ZO+YsSD46AI+AIOAKOgCPgCDgCOwRwbPGW6eo8b7DH/8qxdGoWJnEeHAFHwBFwBBwBR8ARuEoEwpIpW8pG/0GqBJgmnjcKlkX5MjDg/9tZUhOexxFwBBwBR8ARcAQuDQH7G1AcXM1CM+MtcIQB90CWJhvzPDgCjoAj4Ag4Ao6AI3CVCARbiP9xfycHF1/laBZaG29mWZql2YxRJyOLbO4AAAFQSURBVOQIOAKOgCPgCDgCjsCGELD3AH5szXNT4y1YljD51L1vravK6TkCjoAj4Ag4Ao7AhhDAkfWTbKPub0Rb8t3UeIOxsPcN96B731rWlNNyBBwBR8ARcAQcgU0gIAfWUzHKFrLmXjcAaG68QVThkQ68b3e7O/9xBBwBR8ARcAQcAUfgChAIK48vJOqT1nvdDL5JjDcxy7dMOF5ZQX52BBwBR8ARcAQcAUfgChB4LRnfyxZq+nmQGLdJjDcKENPf63RXFqht2IvL9WtHwBFwBBwBR8ARcAQuCgHZPHzrlgMbaLIQG29/qNC/eof97VUpAzD/QjT90yGlCHo+R8ARcAQcAUfAEdgKAqw4fl+7XCq76cgmE939auYt3fAWxNiHdavekBDzH8QABtwdHU2/cSJ6HhwBR8ARcAQcAUfAEVgFAsFRxT639w0YYs/caPh/R4t56WzYdLkAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[ \\left( \\frac{\\alpha y}{p_{1}}, \\  \\frac{\\beta y}{p_{2}}, \\  - \\frac{y \\left(\\alpha + \\beta - 1\\right)}{p_{3}}, \\  \\frac{\\left(\\frac{\\alpha y}{p_{1}}\\right)^{\\alpha} \\left(\\frac{\\beta y}{p_{2}}\\right)^{\\beta} \\left(- \\frac{y \\left(\\alpha + \\beta - 1\\right)}{p_{3}}\\right)^{- \\alpha - \\beta + 1}}{y}\\right)\\right]$"
      ],
      "text/plain": [
       "⎡⎛                                α      β                  -α - β + 1⎞⎤\n",
       "⎢⎜                           ⎛α⋅y⎞  ⎛β⋅y⎞  ⎛-y⋅(α + β - 1) ⎞          ⎟⎥\n",
       "⎢⎜                           ⎜───⎟ ⋅⎜───⎟ ⋅⎜───────────────⎟          ⎟⎥\n",
       "⎢⎜α⋅y  β⋅y  -y⋅(α + β - 1)   ⎝ p₁⎠  ⎝ p₂⎠  ⎝       p₃      ⎠          ⎟⎥\n",
       "⎢⎜───, ───, ───────────────, ─────────────────────────────────────────⎟⎥\n",
       "⎣⎝ p₁   p₂         p₃                            y                    ⎠⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "foc_sol = solve( ( diff(lagr,q1), diff(lagr,q2), diff(lagr,q3), diff(lagr,lam) ), q1,q2,q3,lam)\n",
    "foc_sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is a list with one element, because our solution is unique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAABZCAYAAAD4mbreAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2d39XUNrfGw7u4OleEVBDSAYEKEjogoQKgA7KogAUdECogoYKTcHfuEqjgkFNByNtBzvPzSP40HntG8n/PPFrLY1uWtvZ+tCVtbcmeG//+++8XDkbACBgBI3C+CNy4ceOOpPs+kfCW+v5Xyb0vjcDqEJhCb0XzlgT9UcdTtYFvxxT6Rq5BJSY+qOCflP73MRkwLSNgBIyAEZgOAfXdd0X9nvrun2MpivtH118r7jrG+WwE1oRAH71VntfS6addcgSacWKBQfXNkbRMQn7T8W1uO7nqIpbGw6Tu/xJRG1MpML42AkbACKwfgUfqu1NjigHlz9xBYv3imcMzRaCP3mIEdQbp/EcdeGb/6kwUHigdad7peHMqbXx+M150nWVMPdQz3GNfd6VxvBEwAkZgrQioD3sp3j7pwNX/lTrKn+BV8Sx7nbWHRjIywPyh8zNk10F/fi25R13qEE2HC0Ngyna1Fr2lrxAv/+h4out6UtJV1Uc9VCJCB4R1xlLfWXc8XQA53ggYge0ioD4MQ+otnaEOZqZ3FIdxQXiyO531L94oVhaYaYNFtRwiDI7O5JXOwQh0IjBDu1qT3v4gIF7ntJmjBpWIYEyx1HfSMutE3g+MgBE4OwTUudDBzGKQqJxbOj5wLgFS6dmqgIv/Y5KPPRGP9IwOO41PkuRdQl/H5BiojF7yBym+kfx4pKp+XGeMq7c6OveZ5Eg/l+w5vDjNDoGBepINI3WvxJO1q8DIJHqbLWSSMLQZJiTIfTR0LvkJNDYy4h5+cJSCHxoBI3BRCIQO9bY6mlkmWirnWmU+Fsjvdf6O+0zA2apAvjSwLwLvzAPRqZb+0oe513NiMED+LnFY+isyTlNCc8qeluvr4wj00RPVJeM8S+KnAhu4476j7HbVk34XL3t6G/Sw6Wm9p3gmTWn4VbwP7avoKz6J9kPRwrhqDZ0GlVLjncIK9Ub0VugcaQQuDwF1KHhkeGNs1j04Ku9j6EB/Vfm5kzyMhqYX6rPiiH+ho1dYAoM+8otP5Py7RUi8c72MySVkb+HfUR0IlOoJ6UUqtz3FUrPbVR/6uXor2gdeVuX9TfGl8kS5Os+i+ZdoYwthfHYaVK1LfsqIZwrLtVejUz4HI2AEzgwB9QtxNvvdEqKpU6tmmeIj7oE6xQaDRT2DVT4GgkdkEi28XtwXhSUx6CH/PQm391p4wI43/IonykvK3lVJ4umhjrqOu9JtPT7UW5YYPfQki26SaPR2ldDmclS9bdAecos9lO7BPKDV+h0qVR7fnOINmL3GeJDbEUbACFwMAqFf+F39wmITrTB4srmaPRZxCaK1DkJaeKU/wzN1R3leKZ69ENBgX1HnbFPPD8LSGBTKj+GJ4cQARcCALJa5yqmfpWWPfMSz+MFbijzFxmGkMeQc6gJvXwyMmb0/lip61A/Laa0fnNTzo99YikxwDrxltZM0X851oD1qu0rLFf3eequ8Rz1UgXc8W9QbE0QmaZ9y6w36So+H/kudD8KBQaUMFBI/4tlbOQ5KcoQRMAKbRUD9Al5rltu+VGdyvaQg4gWDCONodNf+Mbn6YKA8eE9GG3ThL1d+pXuWO1AckzuUuZr6D/xgfLyRfLyBNVoAW9E8WEpqFqB0jJOjfSw10It60vrBSaXhefbeP2RR+tnbSROr0nvx3VtvlfeoQVXKSzO96Md28IP05GAydrOZQffPQ9zQTVwtpB1lBIzARhFg78ArdSKLGlMBO3hhgyiDxVEv1chYF2Eg/loHXcXzCYchOC4hf5HsI+PeRg5+MBjGDrnLh3x0svbUqk4xdnp/LFW0WEZjnyADdmtQmt/1/KUOPGE5+rOEnrTyPmPkpP2BcH8n/BEHO+nAoLpKBaWidE+F4tbPqbA0u6+NgBE4QwTUL2AYMND03sg9Jizqm+g0WeY56UkYq9yeGPT50vNJlnPkD/wySA8OPWXvVa7K+l7HMx14Q+ug++abW98Lh86lPqXH8HiiA1oYFlXQNWPcoCAatIXqY6nQ1sHS2kvxM4fHFJmjJ+uoHDl6cpTAAg+FJX1Nb72VzHP0CRhSdwOveyjtGVR6wvotYQrLf0fZv0bACGwNAWZjRZOsMNAcHdCUJnfwbMOLAXfyb0AlBRdhINlODrpTyq+BZcw3tItkTzDrc8mSFltN2GxeGT86V56+SCxg2znh13MMnCk/5opBg2HDwEpZ1SAe+NLtpOEPUS8x3OZuJ4OEH1lvB/FyJHO0jx410zQNqrge3Wn5Nwn43ggYgbNHIO4byBK0YEA7OXgeKZA+ig9edi6RHMnb51ERBiogZ9DdivylsvfBN26kfqE6BTs2m0ejifs/E6IYq7xkcBCUl8HunD86CSbInxvmbie5fG02nfQy2kcHE7qbUSopIrOBpiLHxz4bASNwgQioX8A7QIidyO6u4/fIgPZUz6p9ImTVNYNCzuDZWpI6Nb4Lw9Ifs/WDvQytmXpGqpwiDEIx1Zeedc0AGD/38FbXeDP4f7BNyN9TdolYHqhTcqlMMErr9L7u0yU/xipwbQussjxuPIAueO9t6FY5GF9N46TvhyHrj04GzOplxgYv6W36scw0/th1qyHZlWHOdtLFw5nG0x/iYb8rjOslytqg0sO43Jcq8pliYbGMgBHIROAe6eJgl5Ena0CL9DIGz2NF0qlV/B1LNMKzUgy6iqwH3Q3JP5bsXZi0xWPAxmUVnjPRf8FFCBhIGFVtgfh6gAsJMEKIT2mg0xhue0H6ePQtMT2HztGPpYou5Zcsy+3xcOIGA7AyPE+kSx/P1U7SMs/9GgMfvWTZr9a3q0TqqADpTCB57EsjYAQuEAH6hbrDyJA/e0ALtBg80z6HTirLG6Z07F+J3iNdThaKMDgx6O5tthbHa5e/SPYRawAjCG8V+sBbbbUOhuumZykWTbr6WaiLaq+L8vX6mGskHM4YmHvfZ1QZzxTX62OpDdo5t8iG3peEudpJCU9bTxv7KPSzDjfrq521hdUeEyaPfGkEjMCFIlBi4ABRHNDi8g0G1t6AxsDWwLJz8Gyka97GMvbc7s1EI9yXYlA66K5Z/lLZR4C7+ocOloj5HiL7etGpZmDJt+2zGaSPS6rgShru2W+H4YPODFmFqbxnohX3z6Df7Pca9M1GZBEdPGbgzde48dC1fXASA/fAs6a4Y2GudnKMh7N6pvrmExfItDehu0mMHqAUlWJw72AEjIARCAiU9gulA9pPKufU4NlVGdVAoYcMRm2Dble+0vhSDEoG3bXLXyp7KbZ76cNYxKbyyujRPfqULv/F9OxRwrAAvzooHzpxYHAo/iCuzlR4wWCqLKPqW+B7T5YmWwGbkuX3SGKudhLLu5QzOsDnE+pPeFQGlSKxigkR+N2df41ACwJSoLaZYUtKRw1BYEU4Zy8xlAxoYYDIGTy7YIz91e2uBCPGZ2NAmTmD7obkL5J9IOZvlJ9PA/AXQYxLt4XlwUemFcdHLn/QMXZfFHVqoBiTZH8uqhiYpSHKNEc7KeVty+n/FPNMnjiqlb2rIA1uRMKoVveOpH/PCQF1YLjNUSCH6RE4+kecUxfPYBXK+DxRWQyeD6GtsjoHz66yNaheh2d4USYJE2Owavknlr2rvngTkr1O9DO8lfdtV0I9w+uEd3O0+g80W4tUOfR7i4yRKpslRjbMR51v5bEtMskzGk5t5VxgHEvShPu70xdfRA9V9SaHIpkZXGSQwsY1bN7giJsOf0qU8SJxSYUWRgx+XwmTQfsFUppzX6+tnsUPnTQblb9t6prumYXjUn6o6yH7PvrCPHUHzOB5W/IxeKJXnYPnCQF4e26qUIRBqM/cQXft8hfJPkYFlOq50lf7o8Yo+xQNlZVbr6dI9Xn+i8q/7pMxyTNlO0mKuZhLPFSEuMJXG1RxJhpdg7tkF/KrThD5WQet17AVx4yAGWQfF+vZISc86FyfC6O+g97imKylngOW6BaeHyYzsf0dYCS8Wfr4oKPoS+UHhPpF3A7Zsjpy8fhvv2J2uZQfw6oOkr3a9VlHdF9MOfAXYSAWq1mrZOnmtuNJU36SZWIwlfzZsot3eHivo4QX/mC2NlJEo7f+9MEbfLcU2mTM1I8oZkndxDw+dyMQ7aUa1+ihihExQTeJM3sSOgI8AE2vC4Md3gOHHQJsAm3bHDopPqofrH/2SRzsoygpeE31LFmuxXtlqIuvnCVUcAf/0TbWZmJHGyDE/mF31/Fb2Ll3UFld9CVjkC170OlBk60z1Z/VKbQZGgcBdF79d0VM52ov3xUXkXxoFPH2Us5PJHfTmEL2Rzo+XgoIGXL+KJwGGTUZZbQlYTDPGtDbMidxm63ngPuPaqtj4JBAsu3LBI85N02vBrRLl381FbFyRqwnk1ZQdEJVdtSViopemPhg0tK3QFwKiFcEXB5vgd+peQx4nJ1+bKyewf/Hqeu6Qf863Meln8bjxW8jX5HPKRiKtGNZU5TRl2bkKfLYl05Xvkg3ltOVzvHrRiDWX6zPdXO7Le4ippVBxZJf9FDFB73F0QDFsgQhbux+rdn1gZdH6W4pDXuU2Py99/Gy8OyN8k2+dykMqGz8BQOWVcAg4rG3SXgtPIu/XmEg/9RF1gdfVU6WDvQSomcm8YSBnFXPPYuYIxv48zbunF7C7CWfOQBoKYN+hBD53N2N+xtpx7LGpT6MWuQp8jiM2mHuSDeWc5jCMVtAINZfrM8t8LwVHpno3tVRvch2pYu487832BqweL0btzuvdfJmHEtobPB+zzOdeS2aNNUr0rplczNpfuNaRxqYhcd0afwU17yWy3dwAIU3reAHrxQKGA0rXVZhLTxHfkrPQ/hn4/TRZZVQvyU6UMr/kPQl9TyknCnzgm1TJ6csjw3R1yqAI771Oml5PYhPPvNeOQaTyr9y2Xuoy8VmmVRPLhbVneDRbqr65rE8VLzZ8rMaYO3FoDFqkOW1QgwrNtNiJP3MwKtz/DwDM+7mUhJxNR1dzxLEbz3zF494WXgLq9pkuVaec4EZgX8MzKg4XcVm60AXgTnij9XzHOUPKAP8qYe5A+2zyJCTvuF9pu1gBGKQEWjXL4T/x+punJ/IV3x9eRyqh1TWisEc8hfLfgjf+mNm1NklwJhDT5aQaw1lxv6t6psxqGInfWrAbGVeishSGTRetCSgMeLdIPCdmWpXvM7xmzp4o5r5WJ5pxpF/1CC+4ZmlybYA37jxYvi8Bp4jMz3OQ/lnhnPdVW6pDrTRCTSo+2aoZld6jlHeDHgXjy4NF9Zzk351H2i8101sK63pGpF7r4Q3nvW5XWpgw1jps3crenxZOud/r6BD3JjeLow2/ketUzf1fIywVgzmkL+v7GPgPjeNOXR2bpkobw49WUKuNZVZjVMYVNXFAM7obLu+kUNHx1IfHp/KSIqdn+IYPBmgUs/Q3RA3h4cKQ6+rHPioO+kV8Sy2ysMM/BfpQJsE4rHNYGKpGM8mr6S2vYnZRqoZl13PzYzxPuBHp7RkoJ3WOjkjI3gen6gebkU9Ola20jEb/oU8Ot4pT/RIVX1BmjekZUKG8dlHNvqQSD8lPfb1JBhE3RazX+kAN/acdvVJbTLNIX+R7G1Mrj0uV2eVjnEhOgim8LhOBdUcejIV72unG50yt2D0KuG2uEOTgkGE41inxvM/WjpMPAt4GNJyqfhrxR2jpySjhGpfTQelrmXHpXnuYDc7ui//eEeox4MwUAcO6E0Q0aeeJ2BjMEnw/zyYSjkBjCMCbfNkUNuNHiPaUFzaJ9/eJCUMThgRHMUh6B152fc4dZgEAzHNyzevdLAt4rGObFlmlL9I9qkrYgr6uTqrsnEM/Kn0OAHe6mBbyKrDjHqyahwmZC61XyqDqnWgzGFAihWJsVeiK2AgvWt5SGfIQJ2GLkMmTTPWdesAETp6ZiF0cs2wNM9Nfkrv+/LPQE7egzBQBw7oTRDRp54nYGMwSTxUzfYymOgpAqF+8ZrQNktC0zP4SJnj4Fz9ebBol3hjmmXHeh1Co0mz9X4qDFRY6vWkfcX+tJWPRuQs8g+QvcHuJm6P6qwkqJbxg5FCe5h8a8oIqM2iJyPwuXUS9M97HqrouioVDGv9oLOV0uHyJ1QGWzBUdjG7373BQc+peI7sWVpKrOQ6NIi/Eh6r7IEH1tFpOHv8Bfp7cXPyXCLfkbR9+cdjeP8I3b46cITk8EcD6nl44eNTYPCdvG10sE2byN5HJdwxDmj3lYEQ7vFQtU1SFN0r0Ofg4d7T6V6U8jKNjkGDd7wfeJBzw5zyF8meK8Ca0uXoLMZlSBfbwly6NwSqOfVkCJ9bzRtXDXZ2jqT4Rwc3T6Uw9X6mEumkZHQGBIyyqhONtPQsft/pb8W9qlLpR/GUicsUpSQfm1Wrb1Mp3aSKqrIfhnLp+KPnhX0Mezzqfi8syfMeIz1v+vKvfBi67O/o3FCsNMU6kCOG6FJXvfZQhbzoUlE95/A1NI14Y5CiDTAr5gyfGK4sj9ftRPdVUHq8wCxfTto2QnF7J5UNf/QTlH/SI6T0tGP2xLEsQkBvWv9oPMi19823KseJH+WDn8fip837fSJ3+eOJMaCPpH1R/1lhTvlLZc8SYGWJJONJnVWa6u9FYD3B5EvV2/XKxKnZEZ+ztpO64Au5EL6MT/Tl9f9u8oeUHPw1B5GLHCqfAfnDHOVT1hjlzMnzGPw2aZTwr7QM6HebNKa+V5ko7LM+5SjfKPXcp+wx80gODEI+gLtI26RcsNSR1T6VDgMhC3ulQ6/Y8J4tm9I/0/FPSZ4x0qrM0TEINCv5uc7hU+lml79E9hwZ1pZG8h3VWT1nUiG2d3qqa9okDaNId2P+Oc7ibXY9mUOuNZUhjHE2VDYUfF3pJobouor3k51l1b3W8SEWoGuUlRnCatekt8hzxJfzCPwzmLS+iZeWM8E1s7/VzgAnkLeNJEtl4L9YUGcBD7yxy9LdqUAn88exRNDRQYfP3oOXusZwzg3PlfBxbuKx0k2AAX1gNehJfjrlXAxml79Q9rEgn5POUZ2V/PRBT1VPvPFKndEe8diuuW+aXU/mrLCVlFXbTdKLWzfF1H/r+C8d/zsjgygvM4IYWPrjw6CTu+8RWmWxxFgaFuO5lNGO9IP4V93wUdbfdNRu745yRo1WuSeXmNoKHFDPbeQWiwNvFQ7mSxizTbkxqmir6WbqvTTil4EGnh/ouutzKkz1Wd7iOFjeVFxnEE0GM94knLyv6GBiTAw6cewom4nRkvKflL2L7zXHF+hsry0xS8i+sJ4sIfJSZWJQ/Q+Fq0+6vsGZGwU2Ys/SSamymYnR6RLYX/Gryu41cFYUCn5UNoYFH7qkM88OS/KczeSRhGPwLxoYo9QVGx1XHfrW89qEkhxsRO/7nabRxRE/eFXeSgeKDKExGFHZ9BmUz56r2feSRRmWwmAN8i8le8Te59MIrEFPTnN5HimENR57+iTCl4sYVLuy/btFBEJj/V4D2mZma1vEGZ6FNcvgeHkWMx6a2IknjGo6kFk3yCflsrl9lolfU/Z4n/AyGwZJmYvKn/Axm+wRd59PI5DUz6J6cprT80ghvG1QnUdVWgojsAwC6kTwFPFmy3cybmbZQ6Iy8dThHV2FIT83BmuSf27Zl9HybZa6Jj3ZJoJlXAtvG1RlkDm1ETACTQTCoMobTkVL5006OfcqC6/YPZU1y7aAHJ5IMxcGa5R/Ltlz68Lp6k85rK6dnHPdqB3YoDrnCrZsRsAIGAEjYASMwPQINA2qm9MX6RL6IhBmpu+Vnxl6bqj+HiE3sdMZASNgBIyAETACwxEYZFBpwI9vCA7npIOC3PxsnB8tTM3zmPyKFvtTil+tboI1tczN8nw/DwJj6lrk2LoSkfDZCBiBS0RgSL86yKAaUvBSFbVFnodidYkyD8XsUvNbVy615i23ETACQxG4GkpgzvyaPZcsfc3JWmdZW+S5Uxg/WDUC1rVVV4+ZMwJG4MwR2IxBFQaLf3SOH9FafdVskefVg2oGWxGwrrXC4kgjYASMwGwIbMag0lIE+4n4wOGtMHjMBlLfgrbIc19ZnW9ZBKxry+Lv0o2AETACg/ZQzQ2fBo1vZEzl/oHo3Oy1lrdFnlsFceTqEbCurb6KzKARMAJnjEDqocIDtIVwP8zGt8Br5HGLPEfefd4WAta1bdWXuTUCRuBMEEgNqtVv+A5LfX9vCfst8rwlfM3rfxCwrv0HC18ZASNgBOZEAEcPBtVWPFNg86OYnv1f7gdWyhZ5Hiiysy+EgHVtIeBdrBEwAheJwO1UagyqzyFi70GaaC3XMqZW8ceoJXhskecS+Zx2PQhY19ZTF+bECBiBy0MgXfK7POktsREwAkbACBgBI2AERkCAt/xGWfLT/o3vReuBDv59+akO9mQRR7iv47Fm0KOUVVEc4WeLPI8g9mZJuL42W3Vm3Ai0IhD2/T3Rw290fNIYUW/pCM/eKO6H1syONALLI4CdQ6hsm3TJLz7YPS7/fSDF/ykQ/k3nOzSO0ED4ftSv5SQnz7FFnicHZcUFuL5WXDlmzQj0QOB5GCMYM5438v+o+019JqfBv2/PH4G9rVLpd6i+6iu7ZhJ3lPePkJ/rj2ok7xJ6n3T9LN6HmQeN5anSDf7z30i35FzCc+A3NnbkI6zO47Zj6zx/C+uLOqIjZtbARAEP6Qvp2kedHYyAEVgBAo02zeoGE+80EPd7GuFrI7AyBKIjqtqLfiXmGHQI8cHuruz3c2JAseT3tpEdd+4XGCY6eI6Ll/KGlKnsg0I2zyrlpeT7KRy4nwFvM3+BMwil9WQuqi+x/ZXq62cdLCFg7L9fjyjmxAgYASGQtmkm2M1xgy0jeK4cjMBaEYiOqMqOwqAa/JafBq2KmIyluGeqOavAiLomnY6PYZBrzkZmBayEZzH2JJENPl/quBOMQ+4dJkagsL5eiJ20c0bpF9W3ieExeSOwOQQabZrJdf0Wd+hbiWuOJZuT0wyfNQLoKKH2ULEcR4gPdnf9fiu3bWwokFDDgC6G1i/crzDk8Mwm+z9XyPslsnSyvqR/GO0fE3BY/nud3PvSCBiB9SCA1582e52wxJhRTcCTOF8agbUhEPdQVRN2PFRx5h4fDGGYRtCcUbC8R0Nhw/oaw0me1dBZOkobOwbWX41Be42ynSNPJ+srCi1j/pkOlmZfU4cx3mcjYARWhQB7HuM4FBnz/qmIhM9rRgCHEaFyTKUGFUrdOwRPFEt7NR3Fcc9m7u8aBknvcsbM2IfnIBMej0U2048p/9ZoldaXdI79U9/peKC89UsRW5Pb/BqBM0dgz5hSW2XSxOH9U2de8WcgXnREVTp8UwJFZY6WVl8ZaQCEl8ngxWb0bzWwxTKqBCv6KeJZcmEssn8KmVKP1YpEOmtWiuoLJKgn1Rve0Q86/677dCnwrMGycEZgIwjQPt+ofdK38l+t1UtMOjdXOxTlYARWhUB0IO0MqjDgVBxiMOi+r/ET97bQCLbSELJ5BhvJxdt+5GFvGN431vj74gUZhzIETtaX6oWJwf/pwCsajado/N5TfIwrK9mpjYARmAQBtVPaZ/3xTrVhDCv2VLlvnQRxEx0DgTDWMN4QKl292l3XXqpobYXoohPeg60YUlGwLJ4FHLiwqfmFru9y6Jp9VJ91OMyHwMn6Cp0z9ZJ2xuSj017rixHzIeiSjMCKEFBf+lpH/QkaXTNAse+WN3UdjMCaEYj2UvUFAxi9Gbhl1s7DmCBE551CIyDvyTVvpSUdxgiDHJ8ewFDZ+8sB3U8eSngWMzR4Gnrd8GFQgzdyOMyAQGF9Mdt9rjwsHxDu6/AybQWFf4zAqhBgHEjfwH2je14CSj8MvSqGzYwRCAiw4kGoJ+/RoOLDh702WmvQIh8bzwkMYrfVGDrfqNIzCl/0jb8ePH9ZSeefRRDoUV9MEDgcjIARWDcCjAVMrHlphL1TvJG7tZWOdSNs7qZC4GCv3w0pb9wPhPeFdWu/vTYV/KZrBIyAETACRsAIbB4BTQJYkcPD+oPspsqjWhlUSKaH/3LWgxucHYyAETACRsAIGAEjYAQOEZDNxLen2ML0pewm9uh+ccVPCNU6oBKRwMEIGAEjYASMgBEwAkaggYDsJPZUYyvVG9JJkhpUcRMgLiwHI2AEjIARMAJGwAgYgUMEop20t98vNajiG3rVd5YO8zvGCBgBI2AEjIARMAIXj8CjgMDbFIl6DxWRcmP9w1nrgX6rDSAcjIARMAJGwAgYASOQIBBspVvNPeeph4rkfPjwlhLz4UoHI2AEjIARMAJGwAgYgYBAsI/YQxW3SdXYNA2qX8OT6M6qE/rCCBgBI2AEjIARMAIXjkDcP7W33Acme0t+VYSX/S5cVyy+ETACRsAIGAEj0IaAPFR8s5P/PT7YGtX0UJGf/1Bi2S9aYcQ5GAEjYASMgBEwAkbgYhEIy31siWr9N5gDDxVIKRMf+XwnC6z+B3DiHYyAETACRsAIGAEjcIkIyDZiWxR/t1d/zDPFoc1DxfNXOh4qMxuvHIyAETACRsAIGAEjcLEIBHsIYwpnU/Vl9CYYXQYVy36E+KfHuzv/GgEjYASMgBEwAkbg8hB4EkTmD71bQ+uSHylljfHv3y91tLq2SONgBIyAETACRsAIGIFzR0A2Ed/p/EXeqaddsnZ5qPi4J8t+uLXspepCz/FGwAgYASNgBIzAWSMQHExsger0TgFAp4eqenjjBuuFbML6RgZW9efJxDsYASNgBIyAETACRuDcEQh7p/BOPZUd1Pp2X8TgqEFFIhHDoOIT6/6Pv4iaz0bACBgBI2AEjMDZI1BiA3Uu+UWUZEjx6YQ7Iho3ZMVHPhsBI2AEjIARMAJG4CwRkN3D9zg5sj4hddJDBUoiyoes3uv4WgZW6+uCpHMwAkbACBgBI2AEjMA5ICDb55PkYKnv9xx5sgwqCAVL7S8R9l6qHGSdxggYASNgBIyAEdgkArJ52InI1BkAAAANSURBVIR+L9eYQsj/B/TkeqRXkTA3AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\left( \\frac{\\alpha y}{p_{1}}, \\  \\frac{\\beta y}{p_{2}}, \\  - \\frac{y \\left(\\alpha + \\beta - 1\\right)}{p_{3}}, \\  \\frac{\\left(\\frac{\\alpha y}{p_{1}}\\right)^{\\alpha} \\left(\\frac{\\beta y}{p_{2}}\\right)^{\\beta} \\left(- \\frac{y \\left(\\alpha + \\beta - 1\\right)}{p_{3}}\\right)^{- \\alpha - \\beta + 1}}{y}\\right)$"
      ],
      "text/plain": [
       "⎛                                α      β                  -α - β + 1⎞\n",
       "⎜                           ⎛α⋅y⎞  ⎛β⋅y⎞  ⎛-y⋅(α + β - 1) ⎞          ⎟\n",
       "⎜                           ⎜───⎟ ⋅⎜───⎟ ⋅⎜───────────────⎟          ⎟\n",
       "⎜α⋅y  β⋅y  -y⋅(α + β - 1)   ⎝ p₁⎠  ⎝ p₂⎠  ⎝       p₃      ⎠          ⎟\n",
       "⎜───, ───, ───────────────, ─────────────────────────────────────────⎟\n",
       "⎝ p₁   p₂         p₃                            y                    ⎠"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABkAAAApCAYAAADAk4LOAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACoUlEQVRYCe2X4VFbMQyA83IMkBnoBvQ6QekGMELTDeh1gl66AWWDwgawAYUN2g2gbJB+n2P5/Nz3kvQevzh0p1iWZMmS7FivW6/Xsxq6rlvl+QPjG/AcnftaJ2h0F9BLUL1f6H1rZBfwTmc6EYFDFcHjiqeRP8oqnZNKvsr8E/WCn3k6h1zP5hABdxBXMG+CAf0E/RP8nHkaS3KicFO3mf+B8XemY5C3saUn4ByUWDivMcvusl7aeaaLLjpGe9asK7xIlYzrWilo+NZIuWMxXMmP4fc2yPwo847Um+fimfvB4sIXlN+ywPS1cArjvpHp+AlesjmvhBZ9DFxwNSK0NuP1QBiF/w5toXpAlJ4QwUhmzE1DCz0H6BiFeB2KHTtMNML6fqS0INO5xj0Y8h7glbuQZW7gAtRZ3K10d9DdbEAnz4k4crPpNIbdSQ4wZoTFILRReRLLhdXRVCcelnI/oC/BcpciklIThP8N1Mp/AE+X4P/XJYbLP0bi8jPJSRjZNcYR3qU3Sf5ynHTkYXMbJyVk++LXwm/PTyN9OafrNZKmtNunLyddB22guRGwqbBp+AT62tkYCO/Aj7wZvvf7Q7xeMbIy+ltfOV++8pQqAwebwFg/NPZqQhSHGIn+Vtqmre63UkMOPwH6C3AJ2kePQpuux8qo6frarPSJtUUyhW4i0uh8HIbCQ9vFvf5WPcCGrf1E8J3322S0KemlSysZ0qcAC0uB8+51/iOU9h3HnGis7TrsCnUa3yr7+pgdtJp5x9bjMWTwnH8B39fRhXzX+I8TFhiFsML42YZMPdVbHPSa6yzbOQw5iXqYrjZlOw0OKQzVxEiexXhxWB89mJ53z2q55bW8ptHxnqQOPq+x+S59cU83Jih43r25OnFchmzq+BfoVcsqlbahOQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\frac{\\alpha y}{p_{1}}$"
      ],
      "text/plain": [
       "α⋅y\n",
       "───\n",
       " p₁"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q1star = foc_sol[0][0]\n",
    "q1star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABgAAAAwCAYAAAALiLqjAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADLUlEQVRYCe2Y7VEbQQyGfR4KyEw6MB2EpIJAB7gEkg6SSQUZpwNDBQnpgHQQoAPogI8OnPcRK1l7vvPZHvLPmhGr1ddqJZ1PR7NYLEb/Ew76nDdNM5Hss/BBeFj0viqg50JvtnCDNsoS558yn73wMvM2occyqkCRvxHjVMbnlWA0etT+XYs3vG1HIYsvbR57wSXYJVvH660BHh10q2PRRH/kvE3XpkRn+sURqaCQ87JSD+CjdKPAJZXUhQa4k+wHSkCRXYg3rQos2YwDQQHGIJHfsbqsyE1X/FPhU0uGnViLUW+KJIwiK6KZDC6EliLtudVfIXAivDdq+QfeH9v6ydrQPX0FpgYWEfroJruntp32wYsUiWmFdMO84gCjzIMuNhD5QFIKz1Kan4MTGd1K0AXLK9fSqba3sovia0+gz+4r1wDBCijfRPRe2NWi1KI//xLaDUpb3Wul+gHacygP2FQRtR2hV/GKPjZXCAF7DiSg1VAmIhB4K3yQ4+hv46Y/sqMx6C5s/UeRIA8jIBHaL/uf/a6IH+FNtrcUccK2oOjnwhu3K7ch+u/OYx0XAdfbFsj1z2REqs4V/e/EGzXaoPgoQV+LZv2gFRh183rxe8Qv7cvTG1qlyGn/6uTONdg0kv0Bg5mii+xJG9TcUaF6Ze7oY63Zvshr04Nwn6J9igYzMKhwkDX0EuHlwwzEqMLXDS91eMAH4ZleKnkGMsHaPzKIl7wUfaBlVGHg5UPE5MiEV77fdI0bKPqJHPhAC83Elt+vHMgI6eP5N2gBukD37TwSKeT5kl/YiB4dATeAIG1ztysyviX4RohsOL3KeMm5OXKl4oRpzQZgrciPXS6aW8CrviGQj8Vsg837EkYxy2hDsX8VZRrgum3YufcofJUSw1Q7BTa+ix9pdH1WAekbThEOhFhEt4imZfmgWLl+cY6cBug8PLpICoD3/ExpsY4Rj6HqSM6qSRrl0nlEjzxSiiyAKBzF7O0G1/FVuhQ2vptFc5OJy30N5zAEXLXKvyvmFUdCugqnjgS3kqY4AKGQU6r+z46dlg41QbdCl+fVDsCpkO7BgLX6R0g22Jb+B5+fuk2l/dQLAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle \\frac{\\beta y}{p_{2}}$"
      ],
      "text/plain": [
       "β⋅y\n",
       "───\n",
       " p₂"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q2star = foc_sol[0][1]\n",
    "q2star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIgAAAAwCAYAAADKIzJKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGkUlEQVR4Ae2cjXXUOBDHb+9RAHd0kOsgcBUcdAC5CoAOwksFvFwHCRXw0UGSCjjoILkKCOkg9/9pR0a2ZGezsWPJaN6T19bnaObvmbFW9ur6+vqXHGm1Wu2It4sceVsSTzfJ+dccJyum98XXbo68LZCnHZN3cmrZAUTMPhenj2Q9PiU5LiCTu1LpEMErHVl6OBfrGn9X6Vwp4kFyPoUvlSH3iFY5uRibwJl4ehxxWkiG5rAjVp9qDseeZeW90vkz5b3weVP/mizfaZxLpSdKWOTfxMOVfiNS/S/K/KtbnpsFORSTRxH3E2dIOE9NiXcayZTyPASHdYiS7tVlomilF0qvNfb7DSaG3JF/i3IDyF5CuC2GJ7rA9Ebmd4uxXon/fxLt/lbe10R+Nlkm9z0DecNXNgDhLhZXi3tqsXlhPV42Us/3BPnvhew9CC/uem7ow9/+oXQe3k1W9k55fX6YfBcwDfGhfrwZ/GbjHKnP2e9OA8KpfolBMNf4es6hx+Ix6fvXxdkckf8zpSZ+GtuCHBgoTjTIQWfaIDMZKVs9AqnzTpvmEsErUX6iMd7YOG90fWZKIRKnztAYTX8TnBCEflXiLvyohAywGrguDxSdZk3It8XraAAxJX226YPCrrsgb8hCIEiCuT4iyv4kBTR96Jy78l8lgAIBjqbc5cxwEF/HlrBsWDyeJkog5I8eGhrTxVxKKH7tAmvxthllfUKM0c0Lq/yuCxQekcCHyYbxVHuAiPWBWD9J9kGh9QMfXWJsyon4u4RV6HOLrq7awRsuL0XwN/gEY+3PVI9+NiWeUMZ2rfDasiCjAcQrRpNFAUy08WPKQ0DkbXt3A7hTP4bOQ7rShVuY0m8KQE1dtU8BwC8SsbSfegJp2g+cANC+uTF3eOwlm1cOaz/cKC1eR3Mxwey527jrwoEADc/lQ4gHvQ+Dftyp3V3kD7Wl/HNnzG5XU167+KNngJtca0+zWbKRI26modEsSNPj2kSh7JA2ERKMtcwbHaB0gYTT3gBWZYDPuzfq3jdxA0RklhPrkoN1iPhLZGBBWrqbwoK0BpCQEB6JqH6IsBB/9lTAXQGyFqlvHqkhZ3lMIeucezpqTMa+CHhxI+uaOfM0Q6zQkomrkOcBILf0NPp/MSYwonaE4tcq3NrIkKBMoKxpsIYSkcrD9Y8rKqiui3NU5tcdvinv1nGE2vP0s1UMYm2ZK9bPW8BHOt+KF7UbjcQbAAXAWDEHZP1yI+KOIzmpPlYad/kD0LpA0JMlDYhiv2wyhurB4O4mdcesozEByP42fard4TbtcmujeQBuFjdbWBjVxQiB/LXNeoUjnYNarMfg08W6tjsCpuSTRlBnilMsEulnJtaSvJVu5PCgORvnBL+LufeEq2HRaKMAUvWOBaoTpXvdTaZx+x5R/TySv3YD9K1/JNvkmIm8xRcyj27OUWMQDeR8uQmBWOLjbYVvQqddFJTmJlzxyg3BAiF+vVjSPAhMCaYjKzoqQMaSkBgG0a1NN2P1XftpS0CyJgRgEfJHYBpUyRIgAX/1dGYJjBqkzjyXOvwEEqgAmUCoS+qSNWweRe/8T6J8WZ4v2CxJWzPMpcYgMwi9pCGriylJWzPwWgEyg9BLGrICpCRtzcBrBcgMQi9pyAqQkrQ1A68VIDMIvaQhK0BK0tYMvFaAzCD0koasAClJWzPwOvaGoRmmMO6Q+suAPR7sReF9FjbQ8FcEeRCbql+m9k240gUeqgWJlcqmXbbfsXmGjTRuM7Py2OTLngk2Av80VC1IoGrbqOTfL95RES+Ahdsl2VS9TxOryw46gISVwbq8Vf2id5dpDi2qFqQlDrd90AMCF/O+Xew+awE4AAQbfHkXmD23WBeAxb/ii6JqQQJ1StFYAwDgY47uZmZAw1t8vO3X3anPuzDJbXv0WSpVgKQ1R5B64QFDFbMaAMe9rKWyrivB3USvDdC2ZKouJq09gNC1HmzuxcL4b5G4lgIOn7rkXSDeCmy+aOAKF3CoG4Y6SjRL8V3Z7PR2r14oD9dCfMFnIruWw1sXnm74+lH0SqPyi6UKkI7qBAZcBcr2ayHU4B0fXrHsjTEMRFgSvkcWgYhOSqQag8Ra8/EHLqbrZlxtszL/6SK0KC7AVd4TpcUApMYgTuWtQyr+aFWQhQAMl0qhRaEd+R+UFkPVggSqNMvAAhkrqDcRX1I6UBv/bi4LZaV87vKmuTXlNQYxUUjRxB58upOAFBexyKcSzetW9D9iv54H5qygFwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle - \\frac{y \\left(\\alpha + \\beta - 1\\right)}{p_{3}}$"
      ],
      "text/plain": [
       "-y⋅(α + β - 1) \n",
       "───────────────\n",
       "       p₃      "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q3star = foc_sol[0][2]\n",
    "q3star"
   ]
  },
  {
   "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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAA3CAYAAADgzZrvAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAR5ElEQVR4Ae2dgZXVthKG2T0UQJIKHumAkApCOgC2gkAHUEJOXgdABQQ6ACogoQPoALId5P2fV/KTfW1fSZZk617pHK9tWRrN/B6NRiPt9cW///57q6WGwBiBi4uLu8p74OTfka7817lvlw2BhkBiBFL3O9G7IxYf63iq/vsT7F40ow8MLbkISFHu6f6+dOOlzVfeP7r+j/KubV47NwQaAukQCO13Kv9C/fHpHAeGnnXcMPo/UvZyrkLLP2sErkYGH8X5qxn8s9aJJnx+BEL7HbPx2aT++kkHs/MvbqHb7s2pX2vk+0MyftbBlOcHAfIcmZVP6KJ5sDdYoEgfhckznX/Q8VDHtfDppoa6bqkhcJYI5LQfol2s352Npy9QMfavZbxe6mD0u2sMGwr8hD8tdQjg1b/X8VYHmHXTR6OUum2pIXB+CBSwH8X6nZfRl8BvdMDUJklt39HxjnMMA6r3QvWY6nxy6r/T9ZWeIZeb7xQJvxS9zbACHx3ROBlpfxROePZfdDBAMgC81jEbOwxHqZtdbYaTL7/CEu+rpcwI+OKcSL+DpVG7JexHkX6H8LePISCB36jMO9P5jxXP8lxtX4sPQjEfdP6F+8CGWL3+bVSHOBed+lfR68I8o+fBt1tjlQCnOZkJ80QNuFMEt8ZpiqdxnngkvIWODOKh43LtPgkCzLofSn8Xd4el0G+1wyYFwrzHEguf9t17249I+nO89P1OdBl4xk7IfeXjvLrpjfjuN2C4D+z1otEXQQx+5/HZCludJcgn8fO72v9bR7cKHcALBmvszX9THvnQXJ1isFIdXiIe9FcdVqbnKHcsQ2twEj/gAS/jxGwo5cC4C50aC2nvhQPrGKz5LBohW36P5xy6tUZO8YOxxZ78NNZv3b/nObjrmrDibNLzNXbgFvVF/NfZBqYfeNuPGPqS+2i/E92Dmbbq4YyHynILECYPyU6c+/Pc863yxRMjHqPZJN9T+SrPQPHAPtM1IDPaK6v7PwUWcr3pjcuKTjBWqoPBf+LSMnSCZHPru9eiFYMTxp1tYD0WusfjHeS5z0OuRScYpxD6KcqKR3Tj7xS0tqIh/rPqlq9cBksMPbpIH6SzzfY1U2b2uduuygbrt1s/5NrwldN+RPU78YXR7/vq3LXK4cT0tnzS09cIgtIA6u52bEiwp+Lvsw4M5uI0Rvzb9EgXz41cePh3VZd7YuB2Gr/oYVhC47OhGYQV7YrO1HQW3vCIVqdInGj7hfizC9vwiVe+2uONwWk1CHEEcAZ4n0WT8KHjo5e+Oj3JXwndmmx4IlOyXCubvscOOfrZMd0Gd/A/8GqVN0iR+j2gEXCTzX4YHrL0O9PnwLLTLd2D7+fJUaJ7IPAF7OTzrfPFHyMX/yzk5RXk5DcGK9V5NsWT8vGKknj60FcKwknlJ/ma4jU0T7TZCbRbnbLyiMd/7HXJs9rlXa3Gf46G8pPqVig28KVj0dOHppJ3v1bZIP0O5blUeckR9d5Vz8vTH8txqYqDpNEAL+97FUwSwx0QT3Qj3vDKWWTxWZBJ1OohmZRYiRajMSP+eMH5sGHPnFPEyVP0qGLmHdjFuygae6yUQ7cyygn+j33o70W/fXjNVCZKV6fCOxjS4CmmFMsaYLsoSRyYRZNBUjlCBgwsLFwSZ+pDB+bZK+V1U8JBxcMbBiW2J+I9Rgl/SDI4Jxgr8YtxZ+HKhtCYAnNNGixyGTzWYuWFk9piwDl4Xx1X6/8E4xTb5ErM0Du2qC4mteGl64tEMjwUX966laH5VCTBn8VJXxvkpd+pmEtNZ02/k907Ggab5FcV+xCOCnhNwUZ1MFhM3ccLHUzTiFHq1Bk1Ythcd1N8nbup2YgWBk5Z/+dp6VplWRxKssi41M7UM7UbjBV0lPoQh66RlwODC4b33LZ0nwQr0akOJxeHkOs1mBmcBovrbtt6HqTrbt1j16K9OrxjZactpUXdOsZP6ufix6u/GL6DFtJVZzP9To1TCXqXAsxNjBxv1TDep28CcOr0HpKp/5fybYgIhbbe7UdDmNF87KGT19Mx5ZZOv+vhE+PdLZXL8SwGqwEfwol/fuLAw8Z7fGULSCYMTCqsqsbJYnLsnAAzZqEsps8lb12fI1Aqf0m3SvEQ2Q748x5C0pb6HcLnLsr24R0zzcDQ+IRWOuZVh9VgXhCgjxMG/b7JZM8z/2CFG293yRC3G9djejrOMyQmT3aAwKvpw0STJRNmxmBF86oHVl9nWAEvPH6bviXEqiqcLAAR57WYfa82r6faDdX1KRrkGTro+TjRNs+fjh/onv8mX+yXgbo10USvnx/0MMToPhJvKcOC9APsUEjaRL9DGNxVWb2wLrwjpljdD9q5QHkdkyvIysdz5TnnwS4b3aP0zEH7fF1j8MgbhDgsf3Nnlec/0vo9qHPlUuarvWCsaF8JuSflUz7T3wP8TZ3VWIlONTitfVexmIE/dafaN8+CdX2K1lSe6K8K78C3jiDdmuIjV5548w3vIMdBPzjGl+oU1+9jPO31uRveAWw7YupyORnPAo9gaZTn+UcJP/ae8FrwXtx82r9W3hI9FTlIGGD+jTvUOzggFJARhJVDl598mJNvLrSVCquacHIgi7qMxQwvE50dpJW6PqCV8SZGtzKyE00a/L9F1N5CvyPY3L5KZ/Sl1HjZgM1o6ZUcg80C5FzCiNtwjlsGA00Hc9Oc0XPLTF3/aTIxxNlTDFYOU5M8GpqEwuwaiFOlm+qmwKomnFz5Y65j9Qtjc+A8rNT1GP5j6sToVkw7uesQ5hrru0+bRfXbh6G9lumMvpi7Mgx6e/qmPNuqMNaDJCNGjJ3UeU3GqN3k3PwdvFQ9R2E5vAcdS8x0yGvdH/BhyyQ+R2ElGcHii4NNx5aRHS+F2OgAF8P3IC8Wq1pwMjKvPcVixizs55nGY3V9hly67BW6lY6JdJT4FYAa7EA6iQtT6j6XKKUB5PsyDN+Ftq+67p5ljC9rBN0eWz1joZe8r8rrF1qNkr5SPp2ThU327DNQ8POigw6rvKPJ8M/2UPujZUfrxBYwbQVjpXrEbJENT9J6k/yK3gAb3Q9SSqxqwGkgfORNLGaqh+PBFuBJPdLzYF33FUG00Q90uO8ngXWDdcuX/ppykguHBoeHmWzn+OjM4ErY90BWlSdyQKhq13ZAPFabbhvOeSFssQxOejlTIYmOjp5N7URgUGAgIO7aJdOZiPEHv2hDAiWiw5ZIsVj9LPkIdcGrd0qMVQ04eWMzVzAWM9VjWzE7aFgQPXhPygvW9TkeE+dH6VZiHibJCbO+n08WcDKFe+cMqU4NdsDhvK7LSwHN6GtH4Ozcqz1+0Iv9zl0y7ePlh2zVNLX7U7euIFpZDb/htWasThKnXgt0oXe0Vr/w5iedFbedDNc4QhznnBhU7WwqBoci+h3D2J7qENO3oYYOsALMYZhfO+0Q5uEflKYWfJ1ii5d2lsKCdM4UhZUZLOb25y/xmxqrXeO0BETAs1WYSQ8JTZbeDcbs971pO0DUfm99jG4FtZO7sPoIfWvtr4yW0u/ccGSlfxugTQuxU6pQBhnN6VTs2yV2Sgw1dAF53KblfTIWOy684j4WK0JCMTKmxmrvOK14NX3VFJgRkiAWXWpzQM98xEWsbkU0lbUK63/eoaAZTkrp90zzdWRj9NkiRbKA3dxl+isDv8ajn+RKNLv/9tVDK8tkuQSZln4QVrGDWmqs9o5TgveDx7xavwxOfLch5JsNKdgPpiFeY5yJ4HZyVgBn0efzhKvCW+a9wartpznZrpY2Rp8YNWkV4DckNv9rZcnFiKVfO1ZWjobTDAIyIAzs3S60mSItOxECwjo1zrn1O5Hk25C5VLNsGyR9uzlV+xdDnHuEPwWsGk7Vqnhj3AOBEvrtwcZ+i2D0T2VUZNDKLUtu+iU0peFUAuXWxlYIlNDvrWRL0m4f0w+NpykO1/2CWBIuZoiIp4uZR3PZuT39jv7esKodJ+kSg+kHHSGD6sGvO5bQyTnFa/n5ENihfucTtgBljH5UingRUe2cQqWG1fJbFD7XKsG/369KDedV8LXKZ4IA4Z0ulm+8rdrF7mTJKMSpYNVwyqgkjfTmCOTW780FXMMARh8vi5Q7NHLTSr6/8G9lydWKpV8zVg2nXNrR6O4BgRL6vQc5o3nA6Nv/5rsTTWU/FYP2z0ewfSpYNZwiXn6rUg0CufW7GiCmGHU9/annNeUxaOWe1llPvyZcxrw2nMaItPtTQqCEfleNVx/TlxTVhiyc9YjcRtkOKlVi1XCquq825o8gUFC/j3Cy78cYfTsVurtvVhe5s7zb8Mti4RUPa8eq4bTi5bequ0eglH7vHoglBi+1ze2TKbB6y9xSQ5mf2ZdtjXKW5k4Aq4ZTFs1oRHeCQBH93oms0WxcmpoYS36tr9ZkX7YdwHLKUTNWDaecmtFob41ASf3eWtbo9q3Rx1hawKKJbVjxZ9qWJ57V0zfy1YxVw2lDJW1NZ0egpH5nFyZXA/Y/cvlG7kMWQmQ4gxdDVe+B6vPb4/d08NUhVtDJI/Eifouh29X2+0NbpX5itmasqsHJ77WfZqkd9KdagS2p37VidMt6+vaLM9ZQhwrEh4z5eAUDBkax+8Cz8vjwMd43H6TIktRBmKEwyNBuiVQlVhXiVOJd7rWNzfrTXgE5xtcG+n2Mpd0+74y+jDMhCwz2VSinBuyPph4GmA+cux+y4DOM/WCi8nd0PNHxd2hbM+UtbbfNmaLrs0thZXD6Q2eON+ZgcItN1eAUK+Ap1NN7pg959SfK6nimg/7EGT1htn2Oqah+1wzwbYf5P3X92Ln3vfzmGHkUbvyB8+4ThlJGDBYKbV/OGgPm8kZY6Yt4KBHPt+2WwOoPydR/oFv48Tk5BsrYT0LWhJPF+RzPIf2Jj4ij+8yo+Sg8nyDl10q/4/7M0hb6XSXEnadvOCcEgxdujbKXQFI4ZggonK03jq0zEFxTTgezABvy8aLvUeihymAQS6bsWEkYvDeLKbLRwfHsYj25anBC2HNNIf1JGOFgvXaw4iM/JZ0fp+nNL7fQ782FjmGgN/pSNow1CvMohpDq2JG2GwSgIQOFN4/hwjNOnkSfb2uSUn9u7YbqzN9CWOHl2/WDGU78sivGyU/A0yx1tD9JD3Gi3G3KZ2n4ttLvWtWuN/pGALxJPMyY0AvGfezlY5QZBFjkzZGg+1aK3w80ORqZoZkVK8n0ciQXgwBTebeTz7B2kF0rTgeCnFGGd39SfyWeT+jvBXpzRhhZUbfUb8tDPWcpCfvb+0Oc/6ODeHKfd+xa5Rkk+JLWO1tW14QhoHXP5tmz8vBIPtv7mLPq0ylok51C3rymLKu2s2MFv0pgyYI4W2qDZFWdKnEKlfOUyvOedXj3J6Mj1GEH27NTwuKYLJJ3c/0+xuPenl/AkJvMVIkY+Xd65uVBqw5GnDg3U1Ibc2bBkcHjIMZoyvMsdlGS0BEKzjpBbDhK1delQlix+M374POAXu/DlapWnFwZzu3a9A/v/mTxUT36Hh7/T9KVmBmhJVXNeQ/6XQ1YllEpx4HnqGcojre3r7IYJW/PXWVXefqqb0f3YM93St41eeIlG1aijcF/Y/nTNZ3ae2ajstXiZGU+x7Pe29H+pDJ49oOZtO7RF7y4J+eAm+TcjX7XhPc4pi8cu/Sb/hLbR4l8EuCP4/k+9WLL4AVFeb6xDS7Uy4KVwZ7O/zseHIeuiet/W+Bl/KhmnMaynNP90f4kI8OsD11wZ9LUIz/LxgnR3Vvak37vDZtZfg7CO7akjAx7fq+kXIu/vqly1uPACC/+g5TKMohguFBOjBiLTswQun3Guj6aRIMXzexks7DOmMlMWOHFge0gSe6LQcbMTc04zYh0Ftl6byH9iT50peOrAYefPHkuHXEHAvPotE571O9qEJaCHIR3bJ6EwMDOhnn0jDAN4Q2mlJyzTitFv2vP8rens3jbDVY147Snd1qaF/vedC7Sn0rLl6o9i1MqeudGZ9bTF7Bd0ojKgilx5U23gokPvJpXOn4RL0xhd5f2gFXDaXdq0RhKiEAN+p1Q3Cyk5mL6fWMysOzIeSSwCclsktQ2U95dG3yA2RqrhtMm6tkaLYRALfpdCI7oZv4HM8cOwhKw5/cAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left(\\frac{\\alpha y}{p_{1}}\\right)^{\\alpha} \\left(\\frac{\\beta y}{p_{2}}\\right)^{\\beta} \\left(- \\frac{y \\left(\\alpha + \\beta - 1\\right)}{p_{3}}\\right)^{- \\alpha - \\beta + 1}$"
      ],
      "text/plain": [
       "     α      β                  -α - β + 1\n",
       "⎛α⋅y⎞  ⎛β⋅y⎞  ⎛-y⋅(α + β - 1) ⎞          \n",
       "⎜───⎟ ⋅⎜───⎟ ⋅⎜───────────────⎟          \n",
       "⎝ p₁⎠  ⎝ p₂⎠  ⎝       p₃      ⎠          "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = u.subs([(q1,q1star), (q2,q2star), (q3,q3star)])\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}