多重指数移動平均(DEMA, TEMA)のPythonによる実装比較

[指数移動平均(EMA)、２重指数移動平均(DEMA)、３重指数移動平均(TEMA)の違い](http://qiita.com/toyolab/items/a071d919e4132db71285)

の記事で、DEMA、TEMAを紹介しましたが、今回は、

[Pythonで書いた指数移動平均(EMA)のコードの比較](http://qiita.com/toyolab/items/6872b32d9fa1763345d8)

と同様に、scipyのlfilter使ったり、直接コーディングするけどNumbaで高速化したりとか、何種類かの実装を比較してみます。

# 入力データ
10万サンプルのランダムウォークを入力データとします。

````py
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)のコードの比較](http://qiita.com/toyolab/items/6872b32d9fa1763345d8)
で最速だったEMAのコードです。

````py
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の計算式のまま実装します。

````py
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)の違い](http://qiita.com/toyolab/items/a071d919e4132db71285)

で求めた通りです。

````py
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を使った直接実装

````py
@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を使った実装です。

````py
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を使った実装

````py
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を使った直接実装

````py
@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](https://github.com/toyolab/MT5IndicatorsPy)