11
10

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

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

Posted at

#前提記事
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で実装してみます。処理するデータは前回の記事と同じです。37万個ほどの4本値時系列です。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ミリ秒くらいだったので、2倍ほど遅くなっています。

#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の場合、システム関数の分子は定数、分母は1次多項式なので、システム関数の一般式は次のように書けます。

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

11
10
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
11
10

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?