Automatizando Testes de Raiz Unitária (ADF, PP e KPSS)

Testes de Raiz Unitária

O intuito desse post é mostrar como automatizar os testes de raiz unitária para um conjunto de séries temporais. Para isso, iremos utilizar os testes Dickey-Fuller Aumentado (ADF), Phillips-Perron (PP) e Kwiatkowski-Phillips-Schmidt-Shin (KPSS). A dificuldade em automatizar é que cada um desses testes precisa testar a sua parte determinística como sendo ou não estatisticamente significativa. Além disso, o pacote urca reporta o valor da estatística do teste e os valores críticos. Muitas vezes, ao reportar o teste em um artigo, desejamos ter uma tabela com o valor do teste e um indicativo do p-valor associado.

Dessa forma, iremos criar três funcões para automatizar cada um dos três testes e aplicá-las para um conjunto de dados artificialmente gerado para ser não estacionário.

Geração de Séries Temporais Não Estacionárias

Iremos gerar 10 distintas séries temporais em formato de data.frame (ao invés de ts()). Todavia, antes disso iremos carregar os pacotes dplyr, purrr, tibble, urca e kableExtra. O primeiro pacote é para manipular uma base de dados, o segundo servirá para aplicar uma função em cada uma das colunas de um data.frame, a terceira é para formar um tipo específico de classe de dados, o quarto é para obter os testes de raiz unitária e o último é para gerar uma tabela em html no final do post.

Além disso, vamos fixar a semente aleatória, set.seed(123), para fim de reprodução dos resultados.

Os dados serão gerados como passeios aleatórios, sendo \(t=[1,\ldots,1000]\): \[y_t=y_{t-1}+\varepsilon_{t},\;\;\;\varepsilon_t\sim N(0,1)\]

No intuito de apenas visualizar melhor os dados, transformamos em ts() e plotamos. Posteriormente, todos os testes são feitos na classe de data.frame().

rm(list=ls())

library(dplyr)
library(purrr)
library(tibble)
library(urca)
library(kableExtra)

set.seed(123)

# Geração dos Dados
data <- data.frame(
  serie1  = cumsum(rnorm(1000)),
  serie2  = cumsum(rnorm(1000)),
  serie3  = cumsum(rnorm(1000)),
  serie4  = cumsum(rnorm(1000)),
  serie5  = cumsum(rnorm(1000)),
  serie6  = cumsum(rnorm(1000)),
  serie7  = cumsum(rnorm(1000)),
  serie8  = cumsum(rnorm(1000)),
  serie9  = cumsum(rnorm(1000)),
  serie10 = cumsum(rnorm(1000))
)

# Tranformação dos Dados em Série Temporal apenas para a sua Visualização
data_ts <- ts(data, start = c(1940,1), freq= 12)

# Plot dos Dados em ts().
plot(data_ts, main = "Séries Temporais Geradas por Simulação")

Automatizando o Teste de ADF

Como já mencionado, o problema dos testes é a necessidade da melhor especificação para o termo determinístico. No caso da função ADF, ela pode ser none, drift, trend. O que faremos é começar pela especificação de trend. Caso o p-valor do termo seja maior do que 5%, iremos descartar esse modelo e rodar novamente com uma especificação para drift. Caso o p-valor do termo seja maior que 5%, iremos rodar novamente sem termo determinístico.

Para além disso, faremos comparações entre a estatística obtida pelo teste e o seu valor crítico. Se a comparação for maior do o valor crítico para 10%, deixaremos em branco o valor da variável p_valor e a conclusão é que os dados são Não Estacionários para aquela série testada. O mesmo teste condicional é feito para menor que 10% e se atendida a conclusão é de p_valor = "*" e Estacionária. Repete-se para 5% e 1%. Atribuindo-se p_valor = "**" e p_valor = "***".

O teste é especificado para a escolha automática de defasagens conforme o critério de informação de Akaike.

Ao fim, o teste é reportado com uma tabela, tibble(), onde a estatística do teste é reportada no formado de texto através do mergeamento do número obtido arredondado para 3 casas decimais e a quantidade de asteriscos correspondente ao seu p-valor. Revela-se também o termo determinístico escolhido e conclusão obtida pelo teste. Deixou-se oculto, mas poderia também ser reportado os valores críticos do teste.

