#####(1) Geração computacional de dois processos estocásticos básicos#####(1) Geração computacional de dois processos estocásticos básicos
#####(1) Geração computacional de dois processos estocásticos básicos#####(1) Geração computacional de dois processos estocásticos básicos
#####(1) Geração computacional de dois processos estocásticos básicos#####(1) Geração computacional de dois processos estocásticos básicos


#Nesta primeira parte vamos responder à questão: como é possível gerar,
#computacionalmente, processos estocásticos?
#Antes de responder, cabe esclarecer que, rigorosamente, trata-se de
#simular tais processos, pois, a rigor, eles são gerados pela sociedade
#ou natureza, cabendo ao analista apenas observá-los depois que eles
#ocorrem.
#Por isso a pergunta tem de ser refeita: como é possível simular, a partir
#de equações a diferenças finitas, processos estocásticos?

#Uma forma de entender é procurando descobrir como um comando pré-disponível
#simula tais processos. Vamos usar o comando arima.sim, do pacote "stats"

#Documentação da função arima.sim, a qual simula um modelo
#de séries temporais
?arima.sim

#Documentação da função arima, a qual estima um modelo
#de séries temporais (será usada na segunda parte desta aula)
?arima

#(1.a) Ruído branco automático
rb_sim<-arima.sim(n=100, list(order=c(0,0,0)))

#(1.b) Ruído branco manual (supondo independência ao longo do tempo)
rb_manual<-rnorm(100,0,1)

#Comparando: automático vs manual
par(mfrow=c(1,2))
plot(ts(rb),main="Ruído branco automático")
plot(ts(rb_1),main="Ruído branco manual")

#Ruído branco automático ~ ruído branco manual

###Processo 2: auto-regressivo
#(3.a) Auto-regressivo de ordem 1 (AR(1)), automático
#com coeficiente 0.7
ar_1_sim<-arima.sim(n=100, list(order=c(1,0,0),ar=0.7))

###Processo 3: médias móveis
# Médias móveis de ordem 1 (MA(1)), automático
ma_1_sim<-arima.sim(n=100, list(order=c(0,0,1),ma=0.5))


###Processo 4: ARMA(1,1)
#Completar
arma_1_sim<-arima.sim(n=100, list(order=c(1,0,1),ar=0.7,ma=0.5))


###Visualizando todos os processos simulados
ts_set<-ts(data=cbind(rb_sim,rb_manual,ar_1_sim,ma_1_sim,arma_1_sim),start=1,end=100,frequency=1)
par(mfrow=c(2,2))
plot(ts_set)





############################(2) Metodologia de Box Jenkins, introdução, séries artificiais############################(2) Metodologia de Box Jenkins, introdução, séries artificiais
############################(2) Metodologia de Box Jenkins, introdução, séries artificiais############################(2) Metodologia de Box Jenkins, introdução, séries artificiais
############################(2) Metodologia de Box Jenkins, introdução, séries artificiais############################(2) Metodologia de Box Jenkins, introdução, séries artificiais

####Referências computacionais:
##Kleiber & Zeileis, cap., seção 6.2

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


##########Série de tempo 1: st_1##########Série de tempo 1: st_1##########Série de tempo 1: st_1##########Série de tempo 1: st_1
##########Série de tempo 1: st_1##########Série de tempo 1: st_1##########Série de tempo 1: st_1##########Série de tempo 1: st_1
##########Série de tempo 1: st_1##########Série de tempo 1: st_1##########Série de tempo 1: st_1##########Série de tempo 1: st_1
 
###(1) B.J pt 1.a/b:inspeção da série
#Gráfico
plot(ts(st_1))

#Correlograma automático
par(mfrow=c(2,1))
acf(st_1)
pacf(st_1)

#Nota: o pico na defasagem zero deve ser ignorado pois trata-se
#da variância (conforme visto em sala de aula)

#Vamos usar a opção lag.max dos comandos acf e pacf para alterar
#a defasagem máxima para a regra de bolso de Enders

