状態遷移の式 α(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])
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
(yと推定結果が重なってます)