SymPy - Python Programming for Economics and Finance

16.1Overview¶

Unlike numerical libraries that deal with values, SymPy focuses on manipulating mathematical symbols and expressions directly.

SymPy provides a wide range of features including

symbolic expression
equation solving
simplification
calculus
matrices
discrete math, etc.

These functions make SymPy a popular open-source alternative to other proprietary symbolic computational software such as Mathematica.

In this lecture, we will explore some of the functionality of SymPy and demonstrate how to use basic SymPy functions to solve economic models.

16.2Getting Started¶

Let’s first import the library and initialize the printer for symbolic output

from sympy import *
from sympy.plotting import plot, plot3d_parametric_line, plot3d
from sympy.solvers.inequalities import reduce_rational_inequalities
from sympy.stats import Poisson, Exponential, Binomial, density, moment, E, cdf

import numpy as np
import matplotlib.pyplot as plt

# Enable the mathjax printer
init_printing(use_latex='mathjax')

16.3Symbolic algebra¶

16.3.1Symbols¶

First we initialize some symbols to work with

x, y, z = symbols('x y z')

Symbols are the basic units for symbolic computation in SymPy.

16.3.2Expressions¶

We can now use symbols x, y, and z to build expressions and equations.

Here we build a simple expression first

expr = (x+y) ** 2
expr

We can expand this expression with the expand function

expand_expr = expand(expr)
expand_expr

and factorize it back to the factored form with the factor function

factor(expand_expr)

We can solve this expression

solve(expr)

Note this is equivalent to solving the following equation for x

(x + y)^2 = 0

(1)

16.3.3Equations¶

SymPy provides several functions to manipulate equations.

Let’s develop an equation with the expression we defined before

eq = Eq(expr, 0)
eq

Solving this equation with respect to $x$ gives the same output as solving the expression directly

solve(eq, x)

SymPy can handle equations with multiple solutions

eq = Eq(expr, 1)
solve(eq, x)

solve function can also combine multiple equations together and solve a system of equations

eq2 = Eq(x, y)
eq2

solve([eq, eq2], [x, y])

We can also solve for the value of $y$ by simply substituting $x$ with $y$

expr_sub = expr.subs(x, y)
expr_sub

solve(Eq(expr_sub, 1))

Below is another example equation with the symbol x and functions sin, cos, and tan using the Eq function

# Create an equation
eq = Eq(cos(x) / (tan(x)/sin(x)), 0)
eq

Now we simplify this equation using the simplify function

# Simplify an expression
simplified_expr = simplify(eq)
simplified_expr

Again, we use the solve function to solve this equation

# Solve the equation
sol = solve(eq, x)
sol

SymPy can also handle more complex equations involving trigonometry and complex numbers.

We demonstrate this using Euler’s formula

# 'I' represents the imaginary number i 
euler = cos(x) + I*sin(x)
euler

simplify(euler)

If you are interested, we encourage you to read the lecture on trigonometry and complex numbers.

16.3.3.1Example: fixed point computation¶

Fixed point computation is frequently used in economics and finance.

Here we solve the fixed point of the Solow-Swan growth dynamics:

k_{t+1}=s f\left(k_t\right)+(1-\delta) k_t, \quad t=0,1, \ldots

(2)

where $k_t$ is the capital stock, $f$ is a production function, δ is a rate of depreciation.

We are interested in calculating the fixed point of this dynamics, i.e., the value of $k$ such that $k_{t+1} = k_t$ .

With $f(k) = Ak^\alpha$ , we can show the unique fixed point of the dynamics $k^*$ using pen and paper:

k^*:=\left(\frac{s A}{\delta}\right)^{1 /(1-\alpha)}

(3)

This can be easily computed in SymPy

A, s, k, α, δ = symbols('A s k^* α δ')

Now we solve for the fixed point $k^*$

k^* = sA(k^*)^{\alpha}+(1-\delta) k^*

(4)

# Define Solow-Swan growth dynamics
solow = Eq(s*A*k**α + (1-δ)*k, k)
solow

solve(solow, k)

16.3.4Inequalities and logic¶

SymPy also allows users to define inequalities and set operators and provides a wide range of operations.

reduce_inequalities([2*x + 5*y <= 30, 4*x + 2*y <= 20], [x])

And(2*x + 5*y <= 30, x > 0)

