http://ill-identified.hatenablog.com/entry/2016/02/14/205311

Rstan_VARMA11.R

library(tidyr)    # ver. 0.3.1.
library(dplyr)    # ver. 0.4.3
library(rstan)    # ver. 2.9.0
library(forecast) # ver. 6.2

# 並列化させる設定
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

setwd("~/Documents/blog/20160212_StanVARMA/") # 任意のフォルダに書き換える 
stan.varma11 <- stan_model(file="varma.stan") # compile

N <- 1              # num series
Time <- 400         # span
Time.forecast <- 50 # forecasting span
sample.data <- arima.sim(n = Time, list(ar = c(0.8897), ma = c(-0.2279)),
               sd = sqrt(0.1796))

# plot 80 % & 90 % predition interval of y
plot.stan.pred <- function(data,stanout, Time, Time.forecast, N, span=NULL){
  if (is.null(span) ) span <- 1:(Time + Time.forecast) 
  y_pred <- data.frame()
  for ( i in 1:N){
    m.pred <- rstan::extract(stanout, "y_pred")$y_pred[,,i]
    # 90 % pred interval
    temp <- data.frame(t = Time+1:Time.forecast,
                       series = i,
                       y90_l  = apply(m.pred, 2, quantile, probs=.05),
                       y90_u = apply(m.pred, 2, quantile, probs=.95),
                       y80_l  = apply(m.pred, 2, quantile, probs=.1),
                       y80_u = apply(m.pred, 2, quantile, probs=.9),
                       y_med  = apply(m.pred, 2, median),
                       y_mean = apply(m.pred, 2, mean, rm.na=T)
    )
    y_pred <- rbind(y_pred, temp)
  }
  temp <- data.frame(t=1:Time, data)
  colnames(temp)[1+1:N] <- seq(1:N)
  y_pred <- temp %>% gather(key=series, value=y_mean, -t) %>% mutate(series=as.integer(series)) %>% bind_rows(y_pred)
  y_pred <- y_pred %>% filter(t %in% span)
  ggplot(y_pred) + geom_line(aes(x=t, y=y_mean)) + geom_line(aes(x=t, y_med), linetype=2) +
    geom_ribbon(aes(x=t, ymin=y90_l, ymax=y90_u), alpha=.2, fill="blue") +  geom_ribbon(aes(x=t, ymin=y80_l, ymax=y80_u), alpha=.5, fill="grey") +
    xlim(c(min(span),max(span))) + labs(title="ARMA (1,1) Forecasts by MCMC", y="y") + facet_wrap(~series, nrow = N)
}

# estimate ARMA
res.stan <- sampling(stan.varma11, data=list(T=Time, N=N, y=matrix(sample.data), T_forecast=Time.forecast), chain=1 )
res.auto.arima <- auto.arima(sample.data)
res.arima <- arima(sample.data, order = c(1,0,1))

print(res.stan, pars=c("mu","Psi","Theta","Sigma"))
print(plot.stan.pred(sample.data, res.stan, Time, Time.forecast, N, 350:450))
print(res.auto.arima)
plot(forecast(res.auto.arima, h = Time.forecast), xlim=c(Time-Time.forecast,Time + Time.forecast) )
print(res.arima)
plot(forecast(res.arima, h = Time.forecast), xlim=c(Time-Time.forecast, Time + Time.forecast) )