TL;DR

While exploring the properties of various continuous-time stochastic processes, I had a bit of fun simulating them in python using discrete approximations generated via the Euler-Maruyama method.

In order to check my work when manipulating various stochastic random variables using Itô calculus, I built a Martingale Detector - a crude tool which tests a drift hypothesis via Monte Carlo methods.

When the naive detection scheme yielded unexpected results for certain processes, I discovered that I had implicitly baked in a faulty assumption in my code.

I believe that the bug described below is an instructive example of the consequences of eschewing rigor when designing simulations. This essay is equal parts mathematical candy and postmortem.

While my use for the tool was rather pedestrian (checking my work on some practice problems), testing for signals like Martingale and Markov properties in stochastic time-series is actually of great interest in fields like finance and physics.

A brief review of Martingales

Assume a filtered probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{t \ge 0}, \mathbb{P})$ is given. An $\mathcal{F}_t$-adapted process $Y_t$ is a continuous-time martingale if and only if $\mathbb{E}\left[|Y_t|\right] < \infty$ and

Proving the martingale conditions for simple stochastic processes which can be expressed in closed-form often requires little more than some massaging using the tools of Itô calculus and basic probability theory. For example, for processes $F_t = F(W_t, t)_t$ which can be expressed as functions of a standard Wiener process $W_t$ and time $t$, it is sufficient to prove

Still, I wanted a quick-and-dirty way to check my work, preferrably with a method that makes few assumptions about the underlying dynamics.

The Naive Martingale Detector

The problem domain

I didn’t bother formalizing requirements for the types of processes which the tool should handle. The general form of stochastic differential equations describes far too broad a domain of processes and admits many pathological examples which could evade any numerical drift detector I throw at it ¹.

I essentially restricted $\mu(t, X_t)$ and $\sigma(t, X_t)$ to analytic functions or Itô integrals of analytic functions. I may revisit this topic in the future to further generalize the revised tool presented at the end of this essay.

While fancy non-parametric techniques like bootstrapping might be a more mathematically sound tool for the job, Monte Carlo will essentially give me a population standard deviation estimate ($\hat{\sigma} = s$) to any degree of precision my CPU can stand, thus making a $z$-test a reasonable choice. CLT, take the wheel!

A better name for this tool might be the Drift Detector ², since we will be evaluating a Martingale null hypothesis.

Aiming for simplicity, I ran a $z$-test for each timestep $t = i \, \delta t$ and rejected the null upon discovery of any breaches.

Detection procedure

For the following, let $t_i \equiv i \, \delta t$. ³

My naive (indeed faulty!) scheme for drift detection is essentially a series of hypothesis tests driven by Monte Carlo simulation. The details:

1. Sample $n \times T$ stochastic increments $dW_t \equiv \delta W_{t_i}; \ i \in \{0, \dots, T - 1\}$ from a Gaussian⁴; i.e. $\delta W_{t_i} \sim \mathcal{N}(0, \sqrt{\delta t})$.

2. Perform discrete integration of the $dW_t$ to obtain $W_t$, ensuring that $W_0 = 0$.

3. Shift the $dW_t$ to the left such that any vectorized operations on $W_t$ and $dW_t$ are properly aligned (I may describe this step in further detail later, but it is largely unimportant for this discussion).

4. Construct $Y_t = f(W_t, t)$ via Euler-Maruyama. At this point, path realizations of $Y_t$ can be rendered to get a sense of their dynamics.

5. Drift detection step: Compute $z$-statistics for each timestep ⁵ and test the null.

   $$z_t = \frac{\bar{Y}_t - Y_0}{\hat{\sigma}_t / \sqrt{n}} \approx \frac{\frac{1}{n} \sum_{i=1}^n Y_t^{(i)} - Y_0}{s_t / \sqrt{n}}$$

Below are some examples of the detector working as designed.

Notice the growing critical region boundaries (in black) as the processes diffuse.

For good measure, I included the Wiener process, which passes unsurprisingly. I also cooked up some more exotic test cases, like processes which are locally driftless at the boundaries of the interval. This test will catch these types of processes, up to the resolution of the discretization.

Alas, it was too good to be true!

I quickly noticed the detector making Type 2 errors for certain processes, some of which are shown below.

In order to diagnose the problem, we need to analyze the two non-Martingales shown above ⁶.

1. $X_t = e^{\alpha t} \cos W_t$ (correctly identified)

2. $Y_t = (W_t - t) e^{W_t + k t}$ (false negative)

I set out to build a tool to sanity-check my math; now we must rely on the math to debug the tool!

Non-Martingale examples

Non-Martingale 1: correctly identified

Consider the process

Then $\mathbb{E}\left[X_t | \mathcal{F}_s\right] = e^{\alpha t} \mathbb{E}\left[\cos W_t | \mathcal{F}_s\right]$. To proceed further, define $f(W_t, t) = e^{t/2} \cos W_t$. By Itô’s lemma,

This stochastic differential equation is equivalent to the more flexible integral form:

Rearranging,

Finally, have

So $X_t$ is not a martingale for all $\alpha \ne 1/2$.

Non-Martingale 2: false negative

Consider the following process:

The conditional expectation is relatively easy to calculate. There may be more elegant methods, but I computed the expectation by brute force, taking advantage of the normality of the quantity $W_t - W_s$.

There are then just two functions of Gaussian quantities for which the expectation must be found. Let $x \equiv W_t - Ws$

It is trivial to show that

Finally,

So $Y_t$ is not a martingale for all $k \ne -1/2$.

