状態空間モデルで季節成分を柔軟に表現する方法

はじめに

こんにちは、事業会社で働くデータサイエンティストです。

業務の一環で時系列予測をすることがあるのですが、精度の高い予測結果を出すことはなかなか難しいものです。今回は、予測に使われる状態空間モデルにおいて、季節成分を柔軟にモデリングする方法について書きました（同僚の@Gotoubun_taiwanからアドバイスをいただきました）。

状態空間モデル

状態空間モデルは状態方程式と観測方程式から構成される。

$t$期の状態を$\mu_t$と表すと、

$\mu_t = \mu_{t-1} + w_t, \quad w_t \sim N(0, \sigma_w^2) \quad \cdots ① \quad (状態方程式)$

$t$期の観測値を$y_t$と表すと、

$y_t = \mu_t ~ + v_t, \quad v_t \sim N(0, \sigma_v^2) \quad \cdots ② \quad (観測方程式)$

ただし、$w_t$と$v_t$はそれぞれ状態誤差と観測誤差である。

季節性のモデリング

季節ダミーによるモデリング

季節性をモデリングする手段の1つに、季節をダミー変数で表現する方法がある。これは、「季節成分の合計が平均的に0になる」という仮定をおくものである。

具体的に考えるため、ここでは四半期データを扱おう。春夏秋冬で季節性をモデリングしたいような場面がこれにあたる。

4つの季節成分を$\gamma_1$, $\gamma_2$, $\gamma_3$, $\gamma_4$と表すと、

$\gamma_1 + \gamma_2 + \gamma_3 + \gamma_4 = s_t, \quad s_t \sim N(0,\sigma_s^2)$

となるような確率変数$s_t$を考えてやればよい。

このとき、例えば4つ目の季節成分は

$\gamma_4 = -(\gamma_1 + \gamma_2$ + $\gamma_3) + s_t$
$~~~~= - \sum_{t=1}^3\gamma_t + s_t$

となる。

一般に、$j$期ある季節成分は、

$\gamma_t = -\sum_{i=t-(j-1)}^{t-1}\gamma_i + s_t$
$~~~~= - \sum_{t=2}^4\gamma_t + s_t, \quad s_t \sim N(0, \sigma_s^2)$

と表現できる。

よって②に$\gamma_t$を追加してやることで、

