ARIMAモデルを導出したい

## はじめに
ARIMAモデルは，メジャーな時系列モデルであるにも関わらず，モデル式を数式で説明している日本語の文献は少ない．
多くの文献では，ARモデルやMAモデル，ARMAモデルまでの説明はあるが，ARIMAモデルになると途端に説明が少なくなる．
SARIMAモデルになるとさらに少なくなる．Wikipediaですら，ARMAモデルまでは詳しく書いてあるが，ARIMAモデルから説明が少なくなる．
詳しい導出はどこかの本に書いているのかしれないが，インターネットで検索する範囲では簡単には見つからなかった．
インターネットで公開されているブログ記事では誤った情報も散見された．
英語の論文や海外の大学の講義資料などに断片的な情報はあるが，「これだけ見れば良い」という文献はあまりなかった．
そこで，この分野は全くの素人ではあるが，自分で導出してみることにした．（私も間違えた情報を書いているかもしれない．）
もし有用な本や記事があれば頂きたく思います．また，記事の査読や誤植などに対するコメントを頂ければ幸いです．

続編（[SARIMAモデルを導出したい](https://qiita.com/suzuki_/items/14998f8120b33a2e5668)）を書きました．

## 結論
たくさん予防線を張ったところで，結論を書く．
この記事のゴールは，ある時刻$t$の時系列データ$y_t$を，過去の時系列データ$y_{t-1}, y_{t-2}, \cdots, y_{t-N}$を用いて書くことである．
結論として，ARIMA$(p,d,q)$モデルは以下の式で書ける．（自分で導出したので，あっている確証はない．）

```math
\begin{eqnarray}
    y_{t} = \sum_{i=1}^p \phi_i y_{t-i} - \sum_{k=1}^d (-1)^{k} \binom dk y_{t-k} + \sum_{i=1}^p \sum_{k=1}^d (-1)^{k} \binom dk \phi_i y_{t-i-k} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

【追記】以下のようにシンプルにできるとコメントを頂きました（ありがとうございます）

```math
\begin{eqnarray}
    y_{t} = \sum_{i=1}^p \sum_{k=0}^d (-1)^{k} \binom dk \phi_i y_{t-i-k} - \sum_{k=1}^d (-1)^{k} \binom dk y_{t-k} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t} 
\end{eqnarray}
```

この記事では一貫して，ベースライン$c=0$とする．
導出してみて分かったが，なぜARIMAモデルの説明が少ないか考えてみると，説明が面倒だからだと思われる．
あとで記載するオペレータを用いた形式であれば，教科書などの紙面を抑えることができるので，わざわざ$y_t = ...$の形式で書く必要もないのだろう．
この式の説明をするためには，ARモデルから順に説明する必要がある．

## ARモデル
- AR$(p)$の式を以下に示す．

```math
\begin{eqnarray}
y_{t} = \sum_{i=1}^{p} \phi_{i} y_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- ここで，$p$はorder，$\phi_i$はARモデルを決めるパラメータ，$\varepsilon_{t}$は$y_{t}$に依存しないノイズとする．例えば，平均0，分散1の正規分布を用いる．

## MAモデル
- MA$(q)$の式を以下に示す．

```math
\begin{eqnarray}
y_{t} = \varepsilon_{t} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i}
\end{eqnarray}
```

- ここで，$q$はorder，$\theta_i$はMAモデルを決めるパラメータとする．
- この記事では，$\theta_i > 0$とするため，$(-\theta_{i})$で記載しているが，$\theta_i < 0$とするならば，$(\theta_{i})$で書ける．

## ARMAモデル
- ARMA$(p,q)$の式を以下に示す．

```math
\begin{eqnarray}
	y_{t} = \sum_{i=1}^{p} \phi_{i} y_{t-i} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

## ARIMAモデル
- ARIMA$(p,d,q)$モデルの式は，ARMA$(p,q)$と$y_{t}$の差分を用いて定義できる．

```math
\begin{eqnarray}
	y^{(d)}_{t} = \sum_{i=1}^{p} \phi_{i} y^{(d)}_{t-i} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- ここで，$d$次の$y_{t}$の差分を以下で定義する．

```math
\begin{eqnarray}
	y^{(d)}_{t} := y^{(d-1)}_{t} - y^{(d-1)}_{t-1}
\end{eqnarray}
```

- 1次と2次の差分を以下で定義する．

```math
\begin{eqnarray}
	y^{(1)}_{t} &:=& y_{t} - y_{t-1} \\
	y^{(2)}_{t} &:=& y^{(1)}_{t} - y^{(1)}_{t-1} = (y_{t} - y_{t-1}) - (y_{t-1} - y_{t-2}) = y_{t} - 2y_{t-1} + y_{t-2}
\end{eqnarray}
```

