IIRフィルタ型の移動平均をpandasとscipyで比較してみた

前提記事

Pythonで書いた移動平均の計算時間を比較してみた

移動平均計算させるのに、Pythonでfor文使うのは論外で、pandas、scipyの関数使った方がよいという結果。ただ、上記の記事ではSMA、LWMAというFIRフィルタ型の移動平均だったので、今回はEMA、SMMAといったIIRフィルタ型の移動平均について調べてみました。

EMA

EMAとは、Exponential Moving Averageの略で、次のような式で表します。

$$y(n)=\alpha x(n)+(1-\alpha)y(n-1)$$

ここで、$\alpha$は0～1の実数パラメータです。この式では期間を表すパラメータは入っていないのですが、SMA、LWMAでは期間を表すパラメータを使うので、それに合わせてEMAでも期間パラメータを使うことが多いです。期間を$p$とすると、$\alpha=2/(p+1)$としたのがEMAで、$\alpha=1/p$としたのがSMMAです。個人的にはEMAもSMMAも同じ式なので区別する必要はないと思うのですが、MetaTraderで区別して使われているので一応言及しておきます。

pandasで実装

まずはpandasで実装してみます。処理するデータは前回の記事と同じです。３７万個ほどの４本値時系列です。pandasは時系列データ向けなので、EMAもewm()とmean()関数で簡単に書くことができます。

import numpy as np
import pandas as pd
dataM1 = pd.read_csv('DAT_ASCII_EURUSD_M1_2015.csv', sep=';',
                     names=('Time','Open','High','Low','Close', ''),
                     index_col='Time', parse_dates=True)

def EMA(s, ma_period):
    return s.ewm(span=ma_period).mean()

%timeit MA = EMA(dataM1['Close'], 10)

今回も実行時間を比較するので、計測結果を示します。

10 loops, best of 3: 32.2 ms per loop

SMAの場合、16ミリ秒くらいだったので、２倍ほど遅くなっています。

scipyで実装する前にちょっと変換

scipyのlfilter()を使う場合、ただ$\alpha$を入れるだけではダメで、IIRフィルタの形式に変えて、その係数を入れる必要があります。ということで、ちょっと変換します。（詳しい理論は省略します。ディジタル信号処理の基本です。）

EMAの式の両辺を$z$変換します。

$$Y(z)=\alpha X(z)+(1-\alpha)z^{-1}Y(z)$$

$Y(z)$を左辺にまとめると

$$\{1-(1-\alpha)z^{-1}\}Y(z)=\alpha X(z)$$

となり、$Y(z)/X(z)$を$H(z)$とすると、

$$H(z)=\frac{\alpha}{1-(1-\alpha)z^{-1}}$$

と書けます。これがIIRフィルタのシステム関数となります。lfilter()の引数に渡すのは、このシステム関数の分子多項式の係数、分母多項式の係数です。

EMAの場合、システム関数の分子は定数、分母は１次多項式なので、システム関数の一般式は次のように書けます。

$$H(z)=\frac{b_0}{a_0+a_1z^{-1}}$$

係数を比較すると、$b$、$a$の係数は以下のようになることがわかります。

$$b_0=\alpha, \ a_0=1, \ a_1=\alpha-1$$

scipyで実装

ということで、EMAをscipyのlflter()を使って実装すると以下のように書けます。上の$b$、$a$をリストの形で引数に入れます。

from scipy.signal import lfilter
def EMAnew(s, ma_period):
    alpha = 2/(ma_period+1)
    y = lfilter([alpha], [1,alpha-1], s)
    return pd.Series(y, index=s.index)

%timeit MA = EMAnew(dataM1['Close'], 10)

結果は

100 loops, best of 3: 3.08 ms per loop

となり、SMAの場合とほぼ同じ速さとなりました。やはり、lfilter()速い。

lfilter()の初期条件の調整

今回もlfilter()が速いという結果になったのですが、若干処理結果に問題があります。

EMAだと、最初の出力$y(0)$を計算するさいに、データの存在しない$y(-1)$を使うのですが、pandasの場合、時系列データ向きということで、$y(0)=x(0)$として、$y(-1)$を使わないように処理されています。

pd.DataFrame({'Close':dataM1['Close'],'EMA':MA}).head(10)

	Close	EMA
Time
2015-01-01 13:00:00	1.20962	1.209620
2015-01-01 13:01:00	1.20962	1.209620
2015-01-01 13:02:00	1.20961	1.209616
2015-01-01 13:04:00	1.20983	1.209686
2015-01-01 13:05:00	1.20988	1.209742
2015-01-01 13:06:00	1.20982	1.209762
2015-01-01 13:07:00	1.20987	1.209788
2015-01-01 13:08:00	1.21008	1.209855
2015-01-01 13:09:00	1.20996	1.209878
2015-01-01 13:10:00	1.20977	1.209855

この場合、EMAの結果も入力の時系列とそう外れてはいないのですが、lfilter()の場合、$y(-1)=0$として計算してしまうので、最初の方のEMAの値が入力からかなりずれてしまうことになります。

	Close	EMA
Time
2015-01-01 13:00:00	1.20962	0.219931
2015-01-01 13:01:00	1.20962	0.399874
2015-01-01 13:02:00	1.20961	0.547099
2015-01-01 13:04:00	1.20983	0.667596
2015-01-01 13:05:00	1.20988	0.766193
2015-01-01 13:06:00	1.20982	0.846852
2015-01-01 13:07:00	1.20987	0.912855
2015-01-01 13:08:00	1.21008	0.966896
2015-01-01 13:09:00	1.20996	1.011090
2015-01-01 13:10:00	1.20977	1.047213

この問題はlfilter()のオプショナル引数で解決することができるようです。次のように書くことで、pandasとほぼ同じ結果が得られました。

def EMAnew(s, ma_period):
    alpha = 2/(ma_period+1)
    y,zf = lfilter([alpha], [1,alpha-1], s, zi=[s[0]*(1-alpha)])
    return pd.Series(y, index=s.index)

ここで、ziは状態変数の初期値ということなので、単なる入力、出力の初期値のことではないのですが、ここでは、$y(0)=\alpha x(0)+zi=x(0)$になるようなziを入れておくと、それらしい結果になるようです。