Introduction

In the previous post, we introduced the Feynman-Kac equation by starting with a space of functions and equipping it with the Markov semigroup. From there, the infinitesimal generator and its adjoint gave us two (dual) ways to describe the same evolution of random information: we could evolve a function backward from a future observation (terminal boundary condition), or evolve some initial distribution forward. Then Feynman-Kac appeared where those two descriptions met.

This time, we will make it concrete in a SageMath Jupyter Notebook.

If the expectation is difficult but the PDE is tractable, then solve the PDE (usually numerically).
If the PDE is difficult, simulate paths of X and estimate the expectation by Monte Carlo.

We will use geometric Brownian motion and a simple function f(x) = x^2 so that we can follow every step without losing the main idea in all the algebra. First, we will use the generator that we built in Python to verify the backward Kolmogorov equation for the expected future value of f. Then we will use the adjoint to verify the forward Kolmogorov equation for the transition density. Finally, we will bring the two sides together: by integrating f against that density, we will recover exactly the function that solves the backward Kolmogorov equation.

The aim here is to make this connection visible and practical, using only Python objects that we have already built. By the end, the backward PDE solution, the forward probability distribution, and the conditional expectation should really feel like three ways of looking at the same object.

Recall from the previous blog we ended with this picture:

Flowchart illustrating the relationship between observables and distributions in stochastic processes, featuring concepts like Markov semigroups, generators, and Kolmogorov equations.

Notebook Imports

Let’s first define these imports and helper functions.

# noinspection PyUnresolvedReferences
from sage.all import var, SR, exp, sqrt, latex, log, pi, assume
from Sde.sde import SDE1D   # type: ignore
from Ito.ito_lemma import ito_lemma_1d
from Ito.kolmogorov import backward_kolmogorov
from IPython.display import display, Math
def show_eq(lhs, rhs):
    display(Math(rf"{lhs} = {latex(rhs)}"))
def show_text_math(text, expr):
    display(Math(rf"\text{{{text}}}\quad {latex(expr)}"))

Define the Markov Diffusion

We begin with geometric Brownian motion

\displaystyle dX_t = \mu X_t dt + \sigma X_t dW_t

Its future evolution depends on its current value, rather than the path it took to get there. This makes it a Markov diffusion. To see why, suppose we know that X_t=x. Over the next small interval, its drift is \mu x and its diffusion coefficient is \sigma x. For geometric Brownian motion, the same point holds over any future interval, in other words, once we know X_t, the earlier path gives us no further information about the distribution of X_T at the terminal boundary. The current state is enough to determine the probabilities of future states.

Recall from the previous post, that this lets us define a Markov operator P_s which takes a function f and returns its expected value after an elapsed time s:

\displaystyle P_s[f](x) = \mathbb{E}\left[f(X_{t+s})\mid X_t=x\right]

The infinitesimal generator \mathcal L measures how quickly that expectation changes as s starts to increase from zero. For our geometric Brownian motion, it is

\displaystyle\mathcal{L}[f](x) = \mu x f_x(x) + \frac12 \sigma^2 x^2 f_{xx}(x)

We represent the diffusion using our SDE1D class that we built earlier in this series:

t, T, x, mu, sigma = var("t T x mu sigma")
gbm = SDE1D(
    t=t,
    x=x,
    drift=mu*x,
    diffusion=sigma*x,
    name="GBM"
)

Now we choose the function we will follow for the rest of the notebook, f(x) = x^2. We can ask the SDE1D object to apply its spatial generator directly:

f = x**2
Lf = gbm.generator(f)
show_eq(r"f(x)", f)
show_eq(r"\mathcal{L}[f]", Lf)

The result is

\displaystyle\mathcal{L}[f](x)=(2\mu+\sigma^2)x^2

The generator has returned the instantaneous evolution of this observable function. The 2\mu term comes from the drift, while the extra \sigma^2 term comes from the diffusion.

Next, we will turn this instantaneous description into a function that tells us the expected value of X_T^2 at a future time T.

Code snippet defining a geometric Brownian motion simulation in a programming environment, showcasing the application of the SDE1D object with a test function.

Compute the Conditional Expectation

Now that we have seen what the generator does to f(x)=x^2 over an infinitesimal time interval, let us now calculate the expected value after a finite time interval [t,T].

Luckily for us, geometric Brownian motion has an exact transition formula. Conditional on X_t = x, we have

