## Archive for the ‘**Stochastic Processes**’ Category

## Fluctuations of the time-weighted average price

Motivated by an interview question posed to a friend of mine recently, today I’d like to talk about the **time-weighted average price** (**TWAP**) of a financial asset (e.g. a stock).

Let’s consider the **price **of our stock a short time scale (e.g. **intraday**) corresponding to the interval . Its price will start at and evolve through to ; the time-weighted average price is

defined in this way is interesting, since this is the **effective price** obtained when executing a large order in the market by splitting it into smaller chunks and executing them throughout the day at a constant rate. This TWAP algorithm is one of the basic algorithmic execution tools used e.g. by asset managers to minimize market impact costs. More sophisticated ones include volume-weighted average price (VWAP) where one would adjust the volume executed at each time during the day proportionally to the typical trading volume at that time.

With this motivation, let us look at the **statistics of the TWAP** , and especially its fluctuations. Usually, one approximates the logarithm of the stock price by a brownian motion (corresponding to the assumption, that *relative* price changes are independent and identically distributed). However, on small time scales (when fluctuations are small), the compounding effect is negligible. For simplicity, I will hence assume that is a brownian motion with drift and variance . This means that *absolute* price changes are i.i.d..

## Moments of the TWAP

The Gaussian process is fully determined by its first two moments,

From this, we obtain the first moments of the TWAP (for simplicity I subtracted the trivial component):

In other words, the **variance of the TWAP relates to the variance of the intraday stock price change** as

( denotes connected expectation values here).

This procedure can easily be continued to calculate higher moments of the TWAP in terms of and . In fact, given that the linear combination of Gaussians is again a Gaussian, we can directly infer that the TWAP has a Gaussian distribution with mean and variance .

## Distribution of the TWAP via the MSR formalism

There is another way to obtain the full distribution of the TWAP directly, via the Martin-Siggia-Rose formalism. Using the MSR approach, we know that the generating function of any linear functional of the Brownian is

where is the solution of

,

Applying this to a constant , we get . From this follows the **generating function of the TWAP **

This is clearly the generating function of a Gaussian distribution with the same moments as computed above.

## Stochastic Logistic Growth

A classic model for the growth of a population with a finite carrying capacity is **logistic growth**. It is usually formulated as a differential equation,

. (1)

Here is the size of the population at time , is the **growth rate** and is the **carrying capacity**.

The dynamics of eq. (1) is **deterministic**: The initial population grows (or decays) towards the constant carrying capacity , which is a fixed point of eq. (1). This is seen in the solid trajectories in the figure below:

To make this model more realistic, let’s see how we can extend it to include **stochastic fluctuations **(semi-transparent trajectories in the figure above). I’ll look in the following at a simple stochastic logistic growth model (motivated by some discussions with a friend), where the **steady state can be calculated exactly**. The effect of stochasticity on the steady state is twofold:

- The population size for long times is not fixed, but fluctuates on a scale around the carrying capacity.
- The average population size is increased above the carrying capacity , but the shift goes to 0 as increases (i.e. the deterministic model is recovered for large ).

Now let’s look at the calculation in detail…

## A stochastic, individual-based logistic growth model

In eq. (1), the population is described by a real-valued function . In reality, populations consist of **discrete individuals** and the population size doesn’t change continuously. So, a more realistic approach is to describe the population dynamics as a **birth-death process** with the following reactions:

(2)

In other words, we assume that during a time interval two events may happen:

**Birth**: With a probability , any individual may give birth to offspring, thus increasing the population size by 1.**Death due to competition**: With a probability , out of any two individuals one may die due to the competition for common resources. Thus the population size decreases by 1. Note that is related to the carrying capacity in eq. (1) by .

Note that the stochasticity in this model is not due to random external influences but due to the discreteness of the population (**demographic noise**).

## Solution of the stochastic model

The system of reaction equations (2) translates into the following **master equation** for the probabilities of the population size at time being equal to :

(3)

This looks daunting. However, it can be simplified a lot by introducing the **generating function**

After some algebra we obtain

.

This looks simpler than eq. (3) but still finding the full time-dependent solution does not seem feasible. Let’s focus on the **steady state** where . Together with the boundary conditions and , we obtain the steady-state solution

