# Characteristics of Autonomous Market Makers

Date: 2021-03-26 

(Date started: December 2020 Christmas holidays)

Author: [@mattmcd](https://twitter.com/mattmcd)

This notebook describes an approach using [method of characteristics](https://en.wikipedia.org/wiki/Method_of_characteristics) solutions of partial differential equations (PDEs) to examine existing AMM invariants such as Uniswap, Balancer and Curve (a.k.a. Stableswap).  


```python
from IPython.display import HTML
# Hide code cells https://gist.github.com/uolter/970adfedf44962b47d32347d262fe9be  
def hide_code():
    return HTML('''
    The raw code for this IPython notebook is by default hidden for easier reading.
    To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.''')
```


```python
import sympy as sp
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

%matplotlib inline
```


```python
hide_code()
```





    The raw code for this IPython notebook is by default hidden for easier reading.
    To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.



From the [Balancer whitepaper](https://balancer.finance/whitepaper/):

> The bedrock of Balancer’s exchange functions is a surface defined by constraining a value function $V$
 — a function of the pool’s weights and balances — to a constant. We will prove that this surface implies a spot price at each point such that, no matter what exchanges are carried out, the share of value of each token in the pool remains constant.



The [Balancer whitepaper](https://balancer.finance/whitepaper/) shows that the value function
$$V = \sum_{i=1}^{n} x_{i}^{w_{i}}$$
is related to the token spot prices by the ratio of partial derivative.  

# Starting from the Constraints
The idea of constant level sets of a value function creating constraints on system state (including prices) is discussed in ['From Curved Bonding to Configuration Spaces'](https://epub.wu.ac.at/7385/) by Zargham, Shorish, and Paruch (2020).  

The existing Balancer value function implicit state constraint is 'the share of value of each token in the pool remains constant'.  

In this section we look at starting from a set of constraints and see if it is possible to derive the corresponding value function.  We can then use this value function to determine allowed state changes e.g. for a swap the number of output tokens for an initial state and given number of input tokens.

This approach feels familiar to the Lagrangian dynamics approach in classical physics (the author's background).  In economics the standard approach seems to be start from a value function (a.k.a. utility function) and derive substitution functions that give prices.  Here we attempt to solve the inverse problem.

As a starting point, we consider deriving the Balancer value function (a [Cobb-Douglas Utility Function](https://www.wikiwand.com/en/Cobb%E2%80%93Douglas_production_function)) from the set of constraints for swaps 'the share of value of each token in the pool remains constant'.  

We consider below two and three asset pools with tokens $X, Y, Z$.  The state of the system can be defined by three token balances $x$, $y$, $z$, and three weights $w_x$, $w_y$, $w_z$.

## Uniswap
The two asset case is a generalized form of [Uniswap](https://uniswap.org/docs/v2/protocol-overview/how-uniswap-works/) having token weights $w_x$ and $w_y$ summing to 1.  Uniswap uses $w_x = w_y = \frac{1}{2}$.  Here we use $x$ and $y$ to represent the token X and token Y balances.

The total value of the pool in tems of token X is  
$$v_{x}(x,y) = x + p_{x}^{y}y = x - \frac{\partial{x}}{\partial{y}}y$$
i.e. number of tokens X plus the spot price of converting Y tokens into X tokens.

The constant share of value constraint is hence represented by the equations: 
$$w_{x} = \frac{x}{x - \frac{\partial{x}}{\partial{y}}y} \\
w_y = \frac{- \frac{\partial{x}}{\partial{y}}y}{x - \frac{\partial{x}}{\partial{y}}y}$$
where in the two asset case the second equation is redundant since the weights sum to 1. 


```python
x, y, w_x, w_y = sp.symbols('x y w_x w_y', positive=True)
k = sp.symbols('k', real=True)
X, Y = map(sp.Function, 'XY')
V = sp.Function('V')
```

We can specify that swaps happen on some invariant surface $V(x,y)$ which allows us to replace the spot price $-\frac{\partial{x}}{\partial{y}}$ in  $x - y\frac{\partial{x}}{\partial{y}} = w_x$, subsitituting $-\frac{\partial{x}}{\partial{y}} = \frac{\partial{V}}{\partial{y}}/\frac{\partial{V}}{\partial{x}}$ via the [implicit function theorem](https://en.wikipedia.org/wiki/Implicit_function_theorem).


```python
sp.Eq(x/(x + y*V(x,y).diff(y)/V(x,y).diff(x)), w_x)
```




$\displaystyle \frac{x}{x + \frac{y \frac{\partial}{\partial y} V{\left(x,y \right)}}{\frac{\partial}{\partial x} V{\left(x,y \right)}}} = w_{x}$



SymPy's PDE solver balks at this equation as written, so multiply through to make things easier.


```python
const_share_eq = sp.Eq(x*V(x,y).diff(x), w_x*(x*V(x,y).diff(x) + y*V(x,y).diff(y)))
```


```python
const_share_eq
```




$\displaystyle x \frac{\partial}{\partial x} V{\left(x,y \right)} = w_{x} \left(x \frac{\partial}{\partial x} V{\left(x,y \right)} + y \frac{\partial}{\partial y} V{\left(x,y \right)}\right)$



It turns out that SymPy is capable of solving the PDE directly to give a general solution.


```python
V_sol = sp.pdsolve(const_share_eq).subs({(1-w_x): w_y})
```


```python
V_sol
```




$\displaystyle V{\left(x,y \right)} = F{\left(- \frac{\log{\left(x^{- w_{x}} y^{- w_{y}} \right)}}{w_{y}} \right)}$



We can simplify the solution:


```python
sp.Eq(V(x,y), V_sol.rhs.simplify())
```




$\displaystyle V{\left(x,y \right)} = F{\left(\frac{w_{x} \log{\left(x \right)}}{w_{y}} + \log{\left(y \right)} \right)}$



We show below that the spot price is as expected regardless of the exact form of $F$ and if a specific form is chosen we achieve the desired Cobb-Douglas form.  


```python
sp.Eq((V(x,y).diff(y)/V(x,y).diff(x)), V_sol.rhs.diff(y)/V_sol.rhs.diff(x))
```




$\displaystyle \frac{\frac{\partial}{\partial y} V{\left(x,y \right)}}{\frac{\partial}{\partial x} V{\left(x,y \right)}} = \frac{w_{y} x}{w_{x} y}$



By taking the exponential of the constant we obtain the general Uniswap invariant.


```python
sp.Eq(V(x,y), sp.exp(w_y*V_sol.rhs.args[0]).simplify())
```




$\displaystyle V{\left(x,y \right)} = x^{w_{x}} y^{w_{y}}$



Interestingly, the general solution offers the opportunity to use different functional forms to achieve the same constant share of value contraint.  


```python
sp.Eq(V(x,y), (w_y*V_sol.rhs.args[0]).expand())
```




$\displaystyle V{\left(x,y \right)} = w_{x} \log{\left(x \right)} + w_{y} \log{\left(y \right)}$



We could derive the swap formulae for each form.  Would the formulae from log form of invariant be easier to implement on the Ethererum blockchain?  This is left as an exercise for the reader.

## Balancer
We now consider the general case of a Balancer pool with $n$ assets and weights summing to 1.  Here we examine the three asset case and identify the geometric constraints imposed by the share of value conditions
$$
w_x = \frac{x}{x - \frac{\partial{x}}{\partial{y}}y - \frac{\partial{x}}{\partial{z}}z}\\
w_y = \frac{-\frac{\partial{x}}{\partial{y}}y}{x - \frac{\partial{x}}{\partial{y}}y - \frac{\partial{x}}{\partial{z}}z}\\
w_z = \frac{-\frac{\partial{x}}{\partial{z}}z}{x - \frac{\partial{x}}{\partial{y}}y - \frac{\partial{x}}{\partial{z}}z}
$$
for tokens X, Y and Z having total value $v_x(x,y,z) = x + p_{x}^{y}y + p_{x}^{z}z = x - \frac{\partial{x}}{\partial{y}}y - \frac{\partial{x}}{\partial{z}}z$.

As before we can replace the spot prices given by the partial derivatives with expressions for the invariant surface V.  Below we show the $w_x$ condition, the others are similar.  The weight sum to 1 condition again allows us to eliminate one of the constraint equations.


```python
x, y, z, w_x, w_y, w_z = sp.symbols('x y z w_x w_y w_z', positive=True)
u, v, w = sp.symbols('u v w', positive=True)
xi = sp.Function('xi')
eta = sp.Function('eta')
zeta = sp.Function('zeta')
V = sp.Function('V')
V_d = sp.Function('\mathcal{V}')
```


```python
const_share_eq_3_1 = sp.Eq(x*V(x,y,z).diff(x), w_x*(x*V(x,y,z).diff(x) + y*V(x,y,z).diff(y) + z*V(x,y,z).diff(z)))
```


```python
const_share_eq_3_1
```




$\displaystyle x \frac{\partial}{\partial x} V{\left(x,y,z \right)} = w_{x} \left(x \frac{\partial}{\partial x} V{\left(x,y,z \right)} + y \frac{\partial}{\partial y} V{\left(x,y,z \right)} + z \frac{\partial}{\partial z} V{\left(x,y,z \right)}\right)$



We can simplify this to a constant coefficient first order PDE by the change of variables 
$$
u = \xi(x) = \log{x}\\
v = \eta(y) = \log{y}\\
w = \zeta(z) = \log{z}
$$

The general form for the the change of variables is given by


```python
const_share_eq_3_1.subs({
    V(x,y,z): V(xi(x),eta(y),zeta(z)),
}).simplify().subs({V(xi(x),eta(y),zeta(z)): V_d })
```




$\displaystyle x \frac{d}{d \xi{\left(x \right)}} \mathcal{V} \frac{d}{d x} \xi{\left(x \right)} = w_{x} \left(x \frac{d}{d \xi{\left(x \right)}} \mathcal{V} \frac{d}{d x} \xi{\left(x \right)} + y \frac{d}{d \eta{\left(y \right)}} \mathcal{V} \frac{d}{d y} \eta{\left(y \right)} + z \frac{d}{d \zeta{\left(z \right)}} \mathcal{V} \frac{d}{d z} \zeta{\left(z \right)}\right)$



where $\mathcal{V} = V(\xi(x), \eta(y), \zeta(z))$

Substituting the logarithmic functional form for the transformed variables gives


```python
const_share_eq_3_1.subs({
    V(x,y,z): V(xi(x),eta(y),zeta(z)),
    xi(x): sp.log(x),
    eta(y): sp.log(y), 
    zeta(z): sp.log(z)
}).simplify()
```




$\displaystyle \left. \frac{\partial}{\partial \xi_{1}} V{\left(\xi_{1},\log{\left(y \right)},\log{\left(z \right)} \right)} \right|_{\substack{ \xi_{1}=\log{\left(x \right)} }} = w_{x} \left(\left. \frac{\partial}{\partial \xi_{1}} V{\left(\xi_{1},\log{\left(y \right)},\log{\left(z \right)} \right)} \right|_{\substack{ \xi_{1}=\log{\left(x \right)} }} + \left. \frac{\partial}{\partial \xi_{2}} V{\left(\log{\left(x \right)},\xi_{2},\log{\left(z \right)} \right)} \right|_{\substack{ \xi_{2}=\log{\left(y \right)} }} + \left. \frac{\partial}{\partial \xi_{3}} V{\left(\log{\left(x \right)},\log{\left(y \right)},\xi_{3} \right)} \right|_{\substack{ \xi_{3}=\log{\left(z \right)} }}\right)$



The share of value weight constraints can written as a matrix equation:


```python
balancer_constraints = sp.ImmutableMatrix([
    [1-w_x, -w_x, -w_x],
    [-w_y, 1-w_y, -w_y],
    [-w_z, -w_z, 1-w_z]
])
v_grad = sp.ImmutableMatrix([V(u,v,w).diff(i) for i in [u,v,w]])
```


```python
sp.Eq(sp.MatMul(balancer_constraints, sp.UnevaluatedExpr(v_grad), evaluate=False).subs({V(u,v,w): V_d }),
      (balancer_constraints*v_grad).subs({V(u,v,w): V_d }), 
     evaluate=False)
```




$\displaystyle \left[\begin{matrix}1 - w_{x} & - w_{x} & - w_{x}\\- w_{y} & 1 - w_{y} & - w_{y}\\- w_{z} & - w_{z} & 1 - w_{z}\end{matrix}\right] \left[\begin{matrix}\frac{d}{d u} \mathcal{V}\\\frac{d}{d v} \mathcal{V}\\\frac{d}{d w} \mathcal{V}\end{matrix}\right] = \left[\begin{matrix}- w_{x} \frac{d}{d v} \mathcal{V} - w_{x} \frac{d}{d w} \mathcal{V} + \left(1 - w_{x}\right) \frac{d}{d u} \mathcal{V}\\- w_{y} \frac{d}{d u} \mathcal{V} - w_{y} \frac{d}{d w} \mathcal{V} + \left(1 - w_{y}\right) \frac{d}{d v} \mathcal{V}\\- w_{z} \frac{d}{d u} \mathcal{V} - w_{z} \frac{d}{d v} \mathcal{V} + \left(1 - w_{z}\right) \frac{d}{d w} \mathcal{V}\end{matrix}\right]$



We can use the weights summing to 1 condition to eliminate one of these equations.  The nullspace of the constraint matrix then defines a plane of constant value of the invariant in $u,v,w$ space which we can then map back to the original X, Y, Z token balances.


```python
(sp.simplify(sp.ImmutableMatrix([
#     [1-w_x, -w_x, -w_x],
    [-w_y, 1-w_y, -w_y],
    [-w_z, -w_z, 1-w_z]
]).nullspace()[0]).subs({(1-w_y-w_z): w_x})*w_z
).T*sp.ImmutableMatrix([sp.log(x), sp.log(y), sp.log(z)])
```




$\displaystyle \left[\begin{matrix}w_{x} \log{\left(x \right)} + w_{y} \log{\left(y \right)} + w_{z} \log{\left(z \right)}\end{matrix}\right]$



It can be seen by inspection that exponentiating this invariant results in the original Balancer form.  As for the two asset case it is also possible to retain the value function in this form and derive new forms for the trading functions, which is again an exercise left to the reader.

## Curve 
Curve uses the [StableSwap](https://www.curve.fi/stableswap-paper.pdf) invariant to reduce slippage for stablecoins all having equivalent value.  For example a pool could consist of USDC, USDT, BUSD and DAI, all of which are designed to track USD.  Pools tracking other assets are possible e.g. a BTC pool backed by sBTC, renBTC, wBTC.  The market maker token pool ideally consists of a balanced mix of each token type.  We consider a two asset pool below, this analysis extends to an arbitrary number of assets.

The StableSwap invariant is designed to act as a constant sum market maker $x+y=1$ for small imbalances, and a constant product Uniswap market maker $xy=k$ as the pool becomes more imbalanced.  These are the price constraints defining the system i.e.

* at small imbalance $-\frac{\partial{x}}{\partial{y}}=1$ and tokens are freely interchangeable
* at larger imbalance $-\frac{\partial{x}}{\partial{y}}=x/y$ as for Uniswap (or in general an equal weight Balancer pool)

The StableSwap invariant can be written as $V{\left(x,y \right)} = s \left(x + y\right) + x^{w_{x}} y^{w_{y}}$ where $s$ is an amplification parameter that determines the transition between constant sum and constant product behaviour.



```python
s = sp.symbols('s', positive=True)
```


```python
V_ss = s*(x+y) + x**w_x*y**w_y
```


```python
sp.Eq(V(x,y), V_ss)
```




$\displaystyle V{\left(x,y \right)} = s \left(x + y\right) + x^{w_{x}} y^{w_{y}}$



Below we show the spot price of Y tokens in terms of X tokens that results from this invariant function.  We can see that the limit $s \rightarrow \infty$ gives the constant sum behaviour while the $s \rightarrow 0$ limit gives constant product behaviour.


```python
V_x = (V_ss).diff(x).subs({w_x: sp.Rational(1,2), w_y: sp.Rational(1,2)}).simplify()
V_y = (V_ss).diff(y).subs({w_x: sp.Rational(1,2), w_y: sp.Rational(1,2)}).simplify()
ss_spot = V_y/V_x.simplify()
ss_spot_eq = sp.Eq(V(x,y).diff(y)/V(x,y).diff(x), ss_spot)
ss_spot_eq
```




$\displaystyle \frac{\frac{\partial}{\partial y} V{\left(x,y \right)}}{\frac{\partial}{\partial x} V{\left(x,y \right)}} = \frac{s + \frac{\sqrt{x}}{2 \sqrt{y}}}{s + \frac{\sqrt{y}}{2 \sqrt{x}}}$




```python
sp.Eq(sp.Limit(ss_spot, s, sp.oo), (ss_spot).limit(s, sp.oo))
```




$\displaystyle \lim_{s \to \infty}\left(\frac{s + \frac{\sqrt{x}}{2 \sqrt{y}}}{s + \frac{\sqrt{y}}{2 \sqrt{x}}}\right) = 1$




```python
sp.Eq(sp.Limit(ss_spot, s, 0), (ss_spot).limit(s, 0))
```




$\displaystyle \lim_{s \to 0^+}\left(\frac{s + \frac{\sqrt{x}}{2 \sqrt{y}}}{s + \frac{\sqrt{y}}{2 \sqrt{x}}}\right) = \frac{x}{y}$



Following the previous procedure we'd hope to be able to solve the PDE for $V(x,y)$ from the spot price constraint:
$$\left(s + \frac{\sqrt{y}}{2 \sqrt{x}}\right) \frac{\partial}{\partial y} V{\left(x,y \right)} = \left(s + \frac{\sqrt{x}}{2 \sqrt{y}}\right) \frac{\partial}{\partial x} V{\left(x,y \right)}$$
Attempting this with the actual StableSwap spot price above doesn't immediately work although it should be possible to solve numerically.  


```python
ss_spot_denom = sp.denom(ss_spot_eq.rhs)* sp.denom(ss_spot_eq.lhs)
sp.Eq(ss_spot_eq.lhs * ss_spot_denom, ss_spot_eq.rhs * ss_spot_denom)
```




$\displaystyle \left(s + \frac{\sqrt{y}}{2 \sqrt{x}}\right) \frac{\partial}{\partial y} V{\left(x,y \right)} = \left(s + \frac{\sqrt{x}}{2 \sqrt{y}}\right) \frac{\partial}{\partial x} V{\left(x,y \right)}$



## A new Curve
We can look at the form of the PDE from the StableSwap constraint and explore related functional forms.  It looks like the same limiting behaviour could be achieved with any ratio of $x/y$ i.e. without requiring the $\surd$.  We hence try $$\left(s x y + y\right) \frac{\partial}{\partial y} V{\left(x,y \right)} = \left(s x y + x\right) \frac{\partial}{\partial x} V{\left(x,y \right)}$$
which is easily solvable by SymPy and results in a new Curve invariant:


```python
sol2 = sp.pdsolve((s*x*y + y)*V(x,y).diff(y) - (s*x*y + x)*V(x,y).diff(x), V(x,y)).rhs
V_ss_new = sp.log(sol2.args[0]).expand()
sp.Eq(V(x,y),V_ss_new)
```




$\displaystyle V{\left(x,y \right)} = s x + s y + \log{\left(x \right)} + \log{\left(y \right)}$



We see that the spot price has the same desired limiting behaviour: 


```python
ss_spot_new = (V_ss_new.diff(y)/V_ss_new.diff(x)).simplify()
```


```python
ss_spot_new
```




$\displaystyle \frac{x \left(s y + 1\right)}{y \left(s x + 1\right)}$




```python
sp.Eq(sp.Limit(ss_spot_new, s, sp.oo), (ss_spot_new).limit(s, sp.oo))
```




$\displaystyle \lim_{s \to \infty}\left(\frac{x \left(s y + 1\right)}{y \left(s x + 1\right)}\right) = 1$




```python
sp.Eq(sp.Limit(ss_spot_new, s, 0), (ss_spot_new).limit(s, 0))
```




$\displaystyle \lim_{s \to 0^+}\left(\frac{x \left(s y + 1\right)}{y \left(s x + 1\right)}\right) = \frac{x}{y}$



Comparing to the characteristic surfaces obtained for Uniswap and Balancer using this approach, we see that the functional form is quite similar i.e. a plane in $\log{x}-\log{y}$ space now with the addition of an exponential surface for the imbalance part.  

A nice feature of this new form is the possibility of investigating weighted Curve pools.

# Conclusion
This notebook has demonstrated using the method of characteristics to solve the inverse problem in autonomous market maker design i.e. deriving the invariant function that satisfies some constraints.  

Some possibly original results for alternative Balancer and Curve invariants are presented.

It is hoped that the wider DeFi community finds this technique and the use of a computer algebra system such as [SymPy](https://www.sympy.org/en/index.html) useful.  It offers a complementary approach to simulation based methods such as [cadCAD](https://twitter.com/cadcad_org).

# References and Related Work


## DeFi 
The DeFi space moves quite rapidly so as well as the references below also see Medium for Balancer, Uniswap, Curve etc.  Also [The Defiant](https://thedefiant.substack.com/) newsletter to keep up with things.

[Balancer whitepaper](https://balancer.finance/whitepaper/) by [@fcmartinelli](https://twitter.com/fcmartinelli).  

[From Curved Bonding to Configuration Spaces](https://epub.wu.ac.at/7385/) by [@mZargham](https://twitter.com/mZargham) 

Lots of papers by [@alexhevans](https://twitter.com/alexhevans), [@tarunchitra](https://twitter.com/tarunchitra) and [@GuilleAngeris](https://twitter.com/GuilleAngeris).  A couple of highly relevant entry points below, which annoyingly the author only got around to reading _after_ doing most of the work in this notebook.
* [When does the tail wag the dog? Curvature and Market Making](https://arxiv.org/abs/2012.08040)
* [Liquidity Provider Returns in Geometric Mean Markets](https://arxiv.org/abs/2006.08806)
* [Optimal Fees for Geometric Mean Market Makers](https://stanford.edu/~guillean/papers/g3m-optimal-fee.pdf)

## Mathematics
The author first encountered the method of characterstics last millenium during their undergraduate degree at University of Western Australia.  "Partial Differential Equations: An Introduction" by Strauss was the text but any reasonably advanced undergraduate engineering mathematics book should have something relevant.

Comparing the Balancer or Uniswap invariants to the Cobb-Douglas utility or production function isn't something the author has seen that widely and needs to learn more microeconomics to see how it could fit in.  Seems related to  [Optimization Methods in Finance](https://www.cambridge.org/core/books/optimization-methods-in-finance/8A4996C5DB2006224E4D983B5BC95E3B) and the [Production Possibility Frontier](https://www.wikiwand.com/en/Production%E2%80%93possibility_frontier) e.g. this [handout by Kyle Woodward](https://kylewoodward.com/blog_data/pdfs/handout_micro_equilibrium.pdf).


```python

```

a7aa78056fc4b7252aaa64a35b9dd8b5442056fb
ipynb
Jupyter Notebook
scripts/amm_math_20210308.ipynb
mattmcd/PyBayes
c3931cd495f1e96eea7e8f6b5c527403726de7d7
Apache-2.0
scripts/amm_math_20210308.ipynb
mattmcd/PyBayes
c3931cd495f1e96eea7e8f6b5c527403726de7d7
Apache-2.0
scripts/amm_math_20210308.ipynb
mattmcd/PyBayes
c3931cd495f1e96eea7e8f6b5c527403726de7d7
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES