SARIMAモデルを導出したい 

## はじめに
前回の記事（[ARIMAモデルを導出したい](https://qiita.com/suzuki_/items/d0322bf7560ea115c1ef)）の続編です．
前回の記事の内容の続きから書きます．
記事の査読や誤植などに対するコメントを頂ければ幸いです．

この分野は全くの素⼈ではあるが，⾃分でSARIMAモデルの式を導出してみることにした．
この記事のゴールは，ある時刻$t$の時系列データ$y_t$を，過去の時系列データ$y_{t-1}, y_{t-2}, \cdots, y_{t-N}$を用いて書くことである．

## 結論
SARIMA$(p,d,q) \times (P,D,Q)_m$のモデル式は，以下の式で書ける．（⾃分で導出したので，あっている確証はない．）

```math
\begin{eqnarray}
y_{t} &=& - \sum_{i=1}^p \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
&& + \sum_{i=1}^p \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i y_{t-i-k-ml} \\
&& + \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& - \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} y_{t-k-ml} \\
&& + \sum_{i=1}^q \sum_{j=1}^Q \theta_i \Theta_j \varepsilon_{t-i-mj} - \sum_{i=1}^q \theta_i \varepsilon_{t-i} - \sum_{j=1}^Q \Theta_j \varepsilon_{t-mj} + \varepsilon_{t}  
\end{eqnarray}
```

導出してみて分かったが，なぜSARIMAモデルの説明が少ないか考えてみると，説明がものすごく⾯倒だからだと思われる．
ただし，ARIMAモデルを導出するためのスキームを理解した後であれば，ただの計算であることが分かった．
あとSARIMAモデルについて教えて」と⾔われたら，このURLを送ることができるのが良かった．

## 前回のあらすじ
- 1次の差分と$d$次の差分を以下のように定義する．

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

- 後退作⽤素$B$と差分作⽤素$\nabla$を以下のように定義する．

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

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

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

- ARMA$(p,q)$のモデル式はオペレータ形式で以下のように書ける．

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

- ARIMA$(p,d,q)$のモデル式はオペレータ形式で以下のように書ける．

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

- オペレータ$\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}
```

## 準備
- 差分作⽤素の定義を以下のように拡張する．

```math
\begin{eqnarray}
    \nabla^n_m y_t := (1 - B^m)^n y_t
\end{eqnarray}
```

- 以下のように書ける．

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

## SARIMAモデル
- SARIMAモデルは，ARIMAモデルに対して周期的なデータの変動を追加したものである．
- 周期を$m$とすると，SARIMAモデルは7個の次数$(p,d,q) \times (P,D,Q)_m$で指定される．
- 通常のARIMA$(p,d,q)$と$m$時刻ごとに過去のデータを抽出したARIMA$(P,D,Q)$を重ね合わせたモデルである．
- SARIMA$(p,d,q) \times (P,D,Q)_m$のモデル式は，オペレータを⽤いて，以下のように書ける．

```math
\begin{eqnarray}
    \phi_p (B) \Phi_P (B) \nabla^d \nabla^D_m y_{t} = \theta_q (B) \Theta_Q (B) \varepsilon_{t} 
\end{eqnarray}
```

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

- ここで，$\Phi_P (B)$と$\Theta_Q (B)$は，以下のように定義される．

```math
\begin{eqnarray}
    \Phi_P (B) &:=& 1 - \Phi_1 B^m - \Phi_2 B^{2m} - \cdots  - \Phi_P B^{Pm} = \sum_{i=0}^P (-\Phi_i) B^{mi} \\
    \Theta_Q (B) &:=& 1 - \Theta_1 B^m - \Theta_2 B^{2m} - \cdots  - \Theta_Q B^{Qm} = \sum_{i=0}^Q (-\Theta_i) B^{mi}
\end{eqnarray}
```

- 式から分かる通り，周期$m$で時系列データを取り出したときに，その時系列データがARIMAモデルに従うようなオペレータを$y_t$に作⽤させていることに相当する．
- ただし，$\Phi_0 = -1$，$\Theta_0 = -1$とする．

## 簡単な例
- SARIMAモデルを導出する前に，モデルの意味をイメージしやすいように，簡単な例をあげる．

**$(0,0,0) \times (P,0,Q)_m$ order SARIMA model**

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

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

- $m$周期ごとに過去のデータを抽出した時系列データがARMAモデルに従っていることが分かる．

**$(p,0,q) \times (0,1,0)_m$ order SARIMA model**

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

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

- $y_{t-m}$の項があることで，$m$時刻前の時系列データに強く関係していることが分かる．
- これにより，周期性が発⽣する．
- そのため，周期性のある時系列データを⽣成するためには，$D>0$に設定する必要がある．

## SARIMAモデルの式を導出する
- SARIMA$(p,d,q) \times (P,D,Q)_m$のモデル式は，オペレータを⽤いて，以下のように書ける．

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

- 基本的には，ARIMAモデルの導出と同じである．
- 左辺を計算する．

```math
\begin{eqnarray}
    ({\rm left}) &=& \phi_p (B) \Phi_P (B) (1 - B)^d (1 - B^m)^D y_{t} \\
    &=& \left\{\sum_{i=0}^p \sum_{j=0}^P (-\phi_i) (-\Phi_j) B^{i+j} \right\} \left\{ \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} B^{k+ml} \right\} y_t \\
    &=& \sum_{i=0}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} 
