rm(list=ls())

riw = function(v,V){
	q  = nrow(V)
  iV = solve(V)
	U  = chol(iV)
	T  = matrix(0,q,q)
	for (i in 2:q)
    for (j in 1:(i-1))
      T[i,j] = rnorm(1)
  for (i in 1:q)
    T[i,i] = sqrt(rgamma(1,(v-i+1)/2,1/2))
	C = solve(t(T)%*%U)
	return(C%*%t(C))
}


names = c("Comercio Varejista  - índice",
"Consumo de energia - residencial (GWh)",
"Consumo de energia - industrial (GWh)",
"Carga energia elétrica (Mwmed)",
"Produção industrial - índice",
"UCI (%)",
"Desemprego (%)",
"VIX Adj Close",
"Brasil Embi Plus",
"Câmbio",
"PPI",
"ÍNDICE DE PREÇO DAS EXPORTAÇÕES",
"ÍNDICE DE PREÇO DAS IMPORTAÇÕES",
"ÍNDICE DE QUANTUM DAS EXPORTAÇÕES",
"ÍNDICE DE QUANTUM DAS IMPORTAÇÕES",
"Selic",
"Selic real ex-post 3 meses IPCA",
"Selic real ex-post 12 meses IPCA",
"Selic real ex-post 3 meses IGP-DI",
"Selic real ex-post 12 meses IGP-DI",
"Spread_bndes",
"Spread_geral",
"Spread_geral_PJ",
"Spread_geral_PF",
  "IPCA livres Índice",
"IPCA adm. Índice",
"IGP-DI índice",
"IPC-Br índice",
"IPC-Fipe índice",
"M1 (R$ mil)",
"M2 (R$ mil)",
"M3 (R$ mil)",
"M4 (R$ mil)",
"PMPP (R$ mil)",
"Base (R$ mil)",
"DV (R$ mil)",
"CRB Spot Index",
"Petróleo WTI Spot Price FOB (Dollars per Barrel)",
"IPA_IPC ln(ipa-di/ipc_di)",
"Óleo lubrificante (D% mensal no IPCA)",
"Gasolina (D% mensal no IPCA)",
"Energia elétrica residencial (D% mensal no IPCA)")

data  = read.table("economiabrasileira.txt",header=FALSE)
q     = ncol(data)-1
n     = nrow(data)
date  = data[2:n,1]
data  = data[2:n,2:(q+1)]-data[1:(n-1),2:(q+1)]
m     = apply(data,2,mean)
d     = diag(sqrt(apply(data,2,var)))
for (i in 1:q)
  data[,i] = (data[,i]-m[i])/d[i,i]
n     = nrow(data)
ind   = trunc(seq(1,n,length=6))
dat   = date[ind]

pdf(file="economiabrasileira.pdf",width=30,height=20)
par(mfrow=c(6,7))
for (i in 1:q){
  plot(data[,i],main=names[i],axes=FALSE,xlab="Month",ylab="",pch=16,type="b",ylim=range(data))
  axis(2);axis(1,at=ind,lab=dat);box()
}


