#Limpando área de trabalho
rm(list=ls())
#Importando dados
setwd("D:/UFABC_Ensino/Econometria_III/labs/dados")
load("dados_lab_3__08_10_2026.Rdata")
ls()

############################(2) Metodologia de Box Jenkins, série real############################(2) Metodologia de Box Jenkins, série real
############################(2) Metodologia de Box Jenkins, série real############################(2) Metodologia de Box Jenkins, série real
############################(2) Metodologia de Box Jenkins, série real############################(2) Metodologia de Box Jenkins, série real

ls()
#Resumindo os dados
summary(dados_cn_2)

#Criando a série de PIB per capita não corrigida pela inflação
dados_cn_2$pib_pc<-dados_cn_2$pib/dados_cn_2$pop

#Visualizando: PIB real per capita em nível e primeira diferença
par(mfrow=c(1,2))
plot(dados_cn_2$ano,dados_cn_2$pib_pc,type="l")
abline(v=2020,col="red")
plot(dados_cn_2$ano,c(NA,diff(dados_cn_2$pib_pc)),type="l")
abline(v=2020,col="red")

#Há uma quebra de tendência em 2020.

#Criando a série com a primeira diferença do PIB per capita
dados_cn_2$pib_pc_diff<-c(NA,diff(dados_cn_2$pib_pc))

#Objeto com o tamanho da série (T)
T_obs<-nrow(dados_cn_2)-1
#OBS1: é perdida uma observação por conta da subtração
#OBS2: o critério de Enders para a ordem máxima de defasagens
#T/4, gera, neste caso, nível muito baixo, de 23/4 = 5.75 ~ 6
#Por isso vamos usar a ordem máxima computável, ou seja, 
#T_obs - 1, notando-se que há perda de uma defasagem pois se
#a série tem tamanho T_obs, a defasagem máxima possível entre
#duas observações é T_obs - 1 (p.ex., tomando a observação t = 1
#como referência, a observação mais afastada temporalmente é
#T_obs, existindo, pois, entre as duas uma diferença temporal de
#T_obs - 1 (pense em 1 e 23, defasagem máxima = 22)
#Defindo a ordem máxima de defasagens
OM_def<-T_obs - 1
#Primeiro estágio: inspeção da série

############################(3) Previsão dentro da amostra############################(3) Previsão dentro da amostra
############################(3) Previsão dentro da amostra############################(3) Previsão dentro da amostra
############################(3) Previsão dentro da amostra############################(3) Previsão dentro da amostra


#O modelo selecionado no passo anterior, AR(1), pode ter sua performance
#avaliada a partir da previsão dentro da amostra.
#Vamos tomar: 
#(1) Sub-amostra de estimação: T1 = 2/3 . T_obs ~ 15
#(2) Sub-amostra de avaliação: T2 = 1/3 . T_obs ~ 8
#O que significa: vamos tomar as 15 primeiras observações
#para estimar o modelo, e com base nisso prever as 8 últimas
#observações. Ou seja, realizar previsões a um número de passos
#adiante iguais a um passo, dois passos, três passos, ..., oito passos

#O comando que gera previsões é o predict(obj,n.ahead=h,se.fit=T)
#Em que: obj = objeto armazenando os resultados da estimação da especificação ARMA
#n.ahead = número de passos à frente
#se.fit=T --> reportar o erro padrão da previsão

####(A) previsão dinâmica
#Estimação com sub-amostra (T=T1)
ar_1_est_T1<-arima(dados_cn_2$pib_pc_diff[2:16],order=c(1,0,0))
#Notar a restrição ao tamanho da série: "2:16", anos: 2001 a 2015
#Atenção: por conta da perda da primeira observação devido à primeira diferença
#A série começa em t = 2, e tem-se de considerar até t = 16 para ter-se 
#15 observações

#Gerando a previsão dentro da amostra de avaliação
#Gerando em massa as previsões para todos os passos adiante
vec_prev<-predict(ar_1_est_T1,n.ahead=8,se.fit=T)$pred
vec_ep<-predict(ar_1_est_T1,n.ahead=8,se.fit=T)$se
res_prev<-as.data.frame(cbind(vec_prev,vec_ep))
names(res_prev)<-c("prev","ep")
#Anos previstos: 2016 a 2023

####(b) previsão estática
vec_prev_estat_e<-numeric(0)
##(b.1) para T+1 (2016)
ar_1_est_T1<-arima(dados_cn_2$pib_pc_diff[2:16],order=c(1,0,0))
vec_prev_estat_f<-predict(ar_1_est_T1,n.ahead=1,se.fit=T)
#armazenando
vec_prev_estat_e<-rbind(vec_prev_estat_e,c(as.numeric(vec_prev_estat_f[1]),as.numeric(vec_prev_estat_f[2])))
vec_prev_estat_e

