# Taxa de desemprego na RMSP
# Frequência: Mensal de 1984.12 até 2018.03
# Fonte: Fundação Sistema Estadual de Análise de Dados, 
# Pesquisa de Emprego e Desemprego (Seade/PED)
# Unidade: (%) 
# Comentário: Quadro: Taxas de desemprego, por tipo - Região 
# Metropolitana de São Paulo, 
# município de São Paulo e demais municípios da Região 
# Metropolitana de São Paulo. In: http://www.seade.gov.br . 
# Obs.: Região Metropolitana de São Paulo (RMSP). 
# Compreende desemprego oculto (trabalho precário e desemprego por 
# desalento) e desemprego aberto.
# Atualizado em: 04/05/2018

data = read.table("unemployment-RMSP.txt",header=TRUE)
y = data[,3]
date = data[,1]+data[,2]/12-1/24
plot(date,y,xlab="Months",ylab="Unemployment rate",type="l")


install.packages("forecast") 
library(forecast) 

un = ts(y,frequency=12,start=c(1984,12))
plot(un)

sarima.fit = Arima(window(un,end=c(2010,12)), order=c(1,1,0),seasonal=list(order=c(1,1,1),period=12), lambda=0) 

#plot(un)
#lines(date[1:373],sarima.fit$fitted,col=2)

plot(forecast(sarima.fit,h=87))
lines(un)

pdf("unemployment-RMSP-graph1.pdf",width=20,height=15)
par(mfrow=c(4,4))
for (i in 0:1)
for (j in 0:1)
for (k in 0:1)
for (l in 0:1){
sarima.fit = Arima(window(un,end=c(2015,12)), order=c(i,1,j),seasonal=list(order=c(k,1,l),period=12), lambda=0) 
plot(forecast(sarima.fit,h=27))
lines(un)
}
dev.off()



pdf("unemployment-RMSP-graph2.pdf",width=20,height=10)
l = 1
bic = matrix(0,8,7)
par(mfrow=c(2,4))
for (i in 0:1)
for (j in 0:1)
for (k in 0:1){
sarima.fit = Arima(window(un,end=c(2015,12)), order=c(i,1,j),seasonal=list(order=c(k,1,1),period=12), lambda=0) 
plot(forecast(sarima.fit,h=27))
lines(un)
bic[l,] = c(i,1,j,k,1,1,sarima.fit$bic)
l = l + 1
} 
dev.off()




sarima.fit = Arima(window(un,end=c(2017,12)), order=c(1,1,0),seasonal=list(order=c(1,1,1),period=12), lambda=0) 

plot(forecast(sarima.fit,h=48))
lines(un)