16.3.5Series¶

Series are widely used in economics and statistics, from asset pricing to the expectation of discrete random variables.

We can construct a simple series of summations using Sum function and Indexed symbols

x, y, i, j = symbols("x y i j")
sum_xy = Sum(Indexed('x', i)*Indexed('y', j), 
            (i, 0, 3),
            (j, 0, 3))
sum_xy

To evaluate the sum, we can lambdify the formula.

The lambdified expression can take numeric values as input for $x$ and $y$ and compute the result

sum_xy = lambdify([x, y], sum_xy)
grid = np.arange(0, 4, 1)
sum_xy(grid, grid)

np.int64(36)

16.3.5.1Example: bank deposits¶

Imagine a bank with $D_0$ as the deposit at time $t$ .

It loans $(1-r)$ of its deposits and keeps a fraction $r$ as cash reserves.

Its deposits over an infinite time horizon can be written as

\sum_{i=0}^\infty (1-r)^i D_0

(5)

Let’s compute the deposits at time $t$

D = symbols('D_0')
r = Symbol('r', positive=True)
Dt = Sum('(1 - r)^i * D_0', (i, 0, oo))
Dt

We can call the doit method to evaluate the series

Dt.doit()

Simplifying the expression above gives

simplify(Dt.doit())

This is consistent with the solution in the lecture on geometric series.

16.3.5.2Example: discrete random variable¶

In the following example, we compute the expectation of a discrete random variable.

Let’s define a discrete random variable $X$ following a Poisson distribution:

f(x) = \frac{\lambda^x e^{-\lambda}}{x!}, \quad x = 0, 1, 2, \ldots

(6)

λ = symbols('lambda')

# We refine the symbol x to positive integers
x = Symbol('x', integer=True, positive=True)
pmf = λ**x * exp(-λ) / factorial(x)
pmf

We can verify if the sum of probabilities for all possible values equals 1:

\sum_{x=0}^{\infty} f(x) = 1

(7)

sum_pmf = Sum(pmf, (x, 0, oo))
sum_pmf.doit()

The expectation of the distribution is:

E(X) = \sum_{x=0}^{\infty} x f(x)

(8)

fx = Sum(x*pmf, (x, 0, oo))
fx.doit()

SymPy includes a statistics submodule called Stats.

Stats offers built-in distributions and functions on probability distributions.

The computation above can also be condensed into one line using the expectation function E in the Stats module

λ = Symbol("λ", positive = True)

# Using sympy.stats.Poisson() method
X = Poisson("x", λ)
E(X)

16.4Symbolic Calculus¶

SymPy allows us to perform various calculus operations, such as limits, differentiation, and integration.

16.4.1Limits¶

We can compute limits for a given expression using the limit function

# Define an expression
f = x**2 / (x-1)

# Compute the limit
lim = limit(f, x, 0)
lim

16.4.2Derivatives¶

We can differentiate any SymPy expression using the diff function

# Differentiate a function with respect to x
df = diff(f, x)
df

16.4.3Integrals¶

We can compute definite and indefinite integrals using the integrate function

# Calculate the indefinite integral
indef_int = integrate(df, x)
indef_int

Let’s use this function to compute the moment-generating function of exponential distribution with the probability density function:

f(x) = \lambda e^{-\lambda x}, \quad x \ge 0

(9)

λ = Symbol('lambda', positive=True)
x = Symbol('x', positive=True)
pdf = λ * exp(-λ*x)
pdf

t = Symbol('t', positive=True)
moment_t = integrate(exp(t*x) * pdf, (x, 0, oo))
simplify(moment_t)

Note that we can also use Stats module to compute the moment

X = Exponential(x, λ)

moment(X, 1)

E(X**t)

Using the integrate function, we can derive the cumulative density function of the exponential distribution with $\lambda = 0.5$

λ_pdf = pdf.subs(λ, 1/2)
λ_pdf

integrate(λ_pdf, (x, 0, 4))

Using cdf in Stats module gives the same solution

cdf(X, 1/2)

# Plug in a value for z 
λ_cdf = cdf(X, 1/2)(4)
λ_cdf

# Substitute λ
λ_cdf.subs({λ: 1/2})

16.5Plotting¶

SymPy provides a powerful plotting feature.

