2種類の季節変動をstanで推定

状態遷移の式 α(t+1)=T*αt+ε  (記号はKFASから) でTの中に2種類のQが入るだけなので、何も問題がないはずです。

nn=3
a=rep(c(1,2,3),5*nn)
b=rep(c(5,5,0,0,0),3*nn)
d=10+a+b
ts.plot(d[1:30])

f:id:biones:20151019181030j:plain

scode="
data{
  int N;
  real y[N];
  
}
parameters{
  real sigma;
  real trend[N];
  real s_trend;

  real season[N];
  real s_season;

  real season2[N];  #2つ目の季節変動
  real s_season2;
}

model{
  real ssum;
  real mu[N];
  
  for(t in 2:N){
    trend[t]~normal(trend[t-1],s_trend);
  }

  for(t in 5:N){
    ssum<-0;
    for(k in 1:4){
      ssum<-ssum+season[t-k];
    }
    season[t]~normal(-ssum,s_season);
  }

  for(t in 3:N){
    ssum<-0;
    for(k in 1:2){
      ssum<-ssum+season2[t-k];
    }
    season2[t]~normal(-ssum,s_season2);
  }

  for(t in 1:N){
    mu[t]<-season[t]+trend[t]+season2[t];
    y[t]~normal(mu[t],sigma);
  }

}

"
library(rstan)
N=length(d)
slist=list(y=d,N=N)
fit=stan(model_code = scode,data=slist ,iter=1000,warmup=300,chains=1)

#stanの後処理
res=stan.ss(fit,y=d,varnames = c("season","season2"))

r=d-(res$trend+res$season+res$season2)
#> mean(r^2)
#[1] 3.102593e-21

f:id:biones:20151019180832j:plain (yと推定結果が重なってます)

f:id:biones:20151019180812j:plain