Building Ito Calculus with SageMath: Deriving Black-Scholes

Introduction

In the previous post I shared my journey of getting SageMath installed in Windows and writing my first unit tests.

In this blog we will build our own symbolic Ito calculus layer on top of SageMath and use it to derive the Black-Scholes partial differential equation (PDE) directly from the stochastic differential equation (SDE) of the underlying asset.

Rather than manipulating the formulas by hand, we let the SageMath algebra engine apply our own custom-built Ito lemma, symbolically and verify the structure of the PDE step-by-step.

Setting Up Geometric Brownian Motion in Python

We begin with the standard geometric Brownian motion model for the underlying asset price:

$\displaystyle dS_t =\mu S_t dt + \sigma S_t dW_t$

where $\mu$ is the drift, $\sigma$ is the volatility and $W_t$ is Brownian motion. This is implemented as

Fig 1 – Symbolic representation of geometric Brownian motion in SageMath python.

with variables already declared as SageMath var types:

Fig 2 – Using SageMath var types, we first declare the variables that we will use.

Now, let’s assume the derivative price is a smooth function $V = V(t,s)$ , which is also declared in Python as

Fig 3 – Declaring our derivative price functional with the symbolic arguments that it needs.

Our goal is to compute $dV$ using Ito’s lemma.

Spoiler alert, it will be as simple as executing the following in python:

Code snippet demonstrating a Python function using the Ito calculus, importing necessary libraries and defining variables for stochastic processes.

with console output:

A code snippet showing formulas for Drift and Diffusion, including variables and mathematical expressions, displayed in a programming environment. — **Fig 4** – The output from running the above code.

Let’s get into it.

Encoding the Ito Algebra

Since SageMath does not natively implement stochastic calculus rules, such as $(dW_t)^2 = dt$ , we must first built a small algebra layer where differentials are represented as

$a + bdt + \displaystyle\sum_i c_i dW_i$

and multiplication rules enforce:

$dt^2 = 0$
$dt dW = 0$
$dW_i dW_j = \delta_{ij} dt$

and we must allow symbolic quadratic variation to emerge automatically.

from dataclasses import dataclass
from typing import Dict
from sage.all import SR, var, diff  # type: ignore
@dataclass(frozen=True)
class Ito:
    a: object
    b: object
    c: Dict[int, object]
    def __add__(self, other):
        if not isinstance(other, Ito):
            other = Ito.scalar(other)
        c_new = dict(self.c)
        for i, ci in other.c.items():
            c_new[i] = c_new.get(i, 0) + ci
            if c_new[i] == 0:
                del c_new[i]
        return Ito(self.a + other.a, self.b + other.b, c_new)
    def __radd__(self, other):
        # allows: 0 + Ito(...) and scalar + Ito(...)
        return self.__add__(other)
    def __neg__(self):
        return Ito(-self.a, -self.b, {i: -ci for i, ci in self.c.items()})
    def __sub__(self, other):
        return self.__add__(-other)
    def __rsub__(self, other):
        return Ito.scalar(other).__sub__(self)
    @staticmethod
    def scalar(x):
        return Ito(SR(x), 0, {})
    @staticmethod
    def dt(coeff=1):
        return Ito(0, SR(coeff), {})
    @staticmethod
    def dW(i, coeff=1):
        return Ito(0, 0, {i: SR(coeff)})
    def __mul__(self, other):
        if not isinstance(other, Ito):
            other = Ito.scalar(other)
        a1, b1, c1 = self.a, self.b, self.c
        a2, b2, c2 = other.a, other.b, other.c
        a = a1 * a2
        c = {}
        for i in set(c1) | set(c2):
            ci = a1 * c2.get(i, 0) + a2 * c1.get(i, 0)
            if ci != 0:
                c[i] = ci
        b = a1 * b2 + a2 * b1
        for i in set(c1) & set(c2):
            b += c1[i] * c2[i]
        return Ito(a, b, c)
    def __repr__(self):
        parts = [str(self.a)]
        if self.b != 0:
            parts.append(f"{self.b}*dt")
        for i in self.c:
            parts.append(f"{self.c[i]}*dW{i}")
        return " + ".join(parts)

The core idea of this class is simple: it represents stochastic differentials as truncated objects of the form

$a + bdt + \displaystyle\sum_i c_i dW_i$

where

$a$ is a scalar,
$b$ is the coefficient of $dt$ , and
$c_i$ are the coefficients of the Brownian increments.