\end{eqnarray}
```

- 右辺を計算する．

```math
\begin{eqnarray}
    ({\rm right}) = \theta_q (B) \Theta_Q (B) \varepsilon_{t} = \sum_{i=0}^q \sum_{j=0}^Q \theta_i \Theta_j \varepsilon_{t-i-mj}
\end{eqnarray}
```

- 結果として以下の式で書ける（もうここでゴールでいいんじゃないか）

```math
\begin{eqnarray}
    \sum_{i=0}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} = \sum_{i=0}^q \sum_{j=0}^Q \theta_i \Theta_j \varepsilon_{t-i-mj}
\end{eqnarray}
```

- 上式を変形するにあたり，以下のテクニックを⽤いる．

```math
\begin{eqnarray}
    \sum_{i=0}^p f(i) = f(0) + \sum_{i=1}^p f(i) 
\end{eqnarray}
```

```math
\begin{eqnarray}
    \sum_{i=0}^p \sum_{j=0}^P \sum_{k=0}^d f(i,j,k) = \sum_{j=0}^P \sum_{k=0}^d f(0,j,k) + \sum_{i=1}^p \sum_{j=0}^P \sum_{k=0}^d f(i,j,k)
\end{eqnarray}
```

- 左辺を計算する．

```math
\begin{eqnarray}
({\rm left}) &=& \sum_{i=0}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
%
&=& - \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& + \sum_{i=1}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
%
&=& \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} y_{t-k-ml} \\
&& - \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& + \sum_{i=1}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
%
&=& \sum_{l=0}^D \binom Dl (-1)^{l} y_{t-ml} \\
&& + \sum_{k=1}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} y_{t-k-ml} \\
&& - \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& + \sum_{i=1}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
%
&=& y_{t} + \sum_{l=1}^D \binom Dl (-1)^{l} y_{t-ml} \\
&& + \sum_{k=1}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} y_{t-k-ml} \\
&& - \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& + \sum_{i=1}^p \sum_{j=0}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
%
&=& y_{t} + \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} y_{t-k-ml} \\
&& - \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& - \sum_{i=1}^p \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i y_{t-i-k-ml} \\
&& + \sum_{i=1}^p \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml}
\end{eqnarray}
```

- 右辺を計算する．

```math
\begin{eqnarray}
    ({\rm right}) &=& \sum_{i=0}^q \sum_{j=0}^Q \theta_i \Theta_j \varepsilon_{t-i-mj} \\
    &=& - \sum_{j=0}^Q \Theta_j \varepsilon_{t-mj} + \sum_{i=1}^q \sum_{j=0}^Q \theta_i \Theta_j \varepsilon_{t-i-j} \\
    &=& \varepsilon_{t} - \sum_{j=1}^Q \Theta_j \varepsilon_{t-mj} + \sum_{i=1}^q \sum_{j=0}^Q \theta_i \Theta_j \varepsilon_{t-i-j} \\
    &=& \varepsilon_{t} - \sum_{j=1}^Q \Theta_j \varepsilon_{t-mj} - \sum_{i=1}^q \theta_i \varepsilon_{t-i} + \sum_{i=1}^q \sum_{j=1}^Q \theta_i \Theta_j \varepsilon_{t-i-mj} 
\end{eqnarray}
```

- 結果として以下の式が導ける．

```math
\begin{eqnarray}
y_{t} &=& - \sum_{i=1}^p \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i \Phi_j y_{t-i-mj-k-ml} \\
&& + \sum_{i=1}^p \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \phi_i y_{t-i-k-ml} \\
&& + \sum_{j=1}^P \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} \Phi_j y_{t-mj-k-ml} \\
&& - \sum_{k=0}^d \sum_{l=0}^D \binom dk \binom Dl (-1)^{k+l} y_{t-k-ml} \\
&& + \sum_{i=1}^q \sum_{j=1}^Q \theta_i \Theta_j \varepsilon_{t-i-mj} - \sum_{i=1}^q \theta_i \varepsilon_{t-i} - \sum_{j=1}^Q \Theta_j \varepsilon_{t-mj} + \varepsilon_{t}  
\end{eqnarray}
```

- ただし，各orderがゼロでない場合に限る．それ以外の場合は関連する項を0にする．

## 検算
- したくない