# Função para aplicar o teste ADF e extrair estatísticas
teste_adf <- function(x) {
  adf <- ur.df(x, type = c("trend"), selectlags = c("AIC"))
  
  if (adf@testreg$coefficients["tt",4] > 0.05) {
  adf <- ur.df(x, type = c("drift"), selectlags = c("AIC"))  
  }

  if (adf@testreg$coefficients["(Intercept)",4] > 0.05) {
    adf <- ur.df(x, type = c("none"), selectlags = c("AIC"))  
  }
  
  if(adf@teststat[1] > adf@cval[1,3]) {
    p_valor = ""
    Conclusao = "Não Estacionário"
  }
  if(adf@teststat[1] <= adf@cval[1,3]) {
    p_valor = "*"
    Conclusao = "Estacionário"
  }
  if(adf@teststat[1] <= adf@cval[1,2]) {
    p_valor = "**"
    Conclusao = "Estacionário"
  }
  if(adf@teststat[1] <= adf@cval[1,1]) {
    p_valor = "***"
    Conclusao = "Estacionário"
  }

  tibble(
    ADF = paste(round(adf@teststat[1],3),p_valor),
    #c_val_1    = adf@cval[1,1],
    #c_val_5    = adf@cval[1,2],
    #c_val_10   = adf@cval[1,3],
    Model_ADF    = adf@model,
    Conclusao_ADF = Conclusao
  )
}

Automatizando o Teste PP

Tal como no teste ADF, no teste PP temos a necessidade da melhor especificação para o termo determinístico. No caso do teste PP, ele pode ser constant ou trend. O que faremos é começar pela especificação de trend. Caso o p-valor do termo seja maior do que 5%, iremos descartar esse modelo e rodar novamente com uma especificação para constant.

Tal como no teste ADF, faremos comparações entre a estatística obtida pelo teste e o seu valor crítico. Se a comparação for maior do o valor crítico para 10%, deixaremos em branco o valor da variável p_valor e a conclusão é que os dados são Não Estacionários para aquela série testada. O mesmo teste condicional é feito para menor que 10% e se atendida a conclusão é de p_valor = "*" e Estacionária. Repete-se para 5% e 1%. Atribuindo-se p_valor = "**" e p_valor = "***".

O teste é especificado para a escolha de desagens curtas, através do comando lags = "short". A escolha do type = "Z-tau" é justamente para termos o report dos termos críticos do teste.

Ao fim, o teste é reportado com uma tabela, tibble(), onde a estatística do teste é reportada no formado de texto através do mergeamento do número obtido arredondado para 3 casas decimais e a quantidade de asteriscos correspondente ao seu p-valor. Revela-se também o termo determinístico escolhido e conclusão obtida pelo teste. Deixou-se oculto, mas poderia também ser reportado os valores críticos do teste.

# Função para aplicar o teste PP e extrair estatísticas
teste_pp <- function(x) {
  pp <- ur.pp(x, model = c("trend"), type = "Z-tau",lags = "short")
  
  if (pp@testreg$coefficients["trend",4] > 0.05) {
    pp <- ur.pp(x, model = c("constant"), type = "Z-tau", lags = "short")  
  }

  
  if(pp@teststat[1] > pp@cval[1,3]) {
    p_valor = ""
    Conclusao = "Não Estacionário"
  }
  if(pp@teststat[1] <= pp@cval[1,3]) {
    p_valor = "*"
    Conclusao = "Estacionário"
  }
  if(pp@teststat[1] <= pp@cval[1,2]) {
    p_valor = "**"
    Conclusao = "Estacionário"
  }
  if(pp@teststat[1] <= pp@cval[1,1]) {
    p_valor = "***"
    Conclusao = "Estacionário"
  }
  
  tibble(
    PP = paste(round(pp@teststat[1],3),p_valor),
    #c_val_1    = pp@cval[1,1],
    #c_val_5    = pp@cval[1,2],
    #c_val_10   = pp@cval[1,3],
    Model_PP    = pp@model,
    Conclusao_PP= Conclusao
  )
}