First we plot a simple function using the plot function

f = sin(2 * sin(2 * sin(2 * sin(x))))
p = plot(f, (x, -10, 10), show=False)
p.title = 'A Simple Plot'
p.show()

Similar to Matplotlib, SymPy provides an interface to customize the graph

plot_f = plot(f, (x, -10, 10), 
              xlabel='', ylabel='', 
              legend = True, show = False)
plot_f[0].label = 'f(x)'
df = diff(f)
plot_df = plot(df, (x, -10, 10), 
            legend = True, show = False)
plot_df[0].label = 'f\'(x)'
plot_f.append(plot_df[0])
plot_f.show()

It also supports plotting implicit functions and visualizing inequalities

p = plot_implicit(Eq((1/x + 1/y)**2, 1))

p = plot_implicit(And(2*x + 5*y <= 30, 4*x + 2*y >= 20),
                     (x, -1, 10), (y, -10, 10))

and visualizations in three-dimensional space

p = plot3d(cos(2*x + y), zlabel='')

16.6Application: Two-person Exchange Economy¶

Imagine a pure exchange economy with two people ( $a$ and $b$ ) and two goods recorded as proportions ( $x$ and $y$ ).

They can trade goods with each other according to their preferences.

Assume that the utility functions of the consumers are given by

u_a(x, y) = x^{\alpha} y^{1-\alpha}

(10)

u_b(x, y) = (1 - x)^{\beta} (1 - y)^{1-\beta}

(11)

where $\alpha, \beta \in (0, 1)$ .

First we define the symbols and utility functions

# Define symbols and utility functions
x, y, α, β = symbols('x, y, α, β')
u_a = x**α * y**(1-α)
u_b = (1 - x)**β * (1 - y)**(1 - β)

u_a

u_b

We are interested in the Pareto optimal allocation of goods $x$ and $y$ .

Note that a point is Pareto efficient when the allocation is optimal for one person given the allocation for the other person.

In terms of marginal utility:

\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}}

(12)

# A point is Pareto efficient when the allocation is optimal 
# for one person given the allocation for the other person

pareto = Eq(diff(u_a, x)/diff(u_a, y), 
            diff(u_b, x)/diff(u_b, y))
pareto

# Solve the equation
sol = solve(pareto, y)[0]
sol

Let’s compute the Pareto optimal allocations of the economy (contract curves) with $\alpha = \beta = 0.5$ using SymPy

# Substitute α = 0.5 and β = 0.5
sol.subs({α: 0.5, β: 0.5})

We can use this result to visualize more contract curves under different parameters

# Plot a range of αs and βs
params = [{α: 0.5, β: 0.5}, 
          {α: 0.1, β: 0.9},
          {α: 0.1, β: 0.8},
          {α: 0.8, β: 0.9},
          {α: 0.4, β: 0.8}, 
          {α: 0.8, β: 0.1},
          {α: 0.9, β: 0.8},
          {α: 0.8, β: 0.4},
          {α: 0.9, β: 0.1}]

p = plot(xlabel='x', ylabel='y', show=False)

for param in params:
    p_add = plot(sol.subs(param), (x, 0, 1), 
                 show=False)
    p.append(p_add[0])
p.show()

We invite you to play with the parameters and see how the contract curves change and think about the following two questions:

Can you think of a way to draw the same graph using numpy?
How difficult will it be to write a numpy implementation?

16.7Exercises¶

Solution to Exercise 1

Let’s define the function first

f_upper = y**x - 1
f_lower = x
f = f_upper/f_lower
f

Sympy is smart enough to solve this limit

lim = limit(f, x, 0)
lim

We compare the result suggested by L’Hôpital’s rule

lim = limit(diff(f_upper, x)/
            diff(f_lower, x), x, 0)
lim

Solution to Exercise 2

First, we define the binomial distribution

n, x, θ = symbols('n x θ')

binomial_factor = (factorial(n)) / (factorial(x)*factorial(n-r))
binomial_factor

bino_dist = binomial_factor * ((θ**x)*(1-θ)**(n-x))
bino_dist

Now we compute the log-likelihood function and solve for the result

log_bino_dist = log(bino_dist)

log_bino_diff = simplify(diff(log_bino_dist, θ))
log_bino_diff

solve(Eq(log_bino_diff, 0), θ)[0]