#Correlograma customizado
par(mfrow=c(2,1))
acf(st_1,lag.max=100/4)
pacf(st_1,lag.max=100/4)

#Critério 1, tabela:
##linha 1 da tabela, FAC: 
#esta indica decaimento exponencial, sugerindo uma especificação AR(p)
#linha 2 da tabela, FACP:
#a FACP aponta para truncamento na ordem 1, mas também na ordem 9.
#Quanto à ultima, ela é indicada pela existência de significância na ordem 8.
#Há, pois, duas especificações candidatas: AR(1) e AR(9)

#Critério 2, ordem máxima de significância
#gráfico 1, FAC: a ordem máxima de significância é..., sugerindo ...
#gráfico 2, FACP: a ordem máxima de significância é..., sugerindo ...

#Combinação dos dois critérios:
#A lista inicial, a qual deve compreender todas as especificações sugeridas
#pelos dois critérios é: [completar]

###(1) B.J pt 3.a: estimação e diagnóstico
#Para estimar as especificações da lista inicial
#utiliza-se o comando "ARIMA(série,p,0,q)"
#A sintaxe é tal que o primeiro argumento é a série
#sob estudo.
#Abaixo são estimadas todas as especificações [[T: deixar apenas a primeira]]
ar_1_est<-arima(st_1,order=c(1,0,0))
ar_9_est<-[completar]
ar_20_est<-[completar]
ma_20_est<-[completar]
#arma_9_20_est<-arima(st_1,order=c(9,0,20))
#Não foi possível estimar, não houve convergência: reduzir a ordem maior
#arma_9_19_est<-arima(st_1,order=c(9,0,19))
#arma_9_18_est<-arima(st_1,order=c(9,0,18))
#arma_9_17_est<-arima(st_1,order=c(9,0,17))
#arma_9_16_est<-arima(st_1,order=c(9,0,16))
#arma_9_15_est<-arima(st_1,order=c(9,0,15))
#arma_9_14_est<-arima(st_1,order=c(9,0,14))
#arma_9_13_est<-arima(st_1,order=c(9,0,13))
#arma_9_12_est<-arima(st_1,order=c(9,0,12))
#arma_9_11_est<-arima(st_1,order=c(9,0,11))
arma_9_10_est<-arima(st_1,order=c(9,0,10))
#Apenas o ARMA(9,10) converge

#Nota: usa-se o sufixo "est" para denotar que se trata de
#modelo estimado;

#A estimação do modelo geral diversos resultados, os quais ficam
#armazenados em sub-objetos. Esses últimos são listas com o comando
#"names":
names(ar_1_est)
#O objeto relevante no momento é "residuals", pois agora é o momento
#de inspecionar FAC e FACP dos resíduos da especificação candidata.
#O objeto da inspeção é verificar se há comportamento ruído branco, i.e.:
#Todas as defasagens apresentam, individualmente, segundo o teste de 
#significância individual e conjuntamente, segundo o teste de Ljung-Box
#não rejeitam a hiótese nula.

#Inspeção 1: testes de sig individual contidos no correlograma

#Candidata 1: AR 1
#em caso de bug: par(mar=c(1, 1, 1, 1))
par(mfrow=c(2,1))
#[completar aqui com FAC e FACP dos resídus do AR1]
#[completar aqui com FAC e FACP dos resídus do AR1]

#Observar o gráfico e responder às perguntas:
#Há autocorrelações significantes para alguma defasagem?
#Se houver, quais são as defasagens?

#Na FAC, ...
#Na FACP, ...

#Inspeção 2: Teste de Lung-Box
#O comando que aplica o teste de Ljung-Box é "Box.test"
#Vamos ler a documentação para aprender sobre a sintaxe
?Box.test

#Com base na documentação, os argumentos a incluir no comando são:
#primeiro argumento = resíduos do modelo estimado
#segundo argumento: lag = ordem de defasagem
#terceiro argumento: type = Ljung-Box
#quarto argumento: fitdf = número de graus de liberdade perdidos =
#número de parâmetros estimados = p + q (= 6, no caso)
#(sempre há perda quando o teste se aplicar aos resíduos)

