求geometric brownian motion的概率密度函数

若已知几何布朗运动：

$\frac{dS_t}{S_t}=\mu dt+\sigma dW_t$

利用伊藤积分，求解上式可以得到

$S_t=S_0exp((\mu-\frac{\sigma^2}{2})t+\sigma W_t)$

对任意t，上面是一个对数正态分布，故有

$E(S_t)=S_0e^{\mu t}$

$Var(S_t)=S_0^2e^{2\mu t}(e^{\sigma^2t}-1)$
则可以得到 $S_t$ 的概率密度函数

$f_{S_t}(s;\mu,\sigma,t)=\frac1{\sqrt{2\pi}}\frac1{s\sigma\sqrt t}exp(-\frac{(ln s-lnS_0-(\mu-\frac12\sigma^2)t)^2}{2\sigma^2t})$

这样得到了每个时间点各个价格的概率。

交易策略

假设已知所有参数。

先对一般布朗运动做monte Carlo simulation，然后按定义模拟几何布朗运动：

from scipy.stats import norm
import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt
def generate_gbm(mu, sigma, T, S0,dt):
    N = round(T/dt) 
    t = np.linspace(0, T, N) 
    W = np.random.standard_normal(size = N) 
    W = np.cumsum(W)*np.sqrt(dt) ### standard brownian motion ### 
    X = (mu-0.5*sigma**2)*t + sigma*W 
    S = S0*np.exp(X) ### geometric brownian motion ### 
    return S

参数设置如下，模拟效果如图

T = 2 
mu = 0.1 
sigma = 0.01
S0 = 200
dt = 0.01
S = generate_gbm(mu,sigma,T,S0,dt)
plt.plot(S)
plt.show()

定义上面算出来的概率密度函数

def f(s,mu,sigma,t,S0):
    return (2*math.pi)**(-0.5)*(s*sigma*t**0.5)**(-1)*math.e**(-(math.log(s)-math.log(S0)-(mu-0.5*sigma**2)*t)**2/(2*sigma**2*t))

针对每个点，无论是卖单还是买单，期望的收益为

$profit=|S-S_\mu|Q-KQ^n$

可以得到 $p_{best}(S),Q_{best}(S)$ 这里以n=2为例

def profit(s,s_mu,K,Q,n):
    return abs(s-s_mu)*Q-K*Q**n
def best(s,s_mu,K,n=2):
    Qbest=abs(s-s_mu)/(2*K)
    pbest=profit(s,s_mu,K,Qbest,n=2)
    if pbest>0:
        return Qbest,pbest
    else:
        return 0,1

期望收益： $eprofit$

s_mu=s.mean()
K=0.001
m=[]
for s in np.arange(int(0.9*S0),int(1.1*S0),0.001*S0):
    sum=0
    for t in range(int(T/dt)):
        sum+=best(s,s_mu,K,2)[0]*f(s,mu,sigma,t,S0)
    m.append(sum)    
m=np.array(m)
eprofit=m.mean()

该参数下期望收益eprofit=252.49622608022156

故对每个时间点的价格s，比较s与上面模拟的平均价格 $s_mu$ ,如果价格高卖出，价格低买入，

然后算得profit与期望收益比较，高于期望收益就进行交易，单量Q=best(s,s_mu,K,n=2)

def trade(s):
    if s<s_mu:
        print ('buy')
    else:
        print ('sell')
    if best(s,s_mu,K,n)[1]>eprofit:
        print ('size=',best(s,s_mu,K,n)[0])
    else:
        print ('zero')
    return

下面是几个测试结果：