. (4)

Here we set to connect to the notation of eq. (1). Correspondingly, the steady-state probabilities for population size are

. (5)

Of course, these can equivalently also be obtained by solving the recursion relation (3) with . This result can be easily checked against simulations, see figure below.

## Steady state with stochasticity

Let’s try to understand what the steady state of our stochastic model looks like. From eq. (4) we easily obtain the first moments of the distribution of the population size. The **mean population size** is:

So for large carrying capacities , we recover the same result as in the deterministic model (which is reasonable since for large population sizes demographic noise is less pronounced!). However, for smaller **fluctuations increase the average population size**. E.g. for , the average population size in our stochastic model is .

Now let us look at the **variance**,

.

For large , we have . So we see that stochastic fluctuations **spread out the steady state **to a width around the carrying capacity.

## Estimating expected growth

Let’s take some fluctuating time series — say the value of some financial asset, like a stock price. **What is its average growth rate?** This seemingly trivial question came up in a recent discussion I had with a friend; obviously, it is relevant for quantitative finance and many other applications. Looking at it in more detail, it turns out that a precise answer is actually not that simple, and **depends on the time range** for which one would like to estimate the expected growth. So I’d like to share here some insights on this problem and its interesting connections to stochastic processes. Suprisingly, this was only studied in the literature quite recently!

## The setup

Let’s consider the price of an asset at discrete times . The corresponding **growth rates** are defined by

(1)

I.e. if , the growth rate is . Let’s also assume that our growth process is **stationary**, i.e. that the distribution of the growth rates does not change in time. Say for concreteness that the growth rates are i.i.d. (independent identically distributed) random variables.

## The arithmetic average

The most immediate idea for computing the average growth rate is just to take the **arithmetic average**, i.e.

(2)

What does this give us? Obviously, by the assumption of stationarity made above, with increasing sample size the expectation value for the growth rate at the next time step (or any other future time step) is approximated better and better by :

This seems to be what we’d expect for an average growth rate, so what’s missing? Let’s go back to our sample of asset prices and take the total growth

By the relationship between the arithmetic and the geometric mean, we know that

(3)

The left-hand side is the total growth after time periods and the right-hand side is its estimate from the arithmetic mean in eq. (2). Equality only holds in eq. (3) when all the are equal. When the fluctuate the total growth rate will always be strictly less than the estimate from the arithmetic mean, even as . **So, by taking the arithmetic mean to obtain the total growth in the price of our asset at the end of the observation period, we systematically overestimate it.** The cause for this is the **compounding** performed to obtain the total growth over more than one time step. This makes our observable a **nonlinear function of the growth rate in a single time step**. Thus, fluctuations in the growth rate don’t average out, and yield a net shift in its expectation value.

The first explicit observation of this effect I’ve found is in a 1974 paper by M. E. Blume, *Unbiased Estimators of Long-Run Expected Rates of Return*. Quantitative estimates from that paper and from a later study by Jacquier, Kane and Marcus show that in typical situations this overestimation may easily be **as large as 25%-100%**.

## The geometric average

So we see that the arithmetic mean does not provide an unbiased estimate of the compounded growth over a longer time period. Another natural way to obtain the average growth, which intuitively seems more adapted to capture that effect, is to take the **geometric average of the growth rates**:

(4)

Now, by construction, the total growth at the end of our observation period, i.e. after time periods, is correctly captured:

But this only solves the problem which we observed after eq. (3) for this specific case, when estimating growth for a **time period whose length is exactly the same as the observation time span** from which we obtain the average. For shorter time periods (in particular, when estimating the growth rate for a single time step ), the geometric mean will now **underestimate** the growth rate. On the other hand, for time periods even longer than the observation period, it will still overestimate it like the arithmetic mean (see again the paper by Blume for a more detailed discussion).

## Unbiased estimators

Considering the results above, the issue becomes clearer:** A reasonable (i.e. unbiased) estimate for the compounded growth over a time period requires a formula that takes into account both the number of observed time steps , and the number of time steps over which we’d like to estimate the compounded growth . ** For we can use the arithmetic mean, for we can use the geometric mean, and for general we need another estimator altogether. For the general case, Blume proposes the following (approximately) **unbiased estimator** :

