注）本記事は途中で詰まったので結局モデル立てただけで結果は出てません。

状態空間モデルの勉強をしたものの、手を動かしていないので練習。

データ

kaggleのrossmann-storeコンペ。

Rossmann Store Sales | Kaggle

カラム名の意味とかはここにまとまっている。

https://adtech.cyberagent.io/techblog/archives/259

アプローチ

練習なので、いったんトレンドと曜日効果のみで考える。

また、これはパネルデータなので、簡易化のためいったんある1店舗について考える。
②で、店舗タイプで縛った階層ベイズをやってみる予定。

モデル

以下のモデルで考える。

トレンド項

2階差分のトレンド項で考える。

$Trend_t = 2\mu_{t-1} - \mu_{t-2} + \sigma_{trend}$

これは、以下の式を変形することで求まる。これは、「当期 - 1期前の差分」と、「1期前 - 2期前の差分」は似ていることから以下のように表すことができるため。
似ているが微小に違う $\sigma_{trend}$ があるため、trendは滑らかに変化するという状態を表すことができている。

$\mu_t - \mu_{t-1} = \mu_{t-1} - \mu_{t-2} + \epsilon_{t-2}$

状態空間モデル

状態方程式

下記のような状態方程式を仮定する。

$SalesState_t \sim Normal(2\mu_{t-1} - \mu_{t-2}, \sigma_{SalesState})$

観測方程式

曜日成分

曜日は7日間の合計が0となる周期で表現する。

$\sum_{l=1}^{7} s_{t-1} = \sigma_{s}$
$∴ s_t = -\sum_{l=1}^{6} s_{t-1} + \sigma_{s}$

また、上式の $s_t$ を用いて以下のような表現で祝日効果、祝前日効果を付与することもできる。

$b_1 (s_{日,t} - s_t)$ ・・・祝日効果:祝日は直近の日曜日と類似しているという仮定。
$b_2 (s_{金,t} - s_t) + b_3 (s_{土,t} - s_t)$ ・・・祝日前効果:祝日前日は直近の金土曜日と類似しているという仮定。
ただしこのとき、 $0 \leq b_1$ 、 $0 \leq (b_2 + b_3) \leq 1 かつ 0 \leq b_2 かつ 0 \leq b_3$

これらをまとめて、祝日要素も入れた曜日効果として、

$Week_t = s_t + d_{1,t} b_1 (s_{日,t} - s_t) + d_{2,t} (b_2 (s_{金,t} - s_t) + b_3 (s_{土,t} - s_t))$
※ $d_{1,t}$ はt日目が月〜金の祝日のとき1、それ以外で0。 $d_{2,t}$ はt日目が祝日でない月〜金かつ翌日が祝日のとき1、それ以外で0となる。

柔軟に曜日の効果を取り入れられる反面、パラメータが多いため時間がかかったり収束がしづらいという欠点がある。今回は簡略化のため1つ目の曜日7日間の周期効果のみ使用する。

プロモーション成分

以下のように定義する。もうちょっと細かく考えると、プロモーション効果は減衰していきそうなのでその点を踏まえたモデルにするとより良さそうだが、いったん定数で考える。 $PromotionEffect_t = C_{promo} * PromotionFlg_t$

閉店日の影響

