#前提記事
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を入れておくと、それらしい結果になるようです。