【状態空間モデル】KFAS::SSMcustom で数式と対応させたコーディング

2022/01/222023/10/27R,データサイエンス,時系列分析,統計

R で状態空間モデルを使えるパッケージの一つとして KFAS があります。こちらでは様々なモデルが想定されており、例えばローカル線形トレンドモデルであれば SSModel(y ~ SSMtrend(degree = 2, Q = list(NA, NA)), H = NA) というようにモデリングします。この SSModel() の中の formula（この例だと SSMtrend）をいじれば様々なモデルが使えますが、便利な反面この方法だと手元の数式とプログラムとの対応関係がいまいちわかりにくいです。

本記事では SSMcustom() を使って、モデルの係数や変換方法をそのままプログラムに落とし込むような方法で状態空間モデルを推定していきます。

なお、本記事では「モデリング」「モデルの設計」を以下のように定義します。

• モデリング：入力、出力の関係を数式で表し、それに付随する未知パラメータを推定し、予測値を計算することができるようにするまでの一連の過程
• モデルの設計：入力、出力の関係式を定義すること

状態空間モデルの設計について

まず、状態空間モデルの設計にあたって、状態空間モデルの一般形を確認しながら、モデル中のどこをいじるのかを確認します。

状態空間モデルの一般形は以下の通りです。

観測方程式：

$$
y_t = Z_t \alpha_t + \epsilon_t, \quad \epsilon_t \sim N(0, H_t)
$$

状態方程式：

$$
\alpha_t = T_t \alpha_{t-1} + R_t \eta_t, \quad \eta_t \sim N(0, Q_t) \\
\alpha_1 \sim N(a_1, P_1) \\
t = 1, \ldots, n.
$$

それぞれの次元をまとめると以下の通りです。

このうち、\(y_t\) は観測値で \(\alpha_t\) は推定対象なので、モデルの設計にあたり調整する部分は \(Z_t, H_t, T_t, R_t, Q_t\) となります。なお、本記事では
\(Z_t, T_t, R_t\)
の3つをいじっていきます（これだけでも様々な系列に対応できます）。

KFAS の使い方

KFAS によるモデリングの基本的な流れは以下の通りです。

1. モデルの設計（数式の決定）：SSModel(), SSMcustom()
2. ハイパーパラメータの推定：fitSSM()
3. フィルタリング、スムージング：KFS()
4. 予測：predict()

1については冒頭で述べたとおり様々な関数が用意されていますが、今回は SSMcustom() 縛りでいきます。

help() で SSModel() と SSMcustom() の引数を見てみると

SSModel(formula, data, H, u, distribution, tol = .Machine$double.eps^0.5)
SSMcustom(Z, T, R, Q, a1, P1, P1inf, index, n = 1, state_names = NULL)

となっており、上記の数式では引数の文字と対応する表記にしています。SSModel() の formula は y ~ X という形をとり、チルダの左側に観測値、右側にモデル情報が入ったリスト（SSM~で生成）を入力します。

既出でない文字の説明は以下の通りです。

• u：正規分布を仮定しない場合に設定するパラメータ。H、u のどちらかを決める
• distribution：観測値に応じて仮定する分布（正の整数なら poisson とか）
• tol：一期先予測誤差の分散が0かどうかを決める閾値（あまり気にしなくて○）
• a1, P1：状態変数の初期分布の平均値と分散
• P1inf：散漫な初期分散を設定する状態変数の位置を指定するもの。「散漫な」とは、初期の分布は未知であるため「わからないから分散をめちゃくちゃでかくしとこう」というアイデア
• index：観測値が多次元である時、あるモデルをどの系列に適用するかを指定するインデックス。デフォルトでは全系列に適用するように設定されてる
• state_names：状態変数に新たにつける名前

基本的に、状態変数の初期分散はわからないものとして、P1inf のうち該当する状態変数のインデックスの部分を1として「散漫初期化」をおこないます（デフォルト）。一方、過去の研究やデータ自体の性質から、状態変数の初期値に関して妥当と考えられる値を設定できるのであれば、a1, P1 を手動で設定するのが良いと思います。こういった調整は得られた観測値の時点数（サンプルサイズ）が小さい場合に特に有効です。本記事では状態について散漫初期化を行い、観測誤差分散や状態誤差分散は全てデータから推定するものとします。

以上を踏まえ、2種類のモデリングし、推定されたハイパーパラメータ、フィルタ化推定量、平滑化状態と平滑化状態の信頼区間を確認していきます。

なお、今回使用するデータはこちらの記事で取得した、私のタイピングのスコアデータです。