- ここで$d$次の差分が以下で書けることに瞬時に気づけば，ARIMAモデルの導出は問題ない（無理）

```math
\begin{eqnarray}
	y^{(d)}_{t} = \sum_{k=0}^d \binom dk (-1)^k y_{t-k} 
\end{eqnarray}
```

- これを用いれば，代入して整理すればよい．

```math
\begin{eqnarray}
	\sum_{k=0}^d \binom dk (-1)^k y_{t-k} = \sum_{i=1}^{p} \phi_{i} \sum_{k=0}^d \binom dk (-1)^k y_{t-i-k} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- ただし，問題は両辺に$y_{t-k}$に関する項が出てくるので，式変形が大変．
- 便利な略式記号（作用素・オペレータなどと呼ばれる）を定義することで，式変形を簡単にできる（ただし簡単とは言っていない）

## 後退オペレータ (Backshift operator)
- Backward shift operator $B$を定義する（別名 Lag operator $L$）

```math
\begin{eqnarray}
	B y_{t} &:=& y_{t-1} \\
	B^k y_{t} &:=& y_{t-k}
\end{eqnarray}
```

- このオペレータを使用すると，差分を以下のように書ける．

```math
\begin{eqnarray}
	y^{(1)}_{t} &=& y_{t} - y_{t-1} = y_{t} - B y_{t} = (1-B) y_{t} \\
	y^{(2)}_{t} &=& y_{t} - 2y_{t-1} + y_{t-2} = (1 - 2B + B^2) y_{t} = (1 - B)^2 y_{t}
\end{eqnarray}
```

- $n$次の差分が以下で書ける（証明略）

```math
\begin{eqnarray}
	y^{(n)}_{t} = (1 - B)^n y_{t}
\end{eqnarray}
```

- 以下のような書き方ができる．

```math
\begin{eqnarray}
	(1 - B)(1 - B^m) y_{t} &=& (1 - B - B^m + B^{m+1}) y_{t} \\
         &=& y_{t} - y_{t-1} - y_{t-m} + y_{t-m-1}
\end{eqnarray}
```

## 差分作用素 (Difference operator)
- もう一つオペレータを定義する．
- 差分作用素$\nabla$を以下のように定義する．

```math
\begin{eqnarray}
    \nabla y_t &:=& y_t - y_{t-1} \\
    \nabla^n y_t &:=& y^{(n)}_{t} = y^{(n-1)}_{t} - y^{(n-1)}_{t-1}
\end{eqnarray}
```
- 差分作用素$\nabla$は後退オペレータ$B$で書ける．

```math
\begin{eqnarray}
    \nabla y_t &=& y_t - y_{t-1} = (1 - B) y_t
\end{eqnarray}
```

- 二項定理の展開式を示す．

```math
\begin{eqnarray}
    (x+y)^n = \sum_{k=0}^n \binom nk x^{n-k}y^k \\
    (1-B)^n = \sum_{k=0}^n \binom nk (-B)^k
\end{eqnarray}
```

- 展開式を用いて，差分作用素を式変形する．

```math
\begin{eqnarray}
	\nabla^n y_t = (1 - B)^n y_{t} = \left\{ \sum_{k=0}^n \binom nk (-B)^k \right\} y_t = \sum_{k=0}^n \binom nk (-1)^k y_{t-k}
\end{eqnarray}
```

## オペレータ形式のARMAモデル
- オペレータでARMAモデルを書き直してみる．
- オペレータを用いるとモデル式が簡単に書けることを示す．

#### ARモデル
- AR$(p)$の式を再掲する．

```math
\begin{eqnarray}
	y_{t} = \sum_{i=1}^{p} \phi_{i} y_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- AR$(p)$ operator を以下で定義する．

```math
\begin{eqnarray}
	\phi_p (B) := 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p 
\end{eqnarray}
```

- オペレータに$y_t$を作用させる．

```math
\begin{eqnarray}
	\phi_p (B) y_{t} &=& y_{t} - \phi_1 y_{t-1} - \phi_2 y_{t-2} - \cdots - \phi_p y_{t-p} \\
	&=& y_{t} - \sum_{i=1}^{p} \phi_{i} y_{t-i} \\
	&=& \varepsilon_{t}
\end{eqnarray}
```

- ARモデルはAR$(p)$ operator を用いて以下の式で書ける．

```math
\begin{eqnarray}
	\phi_p (B) y_{t} = \varepsilon_{t}
\end{eqnarray}
```

#### MAモデル
- MA$(q)$の式を再掲する．

```math
\begin{eqnarray}
	y_{t} = \varepsilon_{t} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i}
\end{eqnarray}
```

- MA$(q)$ operator を以下で定義する．

```math
\begin{eqnarray}
	\theta_q (B) := 1 - \theta_1 B - \theta_2 B^2 - \cdots - \theta_q B^q
\end{eqnarray}
```