\displaystyle X_T=x\exp\left[\left(\mu-\frac12\sigma^2\right)(T-t)+\sigma\sqrt{T-t}Z\right]

where Z \sim \mathcal{N}(0,1).

The random variable Z accounts for the Brownian motion accumulated between times t and T. Squaring this expression and taking its expectation gives:

\displaystyle\mathbb{E}\left[X_T^2\mid X_t=x\right] = x^2 e^{(2\mu - \sigma^2)(T-t)}\cdot\mathbb{E}\left[e^{2\sigma\sqrt{T-t}Z}\right]

which simplifies to

\displaystyle\mathbb{E}\left[X_T^2\mid X_t=x\right] = x^2 e^{(2\mu + \sigma^2)(T-t)}

where in this last line, we used the standard Normal identity \mathbb{E}[e^{aZ}] = e^{a^2 / 2}.

Let us call this conditional expectation u(t,x):

tau = T - t
u = x**2 * exp((2*mu + sigma**2) * tau)
show_eq(r"u(t,x)", u)

The variable tau (or sometimes I use s) is the time remaining until T. This matters because u tells us, at time t, what we expect our future observation f(X_T) = X_T^2 to be.

Before checking any PDE, we should first check what happens at the terminal time. When t = T, there is no time left for X to evolve. We therefore know precisely that X_T = x, and therefore that u(T,x) = f(x) = x^2. Hence,

terminal_residual = (u.subs(t=T) - f).simplify_full()
show_eq(r"u(T,x)", u.subs(t=T))
show_eq(r"u(T,x) - f(x)", terminal_residual)
assert terminal_residual == 0
A code snippet in a programming environment, showing calculations involving a variable 'terminal_residual', assertions, and equations related to 'u(T,x)' and 'f(x)'.

SageMath returns zero for the residual!

We have now constructed u from a conditional expectation and confirmed its terminal condition.

The next question is whether this same function also solves the backward Kolmogorov PDE.

Verify the Backward Kolmogorov PDE

We obtained u(t,x) by taking a conditional expectation. Now we will approach the same function from the PDE side.

For our geometric Brownian motion, we showed in a previous blog that the backward Kolmogorov equation is

\displaystyle u_t + \mu x u_x + \frac12\sigma^2 x^2 u_{xx} = 0

with terminal condition u(T,x) = x^2.

The spatial terms are precisely \mathcal{L}[u](x), so we can write the equation more compactly as,

\displaystyle u_t + \mathcal{L}u = 0

Why the plus sign when we are working backward? The terminal time T is fixed. As the current time t moves forward, the remaining time \tau = T-t gets shorter. The time derivative u_t therefore offsets the evolution described by the generator.

Our backward_kolmogorov() method now computes the left-hand side of the above PDE. If the conditional expectation we found really solves it, then the result should simplify to zero:

pde_residual = backward_kolmogorov(
    u=u,
    t=t,
    state_vars=x,
    drift=mu*x,
    diffusion=sigma*x,
    simplify_result=True
)
show_eq(r"u_t+\mathcal{L}[u]", pde_residual)
assert pde_residual == 0
Code snippet showing the implementation of a backward Kolmogorov equation in Python, including the definition of the pde_residual, the show_eq function, and an assertion check to verify its value.

And it does!

We can check the same calculation through our SDE1D object, whose Ito operator applies (\partial_t + \mathcal{L}) to a function:

sde_operator_residual = gbm.ito_operator(u)
show_eq(
    r"(\partial_t+\mathcal{L})[u]",
    sde_operator_residual
)
assert sde_operator_residual == 0
A code snippet showing the verification of an SDE operator in Python, with an expression involving partial derivatives and an assertion statement.

And as we can see, both routes give zero!

Thus, together with the terminal condition, this confirms that the function we calculated as \mathbb{E}[X_T^2 \mid X_t = x] solves the backward PDE generated by the GBM diffusion. This is the expectation-to-PDE side of our Feynman-Kac example all done.

Before turning to the dual, forward picture, let us see why the backward PDE produces a conditional expectation.

The Martingale Step: Applying Ito’s Lemma to the PDE Solution

We have checked that u satisfies the backward PDE, but there is still the question: why does the PDE lead to a conditional expectation?

To see this, let’s apply Ito’s lemma to u(t,X_t), where X_t follows our GBM:

\displaystyle\left(u_t + \mathcal{L}[u](x)\right)dt + \sigma X_t u_x dW_t

The expression multiplying dt is exactly the left-hand side of the backward Kolmogorov equation. We have just shown that it is zero, hence,

