###########################################
# Model BVAR z rozkładem apriori Sims-Zha #
###########################################
    
library(BMR)

#library(AER)
#data("USMacroG")
#y <- USMacroG[,c(1,6,8)]; # PKB, CPI i tbill dla USA
#z[,1] <- log(y[,1]); 
#z[,2] <- log(y[,2]);
#z[,3] <- log(1+y[,3]/100);    

data(BMRVARData)
y <- USMacroData[,c(2:4)];
#This is an updated version of the Stock and Watson (2001) dataset, with inflation, unemployment,
#and the Federal Funds rate, from Q2 1954 to Q4 2011.
z[,1] <- y[,2]; 
z[,2] <- y[,1];
z[,3] <- 100*log(1+y[,3]/100);    

# let us start with classical VAR
z1 <- CVAR(z,p=2,constant=TRUE,irf.periods=20,boot = 1000)
t(z1$Beta)

# Bayesian VAR with Minnesota prior.
z2 <- BVARM(z,coefprior=c(1,1,1),p=2,constant=TRUE,irf.periods=20, 
keep=10000,burnin=5000,VType=1,HP1=0.2,HP2=0.2,HP3=1)

# coefprior - vector of length m prior for own first-lags
# p 		- the number of lags 
# keep 	- the number of Gibbs sampling replications to keep 
# burnin 	- burn-in length for the Gibbs sampler
# VType    	- type of prior 
# HP1 - tightness for own lags
# HP2 - tightness for other lags
# HP3 - tightness for constant

t(z2$Beta)

# Bayesian VAR with Normal-Wishart prior

z3 <- BVARW(z,cores=1,coefprior=c(1,1,1),p=2,constant=TRUE,
irf.periods=20,keep=10000,burnin=1000,
XiBeta=1,XiSigma=1,gamma=NULL)

# cores  - the number of CPU cores used for the sampling run.
# coefprior, p, keep ... - see above
# XiBeta  - overall tightness or location of covariance matrix for coefficients
# XiSigma - location matrix of the inverse-Wishart prior.
# gamma - prior degrees of freedom of the error covariance matrix (minimum/default value is m+1)
# look at page 12 of BMR manual

t(z3$Beta)

# Steady-state Bayesian VAR (Vilanni)

z4 <- BVARS(z,psiprior=c(2,5,5),coefprior=c(0.9,0.9,0.9),p=2,
irf.periods=20,keep=10000,burnin=1000,
XiPsi=1,HP1=0.5,HP4=2,gamma=NULL)

# psiprior - vector of length m with steady-state prior mean (Psi)
# XiPsi    - tightness for Psi.
# H_1      - overall tightness
# HP4      - decay hyperparameter
# gamma - prior degrees of freedom of the error covariance matrix (minimum/default value is m+1)

z4$Psi
t(z4$Beta)

# plot of posterior for coefficients 
# set save=T and look for figures in My documents folder
plot(z2, type=2, save=F)
plot(z3, type=2, save=F)
plot(z4, type=2, save=F)


########
# IRF	 #
########

# Clasical VAR
IRF(z1,percentiles=c(0.05,0.5,0.95), save=FALSE)

# Clasical Minesotta BVAR
IRF(z2,percentiles=c(0.05,0.5,0.95), save=FALSE)

# Bayesian VAR with Normal-Wishart prior
IRF(z3,percentiles=c(0.05,0.5,0.95), save=FALSE)

# Bayesian VAR with steady-state prior
IRF(z4,percentiles=c(0.05,0.5,0.95), save=FALSE)

###############
# Forecasts   #
###############

# Clasical VAR
# does not take into account uncertainty due to parameter estimates
forecast(z1,periods=20,plot=TRUE,confint=0.95,backdata=20)

# Minessota prior BVAR
# try with shocks=TRUE / FALSE
forecast(z2,periods=20,shocks=TRUE,plot=TRUE,percentiles=c(.05,.50,.95),useMean=FALSE,
backdata=20)

# Bayesian VAR with Normal-Wishart prior
forecast(z3,periods=20,shocks=TRUE,plot=TRUE,percentiles=c(.05,.50,.95),useMean=FALSE,
backdata=20)

# Bayesian VAR with steady-state prior
forecast(z4,periods=200,shocks=TRUE,plot=TRUE,percentiles=c(.05,.50,.95),useMean=FALSE,
backdata=20)