指数平滑移動平均における初期値の決定・調整方法とR、Pythonでの計算方法

指数平滑移動平均は移動平均の一種で、数列のある値 $y_t$ に対応する移動平均値 $S_t$ を、次のような漸化式で定める。

$S_t = \alpha y_t + (1-\alpha)S_{t-1}$

ここで $\alpha$ は新しい値に対する重みの割合を調整する係数で、 $0 \lt \alpha \lt 1$ であり、1に近いほど直近の値の影響が大きくなる。

漸化式を級数展開して、次のような表現をすることもできる。

$S_t = \alpha (y_t + (1-\alpha)y_{t-1} + (1 - \alpha)^{2} y_{t-2} + (1-\alpha)^{3} y_{t-3} + ... )$

$y_t$ に近いほど重みが大きく、遠い値の重みが指数関数的に減衰していくことが分かる。

指数平滑移動平均はExponential Moving Averageの頭文字をとってEMAと呼ばれることが多いので、以下EMAとする。

初期値問題

EMAは漸化式で定義されるので無限の数列を仮定しているが、現実のデータは有限であるので、例えば $y_{t-2}$ 以前の数列が入手できないときに $S_{t-1}$ をどうするかという問題がある。これにはいくつかの方法があるのだが、いくつかの方法があるという前提で解説されたテキストはあまり多くなく、調べるのに苦労した。

$S$ を最初の $y$ で代用する方法

比較的単純な方法として、 $S_{t-1} = y_{t_1}$ としてしまう、つまり利用可能な最も古いデータの値をEMAの初期値とするものがある。ネットで適当に検索して出てくるEMAの解説記事でもよく紹介されている。

RでもPythonでもEMAを計算する標準の関数は無いので何らかのパッケージを利用することになる。

Rの場合、pracmaパッケージのmovavg()関数が指数移動平均に対応しており、定義を抜粋すると次のようになっており、最初の結果が数列の最初の値で初期化されていることが分かる。

# R
# pracma::movavg
    a <- 2/(n + 1)
    y[1] <- x[1]
    for (k in 2:nx) y[k] <- a * x[k] + (1 - a) * y[k - 1]

Pythonの場合はpandas.DataFrame.ewmを利用するのが一般的と思われるが、デフォルトではこの方法の計算が行われず、引数にadjust=Falseを指定する必要がある。つまりデフォルトでは何らかの調整がされているということだが、どのような調整が行われているのか？

その前に、単一の $y$ で $S$ を初期化することの問題を説明する。

この初期化方法は唯一つの値を置き換えただけのように見えるが、実際のところは代入した $y$ は、それより古い無限に続く数列すべての $S$ を $y$ に置き換えたのと同じことになっている。 $S$ は漸化式で定義されているので、初期値を適当に定めると、初期値以前の無限の数列に対して何らかの仮定が入ってしまうということである。

例として $S_{t-2} = y_{t-2}$ とした場合を考える。

$S_t = \alpha y_t + (1-\alpha)S_{t-1}$
$S_{t-1} = \alpha y_{t-1} + (1-\alpha)S_{t-2}$
$S_{t-2} = y_{t-2}$

このとき $S_t$ を計算すると次のようになる。

$\displaystyle{ S _ t = \alpha\left(y _ t + (1-\alpha)y _ {t-1} + \frac{1}{\alpha}(1-\alpha)^ 2y _ {t-2}\right) }$

ここで、 $\frac{1}{\alpha}$ をテイラー展開すると次のようになる。

$\displaystyle{ \frac{1}{\alpha} = 1 + (1-\alpha) + (1-\alpha)^ 2 + ... }$

すなわち、次のように $y _ {t-2}$ に対して無限級数がかかっていることがわかる。

$\displaystyle{ S _ t = \alpha( y _ t + (1-\alpha)y _ {t-1} + (1-\alpha)^ 2y _ {t-2} + (1-\alpha)^ 3y _ {t-2} + (1-\alpha)^ 4y _ {t-2}...) }$

これの何が問題かというと、初期値に設定した $y _ {t-2}$ の影響が、実際のデータ数に比べて相対的に大きくなってしまうという点にある。

そこで、 $y _ {t-2}$ の重みを相対的に小さくして、データ数に見合った程度の影響力に抑えてやりたいと思うかもしれない。そのような調整をしているのがpandas.DataFrame.ewmのデフォルトの方法である。

最初の方を重み付きの有限区間の移動平均とみなす

どのような調整がされているかはpandas.DataFrame.ewmのドキュメント(pandas.DataFrame.ewm — pandas 1.4.3 documentation)を確認すると書いてある。

要するに、 $y _ t$ の重みは1、 $y _ {t-1}$ の重みは $(1-\alpha)$ 、 $y _ {t-2}$ の重みは $(1-\alpha)^{2}$ として、重み付きの平均をとった値を $S _ t$ とするということだ。

