Pricing uniswap v3 with stochastic process, Part2

In the previous article, we encountered a double first-hitting time problem. We obtained the expectation of the stopping time using the following three martingales:

martingales

In this article, we will use the third martingale to derive the density function of the stopping time and verify it using Monte Carlo simulation. We are just one step away from the final pricing.

The third martingale

Following the same approach as the previous article, we will apply the optimal stopping theorem. The following is a rough derivation, skipping the discussion of $E[I_{\{\tau<\infty\}}]=1$.

We know that:

substitute it as:

Note that sigma is not the volatility in this case, it is just an algebraic symbol that holds for any sigma. We set $\sigma = \sqrt{2\alpha}$ to obtain:

Laplace Inverse

In this section, we will handle the Laplace inverse in a very technical way to obtain the distribution of $\tau$. If you are not concerned about the mathematical correctness, you can safely accept the conclusion and skip to the next section for MC verification:

Assuming Wt​ is a standard Brownian motion starting at 0, with a<0<b, we can define a stopping time:

Then, for any 0t>0, the Laplace transform of the stopping time τ is given by

First, let’s discuss two martingales:

By the properties of martingales, it can be observed that their linear combination, for any k∈R, is still a martingale:

According to the martingale stopping theorem, the expectation of the martingale when it touches the boundary should be equal to the initial expectation:

k can be any real number, so we can choose an appropriate k to make the following equation hold for both boundary values:

It means that:

If the above holds, we get:

Then we have the following expectation:

We found that by letting λ=sqrt(2t)​, we obtain the Laplace transform and then perform a hyperbolic trigonometric transformation:

Density Function

According to the previous section, we have

(In order to distinguish it from the time domain variable τ in the Laplace inverse transform, we denote the stopping time in this section as τ(a,b​).

Where.

Solving this Laplace inverse transform is not elementary. By consulting the Tables of Integral Transforms by A. Erdelyi, W. Magnus, F. Oberhettinger, and F.G. Tricomi (Eds.), we can find that

where θ1​ is one of the famous four Jacobi Theta functions. The τ appearing in the Theta function does not represent a stopping time, but a time domain variable:

For convenience of calculation, we can transform θ1​ into a more compact form:

Then we can continue with the operation after the Laplace inverse transform lookup:

We can set:

Finally, we get the density we want:

CDF

In this section, we simply need to integrate the probability density function:

Making the substitution,

The original expression becomes:

Verified by MC

We verify the accuracy of the previously mentioned PDF formula using the Monte Carlo algorithm.
The parameters for the MC algorithm are set as follows:
1. $dS=\sigma dW_t$
2. Time length is 3
3. dt is 3/10000
4. The boundaries (a, b) are set to (-0.5, 0.5)

The following code is used to calculate the hitting time of the paths:

def first_hitting_time(paths, a, b):
    hit_index = np.argmax((paths < a) | (paths > b), axis=1)
    N = paths.shape[1]
    # if the brownian motion never hit the boundary, then the hitting time is N
    hit_index[hit_index == 0] = N
    first_hitting_times = hit_index / N
    return first_hitting_times

In the subsequent discussion of pricing models, we will also use the MC method for comparative verification.

For the MC hitting time, we use the KDE method to calculate the probability density. By comparison, We can find that the KDE method may give a probability of hitting time less than 0, which leads to differences near 0. Other than that, the two distributions are consistent.

Code

import math
import numpy as np
import matplotlib
import matplotlib.pyplot as plt


def iter_k_cc_t(b, a, k_range, t):
    k_s = np.arange(-k_range, k_range+1)
    negative_series = np.where(k_s % 2 == 0, 1, -1)
    terms = k_s*(b-a)+b
    return np.sum(negative_series*terms/(np.sqrt(2*np.pi)*t**(3/2))*np.exp(-terms**2/(2*t)))


# define the brownian motion
def standard_brownian_motion(T, N, M):
    dt = T / N
    dW = np.sqrt(dt) * np.random.randn(M, N)
    W = np.cumsum(dW, axis=1)
    return W


# define the first hitting time of brownian motion,hit the boundary a or b, a<0<b,
# return the hitting time.use numpy broadcasting
def first_hitting_time(paths, a, b):
    hit_index = np.argmax((paths < a) | (paths > b), axis=1)
    N = paths.shape[1]
    # if the brownian motion never hit the boundary, then the hitting time is N
    hit_index[hit_index == 0] = N
    first_hitting_times = hit_index / N
    return first_hitting_times


def main():
    T = 3
    N = 10000
    M = 80000
    paths = standard_brownian_motion(T, N, M)
    a = -0.5
    b = 0.5
    stop_time = first_hitting_time(paths, a, b) * T

    import seaborn as sns
    sns.kdeplot(stop_time, shade=True,color="g", label="MC")
    # plot laplace approximation
    k_range = 1000
    t = np.linspace(0, T, 1000)
    laplace_value = [iter_k_cc_t(b, a, k_range, tt) for tt in t]
    sns.lineplot(x=t, y=laplace_value, label='laplace')
    plt.legend()

    plt.savefig('histplot.png')



if __name__ == '__main__':
    main()