Automatizando o Teste KPSS

De maneira similar aos dois testes anteriores, no teste KPSS temos a necessidade da melhor especificação para o termo determinístico. Agora, ele pode ser mu ou tau, análogos a constante e constante mais tendência.

Uma nova complicação surge. Ocorre que no teste ur.kpss() não conseguimos obter no objeto criado o valor do p-valor para a tendência e a constante. Com isso, como o indicado é a inspeção dos dados pelo pesquisador, utilizaremos como alternativa o próprio valor da significância estatística da tendência no teste PP como critério para seleção de mu ou tau.

Além disso, devemos ter o cuidado que a hipótese nula do teste agora é pela estacionariedade. Portando isso modifica o modo de marcar os termos de asterisco que serão adicionados posteriormente, bem como pela conclusão pela estacionariedade e não estacionariedade.

Não obstante, o teste é especificado para a escolha de desagens curtas, através do comando lags = "short".

Ao fim, o teste é reportado com uma tabela, tibble(), onde a estatística do teste é reportada no formado de texto através do mergeamento do número obtido arredondado para 3 casas decimais e a quantidade de asteriscos correspondente ao seu p-valor. Revela-se também o termo determinístico escolhido e conclusão obtida pelo teste. Deixou-se oculto, mas poderia também ser reportado os valores críticos do teste.

# Função para aplicar o teste PP e extrair estatísticas
teste_kpss <- function(x) {
  
  # Como a função ur.kpss não reporta o p-valor da tendência
  # determinística, iremos utilizar o teste PP para reportar 
  # e selecionar se o melhor modelo KPSS é com "mu"ou "tau".
  kpss <- ur.kpss(x, type = "mu", lags = "short")
  pp <- ur.pp(x, model = c("trend"), type = "Z-tau",lags = "short")
  
  if (pp@testreg$coefficients["trend",4] <= 0.05) {
    kpss <- ur.kpss(x, type = "tau", lags = "short")  
  }
  
  if(kpss@teststat[1] > kpss@cval[1]) {
    p_valor = ""
    Conclusao = "Não Estacionário"
  }
  if(kpss@teststat[1] <= kpss@cval[1]) {
    p_valor = "***"
    Conclusao = "Estacionário"
  }
  if(kpss@teststat[1] <= kpss@cval[2]) {
    p_valor = "**"
    Conclusao = "Estacionário"
  }
  if(kpss@teststat[1] <= kpss@cval[4]) {
    p_valor = "*"
    Conclusao = "Estacionário"
  }
  
  tibble(
    KPSS = paste(round(kpss@teststat,3),p_valor),
    #c_val_1     = kpss@cval[1,1],
    #c_val_5     = kpss@cval[1,2],
    #c_val_10    = kpss@cval[1,3],
    Model_KPSS   = kpss@type,
    Conclusao_KPSS = Conclusao
  )
}

Aplicando os Teste Individualmente e Mergeando os Resultados

Nesta parte, tendo as funções dos testes ADF, PP, KPSS automatizadas já carregadas, faremos uso delas. Para isso, usaremos a função map_dfr() para aplicar em cada série do data.frame. Os resultados são reportados com o comando print() para cada teste no conjunto de dados. A agregação dos três é feita com o uso do pacote dplyr, através da função full_join() ela irá mergear intersectando as séries comuns. Ao fim, fazemos a apresentação da tabela com a função kable() no intuito de criar uma tabela em formato html.

# Olhando os Testes Individualmente nos Dados

# Aplicar o teste em todas as colunas
resultado_adf <- data %>%
  purrr::map_dfr(teste_adf, .id = "serie") 

print(resultado_adf)

## # A tibble: 10 × 4
##    serie   ADF        Model_ADF Conclusao_ADF   
##    <chr>   <chr>      <chr>     <chr>           
##  1 serie1  "-1.018 "  none      Não Estacionário
##  2 serie2  "0.668 "   none      Não Estacionário
##  3 serie3  "-2.347 "  drift     Não Estacionário
##  4 serie4  "-2.38 "   drift     Não Estacionário
##  5 serie5  "-3.362 *" trend     Estacionário    
##  6 serie6  "0.305 "   none      Não Estacionário
##  7 serie7  "0.106 "   none      Não Estacionário
##  8 serie8  "-2.339 "  drift     Não Estacionário
##  9 serie9  "-0.427 "  none      Não Estacionário
## 10 serie10 "0.832 "   none      Não Estacionário