\displaystyle du(t,X_t) = \sigma X_t u_x dW_t

We can see this cancellation using our Ito layer built earlier in this series:

du = ito_lemma_1d(
    f=u,
    t=t,
    x=x,
    dX=gbm.dX(),
    dw_index=gbm.dw_index
)
print("Itô differential of u(t, X_t):")
print(du)
show_eq(r"\text{drift}[du]", du.b)
show_eq(r"\text{diffusion}[du]", du.c[gbm.dw_index])
assert du.b.simplify_full() == 0

And the drift simplifies to zero as expected!

The remaining diffusion coefficient is

\displaystyle 2\sigma x^2 e^{(2\mu+\sigma^2)(T-t)}

The zero drift tells us that u(r, X_r), for (t \leq r \leq T), is a local martingale. To use its conditional expectation property, we also need it to be a true martingale. In this GBM example, u(r, X_r) is proportional to X_r^2, and GBM has a finite fourth moment throughout the finite interval [t,T], thus it is square-integrable and we can write:

\displaystyle u(t,x)=\mathbb{E}\left[u(T, X_T) \mid X_t=x \right]

At the terminal time, u(T,y) = f(y) = y^2. Hence

\displaystyle u(t,x)=\mathbb{E}[X_T^2\mid X_t=x]

This gives us the mechanism behind the expectation-to-PDE connection: the backward equation cancels the drift in Ito’s lemma, and the martingale property carries the terminal value back to time t.

We have now completed the observable side of the picture. Next, we will turn to its dual: evolve the transition density forward using the adjoint \mathcal L^*, then integrate f against that density to recover the same u(t,x).

Following the Dual Side

So far we have followed the observable f(y) = y^2 backward from its terminal value.

This time, we keep the starting state X_t = x fixed and follow the distribution of future states forward.

Let s>0 be the time elapsed since t, and let y>0 denote a possible value of X_{t+s}. The GBM transition formula tells us that

\displaystyle\log\displaystyle\frac{X_{t+s}}{x} \sim \mathcal{N}\left(\left(\mu - \frac12\sigma^2\right)s, \sigma^2 s\right)

Consequently, the transition density is also well-known:

\displaystyle p(s,y\mid x) = \frac{1}{\sigma\sqrt{2\pi s}y}\exp\left(-\frac{\left(\log(y/x) - (\mu - \frac12\sigma^2)s\right)^2}{2\sigma^2 s}\right)

Here x is a parameter specifying where the process started and y is the state variable in the forward equation.

Now let’s create this density in SageMath:

s, y = var("s y")
assume(x > 0, y > 0, s > 0, sigma > 0)
log_return = log(y/x) - (mu - sigma**2/2)*s
p = exp(-log_return**2/(2*sigma**2*s)) / (
    y*sigma*sqrt(2*pi*s)
)
show_eq(r"p(s,y\mid x)", p)

The above generator \mathcal{L} we know acts on functions, and its formal adjoint $\mathcal{L}^*$ acts on densities. For GBM, we already know that the forward Kolmogorov equation is

\displaystyle\mathcal{L}^*[p](x) = -\frac{\partial}{\partial y}(\mu y p) + \frac12\frac{\partial^2}{\partial y^2} (\sigma^2y^2p)

Then, our existing forward_kolmogorov() functions computes the residual $p_s – \mathcal{L}^*[p]$:

from Ito.kolmogorov import forward_kolmogorov
forward_residual = forward_kolmogorov(
    p=p,
    t=s,
    state_vars=y,
    drift=mu*y,
    diffusion=sigma*y,
    simplify_result=True
).simplify_full()
show_eq(r"p_s-\mathcal{L}^{*}p", forward_residual)
assert forward_residual == 0
Code snippet demonstrating the import and use of the forward_kolmogorov function from the Ito.kolmogorov module. It shows the computation of the forward residual and an assertion that it equals zero.

The above zero residual confirms that the density evolves forward according to the adjoint of the same generator we used on the backward side. Its initial condition is understood as a limit, i.e., as s decreases to zero, the distribution concentrates at y=x, becoming a point mass there.

Bridging the Gap with the Handshake

We can now do the calculation promised at the start of the post. At elapsed time s=T-t, average our observable f(y) = y^2 against the density:

\displaystyle \int_0^{\infty} y^2 p(s, y\mid x)dy

To evaluate the integral, use the change of variable

