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

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

# 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
}

# estimate
res.stan.2 <- sampling(stan.varma11, data=list(T=Time, N=N, y=t(sample.data.2), T_forecast=Time.forecast), chain=1 )

# display result
print(res.stan.2, pars=c("mu","Psi","Theta","Sigma"))
print(plot.stan.pred(t(sample.data.2), res.stan.2, Time, Time.forecast, N, 350:(Time+Time.forecast)))