#É possível aplicar o teste com apenas uma linha de código, com o comando
#sapply. Este comando aplica uma função a um vetor de argumentos.
#Dado que quando aplicado ao resíduo, o teste não logra gerar p-valor
#para uma defasagem máxima j* < p+q, deve-se aplicar o teste apenas
#para j* >= p+q+1

#Especificando j* inicial
n_param_p1<-1+1
#Fazendo o teste
ar_1_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(ar_1_est$residuals, lag = x, type = "Ljung-Box", fitdf = 1)$p.value)
#Adicionado valores faltantes para as deasdagens < p + q + 1
ar_1_est_LB<-c(rep(NA,n_param_p1-1),ar_1_est_LB)

#Visualizando o objeto com todos os p-valores
View(ar_1_est_LB)

#Visualizando testes de significância individual e conjunta
#em um mesmo gráfico

par(mfrow=c(3,1))
acf(ar_1_est$residuals,lag.max=100/4)
pacf(ar_1_est$residuals,lag.max=100/4)
plot(1:(100/4),ar_1_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#O teste de Ljung-Box aponta para aucorrelações não-significativas
#em todas as defasagens
#A especificação AR(1), pois, não deve ser descartada.

#Agora vamos aplicar o painel com três gráficos, e, pois, todos
#os testes, para as demais candidatas. Porém, é preciso gerar
#os p-valores do teste de Ljung-Box, o que é feito abaixo 

#Candidata 2: AR 9
##Especificando j* inicial
n_param_p1<-9+1
##Fazendo o teste
ar_9_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(ar_9_est$residuals, lag = x, type = "Ljung-Box", fitdf = 1)$p.value)
##Adicionado valores faltantes para as deasdagens < p + q + 1
ar_9_est_LB<-c(rep(NA,n_param_p1-1),ar_9_est_LB)