This bears a striking resemblance to the result for $\mathbb{E}\left[X_t | \mathcal{F}_s\right]$. So why does the martingale detection method work for $X_t = e^{\alpha t} \cos W_t$, but not for $Y_t = (W_t - t) e^{W_t + k t}$?

So, what went wrong?

The Monte Carlo techniques used by the Naive Martingale Detector work by path averaging, i.e.

where convergence in probability is guaranteed by the Weak Law of Large Numbers.

Spot the mistake? $\bar{Y}_t$ converges on the unconditional expectation of $Y_t$, whereas the martingale property is defined by the expectation of $Y_t$ conditional upon the filtrations $\{\mathcal{F}_s\}_{0 \le s \le t}$!

For stochastic processes whose expectation conditioned on $\mathcal{F}_0$ is always some static value, the Naive Martingale Detector will "report" a false positive!

The originally proposed null hypothesis is invalid because it implies a conditional expectation only on $\mathcal{F}_0$ and no other values. It is not equivalent to the Martingale drift property.

Hopefully it’s clear now why the Naive Martingale Detector failed or succeeded for the processes discussed earlier.

• The process $X_t = e^{\alpha t} \cos W_t$ has unconditional expectation $\mathbb{E}\left[X_t\right] = e^{\left(\alpha - \frac{1}{2}\right) t}$.
  NMD correctly identified the drift.

• The process $Y_t = (W_t - t) e^{W_t + k t}$ has unconditional expectation $\mathbb{E}\left[Y_t\right] = 0 \; \forall t$.
  NMD false negative.

• The process $Z_t = \frac{W_t^2}{t}$ has unconditional expectation $\mathbb{E}\left[Z_t\right] = 1 \; \forall t$.
  NMD false negative.

Fixing the Naive Martingale Detector

There are many ways to fix the Naive Martingale Detector. The first is my solution and admittedly a bit hacky; the second is the test used in most papers on this topic.

Solution 1: Mind the filtration

This solution is a bit crude but allows us to keep most of the code from the Naive Martingale Detector.

Modify the hypothesis:

for some $0 \le s \le T$ and some $x \in \mathbb{R}$.

The conditional probability is the key here; there may be many expressions which work. I found a simple threshold filter to yield good results.

The best results are obtained by setting test parameters $s$ and $x$ manually for each process; luckily, this is not hard to do. A filter at time $s = \left\lfloor \frac{T}{4 \delta t} \right\rfloor$ (about 25% along the $t$ axis) works for all processes, but doesn’t always lead to good visuals. At first, I somewhat arbitrarily used $x = 2 \hat{\sigma}(s)$ but switched to filtering the upper quartile to accommodate processes with highly clustered distributions.

Here, we revisit the false negatives from above, this time armed with the new detector.

Success! Now, we must ensure there is no regression for the first examples in this essay.

Unfortunately, we have indeed regressed here, erroneously rejecting the null hypothesis for the process $Y_t = e^{t/2} \cos W_t$. On the bright side, it’s not due to a flaw in logic, but it does expose a weakness of this approach. I believe sampling error is the culprit in this instance, since we filter out 75% of the samples before testing. It’s trivial to compensate for this with a larger initial sample size.

Solution 2: Autoregressive model

Up to this point, every solution discussed (valid or otherwise) just throws CPU power at the problem quite brutishly. This approach is often tempting (and incredibly powerful), especially when time is of the essence.

But I would be doing the reader a great disservice if I did not mention the right way to solve this problem. This approach is not my own; full credit goes to the many people who have contributed to its development.

Instead of analyzing cross-sections of samples, model the sequence and test for the strength of the drift term in an AR(p) model. For AR(1), model

To make this useful, assume a unit root ($\theta$ = 1) so that

Thus, we have our Martingale null hypothesis:

Clever indeed! For more, see this paper from Phillips and Jin: Phillips and Jin (2014)

Final thoughts

I’d like to document some of my takeaways here, written as advice to my future self.

Just because it works in practice, doesn’t mean it works in theory!

Wikipedia lists the 4 defining properties of the Wiener process $W_t$. I designed the path generation procedure described in this essay to preserve properties 1, 3, and 4 (obviously sacrificing continuity). The path generation itself worked well.

However, the mere existence of a functioning discretization mechanism should not lead you to assume that some continuous analog of that mechanism even exists, much less is the fundamental engine of the object of study!

For a while, I operated with the mental model that a realization of the Wiener process was constructed by some mysterious mathematical supercomputer running my same algorithm but with an uncountably infinite array of $dW_t$. But of course, this is the wrong mathematical framing! Instead, entire paths (on some interval) are indexed by samples $\omega \in \Omega$, and we observe the properties of these objects through the $t$-indexed lenses provided by a filtration (a proxy for real-world ignorance) on the underlying probability space.

Notice that $dW_t$ appears nowhere in the 4 defining properties of the Wiener process. $dW_t$ really comes from Itô calculus ⁷, when Itô generalized the Riemann-Stieltjes integral to stochastic settings in the 1940s, two decades after Wiener’s formalization.

More lessons: lightning round

1. There is a place for rigor when designing simulations, and it can be rewarding to nail down seemingly abstract mathematical objects.

2. When designing and debugging simulations, be explicit about your assumptions; own and challenge them. Understand what you lose when discretizing.

3. Use caution when switching between mathematical contexts. Brownian motion and its cousins are interesting objects of study because they are fundamentally random and fundamentally continuous.

   Simulations always make compromises. But simulations like this have additional (sometimes quite subtle) sources of potential error (e.g. discretization error, sampling error, even accumulating floating point errors).

4. Quality testing of code is critical; focus on edge cases and combinations of edge cases.

References