C++/Python/Matlab执行效率分析

以前一直觉得C++效率最高，速度最快，但是今天做的一个实验结果大大出乎我的意料—Python使用向量处理效率一点都不慢，甚至高于C++在O2优化后的效率。 Matlab效率更高。 这为以后选取语言提供了一个很好的参考。

问题起源与对场内期权MC定价时，一步到最后与按天到最后在计算精度上有无差别，镜像问题是FDM定价一步到最后，与按天到最后计算精度上有无差别。我的观点显然是前者无差，后者影响显著。无奈，领导观点竟然相反，为了说明问题，讲理论肯定是不行了，只能用数据说话。为了验证MC的结果，有了本实验。

实验初期，先采用Python写的，因为简单，就用了最简单低级的代码思路。写完后运行发现速度还可以（轨道数目较少，1k条）。于是萌生了与C++比较下效率的想法。于是速度照搬改成C++版本。为了能体现出差距，将轨道数目都改成10w条，结果发现，在简单代码结构下，python不是一般的慢！ 考虑到当初Matlab在矩阵处理时体现的优势，于是想python应该也有同样的特点，于是将python代码改成向量化的。这一改，就发现了有趣的东西，在向量化代码结构下，python执行效率极高，甚至高过C++在O2优化下的效率。

废话不多说，先列出结果数据。表中数据单位为秒。



图中列出了几个相关实验:

- C++原始代码及经过O2优化后的代码。

- C++通过OMP并行后的代码及经过O2优化过的代码

- Python 串行代码及向量化代码

- Matlab 串行及向量化代码

图中给出的结果是让人惊奇的：Python经过向量化后代码效率极大的提高，这跟Matlab一样，两者都要好于经过O2优化过的C++并行代码。Matlab与Python对for循环效率极低，使用时要慎重！

附件一 C++ omp并行与非并行代码（注释掉pragma那行即可）

#include <random>
#include <cmath>
#include <ctime>
#include <omp.h>
using namespace std;

default_random_engine eng(time(NULL));
normal_distribution<double> norm(0,1);

double payoff(double &s, double &k)
{
return fmax(s-k,0.);
}

double pnode(double &s0, double& u, double& q, double& vol, double& t)
{
return s0*exp((u-q-0.5*vol*vol)*t + sqrt(t)*vol*norm(eng));
}

double mc(int& npath, const int& msteps)
{
double u = 0;
double q = 0;
double vol = 0.2;
double r = 0.035;
double s0 = 1.;
double t = 1.;
double k = 1.;

double val = 0.;
#pragma omp parallel for reduction(+:val)
for(int i = 0; i<npath; ++i)
{
double dt = t / msteps;
double st = s0;
for(int j = 0; j<msteps; ++j)
st=pnode(st,u,q,vol,dt);
double v=payoff(st,k)*exp(-r * t);
val +=v;
}

printf("%d steps to terminal with %d paths : v=%f\n", msteps, npath, val/npath);

}

int main()
{
clock_t t1 = clock();
double ompt1 = omp_get_wtime();
int n= 100000;
int m = 100;
mc(n, 1);
mc(n,m);
double ompt2 = omp_get_wtime();
printf("Time cost %f\n", ompt2-ompt1);
return 0;

}

附件二 串行Python程序

# -*- coding: utf-8 -*-
"""
Created on Sat Aug 12 10:52:40 2017

@author: kt
"""

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

def payoff(s, k):
# Calculate the tereminal payoff
return np.max([s-k,0.])

def pnode(s0, u, q, vol, t):
return s0*np.exp((u-q-0.5*vol**2)*t+vol*np.sqrt(t)*np.random.normal(0,1))

def mc(npath, msteps):
u = 0
q = 0
vol = 0.2
r = 0.035
s0 = 1.0
t = 1
k = 1.0

val = 0.;
for i in range(npath):
dt = t/msteps
st = s0
for j in range(msteps):
st = pnode(st,u,q,vol,dt)
val += payoff(st,k)*np.exp(-r*t)
print('%d Step to Terminal: v=%f' % (msteps,val/npath))

t1 = time.clock()
npath = 100000
msteps=100
mc(npath,1)
mc(npath, msteps)
print('Time cost %f' % (time.clock() - t1))

附件三 向量化Python程序

# -*- coding: utf-8 -*-
"""
Created on Sat Aug 12 10:52:40 2017

@author: kt
"""

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

def payoff(s, k):
# Calculate the tereminal payoff
v = s-k
v[s-k<0] = 0
return v

def TermV(s0, u, q, vol, t, npath, msteps):
dt = t/msteps
s = s0*np.ones(npath)
return s*np.exp((u-q-0.5*vol**2)*t+ \
np.sum(vol*np.sqrt(dt)*\
np.random.normal(size=[npath,msteps]),1))

def mc(npath, msteps):
u = 0
q = 0
vol = 0.2
r = 0.035
s0 = 1.0
t = 1
k = 1.0

st = TermV(s0,u,q,vol,t,npath,msteps)
v= np.sum(payoff(st,k))/npath *np.exp(-r*t)

print('%d Step to Terminal: v=%f' % (msteps,v))

t1 = time.clock()
npath = 100000
msteps=100
mc(npath,1)
mc(npath, msteps)
print('Time cost %f' % (time.clock() - t1))

附件四 串行Matlab代码

function mc()
tic;
npath= 100000;
msteps = 100;
mc_s(npath,1);
mc_s(npath,msteps);
toc
end

function y=payoff(s,k)
y = max(s-k,0);
end

function y=pnode(s0,u,q,vol,t)
y=s0*exp((u-q-0.5*vol^2)*t+vol*sqrt(t)*randn());
end

function mc_s(npath, msteps)
u = 0;
q=0;
vol=0.2;
r=0.035;
s0=1;
t=1;
k=1;

val = 0;
for i = 1:npath
dt = t/msteps;
st = s0;
for j = 1:msteps
st = pnode(st,u,q,vol,dt);
end
val = val + payoff(st,k)*exp(-r*t);
end
fprintf('%d step to Terminal :v=%f\n', msteps, val/npath);
end

附件五 向量化Matlab代码

function mc2()
tic;
npath= 100000;
msteps = 100;
mc_s(npath,1);
mc_s(npath,msteps);
toc
end

function y=payoff(s,k)
y=s-k;
y(s-k<0) = 0;
end

function y=TermV(s0,u,q,vol,t,npath, msteps)
dt = t/msteps;
s = s0*ones(npath,1);
y=s.*exp((u-q-0.5*vol^2)*t+vol*sqrt(dt)*sum(randn(npath,msteps),2));
end

function mc_s(npath, msteps)
u = 0;
q=0;
vol=0.2;
r=0.035;
s0=1;
t=1;
k=1;

st = TermV(s0,u,q,vol,t,npath,msteps);
v= sum(payoff(st,k))/npath *exp(-r*t);
fprintf('%d step to Terminal :v=%f\n', msteps, v);
end

结论

实验结果说明，在处理好for的关系后，Python的效率还是挺高的，当然Matlab效率更高！

Obviously, 不要忘了实验的初心：对于MC方法定价场内期权，一步到位与中间详细刻画路径结果精度是一样的，理论依据就是：MC方法是基于概率分布的方法，而布朗运动增量是独立同分布的，他们和的分布与一步到位终端变量的分布是同分布的。

1 step to Terminal :v=0.076928

100 step to Terminal :v=0.076514

时间已过 0.632972 秒。