(5)

**Eq. (5) is a reasonable approximation for the compounded growth**, not having any further information on the form of the distribution, correlations, etc.

For , i.e. estimating growth in a single time step, this gives just the arithmetic mean which is fine as we saw above. For , this gives the geometric mean which is also correct. For other values of , eq. (5) is a linear combination of the arithmetic and the geometric mean.

To see how the coefficients in eq. (5) arise, let us start with an ansatz of the form

(6)

Let us further split up the growth rates as , where is the “true” average growth rate and are fluctuations. Inserting this as well as the definitions of into eq. (6) we get

Now let us assume that the fluctuations are small, and satisfy . Expanding to second order in (the first order vanishes), and taking the expectation value, we obtain

To obtain an estimator that is unbiased (to second order in ), we now choose and such that the term of order is just the true growth rate , and the term of order vanishes. This gives the system

The solution of this linear system for and yields exactly the coefficients in eq. (5).

Of course, here we make the assumption of small fluctuations and also a specific ansatz for the estimator. If one has more information on the distribution of the growth rates this may not be the most adequate one, but with what we know there’s not much more we can do!

## Outlook

As you can see from the above discussion, estimating the expected (compounded) growth over a series of time steps is more complex than it appears at first sight. I’ve shown some basic results, but didn’t touch on many other important aspects:

- In addition to the growth rate , it is also interesting to consider the discount factor . Blume’s approach is extended to this observable in this paper by Ian Cooper.
- If one assumes the growth factors in eq. (1) to be
**log-normally distributed, the problem can be treated analytically**. Jacquier, Kane and Marcus discuss this case in detail in this paper, and also provide an exact result for the unbiased estimator. - The assumption of independent, identically distributed growth rates is not very realistic. On the one hand, we expect the distribution from which the annual returns are drawn to vary in time (i.e. due to underlying macroeconomic conditions). This is discussed briefly in Cooper’s paper. On the other hand, we also expect some amount of correlation between subsequent time steps, even if the underlying distribution does not change. It would be interesting to see how this modifies the results above.

Let me know if you find these interesting — I’ll be glad to expand on that in a future post!

## Fokker-Planck equation for a jump diffusion process

One of the simplest stochastic processes is a Brownian motion with drift, or a **diffusion process**:

(1)

Here, is Gaussian white noise with mean zero and variance . Its integral is a Brownian motion.

Continuous Itô stochastic processes such as eq. (1) are insufficient for applications where the random variable may jump suddenly (such as in avalanches). A natural extension of (1) for modelling this is a so-called **jump-diffusion process**. Let us suppose that our jump sizes are positive, independent and identically distributed with density . Then, the jump diffusion process is

(2)

where are i.i.d. jump sizes as above, and is the number of jumps encountered up to time . For simplicitly, let us assume that jumps occur independently with rate , i.e. that the probability to have a jump in a time interval is . Then, is a Poisson process with rate .

It is well-known that the diffusion process in (1) is equivalently described by a partial differential equation for the distribution of , the **Fokker-Planck equation** (FPE)

. (3)

This representation is useful e.g. for first-passage problems: they correspond to various boundaries introduced in the PDE (3). **So how does one generalise the Fokker-Planck (3) to the case of the jump-diffusion process in (2)? ** I will explain in the following that the answer is

, (4)

and then discuss a specific example.

## 1. Deriving the jump-diffusion FPE

Let us consider a time step from to . The probability for a jump to occur during this interval is , so

, (5)

where denotes averaging over all realizations of the Brownian motion , and denotes averaging over the distribution of the jump size . Since the jump term is already multiplied by the jump probability , the drift and noise contributions there are of higher order in and were dropped.

The averaging over the jump size in (5) yields a convolution with the jump size distribution :

.

The average over the noise in (5) is the same as for standard diffusion. During the interval , the increment of the noise term in (2) is

,

where the last equality is a definition for .

Since is a Brownian motion, is normally distributed:

.

Thus, the average over in (5) is

.

In the last line we used the fact that the normal distribution has mean zero and variance 1, and dropped all terms of higher order than . Inserting all this into (5), we obtain the jump-diffusion FPE (4).

