Trabalho Econometria II

O objetivo do trabalho é aplicar conceitos sobre Modelos Autorregressivos Integrados de Médias Móveis (ARIMA) ou Modelos Sazonais Autorregressivos Integrados de Médias Móveis (SARIMA).

Para tanto, usamos a série temporal ITUB3 e com ela seguimos o seguinte processo para o modelo ARIMA:

1. Visualizar os dados e identificar observações fora do padrão (outliers) e eliminá-los.
2. Se necessário, transformar os dados para estabilizar a variância (logaritmo dos dados, por exemplo)
3. Testar se os dados são estacionários. Caso tenha raiz unitária é preciso diferenciar os dados até se tornarem estacionários
4. Examinar as funções de autocorrelação parcial (FAC) e autocorrelação parcial (FACP) para determinar as ordens máximas \(P\) e \(Q\) para os componentes AR e MA da série estacionária (diferenciada)
5. Estimar todas as combinações para \(P\) e \(Q\).
6. Escolher o modelo com menor AIC
7. Examinar se os resíduos se comportam como ruído branco. Caso contrário, retornar ao passo 3 ou 4.
8. Uma vez que os resíduos são ruído branco, obter as previsões.

Antes de tudo, é necessário instalar os pacotes que contêm as funções necessárias para o desenvolvimento do estudo.

O primeiro passo foi puxar a série da base de dados do yahoo finance usando a função getSymbols do pacote quantmod.

itau=quantmod::getSymbols("ITUB3.SA", src = "yahoo", auto.assign = FALSE, from = '2016-01-01', to = '2018-05-01')

<U+393C><U+3E31>getSymbols<U+393C><U+3E32> currently uses auto.assign=TRUE by default, but will
use auto.assign=FALSE in 0.5-0. You will still be able to use
<U+393C><U+3E31>loadSymbols<U+393C><U+3E32> to automatically load data. getOption("getSymbols.env")
and getOption("getSymbols.auto.assign") will still be checked for
alternate defaults.

This message is shown once per session and may be disabled by setting 
options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.


WARNING: There have been significant changes to Yahoo Finance data.
Please see the Warning section of <U+393C><U+3E31>?getSymbols.yahoo<U+393C><U+3E32> for details.

This message is shown once per session and may be disabled by setting
options("getSymbols.yahoo.warning"=FALSE).

A série temporal, contudo, não apresentava estacionariedade e, portanto, foi necessário transformar os dados para seu formato logaritmico. Isso foi feito por meio da função Return.calculate do pacote PerformanceAnalytics, que aplicou o logaritmo nos valores de fechamento da série. Logo na sequência, foi aplicado o teste ADF para verificar estacionariedade.

daily_return= PerformanceAnalytics::Return.calculate(itau$ITUB3.SA.Close, method = "log")
dygraphs::dygraph(daily_return, main = "Retorno Diário da ITUB3.SA") %>% dyRangeSelector()

unitRoot = fUnitRoots::adfTest(daily_return,lags=2, type=c("nc"))

the condition has length > 1 and only the first element will be usedp-value smaller than printed p-value

print(unitRoot)

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 2
  STATISTIC:
    Dickey-Fuller: -15.1009
  P VALUE:
    0.01 

Description:
 Sun Jun 17 19:11:07 2018 by user: Telmo Diniz

A partir daí foi feito o cálculo das funções de Autocorrelação (FAC) e Autocorrelação Parcial (FACP) para determinar o número de defasagens estatisticamente relevantes da série que impactam em seu valor presente. Para isso foram usadas as funções acf e pacf.

acf_daily_Return = stats::acf(daily_return, na.action = na.pass, plot = FALSE, lag.max = 15)
plot(acf_daily_Return, main = "", ylab = "", xlab = "Defasagem")
title("Função de Autocorrelação (FAC)", adj = 0.5, line = 1)

pacf_daily_return = stats::pacf(daily_return, na.action = na.pass, plot = FALSE, lag.max = 15)
plot(pacf_daily_return, main = "", ylab = "", xlab = "Defasagem")
title("Função de Autocorrelação Parcial (FACP) dos Retornos", adj = 0.5, line = 1)

Como podemos ver nos gráficos, a FACP apresenta uma defasagem estatisticamente relevante, o que sugere uma ordem P de no máximo 1 e a FAC não apresenta defasagens relevantes, o que sugere uma ordem Q igual a 1.

A partir daí, é necessário calcular a Verossimilhança e AIC de todas as combinações de modelos com ordens P menores ou iguais a 1 e ordens Q menores ou iguais a 1 e escolher o modelo que apresentar o menor valor de AIC.

pars = expand.grid(ar = 0:1, diff = 0, ma = 0:1)
pars = pars[-1,]
modelo = list()
for (i in 1:nrow(pars)) {
  tryCatch({
    modelo[[i]] = arima(daily_return, order = unlist(pars[i, 1:3]), method = "ML", transform.pars = FALSE)},
    error = function(e) NULL)
}
log_verossimilhanca = list()
for (i in 1:length(modelo)) {
  log_verossimilhanca[[i]] = modelo[[i]]$loglik
}
aicarma = list()
for (i in 1:length(modelo)) {
  aicarma[[i]] = modelo[[i]]$aic
}
especificacao = paste0("ARMA",pars$ar,pars$diff,pars$ma)
quant_paramentros = pars$ar+pars$ma
tamanho_amostra = rep(length(daily_return), length(modelo))
resultado = data.frame(especificacao,ln_verossimilhanca = unlist(log_verossimilhanca),quant_paramentros,tamanho_amostra, aic = unlist(aicarma))
knitr::kable(resultado[], format = "pandoc", digits = c(0,2,0,0,2), align = 'c')

 especificacao    ln_verossimilhanca    quant_paramentros    tamanho_amostra      aic    
---------------  --------------------  -------------------  -----------------  ----------
    ARMA100            1549.13                  1                  585          -3092.26 
    ARMA001            1549.31                  1                  585          -3092.61 
    ARMA101            1550.44                  2                  585          -3092.87 

A partir da tabela podemos concluir que o melhor modelo é um ARIMA (1,0,0), lembrando que diferimos a série itau uma vez para daily_return, o que significa que o modelo correto é um ARIMA (1,1,0). Agora é preciso verificar o comportamento do ruído do modelo, de forma a garantir que seja um ruído branco.

Por fim, foi usada a função autoplot do pacote ggplot2 para plottar a previsão do modelo para o próximo período.

autoplot(forecast(modelo[[1]], h=4), main = "Previsão do modelo ARIMA(1,1,0)")