$y_t = \mu_t ~ + \gamma_t + v_t, \quad v_t \sim N(0, \sigma_v^2) \quad \cdots ②^{'}$

という、季節性を考慮した観測方程式をたてることができる。

月日のモデリング

上記のモデリングの欠点は、全ての時点の季節成分が$N(0, \sigma_s^2)$という単一の正規分布に従うことである。よりきめ細かいモデリングが必要なとき、代替案として月や日の単位で季節性をモデリングすることが考えられる。

例えば、1月~12月の各月をモデリングしようと思えば、以下のようなパラメータを考えることができる。

$\beta_{month~of~year[i]} \sim N(0,1) \quad for \quad i = 1$
$\beta_{month~of~year[i]} \sim N(\beta_{month~of~year[i-1]},1) \quad for \quad i = 2,\cdots 12$

同様にして、1日目～31日目の各日もモデリングできる。
$\beta_{day~of~month[i]} \sim N(0,1) \quad for \quad i = 1$
$\beta_{day~of~month[i]} \sim N(\beta_{day~of~month[i-1]},1) \quad for \quad i = 2,\cdots 31$

（注）ここで各パラメータ$\beta_{month~of~year[i]}$と$\beta_{day~of~month[i]}$の分散は1としているが、分散を$\sigma_{\beta_{month~of~year[i]}}^2$や$\sigma_{\beta_{day~of~month[i]}}^2$という新たなパラメータにしてやればより柔軟なモデリングができるかもしれない。

よって$②^{'}$の$\gamma_t$を$\beta_{month~of~year[i]}$や$\beta_{day~of~month[i]}$で置き換えることで、以下のような観測方程式をたてられる。

$y_t = \mu_t ~ + \beta_{month~of~year[i]} + \beta_{day~of~month[i]} + v_t, \quad v_t \sim N(0, \sigma_v^2) \quad \cdots ②^{''}$

Stanで実装

$②^{''}$で表される観測方程式を含んだ状態空間モデルをStanで実装してみよう。ここでは、以下のモデルを考える。

$\mu_t = \mu_{t-1} + \delta_{t-1} + w_t, \quad w_t \sim N(0, \sigma_w^2) \quad \cdots \quad (状態方程式)$

$y_t = \mu_t ~ + \beta_{month~of~year[i]} + \beta_{day~of~month[i]} + v_t, \quad v_t \sim N(0, \sigma_v^2) \quad (観測方程式)$

ただし$\delta_{t-1}$はトレンド成分で、

$\delta_{t} = \delta_{t-1} + z_t, \quad z_t \sim N(0, \sigma_z^2)$

である。

dataブロックは以下のようになる。

data {
  int time_type;       
  vector[time_type] y;
  int<lower=1, upper=12> month_of_year[time_type];
  int<lower=1, upper=31> day_of_month[time_type];
}

続いてparametersブロック。

parameters {
  vector[time_type] mu;       // 水準（状態）成分
  vector[time_type] delta;    // トレンド成分
  vector[4] beta_month_of_year_unnormalized;  //季節成分（月）
  vector[31] beta_day_of_month_unnormalized;  //季節成分（日）
  real<lower=0> s_w;  // 水準成分の標準偏差
  real<lower=0> s_z;  // トレンド成分の標準偏差
  real<lower=0> s_v;  // 観測誤差の標準偏差
}

transformed parametersブロックも必要になる。

transformed parameters{
  vector[time_type] alpha;
  vector[12] beta_month_of_year = beta_month_of_year_unnormalized - mean(beta_month_of_year_unnormalized);
  vector[31] beta_day_of_month = beta_day_of_month_unnormalized - mean(beta_day_of_month_unnormalized);
  
  for (i in 1:time_type) {
    alpha[i] = mu[i] + beta_month_of_year[month_of_year[i]] + beta_day_of_month[day_of_month[i]];
  }
}

最後にmodelブロックを記述する。

model {
  s_w ~ normal(2, 2);
  s_z ~ normal(0.5, 0.5);
  s_v ~ normal(10, 5);
  
  beta_month_of_year_unnormalized[1] ~ normal(0,1);
  for (i in 2:12) {
    beta_month_of_year_unnormalized[i] ~ normal(beta_month_of_year_unnormalized[i-1], 1);
  }
  
  beta_day_of_month_unnormalized[1] ~ normal(0,1);
  for (i in 2:31) {
    beta_day_of_month_unnormalized[i] ~ normal(beta_day_of_month_unnormalized[i-1], 1);
  }
  
  // 状態方程式に従い、状態が遷移する
  for(i in 2:time_type) {
    mu[i] ~ normal(mu[i-1] + delta[i-1], s_w);
    delta[i] ~ normal(delta[i-1], s_z);
  }
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:time_type) {
    y[i] ~ normal(alpha[i], s_v);
  }

}

以上のようにStanファイルを設定した後は、$month~of~year$や$day~of~month$をlist型にまとめられるよう、Rコードでデータを整形してやればよい。例えば、$yy $-$mm$-$dd$表記の日付が入った$date$変数があるならば、

df_sales |>
  dplyr::mutate(date = as.Date(date)) |>
  dplyr::mutate(month = stringr::str_sub(date, start = 6, end = 7),
                month = as.numeric(month)) |>
  dplyr::mutate(day = stringr::str_sub(date, start = 9, end = 10),
                day = as.numeric(day)) 

という風に処理してやればよい。

参考文献

馬場真哉. 2019. RとStanではじめる ベイズ統計モデリングによるデータ分析入門. 講談社.