在gBM模拟的不同时间点绘制直方图

Question

bounty到期在6天内。这个问题的答案有资格获得+ 50声望赏金。Thamu Mnyulwa想要引起更多关注来回答这个问题。

[Geometric Brownian motion (gBM)是随机过程，可以认为是标准布朗运动的扩展。

我正在尝试编写一个函数来模拟gBM的不同路径（ntraj路径），然后在列表tcheck中指定的某些点绘制直方图。一旦绘制了这些图，该函数就应每次将对数正态分布叠加在图上。

输出应该看起来像这样

除了gBM而不是标准的布朗运动过程。到目前为止，我具有生成gBM的多个路径的功能，]

def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
    '''
    DOCSTRING:
    1D geomwtric brownian motion
    INPUTS:
    ntraj = "number of trajectories"
    n = "length of a trajectory"
    T = "last time point, i.e final tradjectory t = {0,1,...,T}"
    sigma= volatility
    mu= percentage drift

    '''
    np.random.seed(52323)
    S_0 = 0

    # Discretize, dt =  time step = $t_{j+1}- t_{j}$
    dt = T/(n)  
    sqrtdt = np.sqrt(dt)

    # Container for different colors for each trajectory
    colors = plt.cm.jet(np.linspace(0,1,nTraj))
    # Container for trajectories
    xtraj=np.zeros(n+1,float)
    ztraj=np.zeros(n+1,float)
    trange=np.linspace(start = 0,stop = T ,num = n+1)

    # Simulation
    # Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1) 
    for j in range(nTraj):
        # Loop over time
        for i in range(n):
            xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
        # Loop again over time in order to make geometric drift
        ztraj = S_0 * np.exp(xtraj) # ztraj[z+1]=  ztraj[0]+ np.exp(xtraj[z])

        plt.plot(trange , xtraj,'b-',alpha=0.2, color=colors[j], lw=3.0,label="$\sigma$={}, $\mu$={:.5f}".format(sigma,mu))
    plt.title("1D Geometric Brownian Motion:\n nTraj={}, T={},\n $\Delta t$={:.3f}, $\sigma$={}, $\mu$={:.3f}".format(nTraj,T,dt,sigma,mu))    
    plt.xlabel(r'$t$')
    plt.ylabel(r'$Z_t$');

oneDGeometricBM(nTraj=5,n=10**3,T=10.0,sigma=0.8,mu=1.1)

我已经看到有关如何绘制gBM的多个路径的问题的许多答案，但是我对如何查看特定时间的直方图然后查看分布感兴趣。下面是到目前为止的功能。它不起作用，但是我无法弄清楚我在做什么错。我还添加了我得到的输出。

import yfinance as yf
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import norm, lognorm
ntraj = 10000
S_0 =0
sigma=1
mu=1
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 1.0
'''
ntraj = 10**4
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 5.0 # limits
'''
n=int(tfinal/dt)
sqrtdt = np.sqrt(dt)

x=np.zeros(shape=[ntraj,n+1], dtype=float)
z=np.zeros(shape=[ntraj,n+1], dtype=float)
zrange=np.arange(start=-xv, stop=xv, step=dt)

# Calculate the number of the bins 
binval = math.ceil(np.sqrt(ntraj))
# Nested for loop to create Drifted BM
for i in range(n):
    for j in range(ntraj):
        x[j,i+1]=x[j,i]+ sqrtdt*np.random.randn()


 #Nested loop to create gBM
for j0 in range(ntraj):
    for i0 in range(n+1):
        z[j0,i0] = 0 + np.exp(x[j0,i0])

# Loop to plot the distribution of gBM tradjectories at different times       
for i1 in range(n):
    # Compute histogram at every tsample , sample at time t
    t=(i1+1)*dt
    if t in tcheck:
       # Plot histogram on sample
       plt.hist(z[:,i1],bins=30,density=False,alpha=0.6,label=['t ={}'.format(t)] )
       # Superimpose each samples mean
       xbar = np.average(z[:,i1])
       plt.axvline(xbar, color='RED', linestyle='dashed', linewidth=2) 
       # Plot theoretic distribution { N(0, sqrt[t] ) }
       #plt.plot(xrange,norm.pdf(xrange,0.0,np.sqrt(t)),'k--')