## 2. Example: Exponential jumps

Now let us consider a simple example, where the jump sizes are distributed exponentially:

.

We will compute the distribution of , starting from .

For this, we solve the jump-diffusion FPE (4) with the initial condition . Let us take a Fourier transform of (4) from to :

.

The different spatial modes decouple, and instead of a partial differential equation we obtain an ordinary differential equation for each value of the Fourier variable .

For the above exponential form of the jump size distribution ,

.

Furthermore, the initial condition gives . Hence, the solution $\tilde{P}$ reads

.

While it does not seem feasible to invert this Fourier transform to obtain a closed expression for (but if anyone has an idea, let me know!), is already enough to determine the moments of . Taking derivatives, we obtain for example

.

Similarly, solving the FPE (4) with an absorbing boundary allows computing first-passage times (or at least their moments) for our jump-diffusion process.

Have fun, and do let me know if you have any comments!

## How does an absorbing boundary modify the propagator of a stochastic process?

A recurring theme throughout my blog are stochastic processes, widely used to model quantities under the influence of both deterministic evolution laws and random fluctuations. A complete solution for a stochastic process is its *propagator*: The probability density for the final state at some time given an initial condition at time . The propagator depends on the boundary conditions imposed on the stochastic process. For example, an absorbing boundary is useful for describing extinction in a stochastic model of population dynamics: When the population of a given species reaches zero, it cannot be re-introduced.

Typically, however, the “free” propagator for a stochastic process — i.e. the propagator without any boundaries — is much easier to determine than the propagator with any given boundary condition.

However, common sense tells us that the free propagator for a stochastic process contains already all there is to know about it. So I have wondered for quite a while if and how one can recover e.g. the propagator with an absorbing boundary, knowing only the free propagator. I found the solution only quite recently in this preprint on arXiv, and think that it is interesting enough to explain it here briefly.

## The result

Consider a one-dimensional, stationary, Markovian stochastic process. Let us denote its *free propagator* by

This is the density for the final value at time , given an initial value at time . Stationarity of the process implies that is a function of the time difference only. Define the Laplace transform of the free propagator by

Now let us add an absorbing boundary at . Then, the propagator

is the density for the final value at time , given an initial value at time , *without touching the line between and *. **The upshot of this post is the following simple relationship between the Laplace transforms of the free propagator and of the propagator with absorbing boundary :**

(1)

In case the denominator on the right-hand side is zero, one may have to perform some limiting procedure in order to make sense of this formula. We can now check how nicely it works for simple cases, like Brownian motion with drift, where the result is known by other means.

## Example: Brownian motion with drift

Let us consider an example, the Brownian motion with drift . Its stochastic differential equation is

where is white noise with .

The free propagator is a simple Gaussian

Its Laplace transform is given by:

To compute , let us assume , . Note this is without loss of generality, since if and are on opposite sides of the absorbing boundary the propagator is zero. Applying the result (1), we have

The inverse Laplace transform of this expression is:

This is, of course, exactly the same result (the so-called **Bachelier-Levy formula**) which I had already discussed earlier in a guest blog post at Ricardo’s blog, using slightly different methods (a variant of the Cameron-Martin-Girsanov theorem). Note that some signs which were wrong in the original post are corrected here.

## The proof

Now let us go back to the case of a general stochastic process, and discuss how formula (1) can be proven. Let us assume that the stochastic process in question has a path-integral representation, i.e. that

for some action .

Now, to add an absorbing boundary at we weigh each crossing of by a factor in the path integral, and then let . This effectively retains only those trajectories in the path integral, which never touch . We thus have:

.

Now, to evaluate the path integral, we expand it in a series in (or, in other words, in the number of crossings of ). We get:

.

Note that here I cancelled the factor , coming from the Taylor expansion of , against the possibilities of ordering the times . This is why I can assume them to be in ascending order above. Also, the functions actually set all .

Taking Laplace transforms, the time integrals (which are a -fold convolution) simplify to a simple product:

.

Now, one easily obtains the claimed result:

.

This computation can, of course, be generalized to more complicated situations like several boundaries, higher dimensions, etc. For details see again the original preprint by C. Grosche. Have fun and let me know if you have questions!