You can see that we have explicitly encoded the Ito multiplication rules, which is where stochastic calculus differs fundamentally from classical calculus.

The Structure of the `Ito` Class

The Ito object stores the aforementioned three components:

Python code snippet defining a frozen dataclass named 'Ito' with three attributes: 'a', 'b', and 'c', where 'c' is a dictionary with integer keys and object values.

So for example,

dX = Ito.dt(mu * S) + Ito.dW(1, sigma * S)

would represent

$dS = \mu S dt + \sigma S dW_1$

Multiplication

The multiplication method __mul__ is important, because it obeys the Ito calculus multiplication rules. This means that when we multiply, say

$(a_1 + b_1 dt + c_1 dW)(a_2 + b_2 dt + c_2 dW)$

the class method will drop all $dt^2$ and $dt dW$ terms, and convert any $dW_i dW_j$ terms into $\delta_{ij} dt$ .

The quadratic variation term is handled explicitly,

A fragment of code displaying a Python loop that iterates over the intersection of two sets and modifies a variable 'b' by summing specific indexed elements from two lists 'c1' and 'c2'.

This line encodes $(dW_i)^2 = dt$ .

Everything else is standard algebra.

Why Does This Matter?

Because once multiplication is handled, quadratic variation should emerge automatically. So, for example,

			
dX = mu dt + sigma dW
dX2 = dX * dX

yields:

$dX^2 = \sigma^2 dt$

without having to write a single special case.

The `Ito Lemma` Layer

With the algebra defined, we implement the (1-dimensional) Ito lemma as a function. For the 1D diffusion $dX = \mu dt + \sigma dW$ , and a smooth function $f(t,x)$ , Ito’s lemma states:

$df = f_t dt + f_x dX + \displaystyle\frac{1}{2}f_{xx}(dX)^2$

where subscripts represent partial derivatives.

Since our Ito algebra already computes $(dX)^2$ , the implementation is direct.

Code snippet for implementing Ito's lemma in Python, demonstrating differentiation and diffusion processes using SageMath. — **Fig 6** – The 1-dimensional Ito Lemma function.

This produces an Ito object containing all the data needed to represent

$df = \text{drift}\cdot dt + \text{diffusion}\cdot dW$

Some Extensions

Multi-factor Extension

For scalar $X$ driven by multiple independent Brownian motions:

$dX = \mu dt + \sum_i \sigma_i dW_i$

Ito’s lemma becomes:

$df = f_t dt + f_x dX + \displaystyle\frac{1}{2} f_{xx}\sum_i \sigma_i^2 dt$

Our implementation allows this, for example, this line

Python code snippet implementing Ito's lemma for scalar state using multiple independent Brownian motions, detailing the differentiation and calculation of drift and output. — **Fig 7** – Implementation of a multi-factor extension of Ito’s lemma.

corresponds to computing

$(dX)^2 = \sum_i \sigma_i^2 dt$

Generator Extraction

We also implement the infinitesimal generator:

$\mathcal{L}[f](t,x) = f_t + \mu f_{x} + \displaystyle\frac{1}{2}\sigma^2 f_{xx}$

which is the drift term of $df$ under Ito’s lemma.

Code snippet illustrating a function definition for an infinitesimal generator in a 1D Itô diffusion process, including derivative calculations. — **Fig 8** – Implementation of the infinitesimal generator.

Matching the drift term against $\mathcal{L}f$ provides a clean consistency check.

Applying Ito’s Lemma Symbolically

Using our custom ito_lemma_1d function, we compute

$dV = V_t dt + V_SdS + \frac{1}{2}V_{SS}(dS)^2$

with the following command

Fig 9 – Applying Ito’s lemma symbolically in Python.

and the symbolic engine produces the drift term,

$\displaystyle\frac{1}{2}\sigma^2 S^2 V_{SS} + \mu S V_{S} + V_t$

and diffusion term,

$\displaystyle\sigma S V_S$

Mathematical expressions for drift and diffusion processes. — **Fig 10** – The drift and diffusion terms.

Conclusion

Our design separates algebra from calculus:

ito_algebra.py encodes the stochastic multiplication rules,
ito_lemma.py encodes the calculus layer,
SageMath handles the symbolic differentiation

and the Black-Scholes PDE drops out automatically!

References

Appendix

Ito_algebra.py

