# --------------------------------------------------------------------------------------------------------------------------------------------
# Sistema IBGE de Recuperação Automática - SIDRA
# Pesquisa Industrial Mensal - Produção Física
# --------------------------------------------------------------------------------------------------------------------------------------------
# Tabela 3650 - Producao Fisica Industrial, por grupos e classes industriais selecionados
# Variavel - indice de base fixa sem ajuste sazonal (Base: media de 2012 = 100) (Numero-indice)
# Unidade Territorial - Brasil
# Mes x Grupos e classes industriais	
# Janeiro 2002 a Fevereiro 2017 (182 observacoes)
# Variavel 1: Fabricacao de celulose e outras pastas para a fabricacao de papel	
# Variavel 2: Fabricacao de biocombustiveis	
# Variavel 3: Fabricacao de cimento
# Variavel 4: Fabricacao de automóveis, camionetas e utilitarios
# URL: https://sidra.ibge.gov.br/tabela/3650
# --------------------------------------------------------------------------------------------------------------------------------------------
data = c(62,15.1,64.9,46.5,57.3,11.4,61.8,45.8,58.3,10.9,73.3,55.9,59.3,37.6,71.6,61.2,56.2,123.6,71.5,56.3,64.7,152.6,67.9,50.2,
67.3,159.4,71.1,50.4,65.7,161,76.1,47.9,60.4,148,74.9,54.1,65.1,147.1,80.3,60.7,70.1,76.2,76.2,55.4,70.9,35.4,67.5,44.6,71.7,
22.6,60.1,49.8,68.4,10.5,58.1,49.6,71.7,4.8,60.2,47.4,75,59.6,54.6,49.1,71.1,122.2,61.6,53,73.3,153.7,58,46.6,75.1,169,63.3,
48.7,73,171.1,64.5,48.8,74.5,176.1,62.6,56.7,78.2,171,67.6,60,72.7,90.3,62.2,57.8,78.4,36.7,58.6,53.7,80.8,18.2,57.4,57.5,
76.7,13.6,54.4,52.8,80.3,11.1,63.6,69.5,73.5,25.4,58.2,60.6,78.4,86,63.7,64.3,76.2,130.4,64.3,65.7,79,151,69,68.8,77.4,
181.4,71.1,71.4,74.1,195.4,69.8,72.6,80.1,177.1,68.5,69.6,78.8,162.2,66.8,72.4,80.4,87.8,66.9,65.9,83.7,31.3,65.2,61.4,74.5,
17.9,60,65.7,85.2,8.2,65.5,77.5,77.3,48.7,67.2,74.8,82,140,69.6,75.7,80.3,155.9,71.9,78.5,85.7,170.7,73.7,75.4,73,189.6,77.3,
78.5,79.6,171.9,74.8,73.3,85.9,174.9,74.8,71,82.8,122,72.5,79.1,86.4,46.8,73.4,76.1,88.7,22.3,71.2,73.8,79.4,8.5,66,72.2,88,6.7,
73.8,82.6,87.5,59,68.1,73.3,86.1,152.4,76.9,86.3,84.2,177.8,73.4,79.3,88,191.5,76.3,78.7,80,200.5,83.7,87.4,82.1,191.4,79.2,71.5,
82.1,173.9,80.7,81.7,83.2,130.7,75.9,78.2,87.7,47.3,75.5,66.8,90.5,21.9,69.8,76.4,81.4,12.7,68.4,70.2,88,9.7,78,89.3,81.8,59.4,
76.5,79,85.1,161.8,80.3,95.4,89.3,205.5,76.8,88.8,92.4,192.7,80.6,95.1,85.3,229.3,83,99.7,78,229.2,80.3,88.9,94.2,233.3,88,104.3,
94.1,176.4,80.9,96.6,98.9,80.8,78.2,78.1,97.3,33.3,78.4,94,93.4,25.6,74.6,88.9,97.1,17.5,84.2,101,89.5,63.5,79.6,106.8,102,168.6,
84.4,104.1,94.6,204.6,88,109,98.8,237.8,94.2,112.7,94.6,229,98.8,110.6,95.1,243.9,95.3,105.7,98.2,232.4,99,103.3,84.9,212.3,
90.1,66.8,102.3,122.1,83.5,33.8,97.6,43.7,79.4,67.4,89,32,73.5,72.1,88.8,35.1,84.4,99.2,89.9,92,78.2,91.7,99.8,173.4,86.7,97.6,
95.5,175.3,85.9,101.3,97.4,176.2,91.2,103.2,95,185.6,93.8,106.3,98.5,159.5,88.4,99.3,100.8,178.7,93.3,114.3,96.1,161.9,93.3,
105.2,102.1,100.3,88,88.4,99.2,42.3,84.3,88.7,89.7,35.7,82.9,88.2,103.5,46.4,95.4,119.2,89.5,118.4,88.8,106.2,97.8,179.2,97.6,
111.2,95.9,199.6,93.5,106.8,104.4,202.2,98.6,110,92.5,221.8,103,119.3,94.2,193,102,109.2,101.9,168.6,103.3,112.5,101.7,122.7,
98.2,113.1,103.2,55.9,97,98.6,101.4,24.2,87,94.9,91.9,16,88.7,111.9,105.4,18.9,92.7,114.2,92.2,52.3,93.6,98.7,94.5,166.7,100.4,
106.3,99.2,168.6,98.3,102.9,99.1,173.4,102.2,108.8,101.5,175,104.7,114.9,93.1,181.7,104.1,89.4,103.3,133.4,106.9,95.6,99.3,
81.8,99.1,97.3,103.2,32.6,95.5,88.5,96.7,30.2,88.3,76.1,100.3,27.8,92.8,75,99.9,23.3,105,107.4,96.6,40.5,96.5,93.4,97,111.6,
103.3,101.5,98.5,110.3,95.5,99.9,104.7,162.4,102.2,109.9,97.9,175.1,109.8,119.1,95.5,162.7,104.2,102.9,102.7,163.3,106.2,
115.8,103,125.1,98.9,109.3,107.3,67.8,97.2,89.5,100.1,27.9,90.2,102.3,88.5,19.5,89.5,81.2,94.4,21.8,102.2,110.9,97.8,85.5,
95.2,116.3,98.9,159.6,100.2,116.1,99.5,139.4,98.2,106.4,101.1,175.6,105.1,106.4,96.6,189.2,111.1,113.8,93.1,168.9,107.2,
107.1,103.1,159.1,110.8,108.4,97.5,138.5,103.5,99.7,98.1,80.9,94.6,78.5,97.9,28.1,92.2,83.9,88.4,20.7,94.2,98.8,96.8,30.6,
100.5,93,95.7,89.7,94.9,91.9,97.4,158.1,102.1,93.7,97.8,173.6,94.4,73,101.2,169.2,101,84.2,99.1,203.2,105.7,90.7,91.9,176.2,
104.9,102.9,96.6,181.8,105.7,101,100,105.4,99.9,90.4,100.7,52.7,87.8,71.8,101.8,30.5,89.2,70.9,91.3,27,81.7,72,94.7,35.6,95.6,
91.1,92.9,95,84.2,79.6,103.6,147.8,94.8,76.3,107.3,178.4,85,65.3,116.9,174.2,88.8,77.5,113.7,213.5,95.2,76.4,109.6,175.9,90.2,
61.1,104.3,184.8,94.2,68.6,101.1,121.7,83.6,58.8,119.8,91.4,78.5,52.3,115.6,30.6,71.7,52.5,109.5,25.1,75.2,49,104.5,48.4,76,70.3,
114.4,129,78.2,60.2,123.9,150.1,75.8,62,109.2,140.5,76,64.4,113.7,183.4,82.5,66.9,113.9,171.3,80.8,64.2,122.6,161.9,64.7,64.8,
122.8,143,74.1,67.1,116.8,100.9,75.2,71.9,128.3,44.6,74,60.3,130.9,21.8,72.1,58.2,100.3,18.1,76.4,65.8)

