sympy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7920edad",
   "metadata": {},
   "source": [
    "\n",
    "<a id='sympy'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d217bfe",
   "metadata": {},
   "source": [
    "# SymPy\n",
    "\n",
    "\n",
    "<a id='index-1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb400117",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [SymPy](#SymPy)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Getting Started](#Getting-Started)  \n",
    "  - [Symbolic algebra](#Symbolic-algebra)  \n",
    "  - [Symbolic Calculus](#Symbolic-Calculus)  \n",
    "  - [Plotting](#Plotting)  \n",
    "  - [Application: Two-person Exchange Economy](#Application:-Two-person-Exchange-Economy)  \n",
    "  - [Exercises](#Exercises)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c5b89f1",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Unlike numerical libraries that deal with values, [SymPy](https://www.sympy.org/en/index.html) focuses on manipulating mathematical symbols and expressions directly.\n",
    "\n",
    "SymPy provides [a wide range of features](https://www.sympy.org/en/features.html) including\n",
    "\n",
    "- symbolic expression  \n",
    "- equation solving  \n",
    "- simplification  \n",
    "- calculus  \n",
    "- matrices  \n",
    "- discrete math, etc.  \n",
    "\n",
    "\n",
    "These functions make SymPy a popular open-source alternative to other proprietary symbolic computational software such as Mathematica.\n",
    "\n",
    "In this lecture, we will explore some of the functionality of SymPy and demonstrate how to use basic SymPy functions to solve economic models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6880a6a7",
   "metadata": {},
   "source": [
    "## Getting Started\n",
    "\n",
    "Let’s first import the library and initialize the printer for symbolic output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "204e618b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.plotting import plot, plot3d_parametric_line, plot3d\n",
    "from sympy.solvers.inequalities import reduce_rational_inequalities\n",
    "from sympy.stats import Poisson, Exponential, Binomial, density, moment, E, cdf\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Enable the mathjax printer\n",
    "init_printing(use_latex='mathjax')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8f50769",
   "metadata": {},
   "source": [
    "## Symbolic algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5f412af",
   "metadata": {},
   "source": [
    "### Symbols\n",
    "\n",
    "First we initialize some symbols to work with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6bb414d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x, y, z = symbols('x y z')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73ca3cab",
   "metadata": {},
   "source": [
    "Symbols are the basic units for symbolic computation in SymPy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1140a2b",
   "metadata": {},
   "source": [
    "### Expressions\n",
    "\n",
    "We can now use symbols `x`, `y`, and `z` to build expressions and equations.\n",
    "\n",
    "Here we build a simple expression first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93046ce6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "expr = (x+y) ** 2\n",
    "expr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69fe51fb",
   "metadata": {},
   "source": [
    "We can expand this expression with the `expand` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2c7fdee",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "expand_expr = expand(expr)\n",
    "expand_expr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0725fb47",
   "metadata": {},
   "source": [
    "and factorize it back to the factored form with the `factor` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afca52c0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "factor(expand_expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47114526",
   "metadata": {},
   "source": [
    "We can solve this expression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64451b72",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve(expr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3fc0061",
   "metadata": {},
   "source": [
    "Note this is equivalent to solving the following equation for `x`\n",
    "\n",
    "$$\n",
    "(x + y)^2 = 0\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">[Solvers](https://docs.sympy.org/latest/modules/solvers/index.html) is an important module with tools to solve different types of equations.\n",
    "\n",
    "There are a variety of solvers available in SymPy depending on the nature of the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bf072d6",
   "metadata": {},
   "source": [
    "### Equations\n",
    "\n",
    "SymPy provides several functions to manipulate equations.\n",
    "\n",
    "Let’s develop an equation with the expression we defined before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdc4df7d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eq = Eq(expr, 0)\n",
    "eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "313de53f",
   "metadata": {},
   "source": [
    "Solving this equation with respect to $ x $ gives the same output as solving the expression directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad12e9a0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve(eq, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f373883",
   "metadata": {},
   "source": [
    "SymPy can handle equations with multiple solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29f0431d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eq = Eq(expr, 1)\n",
    "solve(eq, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e23777bd",
   "metadata": {},
   "source": [
    "`solve` function can also combine multiple equations together and solve a system of equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04c8b9c7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eq2 = Eq(x, y)\n",
    "eq2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32641542",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve([eq, eq2], [x, y])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a087dd5",
   "metadata": {},
   "source": [
    "We can also solve for the value of $ y $ by simply substituting $ x $ with $ y $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1457116c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "expr_sub = expr.subs(x, y)\n",
    "expr_sub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "892a705d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve(Eq(expr_sub, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72a5af46",
   "metadata": {},
   "source": [
    "Below is another example equation with the symbol `x` and functions `sin`, `cos`, and `tan` using the `Eq` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c928a2be",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Create an equation\n",
    "eq = Eq(cos(x) / (tan(x)/sin(x)), 0)\n",
    "eq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab8a1f7b",
   "metadata": {},
   "source": [
    "Now we simplify this equation using the `simplify` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "916c9a09",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Simplify an expression\n",
    "simplified_expr = simplify(eq)\n",
    "simplified_expr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e53af384",
   "metadata": {},
   "source": [
    "Again, we use the `solve` function to solve this equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ae338df",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Solve the equation\n",
    "sol = solve(eq, x)\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b141fdd9",
   "metadata": {},
   "source": [
    "SymPy can also handle more complex equations involving trigonometry and complex numbers.\n",
    "\n",
    "We demonstrate this using [Euler’s formula](https://en.wikipedia.org/wiki/Euler%27s_formula)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d33661b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# 'I' represents the imaginary number i \n",
    "euler = cos(x) + I*sin(x)\n",
    "euler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ec5823f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "simplify(euler)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05abd043",
   "metadata": {},
   "source": [
    "If you are interested, we encourage you to read the lecture on [trigonometry and complex numbers](https://python.quantecon.org/complex_and_trig.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56690dcc",
   "metadata": {},
   "source": [
    "#### Example: fixed point computation\n",
    "\n",
    "Fixed point computation is frequently used in economics and finance.\n",
    "\n",
    "Here we solve the fixed point of the Solow-Swan growth dynamics:\n",
    "\n",
    "$$\n",
    "k_{t+1}=s f\\left(k_t\\right)+(1-\\delta) k_t, \\quad t=0,1, \\ldots\n",
    "$$\n",
    "\n",
    "where $ k_t $ is the capital stock, $ f $ is a production function, $ \\delta $ is a rate of depreciation.\n",
    "\n",
    "We are interested in calculating the fixed point of this dynamics, i.e., the value of $ k $ such that $ k_{t+1} = k_t $.\n",
    "\n",
    "With $ f(k) = Ak^\\alpha $, we can show the unique fixed point of the dynamics $ k^* $ using pen and paper:\n",
    "\n",
    "$$\n",
    "k^*:=\\left(\\frac{s A}{\\delta}\\right)^{1 /(1-\\alpha)}\n",
    "$$\n",
    "\n",
    "This can be easily computed in SymPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f31958ca",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A, s, k, α, δ = symbols('A s k^* α δ')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f23a0da9",
   "metadata": {},
   "source": [
    "Now we solve for the fixed point $ k^* $\n",
    "\n",
    "$$\n",
    "k^* = sA(k^*)^{\\alpha}+(1-\\delta) k^*\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8d9c3e1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Define Solow-Swan growth dynamics\n",
    "solow = Eq(s*A*k**α + (1-δ)*k, k)\n",
    "solow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f67d467",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve(solow, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d43c433b",
   "metadata": {},
   "source": [
    "### Inequalities and logic\n",
    "\n",
    "SymPy also allows users to define inequalities and set operators and provides a wide range of [operations](https://docs.sympy.org/latest/modules/solvers/inequalities.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee948286",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "reduce_inequalities([2*x + 5*y <= 30, 4*x + 2*y <= 20], [x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "974ff5cd",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "And(2*x + 5*y <= 30, x > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e4c672e",
   "metadata": {},
   "source": [
    "### Series\n",
    "\n",
    "Series are widely used in economics and statistics, from asset pricing to the expectation of discrete random variables.\n",
    "\n",
    "We can construct a simple series of summations using `Sum` function and `Indexed` symbols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "830b18f9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x, y, i, j = symbols(\"x y i j\")\n",
    "sum_xy = Sum(Indexed('x', i)*Indexed('y', j), \n",
    "            (i, 0, 3),\n",
    "            (j, 0, 3))\n",
    "sum_xy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d18da53",
   "metadata": {},
   "source": [
    "To evaluate the sum, we can [`lambdify`](https://docs.sympy.org/latest/modules/utilities/lambdify.html#sympy.utilities.lambdify.lambdify) the formula.\n",
    "\n",
    "The lambdified expression can take numeric values as input for $ x $ and $ y $ and compute the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "778b156b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "sum_xy = lambdify([x, y], sum_xy)\n",
    "grid = np.arange(0, 4, 1)\n",
    "sum_xy(grid, grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e13c9e9",
   "metadata": {},
   "source": [
    "#### Example: bank deposits\n",
    "\n",
    "Imagine a bank with $ D_0 $ as the deposit at time $ t $.\n",
    "\n",
    "It loans $ (1-r) $ of its deposits and keeps a fraction $ r $ as cash reserves.\n",
    "\n",
    "Its deposits over an infinite time horizon can be written as\n",
    "\n",
    "$$\n",
    "\\sum_{i=0}^\\infty (1-r)^i D_0\n",
    "$$\n",
    "\n",
    "Let’s compute the deposits at time $ t $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94a34ad3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "D = symbols('D_0')\n",
    "r = Symbol('r', positive=True)\n",
    "Dt = Sum('(1 - r)^i * D_0', (i, 0, oo))\n",
    "Dt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfcee572",
   "metadata": {},
   "source": [
    "We can call the `doit` method to evaluate the series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35a0f0b0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Dt.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c84f6434",
   "metadata": {},
   "source": [
    "Simplifying the expression above gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d582f7b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "simplify(Dt.doit())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32c90ec2",
   "metadata": {},
   "source": [
    "This is consistent with the solution in the lecture on [geometric series](https://python.quantecon.org/geom_series.html#example-the-money-multiplier-in-fractional-reserve-banking)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eddb286c",
   "metadata": {},
   "source": [
    "#### Example: discrete random variable\n",
    "\n",
    "In the following example, we compute the expectation of a discrete random variable.\n",
    "\n",
    "Let’s define a discrete random variable $ X $ following a [Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution):\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{\\lambda^x e^{-\\lambda}}{x!}, \\quad x = 0, 1, 2, \\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "459a49d3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ = symbols('lambda')\n",
    "\n",
    "# We refine the symbol x to positive integers\n",
    "x = Symbol('x', integer=True, positive=True)\n",
    "pmf = λ**x * exp(-λ) / factorial(x)\n",
    "pmf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cb70416",
   "metadata": {},
   "source": [
    "We can verify if the sum of probabilities for all possible values equals $ 1 $:\n",
    "\n",
    "$$\n",
    "\\sum_{x=0}^{\\infty} f(x) = 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "744cfe27",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "sum_pmf = Sum(pmf, (x, 0, oo))\n",
    "sum_pmf.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afec2161",
   "metadata": {},
   "source": [
    "The expectation of the distribution is:\n",
    "\n",
    "$$\n",
    "E(X) = \\sum_{x=0}^{\\infty} x f(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e17d5e02",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fx = Sum(x*pmf, (x, 0, oo))\n",
    "fx.doit()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "511a5b3f",
   "metadata": {},
   "source": [
    "SymPy includes a statistics submodule called [`Stats`](https://docs.sympy.org/latest/modules/stats.html).\n",
    "\n",
    "`Stats` offers built-in distributions and functions on probability distributions.\n",
    "\n",
    "The computation above can also be condensed into one line using the expectation function `E` in the `Stats` module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bd9a140",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ = Symbol(\"λ\", positive = True)\n",
    "\n",
    "# Using sympy.stats.Poisson() method\n",
    "X = Poisson(\"x\", λ)\n",
    "E(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cab068e",
   "metadata": {},
   "source": [
    "## Symbolic Calculus\n",
    "\n",
    "SymPy allows us to perform various calculus operations, such as limits, differentiation, and integration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fac5f6ae",
   "metadata": {},
   "source": [
    "### Limits\n",
    "\n",
    "We can compute limits for a given expression using the `limit` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "388c973d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Define an expression\n",
    "f = x**2 / (x-1)\n",
    "\n",
    "# Compute the limit\n",
    "lim = limit(f, x, 0)\n",
    "lim"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24805542",
   "metadata": {},
   "source": [
    "### Derivatives\n",
    "\n",
    "We can differentiate any SymPy expression using the `diff` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94aa3442",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Differentiate a function with respect to x\n",
    "df = diff(f, x)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd3f2f89",
   "metadata": {},
   "source": [
    "### Integrals\n",
    "\n",
    "We can compute definite and indefinite integrals using the `integrate` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5920c7d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Calculate the indefinite integral\n",
    "indef_int = integrate(df, x)\n",
    "indef_int"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8c27d11",
   "metadata": {},
   "source": [
    "Let’s use this function to compute the moment-generating function of [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution) with the probability density function:\n",
    "\n",
    "$$\n",
    "f(x) = \\lambda e^{-\\lambda x}, \\quad x \\ge 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dd72613",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ = Symbol('lambda', positive=True)\n",
    "x = Symbol('x', positive=True)\n",
    "pdf = λ * exp(-λ*x)\n",
    "pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad32bb0b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "t = Symbol('t', positive=True)\n",
    "moment_t = integrate(exp(t*x) * pdf, (x, 0, oo))\n",
    "simplify(moment_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "438b490c",
   "metadata": {},
   "source": [
    "Note that we can also use `Stats` module to compute the moment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c99ad0df",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = Exponential(x, λ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eace325",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "moment(X, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54367c5a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "E(X**t)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bf3405e",
   "metadata": {},
   "source": [
    "Using the `integrate` function, we can derive the cumulative density function of the exponential distribution with $ \\lambda = 0.5 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cffc091a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ_pdf = pdf.subs(λ, 1/2)\n",
    "λ_pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55ee9523",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "integrate(λ_pdf, (x, 0, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9f976d3",
   "metadata": {},
   "source": [
    "Using `cdf` in `Stats` module gives the same solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ed41df2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "cdf(X, 1/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf8d494c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Plug in a value for z \n",
    "λ_cdf = cdf(X, 1/2)(4)\n",
    "λ_cdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c0674a7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Substitute λ\n",
    "λ_cdf.subs({λ: 1/2})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41310dde",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "\n",
    "SymPy provides a powerful plotting feature.\n",
    "\n",
    "First we plot a simple function using the `plot` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55fe6697",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "f = sin(2 * sin(2 * sin(2 * sin(x))))\n",
    "p = plot(f, (x, -10, 10), show=False)\n",
    "p.title = 'A Simple Plot'\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab9b9107",
   "metadata": {},
   "source": [
    "Similar to Matplotlib, SymPy provides an interface to customize the graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cc1b69d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_f = plot(f, (x, -10, 10), \n",
    "              xlabel='', ylabel='', \n",
    "              legend = True, show = False)\n",
    "plot_f[0].label = 'f(x)'\n",
    "df = diff(f)\n",
    "plot_df = plot(df, (x, -10, 10), \n",
    "            legend = True, show = False)\n",
    "plot_df[0].label = 'f\\'(x)'\n",
    "plot_f.append(plot_df[0])\n",
    "plot_f.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "556c1c6f",
   "metadata": {},
   "source": [
    "It also supports plotting implicit functions and visualizing inequalities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16e2696e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p = plot_implicit(Eq((1/x + 1/y)**2, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a0347d4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p = plot_implicit(And(2*x + 5*y <= 30, 4*x + 2*y >= 20),\n",
    "                     (x, -1, 10), (y, -10, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aa51a44",
   "metadata": {},
   "source": [
    "and visualizations in three-dimensional space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eabac57",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p = plot3d(cos(2*x + y), zlabel='')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "350d7285",
   "metadata": {},
   "source": [
    "## Application: Two-person Exchange Economy\n",
    "\n",
    "Imagine a pure exchange economy with two people ($ a $ and $ b $) and two goods recorded as proportions ($ x $ and $ y $).\n",
    "\n",
    "They can trade goods with each other according to their preferences.\n",
    "\n",
    "Assume that the utility functions of the consumers are given by\n",
    "\n",
    "$$\n",
    "u_a(x, y) = x^{\\alpha} y^{1-\\alpha}\n",
    "$$\n",
    "\n",
    "$$\n",
    "u_b(x, y) = (1 - x)^{\\beta} (1 - y)^{1-\\beta}\n",
    "$$\n",
    "\n",
    "where $ \\alpha, \\beta \\in (0, 1) $.\n",
    "\n",
    "First we define the symbols and utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a14309e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Define symbols and utility functions\n",
    "x, y, α, β = symbols('x, y, α, β')\n",
    "u_a = x**α * y**(1-α)\n",
    "u_b = (1 - x)**β * (1 - y)**(1 - β)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e32a43c2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "u_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67a560c1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "u_b"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a306da7",
   "metadata": {},
   "source": [
    "We are interested in the Pareto optimal allocation of goods $ x $ and $ y $.\n",
    "\n",
    "Note that a point is Pareto efficient when the allocation is optimal for one person given the allocation for the other person.\n",
    "\n",
    "In terms of marginal utility:\n",
    "\n",
    "$$\n",
    "\\frac{\\frac{\\partial u_a}{\\partial x}}{\\frac{\\partial u_a}{\\partial y}} = \\frac{\\frac{\\partial u_b}{\\partial x}}{\\frac{\\partial u_b}{\\partial y}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4be74526",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# A point is Pareto efficient when the allocation is optimal \n",
    "# for one person given the allocation for the other person\n",
    "\n",
    "pareto = Eq(diff(u_a, x)/diff(u_a, y), \n",
    "            diff(u_b, x)/diff(u_b, y))\n",
    "pareto"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd494e5e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Solve the equation\n",
    "sol = solve(pareto, y)[0]\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b39d29a1",
   "metadata": {},
   "source": [
    "Let’s compute the Pareto optimal allocations of the economy (contract curves) with $ \\alpha = \\beta = 0.5 $ using SymPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a901d5da",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Substitute α = 0.5 and β = 0.5\n",
    "sol.subs({α: 0.5, β: 0.5})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddbc55b0",
   "metadata": {},
   "source": [
    "We can use this result to visualize more contract curves under different parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "539ca50c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Plot a range of αs and βs\n",
    "params = [{α: 0.5, β: 0.5}, \n",
    "          {α: 0.1, β: 0.9},\n",
    "          {α: 0.1, β: 0.8},\n",
    "          {α: 0.8, β: 0.9},\n",
    "          {α: 0.4, β: 0.8}, \n",
    "          {α: 0.8, β: 0.1},\n",
    "          {α: 0.9, β: 0.8},\n",
    "          {α: 0.8, β: 0.4},\n",
    "          {α: 0.9, β: 0.1}]\n",
    "\n",
    "p = plot(xlabel='x', ylabel='y', show=False)\n",
    "\n",
    "for param in params:\n",
    "    p_add = plot(sol.subs(param), (x, 0, 1), \n",
    "                 show=False)\n",
    "    p.append(p_add[0])\n",
    "p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e650b983",
   "metadata": {},
   "source": [
    "We invite you to play with the parameters and see how the contract curves change and think about the following two questions:\n",
    "\n",
    "- Can you think of a way to draw the same graph using `numpy`?  \n",
    "- How difficult will it be to write a `numpy` implementation?  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cdb6cb3",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0700481d",
   "metadata": {},
   "source": [
    "## Exercise 15.1\n",
    "\n",
    "L’Hôpital’s rule states that for two functions $ f(x) $ and $ g(x) $, if $ \\lim_{x \\to a} f(x) = \\lim_{x \\to a} g(x) = 0 $ or $ \\pm \\infty $, then\n",
    "\n",
    "$$\n",
    "\\lim_{x \\to a} \\frac{f(x)}{g(x)} = \\lim_{x \\to a} \\frac{f'(x)}{g'(x)}\n",
    "$$\n",
    "\n",
    "Use SymPy to verify L’Hôpital’s rule for the following functions\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{y^x - 1}{x}\n",
    "$$\n",
    "\n",
    "as $ x $ approaches to $ 0 $"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "169c49fc",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 15.1](https://python-programming.quantecon.org/#sympy_ex1)\n",
    "\n",
    "Let’s define the function first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b844393",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "f_upper = y**x - 1\n",
    "f_lower = x\n",
    "f = f_upper/f_lower\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e93e8fb4",
   "metadata": {},
   "source": [
    "Sympy is smart enough to solve this limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a0ae33a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lim = limit(f, x, 0)\n",
    "lim"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b562372",
   "metadata": {},
   "source": [
    "We compare the result suggested by L’Hôpital’s rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4be303c8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lim = limit(diff(f_upper, x)/\n",
    "            diff(f_lower, x), x, 0)\n",
    "lim"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c25e6d8",
   "metadata": {},
   "source": [
    "## Exercise 15.2\n",
    "\n",
    "[Maximum likelihood estimation (MLE)](https://python.quantecon.org/mle.html) is a method to estimate the parameters of a statistical model.\n",
    "\n",
    "It usually involves maximizing a log-likelihood function and solving the first-order derivative.\n",
    "\n",
    "The binomial distribution is given by\n",
    "\n",
    "$$\n",
    "f(x; n, θ) = \\frac{n!}{x!(n-x)!}θ^x(1-θ)^{n-x}\n",
    "$$\n",
    "\n",
    "where $ n $ is the number of trials and $ x $ is the number of successes.\n",
    "\n",
    "Assume we observed a series of binary outcomes with $ x $ successes out of $ n $ trials.\n",
    "\n",
    "Compute the MLE of $ θ $ using SymPy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb2693ac",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 15.2](https://python-programming.quantecon.org/#sympy_ex2)\n",
    "\n",
    "First, we define the binomial distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c39efb61",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n, x, θ = symbols('n x θ')\n",
    "\n",
    "binomial_factor = (factorial(n)) / (factorial(x)*factorial(n-r))\n",
    "binomial_factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f371d506",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "bino_dist = binomial_factor * ((θ**x)*(1-θ)**(n-x))\n",
    "bino_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e471048",
   "metadata": {},
   "source": [
    "Now we compute the log-likelihood function and solve for the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "025d9807",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "log_bino_dist = log(bino_dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4bd89ab",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "log_bino_diff = simplify(diff(log_bino_dist, θ))\n",
    "log_bino_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "636c568e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve(Eq(log_bino_diff, 0), θ)[0]"
   ]
  }
 ],
 "metadata": {
  "date": 1710455933.1206274,
  "filename": "sympy.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "SymPy"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}