par(mfrow=c(3,1))
acf(ar_9_est$residuals,lag.max=100/4)
pacf(ar_9_est$residuals,lag.max=100/4)
plot(1:(100/4),ar_9_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [preencher]
#FACP: [preencher]
#Ljung-Box: [preencher]
#Resultado (eliminar/não eliminar): [preencher]

#Candidata 3: AR 20

#Especificando j* inicial
n_param_p1<-[completar]
#Fazendo o teste
ar_20_est_LB<-[completar]
#Adicionado valores faltantes para as defasagens < p + q + 1
ar_20_est_LB<-c(rep(NA,n_param_p1-1),ar_20_est_LB)

par(mfrow=c(3,1))
acf([completar],lag.max=100/4)
pacf([completar],lag.max=100/4)
plot(1:(100/4),ar_20_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [preencher]
#FACP: [preencher]
#Ljung-Box: [preencher]
#Resultado (eliminar/não eliminar): [preencher]

#Candidata 4: MA 20
#Especificando j* inicial
n_param_p1<-20+1
#Fazendo o teste
ma_20_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(ma_20_est$residuals, lag = x, type = "Ljung-Box", fitdf = 1)$p.value)
#Adicionado valores faltantes para as deasdagens < p + q + 1
ma_20_est_LB<-c(rep(NA,n_param_p1-1),ma_20_est_LB)

par(mfrow=c(3,1))
acf(ma_20_est$residuals,lag.max=100/4)
pacf(ma_20_est$residuals,lag.max=100/4)
plot(1:(100/4),ma_20_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [preencher]
#FACP: [preencher]
#Ljung-Box: [preencher]
#Resultado (eliminar/não eliminar): [preencher]


#Candidata 4: ARMA 9,10
#Especificando j* inicial
n_param_p1<-9+10+1
#Fazendo o teste
arma_9_10_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(arma_9_10_est$residuals, lag = x, type = "Ljung-Box", fitdf = 19)$p.value)
#Adicionado valores faltantes para as defasagens < p + q + 1
arma_9_10_est_LB<-c(rep(NA,n_param_p1-1),arma_9_10_est_LB)

par(mfrow=c(3,1))
acf(arma_9_10_est$residuals,lag.max=100/4)
pacf(arma_9_10_est$residuals,lag.max=100/4)
plot(1:(100/4),arma_9_10_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [preencher]
#FACP: [preencher]
#Ljung-Box: [preencher]
#Resultado (eliminar/não eliminar): [preencher]

###Classificação: agora deve-se obter os critérios de liberdade
#para as candidatas não eliminadas. Vamos calcular manualmente
#cada critério, utilizando as fórmulas abaixo:
#1.	AIC = T.log(SQR) + 2.K
#2.	BIC = T.log(SQR) + log(T).K
#3.	HQ = T.log(SQR) + 2.K.log(log(T))
#Como K = p+q:
#1.	AIC = T.log(SQR) + 2.(p+q)
#2.	BIC = T.log(SQR) + log(T).(p+q)
#3.	HQ = T.log(SQR) + 2.(p+q).log(log(T))

#Para armazenar os resultados, vamos utilizar uma tabela
#de seis colunas: 
#(i) especificação, (ii) p+q, (iii) SQR, (iv) AIC, (v) BIC, (vi) HQ
#Coluna 1:
tabela_CI_c1<-c("AR(1)","AR(9)","AR(20)","MA(20)","ARMA(9,10)")
#Coluna 2:
tabela_CI_c2<-c(1,9,20,20,19)
#Coluna 3: necessário calcular as SQRs de todas as especificações
SQR_ar1<-sum(ar_1_est$residuals^2)
SQR_ar9<-[completar]
SQR_ar20<-[completar]
SQR_ma20<-[completar]
SQR_arma910<-[completar]
tabela_CI_c3<-c(SQR_ar1,SQR_ar9,SQR_ar20,SQR_ma20,SQR_arma910)
#Coluna 4: necessário calcular os AICs (AIC = T.log(SQR) + 2.(p+q))
AIC_ar1<-100*SQR_ar1 + 2*1
AIC_ar9<-100*SQR_ar20 + 2*9
AIC_ar20<-100*[completar] + [completar]
AIC_ma20<-100*[completar] + [completar]
AIC_arma910<-100*[completar] + [completar]
tabela_CI_c4<-c(AIC_ar1,AIC_ar9,AIC_ar20,AIC_ma20,AIC_arma910)
#Coluna 5: necessário calcular os BICs (BIC = T.log(SQR) + log(T).(p+q))
BIC_ar1<-100*log(sum(ar_1_est$residuals^2)) +log(100)*1
BIC_ar9<-100*log(sum(ar_9_est$residuals^2)) + log(100)*[completar]
BIC_ar20<-100*log(sum(ar_20_est$residuals^2)) + log(100)*[completar]
BIC_ma20<-100*log(sum(ma_20_est$residuals^2)) + log(100)*[completar]
BIC_arma910<-100*log(sum(arma_9_10_est$residuals^2)) + log(100)*[completar]
tabela_CI_c5<-c(BIC_ar1,BIC_ar9,BIC_ar20,BIC_ma20,BIC_arma910)
#Coluna 6: necessário calcular os HQs (HQ = T.log(SQR) + 2.(p+q).log(log(T)))
HQ_ar1<-100*log(sum(ar_1_est$residuals^2)) +2*1*log(log(100))
HQ_ar9<-100*log(sum(ar_9_est$residuals^2)) + 2*9*log(log(100))
HQ_ar20<-100*log(sum(ar_20_est$residuals^2)) + 2*20*log(log(100))
HQ_ma20<-100*log(sum(ma_20_est$residuals^2)) + 2*20*log(log(100))
HQ_arma910<-100*log(sum(arma_9_10_est$residuals^2)) + 2*19*log(log(100))
tabela_CI_c6<-c(HQ_ar1,HQ_ar9,HQ_ar20,HQ_ma20,HQ_arma910)

#Formando a tabela completa:
tabela_CI<-cbind(tabela_CI_c1,tabela_CI_c2,tabela_CI_c3,tabela_CI_c4,tabela_CI_c5,tabela_CI_c6)
colnames(tabela_CI)<-c("candidata","p+q","SQR","AIC","BIC","HQ")

#Visualizando a tabela
View(tabela_CI)

#Vamos olhar com calma e detectar qual candidata atinge o 
#menor nível de pelo menos dois critérios de informação

#Será que é possível programar o R para que identificar
#a candidata a selecionar?

#Desafio: é um exercício de minimização numérica

#Transformando a tabela em data frame
tabela_CI<-as.data.frame(tabela_CI)
#Adicionando um asterisco na célula com menor AIC 
[completar]
#Adicionando um asterisco na célula com menor BIC 
[completar]
#Adicionando um asterisco na célula com menor HQ 
[completar]

#Visualizando: qual é a candidata com mais asteriscos?
View(tabela_CI)

#Resultado do diagnóstico: 
#Qual é?
#Está correto? Este é o resultado esperado? (pensar)
#Será que o critério de parcimônia atingiu seu objetivo?


#Vamos visualizar novamente os resultados de estimação
ar_1_est

#Como a candidata selecionada, AR(1), tem ordem 1, não há como
#continuar o algoritmo, pois o próximo passo requereria considerar
#um AR(0), o que não faz sentido.
#Conclui-se, pois, que a especificação AR(1) melhor reproduz a 
#série artificial.


##########Série de tempo 2: st_2##########Série de tempo 2: st_2##########Série de tempo 2: st_2##########Série de tempo 2: st_2
##########Série de tempo 2: st_2##########Série de tempo 2: st_2##########Série de tempo 2: st_2##########Série de tempo 2: st_2
##########Série de tempo 2: st_2##########Série de tempo 2: st_2##########Série de tempo 2: st_2##########Série de tempo 2: st_2##Passo 1: Definição da lista inicial de candidatas a partir da

#inspeção da série original
#Neste segundo exercício vamos incluir o teste de Ljung-Box nesta
#inspeção da série original

#[escrever o código para gerar o painel gráfico]
#[ele deve conter: FAC, FACP e Ljung-Box]

#O que o gráfico revela?
#Abordagem da tabela, linha 1, FAC: [completar]
#A FAC sugere, pois...
#Abordagem da tabela, linha 2, FACP: [completar]
#A FACP sugere ...

#Abordagem da ordem máxima significativa:
#FAC: ordem máxima significativa = [completar] ---> [completar]
#FACP: [completar] --> [completar]

#A lista inicial de candidatas é: [completar]

##Passo 2: Estimação e diagóstico da lista inicial

#Passo 2.a) Estimação e diagnóstico eliminatório

##Candidata 1: AR(3)
###Estimação
ar_3_est<-[completar]

##Resultados
ar_3_est

###Diagnóstico
#Especificando j* inicial
n_param_p1<-3+1
#Fazendo o teste
ar_3_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(ar_3_est$residuals, lag = x, type = "Ljung-Box", fitdf = 3)$p.value)
#Adicionado valores faltantes para as deasdagens < p + q + 1
ar_3_est_LB<-c(rep(NA,n_param_p1-1),ar_3_est_LB)

#Gráfico
par(mfrow=c(3,1))
acf([completar],lag.max=100/4)
pacf([completar],lag.max=100/4)
plot(1:(100/4),ar_3_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [completar]
#FACP: [completar]
#Ljung-Box: [completar]
#Resultado: [completar]

##Candidata 2: MA(22)

###Estimação
ma_22_est<-arima([completar])

ma_22_est

###Diagnóstico
#Especificando j* inicial
n_param_p1<-22+1
#Fazendo o teste
ma_22_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(ma_22_est$residuals, lag = x, type = "Ljung-Box", fitdf = 22)$p.value)
#Adicionado valores faltantes para as deasdagens < p + q + 1
ma_22_est_LB<-c(rep(NA,n_param_p1-1),ma_22_est_LB)


#Gráfico
par(mfrow=c(3,1))
acf([completar],lag.max=100/4)
pacf([completar],lag.max=100/4)
plot(1:(100/4),ma_22_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [completar]
#FACP: [completar]
#Ljung-Box: [completar]
#Resultado: [completar]

##Candidata 3: ARMA(3,22)

###Estimação
arma_3_22_est<-arima(st_2,order=c(3,0,22))

arma_3_22_est

###Diagnóstico
#Especificando j* inicial
n_param_p1<-3+22+1
#Necessário aumentar gráfico para j* = 30
#Fazendo o teste
arma_3_22_est_LB<-sapply(n_param_p1:30,function(x)Box.test(arma_3_22_est$residuals, lag = x, type = "Ljung-Box", fitdf = 3+22)$p.value)
#Adicionado valores faltantes para as deasdagens < p + q + 1
arma_3_22_est_LB<-c(rep(NA,n_param_p1-1),arma_3_22_est_LB)

length(arma_3_22_est_LB)


#Gráfico
par(mfrow=c(3,1))
acf(arma_3_22_est$residuals,lag.max=30)
pacf(arma_3_22_est$residuals,lag.max=30)
plot(1:30,arma_3_22_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: [completar]
#FACP: [completar]
#Ljung-Box: [completar]
#Resultado: [completar]


#Passo 2.b) Classificação e seleção

#Para armazenar os resultados, vamos utilizar uma tabela
#de seis colunas: 
#(i) especificação, (ii) p+q, (iii) SQR, (iv) AIC, (v) BIC, (vi) HQ

candidatas<-c("ar_3_est","ar_22_est")
ordens<-c(3,22)

#O loop simplificado abaixo gera a tabela de classificação
##Início
T<-100
tabela_CI<-numeric(0)
for (j in (1:length(candidatas))) {
fluxo<-get(candidatas[j])
SQR_fluxo<-sum(fluxo$residuals^2)
AIC_fluxo<-T*log(SQR_fluxo) + 2*ordens[j]
BIC_fluxo<-T*log(SQR_fluxo) + log(T)*ordens[j]
HQ_fluxo<-T*log(SQR_fluxo) + 2*ordens[j]*log(log(T))
linha_fluxo<-c(candidatas[j],ordens[j],AIC_fluxo,BIC_fluxo,HQ_fluxo)
tabela_CI<-rbind(tabela_CI,linha_fluxo)
}
colnames(tabela_CI)<-c("candidata","ordem","AIC","BIC","HQ")
##Fim
View(tabela_CI)

###Qual é a melhor candidata?

##Passo 2: Estimação e diagóstico eliminatório da lista inicial
#Seguindo o algoritmo 1, não houve o problema de conjunto vazio.
#Deve-se, pois, reduzir a ordem da especificação selecionada.
#Esta será reduzida para [completar].

###Estimação
ar_2_est<-arima(st_2,order=c(2,0,0))


###Diagnóstico
#Especificando j* inicial
n_param_p1<-2+1
#Necessário aumentar gráfico para j* = 30
#Fazendo o teste
ar_2_est_LB<-sapply(n_param_p1:(100/4),function(x)Box.test(ar_2_est$residuals, lag = x, type = "Ljung-Box", fitdf = 2)$p.value)
#Adicionado valores faltantes para as deasdagens < p + q + 1
ar_2_est_LB<-c(rep(NA,n_param_p1-1),ar_2_est_LB)


#Gráfico
par(mfrow=c(3,1))
acf(ar_2_est$residuals,lag.max=(100/4))
pacf(ar_2_est$residuals,lag.max=(100/4))
plot(1:(100/4),ar_2_est_LB,ylim=c(0,1))
abline(h=0.05,col="red",lwd=2)

#FAC: autocorrelações significativas nas defasagens 2, 5 e 21
#FACP: autocorrelações significativas nas defasagens 1 e 5
#Ljung-Box: todas defasagens são significativas
#cabe examinar os números do p-valor
View(ar_2_est_LB)
#De fato, há significância em todas as defasagens

#A candidata deve ser eliminada

#O algoritmo é interrompido, e conclui-se que a candidata selecionado
#no passo anterior mais bem descreve a série original. 

#Vale a pena observar novamente a estimação da melhor candidata
ar_3_est