- オペレータに$\varepsilon_t$を作用させる．

```math
\begin{eqnarray}
	\theta_q (B) \varepsilon_{t} &=& \varepsilon_{t} - \theta_1 \varepsilon_{t-1} - \theta_2 \varepsilon_{t-2} - \cdots - \theta_p \varepsilon_{t-p} \\
	&=& \varepsilon_{t} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} \\
	&=& y_{t}
\end{eqnarray}
```

- MAモデルはMA$(q)$ operator を用いて以下の式で書ける．

```math
\begin{eqnarray}
    y_{t} = \theta_q (B) \varepsilon_{t}
\end{eqnarray}
```

#### ARMAモデル
- ARMA$(p,q)$の式を再掲する．

```math
\begin{eqnarray}
	y_{t} = \sum_{i=1}^{p} \phi_{i} y_{t-i} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- ARMAモデルは operator を用いて以下の式で書ける．

```math
\begin{eqnarray}
    y_{t} = \left\{ y_{t} - \phi_p (B) y_{t} \right\} + \theta_q (B) \varepsilon_{t} \; \Leftrightarrow \; \phi_p (B) y_{t} = \theta_q (B) \varepsilon_{t}
\end{eqnarray}
```

## 準備
- 後で使う式変形を一般的な形式で書く．
- 2つのオペレータを定義する．

```math
\begin{eqnarray}
	a_N (B) &=& 1 - a_1 B - a_2 B^2 - \cdots - a_N B^N \\
	b_M (B) &=& 1 - b_1 B - b_2 B^2 - \cdots - b_M B^M
\end{eqnarray}
```

- ここで，$a_0 = -1$，$b_0 = -1$とすると，以下のように書ける．

```math
\begin{eqnarray}
    a_N (B) &=& 1 - \sum_{i=1}^N a_i B^i = \sum_{i=0}^N (-a_i) B^i \\
    b_M (B) &=& 1 - \sum_{j=1}^M b_j B^j = \sum_{j=0}^M (-b_j) B^j 
\end{eqnarray}
```

- 2つのオペレータの積は以下のように書ける．

```math
\begin{eqnarray}
    a_N (B) b_M (B) = \sum_{i=0}^N \sum_{j=0}^M (-a_i) (-b_j) B^{i+j} = \sum_{i=0}^N \sum_{j=0}^M a_i b_j B^{i+j} 
\end{eqnarray}
```

- 和のスタートを変える式変更を示す（後で使う）

```math
\begin{eqnarray}
    a_N (B) b_M (B) &=& \sum_{i=0}^N \sum_{j=0}^M (-a_i) (-b_j) B^{i+j} \\
    &=& (-a_0) (-b_0) + \sum_{i=1}^N (-a_i) (-b_0) B^{i} + \sum_{j=1}^M (-a_0) (-b_j) B^{j} \\
    &+& \sum_{i=1}^N \sum_{j=1}^M (-a_i) (-b_j) B^{i+j} \\
    &=& 1 + \sum_{i=1}^N (-a_i) B^{i} + \sum_{j=1}^M (-b_j) B^{j} + \sum_{i=1}^N \sum_{j=1}^M (-a_i) (-b_j) B^{i+j} 
\end{eqnarray}
```

- $(-a_i)$のまま計算した理由は，係数の関係を分かりやすくするためである．
- オペレータ$\alpha_N (B)$，$\beta_M (B)$と関数$f(\cdot)$，$g(\cdot)$に対して，オペレータの積が以下のように書ける．（予想）

```math
\begin{eqnarray}
    \alpha_N (B) &=& \sum_{i=0}^N f(\alpha_i) B^i \\
    \beta_M (B) &=& \sum_{j=0}^M g(\beta_j) B^j \\
    \alpha_N (B) \beta_M (B) &=& \sum_{i=0}^N \sum_{j=0}^M f(\alpha_i) g(\beta_j) B^{i+j}
\end{eqnarray}
```

## ARIMAモデルの式を導出する
- ARIMA$(p,d,q)$モデルの途中式を再掲する．

```math
\begin{eqnarray}
	y^{(d)}_{t} = \sum_{i=1}^{p} \phi_{i} y^{(d)}_{t-i} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- 上記の式は，ARMA$(p,q)$に$y_t$に差分作用素を作用させていることになるので，以下の式と同値である．

```math
\begin{eqnarray}
    \phi_p (B) (1 - B)^d y_{t} = \theta_q (B) \varepsilon_{t}
\end{eqnarray}
```

- 左辺を計算する．