from dataclasses import dataclass
from typing import Dict
from sage.all import SR, var, diff  # type: ignore
@dataclass(frozen=True)
class Ito:
    a: object
    b: object
    c: Dict[int, object]
    def __add__(self, other):
        if not isinstance(other, Ito):
            other = Ito.scalar(other)
        c_new = dict(self.c)
        for i, ci in other.c.items():
            c_new[i] = c_new.get(i, 0) + ci
            if c_new[i] == 0:
                del c_new[i]
        return Ito(self.a + other.a, self.b + other.b, c_new)
    def __radd__(self, other):
        # allows: 0 + Ito(...) and scalar + Ito(...)
        return self.__add__(other)
    def __neg__(self):
        return Ito(-self.a, -self.b, {i: -ci for i, ci in self.c.items()})
    def __sub__(self, other):
        return self.__add__(-other)
    def __rsub__(self, other):
        return Ito.scalar(other).__sub__(self)
    @staticmethod
    def scalar(x):
        return Ito(SR(x), 0, {})
    @staticmethod
    def dt(coeff=1):
        return Ito(0, SR(coeff), {})
    @staticmethod
    def dW(i, coeff=1):
        return Ito(0, 0, {i: SR(coeff)})
    def __mul__(self, other):
        if not isinstance(other, Ito):
            other = Ito.scalar(other)
        a1, b1, c1 = self.a, self.b, self.c
        a2, b2, c2 = other.a, other.b, other.c
        a = a1 * a2
        c = {}
        for i in set(c1) | set(c2):
            ci = a1 * c2.get(i, 0) + a2 * c1.get(i, 0)
            if ci != 0:
                c[i] = ci
        b = a1 * b2 + a2 * b1
        for i in set(c1) & set(c2):
            b += c1[i] * c2[i]
        return Ito(a, b, c)
    def __repr__(self):
        parts = [str(self.a)]
        if self.b != 0:
            parts.append(f"{self.b}*dt")
        for i in self.c:
            parts.append(f"{self.c[i]}*dW{i}")
        return " + ".join(parts)

ito_lemma.py

from __future__ import annotations
from typing import Dict, Tuple
# noinspection PyUnresolvedReferences
from sage.all import var, SR, diff  # type: ignore
from Ito.ito_algebra import Ito
def _as_symbolic(x):
    """Coerce x into Sage's symbolic ring when possible."""
    try:
        return SR(x)
    except Exception:
        return x
def _extract_mu_sigmas(dX: Ito) -> Tuple[object, Dict[int, object]]:
    """
    Given dX = mu*dt + sum_i sigma_i*dW_i (Ito object), return (mu, sigmas_dict).
    We expect dX.a == 0 for a differential. If it's not, we don't hard-fail,
    but you should probably not be passing a nonzero scalar part as a "dX".
    """
    mu = dX.b
    sigmas = dict(dX.c)
    return mu, sigmas
def ito_lemma_1d(f, t, x, dX: Ito, dw_index: int = 1) -> Ito:
    """
    Itô's lemma for f(t, X_t) with 1D state X driven by ONE Brownian component dW_{dw_index}.
    Input:
        f        : Sage expression in (t, x)
        t, x     : Sage symbols
        dX       : Ito object representing dX = mu*dt + sigma*dW_{dw_index}
        dw_index : which Brownian index to use as the 1D driver
    Output:
        df as Ito object: df = drift*dt + diffusion*dW_{dw_index}
    """
    ft = diff(f, t)
    fx = diff(f, x)
    fxx = diff(fx, x)
    mu, sigmas = _extract_mu_sigmas(dX)
    sigma = sigmas.get(dw_index, 0)
    drift = ft + fx * mu + SR(1) / SR(2) * fxx * (sigma ** 2)
    diffusion = fx * sigma
    return Ito.dt(drift) + Ito.dW(dw_index, diffusion)
def ito_lemma_scalar_state_multiW(f, t, x, dX: Ito) -> Ito:
    """
    Itô's lemma for scalar state X driven by multiple independent Brownian motions:
        dX = mu*dt + sum_i sigma_i dW_i
    where dW_i dW_j = delta_ij dt.
    Then:
        df = f_t dt + f_x dX + 1/2 f_xx (sum_i sigma_i^2) dt
           = [f_t + f_x mu + 1/2 f_xx * sum_i sigma_i^2] dt + sum_i f_x sigma_i dW_i
    """
    ft = diff(f, t)
    fx = diff(f, x)
    fxx = diff(fx, x)
    mu, sigmas = _extract_mu_sigmas(dX)
    sig_sq = sum((si ** 2) for si in sigmas.values())
    drift = ft + fx * mu + SR(1) / SR(2) * fxx * sig_sq
    out = Ito.dt(drift)
    for i, si in sigmas.items():
        out = out + Ito.dW(i, fx * si)
    return out