各指標を日毎の平均にしたものからスコアを計算し直して、欠損日のスコアを NA としたものを分析対象データとしておきます。

library(dplyr); library(KFAS)
library(tidyverse)
# 読み込みと名前設定
path_data = 'typing_score.csv'
data_original = read_csv(path_data, local = locale(encoding = 'cp932'))
colnames(data_original) = c('date', 'score', 'level', 'time', 'num_correct', 'num_miss', 'wpm', 'accuracy')

# 日別で平均とって、定義通りスコアを計算
data_mean = data_original %>% 
    group_by(date) %>% 
     summarise(time = mean(time), 
        num_correct = mean(num_correct), 
        num_miss = mean(num_miss)) %>% 
    mutate(score = 60 * (num_correct - num_miss) / time * ((num_correct - num_miss) / num_correct) ^ 2, 
         date = as.Date(date)) %>%
    dplyr::select(date, score) %>%
    arrange(date)

# 欠損してる日付を埋めて、その日のスコアを NA とする
fulldate = seq(min(data_mean$date), max(data_mean$date), by = 'days') %>% 
  as_tibble()
colnames(fulldate) = 'date'
data = full_join(fulldate, data_mean, by = 'date')

ローカル線形トレンドモデル

$$
Z_t = \begin{pmatrix}1 & 0 \end{pmatrix},
\alpha_t = \begin{pmatrix}
\mu_t \\
\beta_t
\end{pmatrix},
H_t = \sigma_o^2, \\
T_t = \begin{pmatrix}
1 & 1 \\
0 & 1
\end{pmatrix},
R_t = I_2 := \begin{pmatrix}1 & 0\\ 0 & 1\end{pmatrix},
\eta_t = \begin{pmatrix}\eta_{t,1} \\ \eta_{t,2} \end{pmatrix},
Q_t = \begin{pmatrix}
\sigma_1^2 && 0 \\
0 && \sigma_2^2
\end{pmatrix}
$$

観測方程式：

$$
\begin{eqnarray} y_t &=& \begin{pmatrix}1 & 0 \end{pmatrix} \begin{pmatrix} \mu_t \\ t_t \end{pmatrix} + \epsilon_t \\
&=& \mu_t + \epsilon_t \end{eqnarray}
$$

状態方程式：

$$
\begin{eqnarray}
\begin{pmatrix} \mu_t \\ \beta_t \end{pmatrix} &=& \begin{pmatrix}1 & 1 \\ 0 & 1\end{pmatrix} \begin{pmatrix}\mu_{t-1} \\ \beta_{t-1} \end{pmatrix} + I_2 \begin{pmatrix}\eta_{t,1} \\ \eta_{t,2} \end{pmatrix}\\
&=& \begin{pmatrix} \mu_{t-1} + \beta_{t-1} + \eta_{t,1} \\ \beta_{t-1} + \eta_{t,2} \end{pmatrix}
\end{eqnarray}
$$

R 実装

# モデルの設計から最尤推定、フィルタリング、平滑化まで
Zt = matrix(c(1,0), nrow = 1)
Ht = matrix(NA)
Tt = matrix(c(1,0,1,1), nrow = 2)
Rt = diag(1, nrow = 2)
Qt = diag(NA, nrow = 2)
P1inf = diag(1, nrow = 2)
# 定数項を入れないため、formula の右辺の先頭に -1 + をつけます
model_LocalTrend = SSModel(
    data$score ~ -1 + SSMcustom(
        Z = Zt, T = Tt, R = Rt, Q = Qt, P1inf = P1inf, 
        state_names = c('level', 'trend')
    ), 
    H = Ht
)

fit_LocalTrend = fitSSM(model_LocalTrend, inits = c(1,1,1), method = "BFGS")
result_LocalTrend = KFS(fit_LocalTrend$model)

# データ整形と可視化
pred_LocalTrend = predict(fit_LocalTrend$model, interval = 'confidence', level = 0.95) %>% 
    as_tibble() %>% 
    mutate(date = data$date, score = data$score, filter = result_LocalTrend$att[,'level']) 

ggplot(data = pred_LocalTrend, aes(x = data$date)) + 
    geom_point(aes(y = score), lwd = 1.5, alpha = 0.6) +
    geom_line(aes(y = fit), col = '4', lwd = 1.5, alpha = 0.6) + 
    geom_line(aes(y = filter), col = '3', lwd = 1.5, alpha = 0.6) + 
    geom_ribbon(aes(ymin = lwr, ymax = upr), col = 4, fill = 4, alpha = 0.15) 

