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

#前提記事
[Pythonで書いた移動平均の計算時間を比較してみた](http://qiita.com/toyolab/items/215a943c5cb0c677fa5a)

移動平均計算させるのに、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()`関数で簡単に書くことができます。

````py
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$をリストの形で引数に入れます。

````py
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)$を使わないように処理されています。

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




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Close</th>
      <th>EMA</th>
    </tr>
    <tr>
      <th>Time</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2015-01-01 13:00:00</th>
      <td>1.20962</td>
      <td>1.209620</td>
    </tr>
    <tr>
      <th>2015-01-01 13:01:00</th>
      <td>1.20962</td>
      <td>1.209620</td>
    </tr>
    <tr>
      <th>2015-01-01 13:02:00</th>
      <td>1.20961</td>
      <td>1.209616</td>
    </tr>
    <tr>
      <th>2015-01-01 13:04:00</th>
      <td>1.20983</td>
      <td>1.209686</td>
    </tr>
    <tr>
      <th>2015-01-01 13:05:00</th>
      <td>1.20988</td>
      <td>1.209742</td>
    </tr>
    <tr>
      <th>2015-01-01 13:06:00</th>
      <td>1.20982</td>
      <td>1.209762</td>
    </tr>
    <tr>
      <th>2015-01-01 13:07:00</th>
      <td>1.20987</td>
      <td>1.209788</td>
    </tr>
    <tr>
      <th>2015-01-01 13:08:00</th>
      <td>1.21008</td>
      <td>1.209855</td>
    </tr>
    <tr>
      <th>2015-01-01 13:09:00</th>
      <td>1.20996</td>
      <td>1.209878</td>
    </tr>
    <tr>
      <th>2015-01-01 13:10:00</th>
      <td>1.20977</td>
      <td>1.209855</td>
    </tr>
  </tbody>
</table>
</div>

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

<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Close</th>
      <th>EMA</th>
    </tr>
    <tr>
      <th>Time</th>
      <th></th>
      <th></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>2015-01-01 13:00:00</th>
      <td>1.20962</td>
      <td>0.219931</td>
    </tr>
    <tr>
      <th>2015-01-01 13:01:00</th>
      <td>1.20962</td>
      <td>0.399874</td>
    </tr>
    <tr>
      <th>2015-01-01 13:02:00</th>
      <td>1.20961</td>
      <td>0.547099</td>
    </tr>
    <tr>
      <th>2015-01-01 13:04:00</th>
      <td>1.20983</td>
      <td>0.667596</td>
    </tr>
    <tr>
      <th>2015-01-01 13:05:00</th>
      <td>1.20988</td>
      <td>0.766193</td>
    </tr>
    <tr>
      <th>2015-01-01 13:06:00</th>
      <td>1.20982</td>
      <td>0.846852</td>
    </tr>
    <tr>
      <th>2015-01-01 13:07:00</th>
      <td>1.20987</td>
      <td>0.912855</td>
    </tr>
    <tr>
      <th>2015-01-01 13:08:00</th>
      <td>1.21008</td>
      <td>0.966896</td>
    </tr>
    <tr>
      <th>2015-01-01 13:09:00</th>
      <td>1.20996</td>
      <td>1.011090</td>
    </tr>
    <tr>
      <th>2015-01-01 13:10:00</th>
      <td>1.20977</td>
      <td>1.047213</td>
    </tr>
  </tbody>
</table>
</div>

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

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