$\displaystyle{S _ t = \frac{y _ t + y _ {t-1}(1 - \alpha) + y _ {t-2}(1-\alpha)^{2}}{1 + (1 - \alpha) + (1 - \alpha)^{2}}}$

元のEMAの定義と比較すると、直近のデータの重みが相対的に大きく、古いデータの重みが相対的に小さくなっていることがわかる。また、データが増えていけば元のEMAの定義に近づいていく（極限においては分母が $1/\alpha$ となるため）。

Rのパッケージでこの計算方法を採用したEMAを計算してくれるものはあるかもしれないが、軽く探した程度では見つからなかった。計算方法はそれほど難しくないので、自前で定義してもそれほど複雑にはならない。pracma::movavgを参考に調整機能を追加した例を示す。

# R
my_ewm <- function(x, a = NULL, adjust = TRUE) {
  nx <- length(x)
  if(is.null(a)) a = 2 / (nx + 1)
  y <- numeric(nx)
  y[1] <- x[1]
  if(adjust){
    a_adj <- 1
    w_sum <- x[1]
  }
  for (k in 2:nx) {
    if(adjust) {
      a_adj <- a_adj + (1 - a)^(k - 1)
      w_sum <- (x[k] + (1 - a) * w_sum)
      y[k] <-  w_sum / a_adj
    } else {
      y[k] <- a * x[k] + (1 - a) * y[k - 1]
    }
  }
  return(y)
}

ただ複雑でないといっても非常にシンプルというほど単純でもなく、実装ミスの危険性もあるのでpandasを使えるならpandasを使ったほうが良い。pandasならdf.ewm(alpha=1/2).mean()などとするだけで済むし、オプションも多い。

また、Pythonの環境が整っているならRからreticulateパッケージを使ってpandasを利用しても良い。標準のデータセットであるAirPassengersを平滑化する例を示す。

# R
library(reticulate)
pd <- import("pandas")
pd$DataFrame$ewm(data.frame(AirPassengers), alpha = 0.5)$mean()

この方法はデータの最初の方で $S$ に対する $y$ の影響が大きくなるので、数列の位置によって $y$ の値の影響の度合いが変わってしまう。そのような調整をしているのだから当然といえば当然だが、これは好ましい場合もあれば好ましくない場合もあるような気がする。

初期値に最初のいくつかの値の平均を使う

EMAはテクニカル分析の文脈でもよく出てくるが、そのような場合にこの定義が使われていることが多いように思う。例えばSMBC日興証券の解説ページ(指数平滑移動平均│初めてでもわかりやすい用語集│ＳＭＢＣ日興証券)はこの説明になっている。

最初の1つだけを無限の数列と仮定してしまうよりは妥当性があるように思われる。

この方法をそのまま実装したRやPythonのパッケージは見当たらなかったが、データ系列の最初に平均値を結合して調整なしのEMAを求めれば同じ結果となる。

最初の方の値を無視する

そもそも調整が必要になるのは数列の最初の方に対応する $S$ を計算する際の $y$ に対する重みが適切ではないというのが原因なので、最初の方の $S$ は無視してしまうというのも手である。英語版Wikipediaでは計算がある程度収束するまでの期間をspin-up期間と呼んでおり、いくつ無視すべきかの目安などが解説されている(Moving average - Wikipedia)。

pandas.DataFrame.ewmのオプションにmin_periodsというものがあり、計算に用いるデータが指定数に満たない場合は結果をnp.nanとするが、おそらくこの目的で利用できる。

ちなみに $\alpha$ が小さい場合を除けば、どのような調整方法をとったとしてもある程度の期間計算すれば結果は一致してくるので、利用できるデータが十分にあれば最初の方をバッサリ切ってしまうこの方法が一番実用的かもしれない。もちろん、 $\alpha$ を小さく設定したい場合には調整方法が重要になるということではあるが。

参考

実践時系列解析 ―統計と機械学習による予測
- p.53あたりから指数平滑法の説明がある。初期値の調整については感覚的な軽い説明がされている。
- 説明は悪くないがRとPythonのコードがちょっと微妙。また、文章に直訳っぽく理解しづらい部分がちょいちょいある。
The correct way to start an Exponential Moving Average (EMA) | David Owen’s blog
- 上記書籍で参考文献として挙げられているブログ記事。初期値の調整について説明された貴重な記事だが、説明が適切かどうかは微妙な気がする。
pandas.DataFrame.ewm — pandas 1.4.3 documentation
- この関数は色々なことができるので、ドキュメントを一読しておくとよい。
Moving average - Wikipedia
- 日本語版の記事は英語版を参考にしていると思われるが、英語版の方が説明が充実している。
【Python】pandasの指数平滑移動平均の値が違った理由 - turtlechanのブログ
- 初期値を最初の数項の平均値とする場合の実装例がある。