################(0) Importação de dados

rm(list=ls())
load("D:/UFABC_Ensino/Econometria_III/labs/dados/dados_lab_4__22_10_26.Rdata")
ls()


################(1) Teste ADF: algoritmo (slides)



################(2) Teste ADF: comando e função de distribuição de probabilidades

##(2.a) Aprendendo ao fazer
##Série: Usar a série artificial st_adf_1
##Converter em objeto do tipo "ts"
adf_aprende<-ts(st_adf_1,start=1,end=100,name="st")
##Fórmula simples do modelo-base: diff(yt,1) = mu + delta.t + lag(yt,-1) + et
##Passo 1: gerando as transformações da série
##(1.a) Defasagem de ordem 1 da série
###Sintaxe do comando lag
?lag
##Aplicando o lag: lag(st,-j), j = ordem da defasagem
##Atenção: o sinal de menos indica que se trata de postergação
##e não antecipação da série
##Visualizando a operação de defasagem:
View(cbind(lag=lag(adf_aprende,-1),ts=adf_aprende))
?ts()
##Executando e armazenando em objeto do tipo ts()
adf_aprende_o<-ts(cbind(adf_aprende,[completar com o lag],start=1,end=100)
#Nota: é sempre necessário especificar start e end para garantir que vão
#coincidir com tempos de início e final da série original

#Importante: notar que a primeira observação não está disponível
#e um ponto de dado foi perdido.

#Atribuindo nomes intuitivos às séries
adf_aprende_o<-ts(adf_aprende_o,start=1,end=100,names=c("st","L1_st"))

##(1.b) Diferença de ordem 1 da série
##(1.b.1) implementação manual (diff(yt,1) = yt - lag(yt,-1)
diff_manual<-[completar]
adf_aprende_o<-ts(cbind(adf_aprende_o,[completar]),start=1,end=100,names=c("st","L1_st","D1_st"))


##(1.b.2) implementação automática

##Sintaxe do operador diferença
?diff

#Calculando a primeira diferença
diff(adf_aprende_o[,1])

#Armazenando
diff_auto<-[completar]
adf_aprende_o<-ts(cbind(adf_aprende_o,[completar],start=1,end=100,names=c("st","L1_st","D1_st","D1a_st"))
#Comparando manual e automático: como saber se duas séries 
#são iguais em todos os pontos?
View(adf_aprende_o)

#Vendo números de colunas das variáveis 
colnames(adf_aprende_o)
#Cálculo manual [completar]
[completar]
#OBS: notar o uso do na.rm=T

##(1.c) Demais diferenças do modelo-base aumentado
##diff(yt,1) = mu + delta.t + lag(yt,-1) + lag(diff(yt,1),-1) + ... + lag(diff(yt,1),-p) + et
##Defasagem da primeira diferença
L1D1_st<-[completar]
##Armazenando
adf_aprende_o<-ts(cbind(adf_aprende_o,[completar]),start=1,end=100,names=c(colnames(adf_aprende_o),"L1D1_st"))

##(1.d) Verificando que: 
#diferença de ordem 1 da defasagem de ordem 1 =
#defasagem de ordem 1 da diferença de ordem 1

D1L1_st<-[completar]
#Armazenando
adf_aprende_o<-ts(cbind(adf_aprende_o,[completar]),start=1,end=100,names=c(colnames(adf_aprende_o),"D1L1_st"))
#Comparando
sum(abs(adf_aprende_o[,6]-adf_aprende_o[,5]),na.rm=T)


##(2.b) Componentes da função: compreendendo
#(1) Transformações da série original (defasagem e diferenças)
#A série original tem de ser transformada de três maneiras:
#(1.a) na série defasada em um instante de tempo (variável para teste de RU)
#(1.b) na primeira diferença da série
#(1.c) nas demais diferenças da série necessárias (diferenças de ordem > 1)
#As transformações (1.a) e (1.b) são feitas no seguinte componente do comando:
#a<-cbind(1,year=ts(year),ck=ts(y),ck=diff(ts(y)),ck_l1=lag(ts(y),-1))

#A transformação (1.c) é feita no seguinte componente do comando:
#for(i in 1:K) {
#eval(parse(text=paste("a<-cbind(a,diff_l",i,"=diff(lag(ts(y),-i)))",sep="")))
#}
#Este loop gera automaticamente as diferenças de todas as defasagens da
#primeira à K-ésima defasagem


#(2) Redação da fórmula do modelo
#(2.a) Uso da função de texto "paste", primeiramente para gerar os nomes das
#diferenças de todas as defasagens:
##P.ex., K = 5
K<-5
paste("diff_l",1:K,sep="")

#(2.b) Uso da função paste para completar a fórmula com as demais variáveis
xnam<-paste("diff_l",1:K,sep="")
paste("y_d ~ 1+ year+y_l1+", paste(xnam, collapse= "+"))
#O que cada componente representa? (lembrar do modelo-base)


#(3) Estimação por mínimo quadrados ordinários
lm(formula,data=a)

##(4) Função completa: acionar

####Selecionar: início
adf<-function(K,year,y,int,trend) {
a<-cbind(1,year=ts(year),ck=ts(y),ck=diff(ts(y)),ck_l1=lag(ts(y),-1))
for(i in 1:K) {
eval(parse(text=paste("a<-cbind(a,diff_l",i,"=diff(lag(ts(y),-i)))",sep="")))
}
colnames(a)<-c("inter","year","y","y_d","y_l1",paste("diff_l",1:K,sep=""))
a<-data.frame(a[(K+2):nrow(a),])
a<-subset(a,is.na(year)==FALSE)
a$uni<-1
xnam <- paste("diff_l",1:K,sep="")
if (int == 1 & trend == 1) {
fmla <- as.formula(paste("y_d ~ 1+ year+y_l1+", paste(xnam, collapse= "+")))
}
if (int == 1 & trend == 0) {
fmla <- as.formula(paste("y_d ~ 1+ y_l1+", paste(xnam, collapse= "+")))
}
if (int == 0 & trend == 1) {
fmla <- as.formula(paste("y_d ~ 0+ year+y_l1+", paste(xnam, collapse= "+")))
}
if (int == 0 & trend == 0) {
fmla <- as.formula(paste("y_d ~ 0+ y_l1+", paste(xnam, collapse= "+")))
}
mqo<<-lm(fmla,data=a)
return(mqo)
}
####Selecionar: fim


#Por que precisamos da função? Quais são os dois desafios de rodar
#a regressão referente ao modelo base do teste ADF?

#Sintaxe da função:
function(K,year,y,int,trend)
#K = p-max = número máximo de defasagens autoregressivas
#year = vetor com os instantes de tempo
#y = vetor com a série de tempo
#trend: 1 se há tendência determinista, 0 caso contrário
#int: 1 se há intercepto, 0 caso contrário



################(3) Implementação 1: taxa de desemprego brasileira (FMI)


###Passo 1: modelo-base com os dois termos deterministas

###Passo 1.a: definição do número máximo de lags
#Vamos usar o critério da nota de aula 4:
#Qual é o número de observações?
dim(desemp)
#Há 35 < 100 observações, usamos o critério de lag máximo estimável
#Critério: pmax = (T-D-2-epilson)
#epilson = 10^-5
#Assumindo de início os dois termos deterministas
#Há necessidade de arrendondar "para baixo", usar a função floor para isso
[completar com o cálculo do pmax]
#p-max = [completar]

###Passo 1.b: estimar a regressão com p = p-max 
adf([completar])

##Verificando os resultados detalhados de estimação
summary(mqo)
mqo

###Passo 1.c: verificar se os resíduos são RBs
par(mfrow=c(2,1))
acf([completar])
pacf([completar])


#Importante: número de observações relevante
#Cada defasagem introduzida representa a perda de uma partícula
#informacional, um instante de tempo. O modelo-base é estimado
#com base nas observações remanescentes e é o número delas que 
#importante e não o tamanho da série empírica.
#Ou seja, é o tamanho efetivamente usado da série.
#Este é igual a tamanho da série - (defasagens auto-regressivas + 1)
#Notar que, como as defasagens autoregressivas são aplicadas à
#primeira diferença da série, uma observação adicional foi perdida
#no cálculo de tal diferença (pois ela contém uma defasagem da série)
T_obs<-(nrow(desemp)-14-1)

#Outra forma de obter o número é recuperando ele do objeto
#com os resultados da regressão linear
#Basta notar que número de graus de liberdade dos
#resíduos = N - #parâmetros estimados
#E também que tamanho do vetor "coefficients" = #parâmetros estimados
#Com isso, N = GL resíduos + #coefficients
mqo$df.residual+length(mqo$coeff)
T_obs

#Teste de Ljung-Box: limite de lags = tamanho da amostra remanescente - 1
#Além disso, #lag > parametros estimados
#Recordar: quantos graus de liberdade perdemos para gerar os resíduos?
Box.test([completar],lag=T_obs-1,fitdf=[completar],type="Ljung-Box")


###Resultado?
###(chegamos a um impasse?)

###Passo 1.c: aplicar o teste para a tendência

#1.c.1: extrair a estatística t para a tendência do objeto
#com os resultados da regressão
#Devemos saber: posição do coeficiente procurado
mqo$coeff
#Também devemos saber: como extrair o erro padrão
#do objeto MQO: usar a função vcov(mqo)

#Estat t:
[completar]
#Conferindo com saída do R
summary(mqo)

#1.c.2: consultar a tabela de Dickey e Fuller, 1981, para 
#obter o valor crítico. Nível de significância: 5%; bicaudal

#Importante: sempre levar em conta o número de observações remanescentes


#1.c.3: resultado do teste: rejeitar ou não a nulidade da tendência?
#Selecionar estatística para tamanho amostral mais próximo do
#disponível

###Passo 2: modelo-base sem tendência

##2.a: número máximo de lags
#Agora há apenas um termo determinista
floor((nrow(desemp)-1-2-10^-5)/2)-1
#p-max = 14

###Passo 2.b: estimar a regressão com p = p-max 
adf([completar])

#Ver resultados detalhados da estimação
summary(mqo)

###Passo 2.c: verificar se os resíduos são RBs
par(mfrow=c(2,1))
acf([completar])
pacf([completar])

#Ljung-Box: lags <= tamanho da amostra remanescente - 1
T_obs<-(nrow(desemp)-14-1)
T_obs

#Teste de Ljung-Box: limite de lags = tamanho da amostra remanescente - 1
#Além disso, #lag > parametros estimados
Box.test(mqo$residuals,lag=T_obs-1,fitdf=length(mqo$coeff),type="Ljung-Box")


###Passo 2.d: aumentar a defasagem em uma unidade
#testar novamente os resíduos
adf([completar])

#Teste de Ljung-Box: limite de lags = tamanho da amostra remanescente - 1
#Além disso, #lag > parametros estimados
T_obs<-([completar])
T_obs

Box.test(mqo$residuals,lag=T_obs-1,fitdf=length(mqo$coeff),type="Ljung-Box")

#(chegamos a um impasse?)

###Passo 2.d: testar os termos deterministas
#Obtendo a estatística t para o intercept
mqo$coeff
#Valor da estatística
[completar]
#Comparando com o VC da tabela de DF(1981)
#Número de observações remanescentes
T_obs<-(nrow(desemp)-15-1)
mqo$df.residual+length(mqo$coeff)
T_obs

	
###Passo 3: modelo sem termos deterministas

##3.a: número máximo de lags
#Agora não há termos deterministas, D =0 
floor([completar]) - 1
#p-max = 15
##Houve mudança do número máximo de lags?

###Passo 3.b: estimar a regressão com p = p-max 
adf([completar])

#Ver resultados detalhados da estimação
summary(mqo)

###Passo 3.c: verificar se os resíduos são RBs
par(mfrow=c(2,1))
acf(mqo$residuals,lag.max=length(mqo$coeff)+1)
pacf(mqo$residuals,lag.max=length(mqo$coeff)+1)

#Teste de Ljung-Box: limite de lags = tamanho da amostra remanescente - 1
#Além disso, #lag > parametros estimados
T_obs<-(nrow(desemp)-15-1)
T_obs

Box.test(mqo$residuals,lag=T_obs-1,fitdf=length(mqo$coeff),type="Ljung-Box")

#Passo 4: teste de RU
##Estatística do teste: alfa^/SE^(alfa^)
#Identificando posição do coeficiente
mqo$coeff
#Posição?
#Cálculo da estatística
estat<-[completar]

#Obtenção do valor crítico: tabela de Mackinnon, 2010
#O artigo adequado é o publicado em 2010, que consta
#no Moodle
#Vamos usar as tabelas 2 e 3 do apêndice, a saber:
#Modelo-base com constante e tendência: tabela 3
#Modelo-base sem constante ou sem tendência: tabela 2
#A tabela não fornece diretamente o valor crítico,
#é necessário fazer o cálculo com a função cúbica
#na página 14:
VC = beta_inf + beta1/T + beta2/T^2 + beta/T^3
#Sempre considerar N = 1, este o número de séries
#de tempo (ou de dimensões)

#Importante: número de observações relevante
#Relembrar discussão anterior: os testes devem ser embasados
#no número de observações remanescente. I.e., N - (defasagens+1)
#Ou como calculado abaixo
T_obs<-[completar]

#O presente caso tem os ingredientes:
#termos deterministas: nenhum --> tabela 2, variant = tau-n-c
#nível de significância = 5% ---> linha 3
#tamanho da amostra
T_obs<-mqo$df.residual+length(mqo$coeff)
#Neste caso:
beta_inf<--1.94100
beta_1<--0.2686 
beta_2<--3.365 
beta_3<-31.223

#Cálculo
summary(mqo)
names(mqo)
T_obs<-nrow(desemp)
VC<-[completar]
VC

#Comparando: estat vs VC
estat
VC
names(mqo)




################(4) Implementação 2: deflator (IBGE)
ls()
par(mfrow=c(1,1))
plot(ts(dados_cn_2$deflator))



################(5) Implementação 3: desmatamento da Amazônia Legal (INPE)
#objeto: desmatamento

plot(desmatamento$year,desmatamento$st,type="l",yaxt="n",ylab="desmatmento",xlab="ano")
axis(1,seq(min(desmatamento$year),max(desmatamento$year),by=1))
abline(v=2004,col="red")
abline(v=2018,col="red")
abline(v=2022,col="red")



################(6) Implementação 4: reservas cambiais do Brasil (Banco Central)
#objeto: reservas