cat(" sigma_o^2: ", result_LocalTrend$model$H, "\n",
    "sigma_1^2: ", result_LocalTrend$model$Q[1,1,1], "\n", 
    "sigma_2^2: ", result_LocalTrend$model$Q[2,2,1], "\n")

 sigma_o^2:  95.33351 
 sigma_1^2:  72.88161 
 sigma_2^2:  0.0006768827 

薄い水色部分は95%信頼区間であり、1月上旬の欠損が続いている部分は区間が大きくなっていることがわかります。

分散を見るとトレンド成分 \(\beta_t\) の分散 \(\sigma_2^2\) の推定値が非常に小さく、トレンド成分を確定的ににしてもそれほどモデルの説明性は変わらなそうです。

季節成分モデル

毎週火曜日に入力文字のお題が変わるため、7日おきに下落→慣れによるスコア上昇がありそうなので、水準＋季節成分で観測値を説明してみます。

$$
Z_t = \begin{pmatrix}1 & 1 & 0 & 0 & 0 & 0 & 0\end{pmatrix},
\alpha_t = \begin{pmatrix}
\mu_t \\
\gamma_{t,1} \\
\gamma_{t,2} \\
\gamma_{t,3} \\
\gamma_{t,4} \\
\gamma_{t,5} \\
\gamma_{t,6}
\end{pmatrix},
H_t = \sigma_o^2,
$$

$$
T_t = \begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & -1 & -1 & -1 & -1 & -1 & -1 \\
0 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 0
\end{pmatrix}
$$

$$
R_t = \begin{pmatrix}1 &0 \\ 0 & 1\\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0\end{pmatrix},
\eta_t = \begin{pmatrix}\eta_{t,1} \\ \eta_{t,2} \end{pmatrix}
Q_t = \begin{pmatrix} \sigma_1^2 & 0 \\ 0 & \sigma_2^2 \end{pmatrix}
$$

観測方程式（tex 辛いので省略）：

$$
\begin{eqnarray} y_t &=& \mu_t + \gamma_{t,1} + \epsilon_t \end{eqnarray}
$$

状態方程式（数式バグの多さにブチギレて別で数式生成してコピペ）：

R 実装

Zt = matrix(c(1,1,0,0,0,0,0), nrow = 1)
Ht = matrix(NA)
Tt = rbind(c(1,rep(0,6)), c(0,rep(-1, 6)), cbind(rep(0,5), diag(5), rep(0, 5)))
Rt = rbind(diag(1,2), matrix(0, nrow = 5, ncol = 2))
Qt = diag(NA, nrow = 2)
P1inf = diag(1, nrow = 7)
model_season = SSModel(
    data$score ~ -1 + SSMcustom(
        Z = Zt, T = Tt, R = Rt, Q = Qt, P1inf = P1inf, 
        state_names = c('level', str_c('season', seq(1,6)))
    ), 
    H = Ht
)

fit_season = fitSSM(model_season, inits = c(1,1,1), method = "BFGS")
result_season = KFS(fit_season$model)

pred_season = predict(fit_season$model, interval = 'confidence', level = 0.95) %>% 
    as_tibble() %>% 
    mutate(date = data$date, score = data$score, filter = result_season$att[,'season1'] + result_season$att[,'level']) 
ggplot(data = pred_season, aes(x = data$date)) + 
    geom_point(aes(y = score), lwd = 1.5, alpha = 0.6) +
    geom_line(aes(y = fit), col = '4', lwd = 1.5, alpha = 0.6) + 
    geom_line(aes(y = filter), col = '3', lwd = 1.5, alpha = 0.6) + 
    geom_ribbon(aes(ymin = lwr, ymax = upr), col = 4, fill = 4, alpha = 0.15) 

ローカル線形トレンドと比べると信頼区間がカクカクしています。

cat(" sigma_o^2: ", result_season$model$H, "\n",
    "sigma_1^2: ", result_season$model$Q[1,1,1], "\n", 
    "sigma_2^2: ", result_season$model$Q[2,2,1], "\n")

 sigma_o^2:  38.34644 
 sigma_1^2:  74.75404 
 sigma_2^2:  2.327765 

# 季節成分
plot(result_season$alphahat[,'season1'])

季節成分を見てみると、「テーマが変わるとスコアが下落し、徐々に慣れてスコアが上がる」という特徴をうまく捉えられているように見えます。

終わりに

この記事を書いてる途中、小野滋氏による上位互換資料を見つけてしまいましたので、より詳しく幅広く学習したい方はこちらをご参照ください。

参考

• KFAS reference manual
• KFAS vignette
• 小野滋. VII. 状態空間モデルの基礎.