データは以下のようになっている(見づらいので一部期間のみ）。

df_sales = read_csv('./rossmann-store-sales/train.csv') %>% 
  mutate(Store = as.factor(Store))

df_sales %>% 
  filter(as.numeric(Store) ==1) %>% 
  filter(Date < '2013-04-01') %>% 
  ggplot(aes(Date, Sales, colour = Store)) +
  geom_line() +
  scale_x_date(labels = date_format("%Y/%m/%d"))

f:id:chito_ng:20190618143027p:plain

0となっている部分は閉店日。そのため、閉店日では0になるようする。

また、0部分を除くと以下のようになる。

df_sales %>% 
  filter(as.numeric(Store) ==1) %>% 
  #filter(Sales>0) %>% 
  filter(Date < '2013-04-01') %>% 
  ggplot(aes(Date, Sales, colour = Store)) +
  geom_line() +
  scale_x_date(labels = date_format("%Y/%m/%d"))

f:id:chito_ng:20190618143145p:plain

時系列では1時点前の値を参照することが多いが、今回1時点前が閉店日(Sales = 0) の場合は不適切となる。上図を見る限り、直近の開店日を参照しても問題がなさそうなのでそうすることにする。

以上を踏まえ、下記のような観測方程式を仮定する。

tが閉店日の場合、 $SalesObserve_t = 0$
それ以外の場合、
$SalesObserve_t = SalesState_t + Week_t + PromotionEffect_t + \sigma_{SalesObserve}$
$SalesState_t = 2\mu_{t-1} - \mu_{t-2} + \sigma_{SalesState}$
$Week_t = s_t + d_{1,t} b_1 (s_{日,t} - s_t) + d_{2,t} (b_2 (s_{金,t} - s_t) + b_3 (s_{土,t} - s_t))$
$PromotionEffect_t = C_{promo} * PromotionFlg_t$

stanコード

data {
   int<lower=1> N;
   real SALES[N];
   real<lower=0, upper=1> PROMO_FLG[N];
   real<lower=0, upper=1> OPEN_FLG[N];
}
parameters {
   real<lower=0, upper=10000> trend[N];
   real dow[N];
   
   real<lower=0, upper=100> sigma_trend;
   real<lower=0, upper=100> sigma_dow;
   real<lower=0, upper=100> sigma_sales;
}

transformed parameters {
  real promo_effect[N];
  real c_promo;
  
  
  //プロモーション効果
  for (i in 1:N) {
    promo_effect[i] = c_promo * PROMO_FLG[i];
  }
  
}

model {
   
  // 状態方程式
  for (i in 3:N) {
      trend[i] ~ normal(2*trend[i-1] - trend[i-2], sigma_trend);
    }
    
    
  // 曜日効果
  for (i in 7:N) {
    dow[i] ~ normal(-dow[i-1] -dow[i-2] -dow[i-3] -dow[i-4] -dow[i-5] -dow[i-6], sigma_dow);
  }
  
  // 観測方程式
  for (i in 3:N) {
    if (OPEN_FLG[i] == 0){ // 閉店日
        SALES[i] ~ normal(0,0);  //0
      } else { 
        SALES[i] ~ normal(trend[i] + dow[i] + promo_effect[i], sigma_sales);
      }
    }
    
    
}


generated quantities{
  real sales_mu_all[N+3]; //3期先まで予測
  real sales_pred[N+3];
  real mu_all[N+3];

  for (i in 1:N){ //実測部分まで
    sales_mu_all[i] = trend[i] + dow[i] + promo_effect[i];
    sales_pred[i] = normal_rng(sales_mu_all[i], sigma_sales);
  }

  for (i in N+1:N+3){ //予測部分
    sales_mu_all[i] = normal_rng(sales_mu_all[i-1], sigma_sales); //状態
    sales_pred[i] = normal_rng(sales_mu_all[i], sigma_sales); // 観測.閉店日の対応はしていないので注意
  }

}

結果

Chain 1: Rejecting initial value: Chain 1: Error evaluating the log probability at the initial value. Chain 1: Exception: normal_lpdf: Location parameter is nan, but must be finite! (in 'model71f5bfe6db5_state_space_model' at line 46)

Chain 1: Chain 1: Initialization between (-2, 2) failed after 100 attempts. Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model. [1] "Error in sampler$call_sampler(args_listi) : Initialization failed."

動かねーｗ
初期値が問題っぽいので、パラメータに制限をかけるべきなんだろうけど皆目検討つかない。てか無情報のdow[1~6], trend[1,2]あたりがなんとなく原因な気がする。