data     = matrix(data,byrow=TRUE,ncol=4)
date     = seq(2002,2017+1/12,by=1/12)
cbind(date,data)
names = c("Celulose",	"Biocombustiveis","Cimento","Automoveis")
n = nrow(data)

par(mfrow=c(2,2))
for (i in 1:4){
  plot(date,data[,i],type="l",main=names[i],ylab="Indice",xlab="")
  abline(v=date[156],lty=2)
}

for (i in 1:4)
acf(data[,i],lag=61,main=names[i])

for (i in 1:4)
pacf(data[,i],lag=61,main=names[i])


install.packages("astsa")
library("astsa")

y = data[,2]
y1 = ts(y[1:n1],frequency=12,start=c(2002,1))
y2 = y[(n1+1):n]
n1 = 156
stat = matrix(0,12,11)
l = 0
for (p1 in 1:6)
for (P1 in 0:1){
	l = l + 1
    fit = sarima(y1,p=p1,d=0,q=0,P=P1,D=0,Q=0,S=12,details=FALSE)
    fit1 =  sarima.for(y1,n.ahead=26,p=p1,d=0,q=0,P=P1,D=0,Q=0,S=12)
    rmse = sqrt(mean((y2-fit1$pred)^2))
    mae = mean(abs(y2-fit1$pred))
    stat[l,] = c(p1,0,0,P1,0,0,fit$AIC,fit$AICc,fit$BIC,rmse,mae)
}
ind = 1:12
AICorder   = ind[stat[,7]==min(stat[,7])]
AICcorder = ind[stat[,8]==min(stat[,8])]
BICorder   = ind[stat[,9]==min(stat[,9])]
RMSEorder  = ind[stat[,10]==min(stat[,10])]
MAEorder =  ind[stat[,11]==min(stat[,11])]
c(AICorder,AICcorder,BICorder,RMSEorder,MAEorder)


sarima(y1,p=1,d=0,q=0,P=1,D=0,Q=0,S=12,details=FALSE)
sarima(y1,p=2,d=0,q=0,P=1,D=0,Q=0,S=12,details=FALSE)

par(mfrow=c(2,1))
sarima.for(y1,n.ahead=26,p=1,d=0,q=0,P=1,D=0,Q=0,S=12)
abline(v=date[156],lty=2)
lines(date[(n1+1):n],y[(n1+1):n],col=4,type="b",pch=15,cex=0.5)
sarima.for(y1,n.ahead=26,p=2,d=0,q=0,P=1,D=0,Q=0,S=12)
abline(v=date[156],lty=2)
lines(date[(n1+1):n],y[(n1+1):n],col=4,type="b",pch=15,cex=0.5)