# Consumo de energia elétrica - industrial
# Fonte: Ministério de Minas e Energia, Balanço Energético Nacional (MME)
# Frequência: Anual de 1968 até 2018
# Unidade: MWh

data = read.table("http://hedibert.org/wp-content/uploads/2020/03/consumo-energia-sudeste-sul.txt",header=TRUE)
ano  = data[,1]
y    = as.matrix(data[,2:8])


# Normalizando pela populacao corrente dos estados
populacao = c(45919049,17264943,21168791,4018650,11377239,7164788,11433957)
for (i in 1:7)
  y[,i] = y[,i]/populacao[i]


# Taxa de crescimento do consumo de energia
ano1 = data[2:nrow(data),1]
y1   = apply(log(y),2,diff)


# Analise exploratoria simples
par(mfrow=c(1,2))
plot(ano,y[,1],ylim=range(y),xlab="Ano",ylab="Consumo de energia",type="l")
for (i in 2:7)
lines(ano,y[,i],col=i)
legend("topleft",legend=c("SP","RJ","MG","ES","RS","SC","PR"),col=1:7,lty=1)

plot(ano1,y1[,1],ylim=range(y1),xlab="Ano",ylab="Taxa de crescimento do consumo",type="l")
for (i in 2:7)
lines(ano1,y1[,i],col=i)


# Analise VAR inicial
# install.packages("MTS")
# library(MTS)

# Escolha aqui y ou y1
##yy   = y[,1:4]
##year = ano
yy   = y1[,1:4]
year = ano1

# Fittinng a VAR(1) for (SP,RJ,MG,ES)
mle.fit = VAR(yy,1,include.mean=T)

# Forecasting with the estimated VAR(1)
# Mean and variance of 1,2,3-step ahead forecasts
# 
# y(t+1) = mu + A*y(t) + e(t+1)
# y(t+2) = mu + A*y(t+1) + e(t+2)
#        = mu + A*[mu + A*y(t) + e(t+1)] + e(t+2)
#        = (I+A)*mu + A*A*y(t) + A*e(t+1) + e(t+2)
# y(t+3) = mu + A*y(t+2) + e(t+3)
#        = mu + A*[mu + A*y(t+1) + e(t+2)] + e(t+3)
#        = (I+A)*mu + A*A*y(t+1) + A*e(t+2) + e(t+3)
#        = (I+A)*mu + A*A*[mu + A*y(t) + e(t+1)] + A*e(t+2) + e(t+3)
#        = (I+A+A*A)*mu + A*A*A*y(t) + A*A*e(t+1) + A*e(t+2) + e(t+3)
#
# E{y(t+1)|D(t)} = (I+0)*mu + A*y(t)
# E{y(t+2)|D(t)} = (I+A)*mu + A*A*y(t)
# E{y(t+3)|D(t)} = (I+A+A*A)*mu + A*A*A*y(t)
#
# V{y(t+1)|D(t)} = S
# V{y(t+2)|D(t)} = A*S*t(A)+S
# V{y(t+3)|D(t)} = A*A*S*t(A*A)+A*S*t(A)+S

q = ncol(yy)
n = nrow(yy)
mle.fit = VAR(yy,1,include.mean=T)
B       = t(mle.fit$coef[2:(q+1),])
BB      = diag(1,q)
hmax    = 20
y.f     = matrix(0,hmax,q)
I1      = diag(0,q)
mu      = mle.fit$coef[1,]
for (h in 1:hmax){
  I1 = I1 + BB
  BB = BB%*%B	
  y.f[h,] = I1%*%mu + BB%*%yy[n,]
}

estado = c("SP","RJ","MG","ES")
par(mfrow=c(2,2))
for (serie in 1:q){
plot(year,yy[,serie],ylim=range(yy[,serie],y.f[,serie]),xlim=c(year[1],year[n]+hmax),
     xlab="Ano",ylab="Consumo de energia",type="l")
points((year[n]+1):(year[n]+hmax),y.f[,serie])
title(estado[serie])
}


# Representacao MA(infinito) para um VAR(2)
#
# Psi(s)  = 0  para s<0
# Psi(0)  = I
#
# Para s > =1
#   Psi(s)  = B1*Psi(s-1)   + B2*Psi(s-2)
#
# Psi(1)  = B1
# Psi(2)  = B1*B1 + B2
# Psi(3)  = B1*B1*B1 + B1*B2 + B2*B1
# Psi(4)  = B1*[B1*B1*B1 + B1*B2 + B2*B1] + B2*[B1*B1 + B2]
#         = B1*B1*B1*B1 + B1*B1*B2 + B1*B2*B1 + B2*B1*B1 + B2*B2
#         .
#         .
#         .
# Representacao MA(infinito) para um VAR(3)
mle.fit  = VAR(yy,3,include.mean=T)
B1       = t(mle.fit$coef[2:(q+1),])
B2       = t(mle.fit$coef[(q+2):(2*q+1),])
B3       = t(mle.fit$coef[(2*q+2):(3*q+1),])
Psi      = array(0,c(hmax,q,q))
Psi[1,,] = diag(1,q)
Psi[2,,] = B1
Psi[3,,] = B1%*%B1 + B2
for (h in 4:hmax){
  Psi[h,,] = B1%*%Psi[h-1,,] + B2%*%Psi[h-2,,] + B3%*%Psi[h-3,,]
}
par(mfrow=c(3,q))
for (i in 1:3)
  for (j in 1:q){
    ts.plot(Psi[,i,j],ylab="",main=paste("Psi(",i,",",j,")",sep=""))
    abline(h=0,lty=2)
  }