def generator_L_1d(f, t, x, mu, sigma):
    """
    Infinitesimal generator for a 1D Itô diffusion:
        dX = mu(t,x) dt + sigma(t,x) dW
    L f = f_t + mu f_x + 1/2 sigma^2 f_xx
    Useful for Feynman–Kac sanity checks.
    """
    ft = diff(f, t)
    fx = diff(f, x)
    fxx = diff(fx, x)
    return ft + mu * fx + SR(1) / SR(2) * (sigma ** 2) * fxx

test_ito_lemma.py

import unittest
# noinspection PyUnresolvedReferences
from sage.all import var, SR, sin  # type: ignore
from Ito.ito_algebra import Ito
from Ito.ito_lemma import ito_lemma_1d, ito_lemma_scalar_state_multiW, generator_L_1d
def assert_sage_expr_zero(testcase: unittest.TestCase, expr, msg: str = ""):
    """
    Robust "expr == 0" for Sage symbolic expressions.
    """
    e = SR(expr).simplify_full()
    testcase.assertTrue(e == 0, msg or f"Expected 0, got {e}")
def assert_sage_expr_equal(testcase: unittest.TestCase, lhs, rhs, msg: str = ""):
    """
    Robust symbolic equality check: lhs == rhs.
    """
    assert_sage_expr_zero(testcase, SR(lhs) - SR(rhs), msg or f"{lhs} != {rhs}")
class TestItoLemma(unittest.TestCase):
    def test_ito_lemma_x_squared_1d(self):
        # f(x) = x^2, dX = mu dt + sigma dW
        t, x, mu, sigma = var("t x mu sigma")
        f = x**2
        dX = Ito.dt(mu) + Ito.dW(1, sigma)
        df = ito_lemma_1d(f, t, x, dX, dw_index=1)
        # Expected: df = (2x*mu + sigma^2) dt + (2x*sigma) dW
        expected_drift = 2 * x * mu + sigma**2
        expected_diff = 2 * x * sigma
        assert_sage_expr_equal(self, df.a, 0, "df should have no scalar term")
        assert_sage_expr_equal(self, df.b, expected_drift, "wrong drift for x^2")
        self.assertIn(1, df.c, "missing dW1 term")
        assert_sage_expr_equal(self, df.c[1], expected_diff, "wrong diffusion for x^2")
    def test_ito_lemma_constant_function_is_zero(self):
        t, x, mu, sigma = var("t x mu sigma")
        f = SR(7)  # constant
        dX = Ito.dt(mu) + Ito.dW(1, sigma)
        df = ito_lemma_1d(f, t, x, dX, dw_index=1)
        assert_sage_expr_equal(self, df.a, 0)
        assert_sage_expr_equal(self, df.b, 0)
        self.assertEqual(df.c, {})
    def test_ito_lemma_linear_function(self):
        # f(t,x) = a*x + b*t
        t, x, a, b, mu, sigma = var("t x a b mu sigma")
        f = a * x + b * t
        dX = Ito.dt(mu) + Ito.dW(1, sigma)
        df = ito_lemma_1d(f, t, x, dX)
        # f_t = b, f_x = a, f_xx = 0
        expected_drift = b + a * mu
        expected_diff = a * sigma
        assert_sage_expr_equal(self, df.b, expected_drift)
        assert_sage_expr_equal(self, df.c[1], expected_diff)
    def test_ito_lemma_multiW_scalar_state(self):
        # Scalar X driven by 2 independent Brownian motions
        t, x, mu, s1, s2 = var("t x mu s1 s2")
        f = x**3 + t * x
        dX = Ito.dt(mu) + Ito.dW(1, s1) + Ito.dW(2, s2)
        df = ito_lemma_scalar_state_multiW(f, t, x, dX)
        # Compute expected manually:
        # f_t = x
        # f_x = 3x^2 + t
        # f_xx = 6x
        fx = 3 * x**2 + t
        ft = x
        fxx = 6 * x
        expected_drift = ft + fx * mu + SR(1) / SR(2) * fxx * (s1**2 + s2**2)
        expected_dW1 = fx * s1
        expected_dW2 = fx * s2
        assert_sage_expr_equal(self, df.a, 0)
        assert_sage_expr_equal(self, df.b, expected_drift)
        self.assertIn(1, df.c)
        self.assertIn(2, df.c)
        assert_sage_expr_equal(self, df.c[1], expected_dW1)
        assert_sage_expr_equal(self, df.c[2], expected_dW2)
    def test_generator_matches_drift_part_1d(self):
        # For 1D: drift from Ito lemma should equal Lf where L is generator,
        # for f(t,x) with dX = mu dt + sigma dW, excluding diffusion terms.
        t, x, mu, sigma = var("t x mu sigma")
        f = sin(x) + t**2
        dX = Ito.dt(mu) + Ito.dW(1, sigma)
        df = ito_lemma_1d(f, t, x, dX)
        Lf = generator_L_1d(f, t, x, mu, sigma)
        assert_sage_expr_equal(self, df.b, Lf, "Ito drift should equal generator Lf")
    def test_ito_lemma_multiW_reduces_to_1d_when_one_factor(self):
        # multiW formula should match 1D formula if only one dW term exists
        t, x, mu, sigma = var("t x mu sigma")
        f = x**4
        dX = Ito.dt(mu) + Ito.dW(1, sigma)
        df1 = ito_lemma_1d(f, t, x, dX)
        dfm = ito_lemma_scalar_state_multiW(f, t, x, dX)
        assert_sage_expr_equal(self, df1.b, dfm.b)
        self.assertEqual(set(df1.c.keys()), set(dfm.c.keys()))
        assert_sage_expr_equal(self, df1.c[1], dfm.c[1])
