仕事でたまーにARIMAモデルを推定することがあるんだけど、そのたびに「いま推定してるモデルって、どんな定数項が入っているんだろう…」とちょっと混乱してしまう。じっくり考えればわかるんだけど、あいにく単変量のARIMAモデルというのは、(自分の仕事のなかでは) 深く考えずにちゃちゃっと推定する場面が多く、考えている時間がないことが多い。時計を睨みながらキーボード叩いているようなときに、いちいちマニュアルなど調べたくないじゃないですか。
　というわけで、このたびイライラが募り、メモを取った。

おさらい
　久しぶりなのでおさらいしておくと…
　ARMA(1,1)モデルは $$ Y_t = c + \phi_1 Y_{t-1} + \varepsilon_t + \theta_1 \varepsilon_{t-1} $$ である(以下、\(\varepsilon \sim WN\)とする)。今期の値は前期の値に係数\(\phi_1\)で影響されるし(AR(1)), 前期の撹乱項にも係数\(\theta_1\)で影響されるわけね(MA(1))。第二項を左に持ってきてラグ演算子\(L\)を使うと $$ (1-\phi_1 L) Y_t = c + (1+\theta_1L) \varepsilon_t$$ と書ける。AR(1)の部分が定常であれば $$ Y_t = \frac{c}{1-\phi_1} + \frac{1+\theta_1 L}{1-\phi_1 L} \varepsilon_t$$ とも書ける。定数\(\frac{c}{1-\phi_1}\)は\(Y_t\)の期待値である。

　ARMA(p,q)モデルに一般化しよう。$$ Y_t = c + \phi_1 Y_{t-1} + \cdots + \phi_p Y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}$$ $$ (1-\phi_1 L – \cdots -\phi_p L^p) Y_t = c + (1+\theta_1 L+ \cdots + \theta_q L^q) \varepsilon_t$$ AR(p)の部分が定常であれば $$ Y_t = \frac{c}{1-\phi_1 – \cdots – \phi_p} + \frac{1+\theta_1 L + \cdots + \theta_q L^q}{1-\phi_1 L – \cdots – \phi_p L^p} \varepsilon_t$$ 以下では第1項を\(\mu\)と略記する。これは\(Y_t\)の期待値である。また、第2項をARMA撹乱項と呼び\(E(p,q)_t\)と略記する。その期待値は0である。

　ARIMA(p,1,q)モデルは、\(Y_t\)の1階差分についてのARIMA(p,q)モデルである。$$ (1-L) Y_t = c + \phi_1 Y_{t-1} + \cdots + \phi_p Y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}$$ $$ (1-\phi_1 L – \cdots -\phi_p L^p) (1-L) Y_t = c + (1+\theta_1 L+ \cdots + \theta_q L^q) \varepsilon_t$$ AR(p)の部分が定常であれば $$ (1-L) Y_t = \frac{c}{1-\phi_1 – \cdots – \phi_p} + \frac{1+\theta_1 L + \cdots + \theta_q L^q}{1-\phi_1 L – \cdots – \phi_p L^p} \varepsilon_t = \mu + E(p,q)_t $$ 左辺に残るラグ演算子をとることを考えよう。毎期 \(\mu = \frac{c}{1-\phi_1 – \cdots – \phi_p}\)づつ増大する時系列\(M_t\)を考える。また前期との差分がARMA(p,q)撹乱項\(E(p,q)_t\)になる項をARIMA撹乱項と呼び\(E(p,1,q)_t\)と略記する。すると$$ (1-L) Y_t = (1-L) M_t + (1-L) E(p,1,q)_t$$ すべての項に\((1-L)\)がついているから、消して $$ Y_t = \mu_0 + \mu t + E(p,1,q)_t $$ と書き換えられる。
　\(\mu = \frac{c}{1-\phi_1 – \cdots – \phi_p}\)は\(Y_t\)の期待値ではなくて\((1-L)Y_t\)の期待値である。知らんうちに\(\mu_0\)なんてのが入ってくるところが気持ち悪いですよね。
　めんどくさいので ARIMA(p,d,q) に一般化するのはやめるけど、最初の式に\(c\)を入れたことで、最後の式にはtのd次多項式がはいることになる。

stats::arima()の引数
　ここからは、ARIMAモデルをRで推定するときの話。

　まずはARIMAモデルのためのRの標準的な関数 stats::arima() について。
　モデルの定数に関する引数として、include.meanという引数が用意されている。ドキュメントによればこういう意味らしい。