```math
\begin{eqnarray}
    \phi_p (B) (1 - B)^d y_{t} &=& \left\{ \sum_{i=0}^p (-\phi_i) B^{i} \right\} \left\{ \sum_{k=0}^d \binom dk (-B)^{k} \right\} y_{t} \\
    &=& \sum_{i=0}^p \sum_{k=0}^d \left\{ (-1)^{k} \binom dk (-\phi_i) \right\} y_{t-i-k}  \\
    &=& y_t + \sum_{i=1}^p (-\phi_i) y_{t-i} + \sum_{k=1}^d (-1)^{k} \binom dk y_{t-k} + \sum_{i=1}^p \sum_{k=1}^d \left\{ (-1)^{k} \binom dk (-\phi_i) \right\} y_{t-i-k} \\
    &=& y_t - \sum_{i=1}^p \phi_i y_{t-i} + \sum_{k=1}^d (-1)^{k} \binom dk y_{t-k} - \sum_{i=1}^p \sum_{k=1}^d \left\{ (-1)^{k} \binom dk \phi_i \right\} y_{t-i-k} 
\end{eqnarray}
```

- 上式の3行目では，和のスタートを変えた．
- スタートを変える理由は式から$y_t$の項を独立させるためである．
- 右辺を計算し，$y_t$以外を移行して整理すると，結論で述べた式が導ける．

```math
\begin{eqnarray}
    y_{t} = \sum_{i=1}^p \phi_i y_{t-i} - \sum_{k=1}^d (-1)^{k} \binom dk y_{t-k} + \sum_{i=1}^p \sum_{k=1}^d \left\{ (-1)^{k} \binom dk \phi_i \right\} y_{t-i-k} + \sum_{i=1}^{q} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

## 具体例で検算する
- 2通りの方法でARIMA$(2,1,1)$を導出し，式が一致することを確かめる．

```math
\begin{eqnarray}
    \phi_2 (B) (1 - B) y_{t} = \theta_1 (B) \varepsilon_{t}
\end{eqnarray}
```

```math
\begin{eqnarray}
    y_{t} = \sum_{i=1}^2 \phi_i y_{t-i} - \sum_{k=1}^1 (-1)^{k} \binom 1k y_{t-k} + \sum_{i=1}^2 \sum_{k=1}^1 \left\{ (-1)^{k} \binom 1k \phi_i \right\} y_{t-i-k} + \sum_{i=1}^{1} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t}
\end{eqnarray}
```

- まずは，オペレータ形式から求める．

```math
\begin{eqnarray}
    && \phi_2 (B) (1 - B) y_{t} = \theta_1 (B) \varepsilon_{t} \\
    & \Leftrightarrow & (1 - \phi_1 B - \phi_2 B^2) (1 - B) y_{t} = (1 - \theta_1 B) \varepsilon_{t} \\
    & \Leftrightarrow &(1 - \phi_1 B - \phi_2 B^2 - B + \phi_1 B^2 + \phi_2 B^3) y_{t} = (1 - \theta_1 B) \varepsilon_{t} \\
    & \Leftrightarrow &\left\{ 1 - (1 + \phi_1) B + (\phi_1 - \phi_2) B^2 + \phi_2 B^3 \right\} y_{t} = (1 - \theta_1 B) \varepsilon_{t} \\
    & \Leftrightarrow &y_{t} = (1 + \phi_1) y_{t-1} - (\phi_1 - \phi_2) y_{t-2} - \phi_2 y_{t-3} + \varepsilon_{t} - \theta_1 \varepsilon_{t-1}
\end{eqnarray}
```

- 次に，導出した式から求める．

```math
\begin{eqnarray}
    && y_{t} = \sum_{i=1}^2 \phi_i y_{t-i} - \sum_{k=1}^1 (-1)^{k} \binom 1k y_{t-k} + \sum_{i=1}^2 \sum_{k=1}^1 \left\{ (-1)^{k} \binom 1k \phi_i \right\} y_{t-i-k} + \sum_{i=1}^{1} (-\theta_{i}) \varepsilon_{t-i} + \varepsilon_{t} \\
    & \Leftrightarrow & y_{t} = \phi_1 y_{t-1} + \phi_2 y_{t-2} - (-1) y_{t-1} + \sum_{i=1}^2 (-1) \phi_i y_{t-i-1} - \theta_{i} \varepsilon_{t-1} + \varepsilon_{t} \\
    & \Leftrightarrow & y_{t} = \phi_1 y_{t-1} + \phi_2 y_{t-2} + y_{t-1} - \phi_1 y_{t-2} - \phi_2 y_{t-3} - \theta_{i} \varepsilon_{t-1} + \varepsilon_{t} \\
    & \Leftrightarrow & y_{t} = (1 + \phi_1) y_{t-1} - (\phi_1 - \phi_2) y_{t-2} - \phi_2 y_{t-3} + \varepsilon_{t} - \theta_{i} \varepsilon_{t-1} 
\end{eqnarray}
```

- 2つの式は一致した．

## 次回
- SARIMAモデル編でお会いしましょう