##(b.2) para T+2 (2017) 
ar_1_est_T2<-arima(dados_cn_2$pib_pc_diff[2:17],order=c(1,0,0))
vec_prev_estat_f<-predict(ar_1_est_T2,n.ahead=1,se.fit=T)
#armazenando
vec_prev_estat_e<-rbind(vec_prev_estat_e,c(as.numeric(vec_prev_estat_f[1]),as.numeric(vec_prev_estat_f[2])))
vec_prev_estat_e

##(b.3) para T+3 (2018)
ar_1_est_T3<-arima(dados_cn_2$pib_pc_diff[2:18],order=c(1,0,0))
vec_prev_estat_f<-predict(ar_1_est_T3,n.ahead=1,se.fit=T)
#armazenando
vec_prev_estat_e<-rbind(vec_prev_estat_e,c(as.numeric(vec_prev_estat_f[1]),as.numeric(vec_prev_estat_f[2])))
vec_prev_estat_e

for (jj in 1:5) {
ar_1_est_Tx<-arima(dados_cn_2$pib_pc_diff[2:18+jj],order=c(1,0,0))
vec_prev_estat_f<-predict(ar_1_est_Tx,n.ahead=1,se.fit=T)
#armazenando
vec_prev_estat_e<-rbind(vec_prev_estat_e,c(as.numeric(vec_prev_estat_f[1]),as.numeric(vec_prev_estat_f[2])))
}
vec_prev_estat_e

res_prev_estat<-as.data.frame(vec_prev_estat_e)
colnames(res_prev_estat)<-c("prev","se")
res_prev_estat

#Comparando previsto e observado 1: visualmente

##Incluindo limites do intervalo de confiança de previsão
#res_prev$ic_li<-res_prev$prev-1.96*res_prev$ep
#res_prev$ic_ls<-res_prev$prev+1.96*res_prev$ep


##Gráfico
par(mfrow=c(1,1))
##Observado
plot(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24],type="l",ylim=range(dados_cn_2$pib_pc_diff[17:24],res_prev_estat$prev))
points(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24],col="black")

##Previsão dinâmica
lines(dados_cn_2$ano[17:24],res_prev$prev,col="blue")
points(dados_cn_2$ano[17:24],res_prev$prev,col="black")

##Previsão estática
lines(dados_cn_2$ano[17:24],res_prev_estat$prev,col="purple",lty=1,lwd=2)
points(dados_cn_2$ano[17:24],res_prev_estat$prev,col="purple",lty=1,lwd=2)

##Visualizando erros de previsão

##Previsão dinâmica
plot(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24]-res_prev$prev,type="l",ylim=range(dados_cn_2$pib_pc_diff[17:24]-res_prev$prev,dados_cn_2$pib_pc_diff[17:24]-res_prev_estat$prev))
points(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24]-res_prev$prev,col="black")

##Previsão estática
lines(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24]-res_prev_estat$prev,col="purple",lty=1,lwd=2)
points(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24]-res_prev_estat$prev,col="purple",lty=1,lwd=2)

#Aparentemente os erros estáticos são maiores com menor frequência

#MSE
mse_din<-sqrt(sum((dados_cn_2$pib_pc_diff[17:24]-res_prev$prev)^2)/8)
mse_estat<-sqrt(sum((dados_cn_2$pib_pc_diff[17:24]-res_prev_estat$prev)^2)/8)
#MAE
mae_din<-sum(abs(dados_cn_2$pib_pc_diff[17:24]-res_prev$prev))/8
mae_estat<-sum(abs(dados_cn_2$pib_pc_diff[17:24]-res_prev_estat$prev))/8
#MAP
map_din<-sum(abs(dados_cn_2$pib_pc_diff[17:24]-res_prev$prev)/dados_cn_2$pib_pc_diff[17:24])/8
map_estat<-sum(abs(dados_cn_2$pib_pc_diff[17:24]-res_prev_estat$prev)/dados_cn_2$pib_pc_diff[17:24])/8

#Qual das três abordagens teve melhor desempenho?

#(Uma razão para o pior desempenho da estática com base na MSE é
#a grande amplificação do erro cometido em 2022)

#Outra razão é uma mudança grande no ARMA reestimado
##ARMA inicial (amostra de estimação até 16o obs)
ar_1_est_T1
##ARMA final (amostra de estimação até 23o obs)
ar_1_est_Tx

#Enfim, o que me parece é que temos um problema de previsão bastante desafiador
#conforme indicado pela quebra brusca de tendência na série empírica

##Gráfico
par(mfrow=c(1,1))
##Observado
plot(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24],type="l",ylim=range(dados_cn_2$pib_pc_diff[17:24],res_prev_estat$prev))
points(dados_cn_2$ano[17:24],dados_cn_2$pib_pc_diff[17:24],col="black")