• ARMAモデルの場合、\(c\)を推定するかどうか。デフォルトはTRUE。
• ARIMAモデルの場合は無視される。

　ここからは私の理解、というか想像になってしまうんだけど…
　ARMAモデルの場合、stats::arima()は\( Y_t = \mu + E(p,q)_t\) を推定する。\(\mu = \frac{c}{1-\phi_1 – \cdots – \phi_p}\)はモデルのパラメータである。
　include.mean = TRUEのとき、stats::arima()が返す係数のなかにinterceptというのがある。これは(\(c\)の推定値じゃなくて) \(\mu\)の推定値だと思う。なお、\(\mu\)は\(Y_{t}\)の期待値の推定値ではあるものの、\(Y_{t}\)の標本平均とぴったり一致するわけではない。そりゃそうですよね、互いに独立ではないARMA撹乱項を想定して最尤推定してるわけだから。
　include.mean = FALSEのとき、stats::arima()は\(\mu\)を0に固定する。

　いっぽうARIMAモデル(d=1)の場合、stats::arima()は \( (1-L) Y_t = E(p,q)_t \)を推定する。だから引数include.meanは意味をなさない。
　ここでちょっと不思議に感じるのは、じゃあstats::arima()が返したARIMAモデルのオブジェクトを使って、(時系列の最後の期を\(T\)として) \(Y_{T+1}\)を予測するときにはどうするの？ という点である。だって、もしオブジェクトがARMA撹乱項のパラメータしか持ってないんなら、\(Y_{T+1}\)の予測なんてできないじゃないですか。値の平均的な大きさについての情報を持っていないと。
　stats::arima()が返すモデルオブジェクトをよーく見ると、出力されるパラメータ推定値のほかにも、なんだかパラメータらしきものが入っている。どうやら、stats::arima()はARIMAモデルを状態空間モデルに変換してカルマンフィルタで推定しており、オブジェクトは状態空間モデルのパラメータと(たぶん第\(T\)期のみの)状態推定値を隠し持っているようである。
　安心してください！持ってますよ！ (名前知らないけど、なんとかさんというタレントさんの真似)

forecast::Arima()の引数
　日頃の仕事ではstats::arima()を使うことはほとんどなくて、Hyndman先生謹製のforecastパッケージに頼っている。便利である反面、引数がややこしい。さあ深呼吸！
　forecast::Arima()には、モデルの定数に関する引数が3つある。まずはドキュメントからメモしておこう。

• include.mean: ARIMAモデルに平均の項を含めるか。デフォルトは、ARMAモデルではTRUE, ARIMAモデルではFALSE。
• include.drift: ARIMAモデルに線形ドリフト項を含めるか(ARIMA誤差つき線形回帰モデルを当てはめるか)。デフォルトはFALSE。
• include.constant: TRUEにすると、ARMAモデルではinclude.mean = TRUEとなり、ARIMAモデルではinclude.drift = TRUEとなる。d>1のとき、この引数をどう指定しようが定数項は含まれない。これは意図的な仕様である(定数項を含めると、二次トレンドやさらに高次な多項トレンドが含まれてしまうことになるから)。

　ぶひー、ややこしい。ここからも私の理解(想像)が入ってしまうけど、Hyndman先生ご自身による解説を手掛かりに、なにをどうするとなにがどうなるのかをまとめておこう。

　include.constant引数はinclude.mean, include.driftの指定を上書きするだけである。また、試してみたところ、ARIMAモデルでは include.meanは単に無視されるようだ。結局、下の組み合わせを考えればよい。

• ARMAモデル(d=0)
  • デフォルト(include.mean = T, include.drift = F) … \(Y_t = \mu + E(p,q)_t\) を推定する。係数として出力される mean は、\(\mu = \frac{c}{1-\phi_1 – \cdots – \phi_p}\)の推定値。
  • include.mean = F, include.drift = F … \( Y_t = E(p,q)_t\) を推定する。
  • include.mean = T, include.drift = T … この辺から自信がなくなってくるんだけど、確定的一次トレンドつきARMAモデル\( Y_t = \beta_0 + \beta t + E(p,q)_t \) を推定しているのではなかろうか。係数として出力される mean と drift は、それぞれ \( \hat{\beta_0}, \hat{\beta} \)なのではないかと思う。
  • include.mean = F, include.drift = T … いったいどんなときに使うモデルなのかわからないが、\( Y_t = \beta t + E(p,q)_t\) を推定しているのではなかろうか。係数として出力される drift は\( \hat{\beta} \) ではないかと思う。