# Aplicar o teste em todas as colunas
resultado_pp <- data %>%
  purrr::map_dfr(teste_pp, .id = "serie") 

print(resultado_pp)

## # A tibble: 10 × 4
##    serie   PP          Model_PP                 Conclusao_PP    
##    <chr>   <chr>       <chr>                    <chr>           
##  1 serie1  "-2.063 "   with intercept           Não Estacionário
##  2 serie2  "-2.596 "   with intercept and trend Não Estacionário
##  3 serie3  "-2.322 "   with intercept           Não Estacionário
##  4 serie4  "-2.495 "   with intercept           Não Estacionário
##  5 serie5  "-3.371 *"  with intercept and trend Estacionário    
##  6 serie6  "-3.433 **" with intercept and trend Estacionário    
##  7 serie7  "-1.416 "   with intercept           Não Estacionário
##  8 serie8  "-2.474 "   with intercept           Não Estacionário
##  9 serie9  "-1.365 "   with intercept           Não Estacionário
## 10 serie10 "-2.734 "   with intercept and trend Não Estacionário

# Aplicar o teste em todas as colunas
resultado_kpss <- data %>%
  purrr::map_dfr(teste_kpss, .id = "serie") 

print(resultado_kpss)

## # A tibble: 10 × 4
##    serie   KPSS      Model_KPSS Conclusao_KPSS  
##    <chr>   <chr>     <chr>      <chr>           
##  1 serie1  "4.446 "  mu         Não Estacionário
##  2 serie2  "1.734 "  tau        Não Estacionário
##  3 serie3  "4.925 "  mu         Não Estacionário
##  4 serie4  "2.342 "  mu         Não Estacionário
##  5 serie5  "0.888 "  tau        Não Estacionário
##  6 serie6  "1.396 "  tau        Não Estacionário
##  7 serie7  "3.479 "  mu         Não Estacionário
##  8 serie8  "3.81 "   mu         Não Estacionário
##  9 serie9  "10.467 " mu         Não Estacionário
## 10 serie10 "1.049 "  tau        Não Estacionário

### Agregação do Resultado dos Testes

Resultados <- resultado_adf %>%
              full_join(resultado_pp, by = "serie") %>%
              full_join(resultado_kpss, by = "serie")

Resultados %>%
  kbl(caption = "Tabela com os Testes de Raiz Unitária") %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                font_size = 14)

Tabela com os Testes de Raiz Unitária

serie	ADF	Model_ADF	Conclusao_ADF	PP	Model_PP	Conclusao_PP	KPSS	Model_KPSS	Conclusao_KPSS
serie1	-1.018	none	Não Estacionário	-2.063	with intercept	Não Estacionário	4.446	mu	Não Estacionário
serie2	0.668	none	Não Estacionário	-2.596	with intercept and trend	Não Estacionário	1.734	tau	Não Estacionário
serie3	-2.347	drift	Não Estacionário	-2.322	with intercept	Não Estacionário	4.925	mu	Não Estacionário
serie4	-2.38	drift	Não Estacionário	-2.495	with intercept	Não Estacionário	2.342	mu	Não Estacionário
serie5	-3.362 *	trend	Estacionário	-3.371 *	with intercept and trend	Estacionário	0.888	tau	Não Estacionário
serie6	0.305	none	Não Estacionário	-3.433 **	with intercept and trend	Estacionário	1.396	tau	Não Estacionário
serie7	0.106	none	Não Estacionário	-1.416	with intercept	Não Estacionário	3.479	mu	Não Estacionário
serie8	-2.339	drift	Não Estacionário	-2.474	with intercept	Não Estacionário	3.81	mu	Não Estacionário
serie9	-0.427	none	Não Estacionário	-1.365	with intercept	Não Estacionário	10.467	mu	Não Estacionário
serie10	0.832	none	Não Estacionário	-2.734	with intercept and trend	Não Estacionário	1.049	tau	Não Estacionário