对离散点进行积分的python程序实现_Stochastic Calculus(Python)(四)

在应用Ito Lemma的过程中，可以看到一个“无法解释”的现象那就是：

• ;
• ;
• ;

多数书中都把这三条做为一个“Rule”，只需要记住，在做伊藤积分的时候遇到带入值计算就可。不过后边这两个貌似还是比较如容易理解，毕竟， 如果把

看作一个时间的很小的增量（

），后两者等于0也不难。不过第一个是为什么呢？我也没能从公式上理解这个，不过可以通过一个模拟来体现。我们把

看作是一个维纳过程的增量，表示为

。

这个

可以作为我们建立一切随机过程模型的“砖块”，每一个随机过程都需要对这个“砖块”进行一定的变形。现在我们把

用伊藤积分“积”出来：

，

为了做一个对比，我们将两个不同的维纳过程（

）做相同的运算，并对其做伊藤积分，由于这两个维纳过程完全独立随机，那么可以知道

：

。

这里可以用一个求和的过程来逼近积分（也是黎曼和的思想）：

，

关于如果建立这个

“砖块”，可以参考之前的文章：

我肯定疯了：Stochastic Calculus(Python)(二)​zhuanlan.zhihu.com

关于模拟积分过程可以参考：

我肯定疯了：一些 Python 常规操作(1.1)​zhuanlan.zhihu.com

那么开始：

T = 1.0
npath = 10
for logstep in range(0,6):
    nstep = 2*10**logstep
    tgrid = np.linspace(0, T, nstep+1)
    dT = T/nstep
    S1 = np.zeros((nstep+1,npath))
    S2 = np.zeros((nstep+1,npath))
    for i in range(0,nstep):
        Z1 = np.random.normal(size=npath)
        Z2 = np.random.normal(size=npath)
        S1[i+1,:] = S1[i,:] + (np.sqrt(dT) * Z1)**2
        S2[i+1,:] = S2[i,:] + (np.sqrt(dT) * Z1) * (np.sqrt(dT) * Z2)
        
    
    plt.plot(tgrid,S1)
    plt.title('S_1 with number of steps = %d'%nstep)
    plt.show()
    plt.plot(tgrid,S2)
    plt.title('S_2 with number of steps = %d'%nstep)
    plt.show()

这里模拟了从从 2 步到 20000 步的过程，我们的预期应该是，随着步数的增加，可以看到

接近于一条

的直线，而

则和之前的多条路径的随机过程一样，只不过这次只模拟了10条。

再来说说Ito Lemma。之前提到的

:

，

写成离散形式：

，

通过Ito Lemma变形：

；或者

。

现在来模拟一下在连续状态下

和离散状态下

的

。同时，如果Ito Lemma没错，那么应该有：

• 和
  具有相同的分布；
• 和
  的 log-return应该相同。两个过程的
  应具有相同的正态分布

假设：

那么开始：

r = 0.00
v = 0.2
T = 1.
npath = 10000
nstep = 10
dt = T/nstep
S0 = 1.
Sa = S0*np.ones((nstep+1,npath))
Sm = S0*np.ones((nstep+1,npath))
lnSm = np.log(Sm)
for i in range(0,nstep):
    Z = np.random.normal(size=npath)
    dSa = r*Sa[i,:]*dt + v*np.sqrt(dt)*Sa[i,:]*Z
    Sa[i+1,:] = Sa[i,:] + dSa
    dSm = (r - v**2/2.)*dt + v*np.sqrt(dt)*Z
    lnSm[i+1,:] = lnSm[i,:] + dSm

Sm = np.exp(lnSm)
Fa = np.mean(Sa)
Fm = np.mean(Sm)

plt.hist(Sa[-1,:], 100)
plt.hist(Sm[-1,:], 100)
plt.show()

lnSaT = np.log(Sa[-1,:])
lnSmT = lnSm[-1,:];
ya, xa = np.histogram(lnSaT, 100)
ym, xm = np.histogram(lnSmT, 100)
mGBM = np.mean(lnSaT)
sGBM = np.std(lnSaT)
yn = scipy.stats.norm.pdf(xa,loc=mGBM,scale=sGBM)
plt.plot(xa[1:],ya/npath/(xa[1]-xa[0]), label='Additive')
plt.plot(xm[1:],ym/npath/(xm[1]-xm[0]), label='Multiplicative')
plt.plot(xa,yn,'b', label='PDF')
plt.legend()
plt.title('PDF of log-returns')
plt.show()

plt.plot(np.diff(np.log(Sa[:,0])), label='Additive')
plt.plot(np.diff(np.log(Sm[:,0])), label='Multiplicative')
plt.legend()
plt.title('Single path log-returns')
plt.show()

DONE！！！

---

欧洲疫情不知何时是头儿，估计实习也凉了。。。