#install.packages("rstan")
library("rstan")
options(mc.cores=4)

# If installing rstan doesn’t work immediately, here are some details 
# on what to do to get c++ configured to work with stan: 
# https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started


#######################################################
# EXAMPLE 1
#######################################################

# Simulating some data
n     = 100
y     = rnorm(n,1.6,0.2)

par(mfrow=c(1,2))
plot(y)
hist(y)

# Running stan code
model = stan_model("example1.stan")

fit = sampling(model,list(n=n,y=y),iter=1000,chains=4)

print(fit)

params = extract(fit)

par(mfrow=c(2,3))
ts.plot(params$mu,xlab="Iterations",ylab="mu")
acf(params$mu)
hist(params$mu,main="",xlab="mu")
ts.plot(params$sigma,xlab="Iterations",ylab="sigma")
acf(params$sigma)
hist(params$sigma,main="",xlab="sigma")



#######################################################
# EXAMPLE 2
#######################################################
# Simulating some data
n     = 100
mu    = 0
sigma = 1
nu    = 30
y     = mu+sigma*rt(n,df=nu)

par(mfrow=c(1,2))
plot(y)
hist(y)

# Running stan code
model = stan_model("example2.stan")

fit = sampling(model,list(n=n,y=y),iter=2000,chains=4)

print(fit)

params = extract(fit)

par(mfrow=c(3,3))
ts.plot(params$mu,xlab="Iterations",ylab="mu")
acf(params$mu)
hist(params$mu,main="",xlab="mu")
ts.plot(params$sigma,xlab="Iterations",ylab="sigma")
acf(params$sigma)
hist(params$sigma,main="",xlab="sigma")
ts.plot(params$nu,xlab="Iterations",ylab="nu")
acf(params$nu)
hist(params$nu,main="",xlab="nu")




#######################################################
# EXAMPLE 3
#######################################################
# Simulating some data
n     = 200
alpha = 0
beta  = 0.80
gamma = 0.15
sigma = 1
nu    = 4
y = rep(0,n)
for (t in 3:n)
  y[t] = alpha+beta*y[t-1]+gamma*y[t-2]+sigma*rt(1,df=nu)

par(mfrow=c(1,1))
ts.plot(y)

# Running stan code
model = stan_model("example3.stan")

fit = sampling(model,list(n=n,y=y),iter=2000,chains=4)

print(fit)

params = extract(fit)

par(mfrow=c(5,2))
ts.plot(params$alpha,xlab="Iterations",ylab="alpha")
#acf(params$alpha)
hist(params$alpha,main="",xlab="alpha")
ts.plot(params$beta,xlab="Iterations",ylab="beta")
#acf(params$beta)
hist(params$beta,main="",xlab="beta")
ts.plot(params$sigma,xlab="Iterations",ylab="sigma")
#acf(params$sigma)
hist(params$sigma,main="",xlab="sigma")
ts.plot(params$nu,xlab="Iterations",ylab="nu")
#acf(params$nu)
hist(params$nu,main="",xlab="nu")


ts.plot(params$gamma,xlab="Iterations",ylab="gamma")
#acf(params$gamma)
hist(params$gamma,main="",xlab="gamma")