http://ill-identified.hatenablog.com/entry/2016/04/02/225155

Rstan_VARMApq.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())

N <- 2              # num series
Time <- 400         # span
Time.forecast <- 50 # forecasting span
psi   <- matrix(c(0.4, 0.3, 0, -0.2), ncol=N)
theta <- matrix(c(0.3, 0.2, 0.1, 0), ncol=N) 
mu    <- matrix(c(1,2), ncol=1)

# ここは VARMA(1,1) と同じ
# generate random numbers
#eps_l <- matrix(mvrnorm(n=1, rep(0,N), 5*diag(N)), nrow=N)
#sample.data.2 <- eps_l + mu
#for( t in 2:Time){
#  eps <- matrix(mvrnorm(n=1, rep(0,N), 2*diag(N)), nrow=N)
#  sample.data.2 <- cbind(sample.data.2, mu + matrix(psi %*% sample.data.2[,t-1] + theta %*% eps_l + eps, nrow=N) )
#  eps_l <- eps
#}

# VARMA(p,q)

# modify plot function
plot.stan.pred.2 <- function(data,stanout, Time, Time.forecast, N, p, q, 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[,Time+1:Time.forecast, 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=paste0("ARMA(", p, ",", q, ") Forecasts by MCMC"), y="y") + facet_wrap(~series, nrow = N)
}

# VARMA(1,1) N=2
N <- 2              # num series
Time <- 400         # span
Time.forecast <- 50 # forecasting span
psi   <- matrix(c(0.4, 0.3, 0, -0.2), ncol=N)
theta <- matrix(c(0.3, 0.2, 0.1, 0), ncol=N) 
mu    <- matrix(c(1,2), ncol=1)

# generate random numbers
eps_l <- matrix(mvrnorm(n=1, rep(0,N), 5*diag(N)), nrow=N)
sample.data.2 <- eps_l + mu
for( t in 2:Time){
  eps <- matrix(mvrnorm(n=1, rep(0,N), 2*diag(N)), nrow=N)
  sample.data.2 <- cbind(sample.data.2, mu + matrix(psi %*% sample.data.2[,t-1] + theta %*% eps_l + eps, nrow=N) )
  eps_l <- eps
}

stan.varmapq <- stan_model(file="varma_pq.stan") # compile
p <- 2
q <- 0
res.stan.3 <- sampling(stan.varmapq, data=list(T=Time, N=N, y=t(sample.data.2), T_forecast=Time.forecast, p=p, q=q), chain=1 )
print(res.stan.3, pars=c("mu","Psi","Theta","Sigma"))
print(plot.stan.pred.2(t(sample.data.2), res.stan.3, Time, Time.forecast, N, p, q, 350:(Time+Time.forecast)))