# ------------------------------------------------------------------------
# AR(1) models
# ------------------------------------------------------------------------
M       = 10000
h       = 24
d       =  1
y       = data[2:(n-h),25]
X       = data[1:(n-h-1),25]
n1      = length(y)
mu      = rep(0,d)
Vb      = diag(1000000,d)
lambda  = 0.000001
nu      = 0.000001
par1    = (nu+n)/2
C1      = solve(t(X)%*%X+solve(Vb))
mu1     = C1%*%(t(X)%*%y+solve(Vb)%*%mu)
par2    = (nu*lambda+t(y)%*%y+t(mu)%*%solve(Vb)%*%mu-t(mu1)%*%solve(C1)%*%mu1)/2
sig2    = 1/rgamma(M,par1,par2)
beta    = rnorm(M,mu1,sqrt(sig2*C1))
yf2     = matrix(0,M,h+1)
yf2[,1] = rep(data[n-h+1,25],M)
for (i in 2:(h+1)){
  for (j in 1:M){
    yf2[j,i] = rnorm(1,beta[j]*yf2[j,i-1],sqrt(sig2[j]))
  }
}
yf2 = yf2[,2:(h+1)]
quants = t(apply(yf2,2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(1,1))
plot(data[,25],ylim=range(data[,25],quants),xlim=c(1,n),main="",xlab="Month",ylab="",pch=16,type="b",axes=FALSE)
axis(2);axis(1,at=ind,lab=dat);box()
lines((n-h+1):n,quants[,2],col=4,type="b",pch=17)
lines((n-h+1):n,quants[,1],type="b",pch=17,col=4)
lines((n-h+1):n,quants[,3],type="b",pch=17,col=4)
abline(v=n-h,lty=2)

# ------------------------------------------------------------------------
# FACTOR ANALYSIS followed by VAR(1) model: IPCA livres + K common factors
# ------------------------------------------------------------------------
nK = 4
X1 = data[1:(n-h),c(1:24,26:q)]
quants2 = array(0,c(nK,q,h,3))
for (K in 1:nK){
fac = factanal(X1,K,"mle",scores="regression")
pairs(cbind(data[1:(n-h),25],fac$scores))
cor(cbind(data[1:(n-h),25],fac$scores))
F = fac$scores
yy = cbind(data[1:(n-h),25],F)
nn = nrow(yy)

# Prior hyperparameters
q   = 1+K
b0  = matrix(0,q,q)
B0  = matrix(1000,q,q)
V0  = diag(0.001,q)
v0  = 0.001

# Other summaries
y1  = yy[2:nn,]
X   = yy[1:(nn-1),]
XtX = t(X)%*%X
Xty = t(X)%*%y1
v1  = v0 + (nn-1)
iB0 = 1/B0
iB0b0 = iB0*b0

# Initial value
b  = B0

# Gibbs sampler
set.seed(132154)
burnin  = 10000
M       = 10000
niter   = burnin + M
draws.b = array(0,c(niter,q,q))
draws.s = array(0,c(niter,q,q))
e       = matrix(0,n-1,q)
for (iter in 1:niter){
  print(iter)
  # Sampling the variance-covariance matrix
  e = y1-X%*%t(b)
  V1 = V0+t(e)%*%e
  s  = riw(v1,V1)
  # Sampling the autoregressive coefficients
  for (i in 1:q){
  	  V     = solve(XtX/s[i,i]+diag(iB0[i,]))
    m     = V%*%(Xty[,i]/s[i,i]+iB0b0[i,])
    b[i,] = m + t(chol(V))%*%rnorm(q)
  }
  # Storing the draws
  draws.b[iter,,] = b
  draws.s[iter,,] = s
}
draws.b = draws.b[(burnin+1):niter,,]
draws.s = draws.s[(burnin+1):niter,,]

# Posterior distributions (marginals)
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
  	  breaks = seq(min(draws.b[,i,j]),max(draws.b[,i,j]),length=20)
    hist(draws.b[,i,j],prob=TRUE,xlab="",main="",breaks=breaks,col=grey(0.9))
    title(paste("B",i,j,sep=""))
  }

par(mfrow=c(q,q))
for (i in 1:q){
  for (j in 1:q){
  	  if (i<=j){
  	    breaks = seq(min(draws.s[,i,j]),max(draws.s[,i,j]),length=20)
      hist(draws.s[,i,j],prob=TRUE,xlab="",main="",breaks=breaks,col=grey(0.9))
      title(paste("S",i,j,sep=""))
  	  	}else{
  	  	  corr = draws.s[,i,j]/sqrt(draws.s[,i,i]*draws.s[,j,j])
  	    breaks = seq(min(corr),max(corr),length=20)
      hist(corr,prob=TRUE,xlab="",main="",xlim=c(-1,1),breaks=breaks,col=grey(0.9))
      title(paste("R",i,j,sep=""))
      abline(v=0,lwd=2)
    }
  }
}

h        = 24
yf2      = array(0,c(M,h+1,q))
yf2[,1,] = matrix(yy[nn,],M,q,byrow=TRUE)
for (i in 2:(h+1)){
  for (j in 1:M){
    yf2[j,i,] = draws.b[j,,]%*%yf2[j,i-1,] + t(chol(draws.s[j,,]))%*%rnorm(q)
  }
}
yf2 = yf2[,2:(h+1),]
for (i in 1:q)
  quants2[K,i,,]= t(apply(yf2[,,i],2,quantile,c(0.025,0.5,0.975)))
}
  
par(mfrow=c(1,1))
plot(data[,25],ylim=range(data[,25],quants,quants2[,1,,]),
     xlim=c(1,n),main="",xlab="Month",ylab="",pch=16,type="b",axes=FALSE)
axis(2);axis(1,at=ind,lab=dat);box()
for (j in 1:3){
  lines((n-h+1):n,quants[,j],lwd=2)
  for (K in 1:nK)
    lines((n-h+1):n,quants2[K,1,,j],col=1+K,lwd=2)
}
abline(v=n-h,lty=2)
legend(n-h,max(data[,25]),legend=c("AR","FAVAR1","FAVAR2","FAVAR3","FAVAR4"),lwd=rep(2,1+K),
       col=1:(1+K),lty=rep(1,1+K),bty="n",cex=2)

dev.off()