\displaystyle z = \frac{\log(y/x) - (\mu - \frac12\sigma^2)s}{\sigma\sqrt{s}}

Under this substitution, p(s, y\mid x)dy becomes the standard normal density \phi(z) dz, while

\displaystyle y^2 = x^2 \exp\left[2\left(\mu-\frac12\sigma^2\right)s + 2\sigma\sqrt{s}z\right]

The integral is therefore a normal exponential moment. Using the known identity: \int_{-\infty}^{\infty} e^{az}\phi(z)dz = e^{a^2/2}, we obtain

\displaystyle \int_{0}^{\infty} y^2p(s, y\mid x)dy = x^2 e^{(2\mu-\sigma^2)s}

We can check that this result agrees with the conditional expectation constructed earlier:

density_moment = (
    x**2
    * exp(2*(mu - sigma**2/2)*s)
    * exp((2*sigma*sqrt(s))**2/2)
).simplify_full()
backward_value = u.subs(t=T-s)
show_eq(r"\int_0^\infty y^2p(s,y\mid x)\,dy", density_moment)
show_eq(r"u(T-s,x)", backward_value)
handshake_residual = (density_moment - backward_value).simplify_full()
show_eq(r"\text{difference}", handshake_residual)
assert handshake_residual == 0
Code snippet showing calculations for density moment, backward value, and handshake residual in a mathematical context.

There is our long-awaited handshake! The difference is zero.

SageMath has confirmed that our analytically evaluated density integral equals the solution of the backward Kolmogorov equation (PDE). The forward equation evolves the density forward from the starting state, whilst the backward equation evolves the observable function backward from its terminal value. Pairing the density with that observable brings the two sides together:

\displaystyle u(t,x) = \int_{0}^{\infty} y^2 p(T-t,y\mid x)dy = \mathbb{E}\left[X_T^2 \mid X_t = x\right] = x^2 e^{(2\mu + \sigma^2)(T-t)}

Conclusion

We started with a question about a future observation: if geometric Brownian motion is at x today, what is the expected value of X_T^2 at the terminal boundary? The exact GBM transition gave us a (solution) function u(t,x), and our SageMath calculations showed that this function satisfies the backward Kolmogorov equation and the terminal condition u(T,x)=x^2.

In the middle, we used Ito’s lemma and explained the link between those calculations. When we evaluate u along a path of the diffusion, the drift coefficient of u(r,X_r) is just the differential operator u_r + \mathcal{L}u. We verified that this expression is zero. In this example u(r,X_r) is a true martingale, so its value at time t is the conditional expectation of its value at the later time T. The terminal condition then turns that future value into X_T^2.

Finally, we also reached the answer from the other direction. Starting at x, the (already known) lognormal transition d0ensity evolves forward according to the adjoint generator \mathcal{L}^*. Averaging y^2 against that density gives exactly the same solution u(t,x)!. The backward PDE follows the observable; the forward PDE follows the probabilities. Their pairing is the conditional expectation we set out to calculate.

That is the Feynman–Kac handshake made concrete in this simple example.

backward_value = u.subs(t=T-s)
handshake_residual = (density_moment - backward_value).simplify_full()
assert handshake_residual == 0

Our Ito and Sde Python Layer

The SageMath calculations are possible because of our Ito and Sde Python layers we built throughout this series. The class object SDE1D holds the GBM coefficients and applies its generator; backward_kolmogorov() and forward_kolmogorov() turn the generator and its adjoint into PDE residuals; and ito_lemma_1d() exposes the drift and diffusion terms when we need to evaluate u(t,X_t) along the process. SageMath then simplifies each expression so we can check the identities. The layers let us carry the same diffusion through the probabilistic and PDE calculations without redefining its dynamics at every step.

Why is it Useful to Have the Dual Representation?

In the end, we often have these types of equations: conditional expectations:

u(t,x) = \mathbb{E}[g(X_t)\mid X_t = x]

and this is often very hard to calculate. But now we know that this is only one way to look at the solution. Feynman-Kac gives us an alternative route: solve the backward PDE with terminal boundary condition u(T,x) = g(x).

But it’s useful in both directions.

If the expectation is difficult but the PDE is tractable, then solve the PDE (usually numerically).
If the PDE is difficult, simulate paths of X and estimate the expectation by Monte Carlo.

In either case, the result is a function u(t,x) that you can evaluate across current states and times. This GBM example was unusually nice and well behaved because both routes have a closed form solution, which makes it super easy to evaluate and use the zero residual metric to show both sides are equivalent.