if __name__ == "__main__":
    unittest.main()

demo_ito_calcs_1.py

# noinspection PyUnresolvedReferences
from sage.all import var, sin, exp, SR
from Ito.ito_algebra import Ito  # type: ignore
from Ito.ito_lemma import ito_lemma_1d, ito_lemma_scalar_state_multiW, generator_L_1d
def show(label, ito_obj: Ito):
    print("\n" + label)
    print("  scalar a =", ito_obj.a)
    print("  dt     b =", ito_obj.b)
    print("  dW     c =", ito_obj.c)
def main():
    t, x, mu, sigma = var("t x mu sigma")
    # 1) Basic differential and quadratic variation
    dX = Ito.dt(mu) + Ito.dW(1, sigma)
    show("dX = mu dt + sigma dW1", dX)
    dX2 = dX * dX
    show("dX^2 (should be sigma^2 dt)", dX2)
if __name__ == "__main__":
    main()

demo_ito_calcs.py

			
# noinspection PyUnresolvedReferences
from sage.all import var, function  # type: ignore
from Ito.ito_algebra import Ito
from Ito.ito_lemma import ito_lemma_1d
def main():
    t, S, mu, sigma = var("t S mu sigma")
    V = function("V")(t, S)
    dS = Ito.dt(mu * S) + Ito.dW(1, sigma * S)
    dV = ito_lemma_1d(V, t, S, dS)
    print("Drift: ")
    print(dV.b.simplify_full())
    print("\nDiffusion: ")
    print(dV.c[1].simplify_full())
if __name__ == "__main__":
    main()

		

	Muhammad Shumail on Hadamard & CNOT
	spiral on A Recipe for Pricing
	Deutsch-Jozsa Algori… on Quantum Oracles
	Raymond Glowner on Fractal Manifolds III –…
	Parvez on A Recipe for Pricing

From First Principles

Mathematics Without Hand Waving

Deriving Black-Scholes with SageMath

Introduction

Setting Up Geometric Brownian Motion in Python

Encoding the Ito Algebra

The Structure of the `Ito` Class

Multiplication

Why Does This Matter?

The `Ito Lemma` Layer

Some Extensions

Multi-factor Extension

Generator Extraction

Applying Ito’s Lemma Symbolically

Conclusion

References

Appendix

Ito_algebra.py

ito_lemma.py

test_ito_lemma.py

demo_ito_calcs_1.py

demo_ito_calcs.py

Leave a comment Cancel reply

Introduction

Setting Up Geometric Brownian Motion in Python

Encoding the Ito Algebra

The Structure of the Ito Class

Multiplication

Why Does This Matter?

The Ito Lemma Layer

Some Extensions

Multi-factor Extension

Generator Extraction

Applying Ito’s Lemma Symbolically

Conclusion

References

Appendix

Ito_algebra.py

ito_lemma.py

test_ito_lemma.py

demo_ito_calcs_1.py

demo_ito_calcs.py

Share this:

Leave a comment Cancel reply

The Structure of the `Ito` Class

The `Ito Lemma` Layer