所以总结我的问题。我正在尝试模拟gBM的多个轨迹，将结果存储在数组中，然后在该数组上循环并使用matplotlib在特定点上绘制直方图，然后最后在我的直方图上叠加对数正态分布。

Answer 1

所以我看了你的问题。我已经编辑了您的函数以停止绘制并返回xtraj，我认为这是您的布朗运动：

def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
    '''
    DOCSTRING:
    1D geomwtric brownian motion
    INPUTS:
    ntraj = "number of trajectories"
    n = "length of a trajectory"
    T = "last time point, i.e final tradjectory t = {0,1,...,T}"
    sigma= volatility
    mu= percentage drift

    '''
    np.random.seed(52323)
    S_0 = 10

    # Discretize, dt =  time step = $t_{j+1}- t_{j}$
    dt = T/(n)  
    sqrtdt = np.sqrt(dt)

    # Container for different colors for each trajectory
    colors = plt.cm.jet(np.linspace(0,1,nTraj))
    # Container for trajectories
    xtraj=np.zeros(n+1,float)
    ztraj=np.zeros(n+1,float)
    trange=np.linspace(start = 0,stop = T ,num = n+1)

    out = []
    # Simulation
    # Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1) 
    for j in range(nTraj):
        # Loop over time
        for i in range(n):
            xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu

        # Loop again over time in order to make geometric drift
    ztraj = S_0 * np.exp(xtraj) # ztraj[z+1]=  ztraj[0]+ np.exp(xtraj[z])

    return ztraj

然后，每个时间步长的位移就是数组xtraj：dx = np.ediff1d(oneDGeometricBM(...))中的差，因此我们计算这些值的直方图：

fig, ax = plt.subplots()

ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=10.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, label='T=10', density=True)
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=1.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, color='k', label='T=1', density=True)
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=5.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, color='r', label='T=5', density=True)

ax.set_xlabel('x')
ax.set_ylabel('$p(x,T)$')

ax.legend()

如示例中，我使用了3个不同的T值。为了标准化直方图，以使y轴现在代表概率p(x,T)，即。所有p*x = 1的总和，我们使用density=True参数。

编辑

我已编辑oneDGeometricBM函数以返回ztraj = S0*np.exp(xtraj)。您的初始S0值为0，因此我将其设为非零。

您可以将ztraj差异绘制为：

fig, ax = plt.subplots()

colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]

for c, T in zip(colors, T):
    ztraj = oneDGeometricBM(nTraj=5,n=10**3,T=T,sigma=0.8,mu=1.1)
    diff = np.ediff1d(ztraj)
    ax.hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-10, 10))

ax.set_xlabel('x')
ax.set_ylabel('$p(x,T)$')

ax.legend()

EDIT2

从更仔细地观察生成的直方图后，我认为您的建模是正确的，随着ztraj对于大T变大，应该调整图的xrange，可以使用range参数限制直方图。因此，我为三个单独的ztraj绘制了d(ztraj)和T。 ztraj的确似乎近似遵循对数正态分布，而ztraj的差异似乎近似遵循洛伦兹分布（必须检验该理论，也许是高斯分布）。复制代码：

fig, ax = plt.subplots(ncols=2, figsize=plt.figaspect(1./2))

colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]

for c, T in zip(colors, T):
    ztraj = oneDGeometricBM(nTraj=5,n=10**4,T=T,sigma=0.8,mu=1.1)

    ax[0].hist(ztraj, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(0, 95))

    diff = np.ediff1d(ztraj)
    ax[1].hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-5, 5))

ax[0].set_xlabel('z')
ax[0].set_ylabel('$p(z,T)$')
ax[0].set_title('Histogram of ztraj positions')

ax[1].set_xlabel('dz')
ax[1].set_ylabel('$p(dz,T)$')
ax[1].set_title('Histogram of d(ztraj) positions\nbetween time steps')

ax[0].legend()
fig.tight_layout()