• ARIMAモデル(d=1)
  • デフォルト(include.drift = F) … \( (1-L)Y_t = E(p,q)_t \) を推定している。
  • include.drift = T … \( (1-L)Y_t = \mu + E(p,q)_t\) を推定している。係数として出力されるdriftは \( \hat{\mu} \)だと思う。
• ARIMAモデル(d>1)
  • デフォルト(include.drift = F) … \( (1-L)^d Y_t = E(p,q)_t \) を推定している。
  • include.drift = T … include.drift はFに置き換えられ、警告が表示される。おまえらに確率的二次トレンドは推定させないぞというHyndman先生の強い意志を感じる。

driftってなんのこと？
　… 上記のメモをとってて思ったんだけど、わたくし、「ドリフト」「トレンド」って言葉が何を指しているのか、いっつも混乱するのであります。

　疑問その1。
　こういうのを「ドリフトとトレンドのついたランダムウォーク」っていいませんか？$$ Y_t = \alpha_1 + \alpha_2 t + Y_{t-1} + \varepsilon_t $$ \(\alpha_1\)をドリフト、\(\alpha_2 t\)をトレンドと呼ぶと思う。
　ラグ演算子で書き換えて $$ (1-L)Y_t = \alpha_1 + \alpha_2 t + \varepsilon_t$$ ホワイトノイズ\(\varepsilon_t = E(0,0)_t\)をARMA撹乱項\(E(p,q)_t\)に一般化すれば、$$(1-L)Y_t = \alpha_1 + \alpha_2 t + E(p,q)_t$$ stats::arima()は \(\alpha_1, \alpha_2\)を0に固定する。forecast::Arima()は\(\alpha_2\)を0に固定し、\(\alpha_1\)のことをdriftと呼ぶ。うん、ここまではわかります。

　いっぽう、こういうモデルを「トレンド付き定常ARMAモデル」っていうと思うんですよね。$$ Y_t = \beta_0 + \beta_1 t + E(p,q)_t $$ \(\beta_0 + \beta_1 t\)をトレンドと呼ぶと思う。
　stats::arima()は\(\beta_1\)を0に固定し、\(\beta_0\)だけ推定する。いっぽうforecast::Arima()は両方を推定できる。さて、ここで頭がこんがらがっちゃうのだが、forecast::Arima()は\(\beta_1\)のことをdriftと呼んでいるのである。えええ？ それってトレンドっていいません？

　というわけで… ARMA/ARIMAモデルのパラメータってなに？ と考える際には、モデルの左辺は\((1-L)^d Y_t\)だ、d=0とd=1についてそれぞれ$$ Y_t = \beta_0 + \beta_1 t + E(p,q)_t $$ $$ (1-L) Y_t = \alpha_1 + E(p,q)_t$$ だ、と考えたほうがわかりやすい。そう考えないと、ARIMAモデルで\(Y_t\)の期待値がパラメータでない理由がよくわからなくなる。
　いっぽう、forecast::Arima()がなにをdriftと呼んでいるかを理解するには、モデルの左辺は\(Y_t\)だ、d=0とd=1についてそれぞれ$$ Y_t = \beta_0 + \beta_1 t + E(p,q)_t $$ $$ Y_t = \alpha_0 + \alpha_1 t + E(p,1,q)_t$$だ、と考えたほうがわかりやすい。forecast::Arima()はこの定式化における \(t\)の係数 \(\beta_1, \alpha_1\)をdriftと呼んでいるわけだ。

　疑問その2。
　ARIMA撹乱項を持つ時系列 $$ Y_t = \alpha_0 + \alpha_1 t + E(p,d,q)_t $$はあてどなくドリフト(漂流)する。いっぽう、ARMA撹乱項を持つ時系列 $$ Y_t = \beta_0 + \beta_1 t + E(p,q)_t $$ は傾き\(\beta_1\)の直線の周りをうろうろするだけである。どちらもドリフトと呼ぶのって、語感としてはちょっと変なのではなかろうか？
　いやまてよ、いくらドリフト走行しても道路を外れるわけじゃなし、ザ・ドリフターズは公開録画の人気番組を毎週毎週収録するくらい勤勉な人たちだったわけだし、別にいいのか…

　全然話はちがうけど、いかりや長介さんリーダーのドリフターズって、アメリカの昔のコーラスグループのThe Drifters にちなんだ名前だったのだろうか、全然関係なかったのだろうか。同じ名前を付けるなんて、いまなら叱られちゃいそうですね。