Skip to content

Instantly share code, notes, and snippets.

@Gedevan-Aleksizde
Last active April 12, 2016 13:30
Show Gist options
  • Save Gedevan-Aleksizde/9e9eb0ff9a1dd53533bf to your computer and use it in GitHub Desktop.
Save Gedevan-Aleksizde/9e9eb0ff9a1dd53533bf to your computer and use it in GitHub Desktop.
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) )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment