Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment