指数移動平均(EMA)、２重指数移動平均(DEMA)、３重指数移動平均(TEMA)の違い
の記事で、DEMA、TEMAを紹介しましたが、今回は、
と同様に、scipyのlfilter使ったり、直接コーディングするけどNumbaで高速化したりとか、何種類かの実装を比較してみます。
入力データ
10万サンプルのランダムウォークを入力データとします。
import numpy as np
dn = np.random.randint(2, size=100000)*2-1
gwalk = np.cumprod(np.exp(dn*0.01))*100
DEMAの実装比較
EMAを使った実装
Pythonで書いた指数移動平均(EMA)のコードの比較
で最速だったEMAのコードです。
from numba import jit
@jit(nopython=True)
def EMA(x, alpha):
y = np.empty_like(x)
y[0] = x[0]
for i in range(1,len(x)):
y[i] = alpha*x[i] + (1-alpha)*y[i-1]
return y
%timeit y = EMA(gwalk, 0.15)
1000 loops, best of 3: 227 µs per loop
これを使って、DEMAの計算式のまま実装します。
def DEMA1(x, period):
alpha = 2/(period+1)
ema = EMA(x, alpha)
ema2 = EMA(ema, alpha)
y = 2*ema-ema2
return y
%timeit y1 = DEMA1(gwalk, 14)
1000 loops, best of 3: 1.19 ms per loop
EMAは２回しか使っていないのですが、EMA単独の５倍くらい遅くなっています。
lfilterを使った実装
scipyのlfilter関数を使った実装です。フィルタ係数は
指数移動平均(EMA)、２重指数移動平均(DEMA)、３重指数移動平均(TEMA)の違い
で求めた通りです。
import scipy.signal as sp
def DEMA2(x, period):
alpha = 2/(period+1)
a = [1, 2*(alpha-1), (1-alpha)**2]
b = [alpha*(2-alpha), 2*alpha*(alpha-1)]
zi = sp.lfilter_zi(b, a)
y,zf = sp.lfilter(b, a, x, zi=zi*x[0])
return y
%timeit y2 = DEMA2(gwalk, 14)
1000 loops, best of 3: 717 µs per loop
少し速くなりました。
Numbaを使った直接実装
@jit(nopython=True)
def DEMA3(x, period):
alpha = 2/(period+1)
a1 = 2*(alpha-1)
a2 = (1-alpha)**2
b0 = alpha*(2-alpha)
b1 = 2*alpha*(alpha-1)
y = np.empty_like(x)
y[0] = x[0]
y[1] = b0*x[1] + b1*x[0] - a1*y[0] - a2*y[0]
for i in range(2,len(x)):
y[i] = b0*x[i] + b1*x[i-1] - a1*y[i-1] - a2*y[i-2]
return y
%timeit y3 = DEMA3(gwalk, 14)
1000 loops, best of 3: 488 µs per loop
さらに速くなり、直接実装が最速となりました。ただ、EMAのときより差は縮まっています。次のTEMAはどうでしょうか。
TEMAの実装比較
EMAを使った実装
まずはEMAを使った実装です。
def TEMA1(x, period):
alpha = 2/(period+1)
ema = EMA(x, alpha)
ema2 = EMA(ema, alpha)
ema3 = EMA(ema2, alpha)
y = 3*ema-3*ema2+ema3
return y
%timeit y1 = TEMA1(gwalk, 14)
100 loops, best of 3: 1.89 ms per loop
まずはこれくらい。
lfilterを使った実装
def TEMA2(x, period):
alpha = 2/(period+1)
a = [1, 3*(alpha-1), 3*(1-alpha)**2, (alpha-1)**3]
b = [3*alpha*(1-alpha)+alpha**3, 3*alpha*(alpha-2)*(1-alpha), 3*alpha*(1-alpha)**2]
zi = sp.lfilter_zi(b, a)
y,zf = sp.lfilter(b, a, x, zi=zi*x[0])
return y
%timeit y2 = TEMA2(gwalk, 14)
1000 loops, best of 3: 718 µs per loop
DEMAとほとんど変わらない速度です。
Numbaを使った直接実装
@jit(nopython=True)
def TEMA3(x, period):
alpha = 2/(period+1)
a1 = 3*(alpha-1)
a2 = 3*(1-alpha)**2
a3 = (alpha-1)**3
b0 = 3*alpha*(1-alpha)+alpha**3
b1 = 3*alpha*(alpha-2)*(1-alpha)
b2 = 3*alpha*(1-alpha)**2
y = np.empty_like(x)
y[0] = x[0]
y[1] = b0*x[1] + b1*x[0] + b2*x[0] - a1*y[0] - a2*y[0] - a3*y[0]
y[2] = b0*x[2] + b1*x[1] + b2*x[0] - a1*y[1] - a2*y[0] - a3*y[0]
for i in range(3,len(x)):
y[i] = b0*x[i] + b1*x[i-1] + b2*x[i-2] - a1*y[i-1] - a2*y[i-2] - a3*y[i-3]
return y
%timeit y3 = TEMA3(gwalk, 14)
1000 loops, best of 3: 604 µs per loop
TEMAの場合も、直接実装が最速をキープしました。ただ、lfilterは次数が上がっても計算時間はほとんど変わらないことから、4次とか5次になると逆転するんじゃないかなと思います。相場の分析に高次のフィルタが必要かどうかは別として。。。
とりあえず、DEMAとTEMAはNumbaを使った直接実装をすることにします。→GitHub