Erros, Warnings e letras vermelhas

Tão importante quanto saber fazer a coisa certa é saber quando temos um problema. Existem dois tipos de erro:

• Error: esses são de fato erros, e o R não consegue proceder. O que você escreveu tem algum problema grave e não pode ser executado.
• Warning: o R conseguiu fazer o que você pediu, mas alguma coisa esquisita aconteceu, e o R está te avisando.

Veja que em ambos os casos a mensagem vai aparecer em vermelho no console. Error e warnings vem das mais diferentes formas (afinal, é possível errar de diferentes formas). Alguns são claros, como tentar somar um número e uma letra:

2 + "a"

## Error in 2 + "a": argumento não-numérico para operador binário

(Apesar de operador binário não ser exatamente óbvio)

Outros são mais misteriosos, e quando encontramos um erro impossível de entender, coloque o output do console no google. Typos geram muitos erros (ex.: esquecer o parêntese, errar o nome de um objeto). Em um código longo, descobrir o erro pode ser difícil. Existem ferramentas para isso, e o próprio source do RStudio indica linhas com algum erro com um x vermelho do lado do número da linha. Em geral, em um código muito longo, não é uma estratégia ruim escrever uma parte do código, testar e debugar (retirar os erros) ela. Funções ajudam a fazer isso: trataremos delas mais a frente.

É importante notar que, apesar de uma linha de código que devolva um warning não está errada, ela pode muito bem fazer uma coisa que você não quer. Assim, um warning surpresa deve ser visto com bastante cuidado e normalmente você não deve ignorar um warning.

Veja que essas mensagens não são as únicas que aparecem em vermelho no R: quando o R instala ou carrega um pacote ele exibe algumas coisas em vermelho¹. Não se preocupe: isso não quer dizer que o R encontrou algum erro ou warning ao instalar o pacote, é apenas o comportamento padrão do R.

install.packages("nome-do-pacote")

library(nome-do-pacote)

install.packages("ctv")
library(ctv)
install.views("Econometrics")

Ajuda

Em 90% do tempo você vai precisar de olhar o help: seja para lembrar que opções estão disponíveis em um comando, ou lembrar exatamente o que o comando faz, ou descobrir qual o comando para fazer alguma coisa específica que você só tem uma palavra chave: Como gera números saídos de uma distribuição normal mesmo?’’²

Se você sabe o comando e quer ler o help deste comando, basta fazer ?comando. Por exemplo, se você quiser saber o que o comando rnorm faz e quais as opções ele oferece, basta digitar no console ?rnorm. Se você não sabe o nome do comando, mas quer saber todos os comandos que estão relacionados a uma determinada palavra chave, use ??palavra. Por exemplo, se você quiser saber todos os comandos que envolvem a distribuição normal, basta usar ??normal. Observe que, se você tiver muitos pacotes, isto pode demorar um pouco, já que exige que o R procure cada pacote que se referencie àquela palavra.

Você também pode estar interessado em todos os comandos de um pacote. Neste caso, a melhor solução é é usar help(package=``nome-do-pacote'').

Comentando

Em geral, quando se escreve um código, é importante explicar o que algumas partes do código fazem. Isso não é exclusivo para o caso em que o código vai ser distribuído: pelo contrário, o você do futuro vai agradecer muito o você do passado se você comentar o seu código. Para comentar o código no R, usa-se o #. Tudo depois do # vira um comentário é ignorado pelo R. Assim, suponha que queremos explicar o que o parâmetro a abaixo é: suponha que ele vem de um outro estudo, de Cleese et. al. (1975). Então:

    a <- 1 #Retirado de Cleese et. al. (1975)

A regra geral para comentários é: eles não devem ser óbvios (ex.: 1 + 1 # somando dois números) mas devem ajudar o leitor - que eventualmente será você mesmo - a entender o que está sendo feito - especialmente em momentos mais obscuros.

Objetos*

Como muitas linguagens de programação, existem vários tipos de objetos no R. Objetos são maneiras diferentes de guardar os dados. Os mais usados e mais comuns são:

• Vetor
• Matriz
• Dataframe
• Lista

Em geral, para dar um nome a um objeto, usamos uma setinha, <-, que é o sinal de menor e o menos. Podemos usar o igual também, mas é preferível a seta. Então, se fizermos a <- 1 e digitarmos a, o R vai mostrar 1. Veja que por padrão o R não mostra o valor de objetos que você acabou de criar, e vai mostrar apenas se você pedir.

a <- 1
a

## [1] 1

Vetores são só coleções unidimensionais de “coisas”. Para criar um vetor, basta usar c() e separar os elementos por vírgula. Suponha que queremos listar os 4 primeiros números primos, então poderíamos fazer:

primos <- c(2,3,5,7)

Assim, ao digitar no console primos:

primos

## [1] 2 3 5 7

Podemos fazer operações com vetores, como somar, subtrair, multiplicar e dividir. Observe que multiplicar um vetor não é multiplicar uma matriz: o R vai multiplicar elemento a elemento. Então c(1,2,3)*c(1,2,3) gera como resultado:

c(1,2,3)*c(1,2,3)

## [1] 1 4 9

Para multiplicar dois vetores como multiplicariamos usualmente, usamos %*%:

c(1,2,3)%*%c(1,2,3)

##      [,1]
## [1,]   14

Podemos agrupar vetores em matrizes, usando os comandos cbind e rbind, que transforma cada vetor em uma coluna ou em uma linha, respectivamente. Assim, se fizermos cbind(c(1,2),c(3,4)), teríamos:

cbind(c(1,2),c(3,4))

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

Se usarmos o rbind(c(1,2),c(3,4)), teríamos:

rbind(c(1,2),c(3,4))

##      [,1] [,2]
## [1,]    1    2
## [2,]    3    4

Existe outra maneira de fazer matrizes, com o comando matrix. Veja que as operações *e %*% funcionam assim como elas funcionam com vetores: *multiplica elemento a elemento a matriz e %*%multiplica a matriz da maneira usual:

A <- matrix(c(1,2,3,4),ncol = 2)
B <- matrix(c(1,1,1,1),ncol = 2)
A*B

##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4

A%*%B

##      [,1] [,2]
## [1,]    4    4
## [2,]    6    6

Observe que matrizes forçam todos os elementos a serem do mesmo tipo. Suponha que você quer fazer uma matriz com nomes e notas de alunos, e quer tirar a média das notas. O R vai acusar um erro, porque as notas não serão do tipo numérico, e sim do tipo character, que é o tipo do nome dos alunos. Para contornar este problema e poder agregar vetores com diferentes tipos - como é o caso do exemplo de notas de alunos e notas é que existe o objeto do tipo dataframe. O comando data.frame funciona como um cbind, que permite diferente juntar vetores de diferentes tipos em uma “matriz”.

Observe que, para formar matrizes e dataframes, os vetores tem que ter o mesmo tamanho, o que nem sempre é possível ou desejável. Neste contexto, existem as listas, que são um anything goes. Você pode ter uma lista de vetores, uma lista de variáveis, uma lista de listas etc. Além disto, você pode dar nomes as coisas dentro dela e chamar pelo nome. Por exemplo, se fizermos:

um.teste <- list(Ola = "Olá", Dia = date())

E depois digitarmos um.teste$Dia, ele deve exibir qual a data atual. Digite um.teste$Ola e ele deve exibir um Olá na sua tela:

um.teste$Dia

## [1] "Fri May  1 11:40:31 2020"

um.teste$Ola

## [1] "Olá"

Observem que eu usei um “Olá”, entre aspas na lista. Aspas são bastante importante. Por exemplo, faça: c("Bom","Dia"). O R vai mostrar na tela as duas palavras. Agora, suponha que você esquecesse as aspas no Dia. Agora, o R dará um erro: ele afirma que o objeto Dia não foi encontrado.

Assim, se quisermos digitar palavras, frases, letras, devemos colocar eles entre aspas. Caso contrário, o R vai buscar o objeto com aquele nome. De maneira bastante grosseira, expressões desse tipo são chamadas de string.

Familiarizados com como instalar, pedir ajuda, e os principais objetos do R, podemos proceder para a primeira etapa de qualquer análise de dados: como inserir os dados no R.

Testes de hipótese

Suponha que queremos testar a significância conjunta de x.2 e x.3. Precisamos fazer um teste F. Uma maneira é estimar um novo modelo, que chamaremos de reg.2, só com o x.1: reg.2 <- lm(y ~ x.1)

Agora, para testar a significância conjunta de x.2 e x.3 basta fazer anova(reg,reg.2) e o R reportará os valores do teste (incluindo o p valor)

Testes mais gerais também estão disponíveis, pelo pacote car⁶. O comando é linearHypotesis com o primeiro argumento sendo o objeto com o modelo. O segundo comando é a hipótese que estamos testando, entre aspas, e entramos ele de maneira extremamente intuitiva: suponha que no exemplo acima queremos testar se o coeficiente de x.2 é igual ao coeficiente de x.4. Nesse caso, o comando se resumiria a linearHypotesis(reg,"x.2 = x.4"). Várias opções podem ser usadas, como usar estimadores de erro padrão robustos a heterocedasticidade. Recomendamos que o leitor olhe o help do comando no R.

Erros padrões robustos

Na presença de heterocedasticidade ou autocorrelação, os erros padrões usuais não são confiáveis. Infelizmente, o comando summary não exibe erros robustos por default e nem permite alterar os erros padrões exibidos. Felizmente, o pacote lmtest traz uma opção para o sumário padrão do R. A primeira coisa a fazer é carregar o pacote sandwich (library(sandwich)). O comando é coeftest e a sintaxe é curiosa. No caso do nosso exemplo, se quisermos obter erros robustos a heterocedasticidade:

coeftest(reg, vcov. = vcovHC(reg))

Uma explicação: O comando coeftest chama o comando vcovHC do pacote sandwich. Por sua vez, o comando vcovHC precisa saber qual o modelo que vai ter erros robustos. Por isso uma função que recebe uma função. Veja que se quisermos erros robustos a heterocedasticidade e autocorrelação, o comando vira vcovHAC.

Logs

Muitas vezes queremos fazer as regressões não com as variáveis em nível, mas com as variáveis em log. Nesse caso, suponha que queremos a variável dependente - y - em log. Para isso, basta fazer log.y <- log(y) e a variável log.y vai ser a versão em log da variável y. Você pode reescrever a regressão como lm(log.y ~ x.1 + x.2 + x.3).

Agora, pode ocorrer de em alguns casos o vetor y ter algum elemento zero. Mas $\log \left(0\right)=-\infty$. O R tem um elemento Inf (e - Inf) para esses casos, mas o comando lm vai acusar um erro ao receber um vetor com algum elemento Inf. A solução é trocar esse valor para alguma coisa, como um NA, que o R vai ignorar⁷. Para fazer isso, suponha que o vetor com os Inf seja o log.y do paragrafo anterior. Precisamos explicar para o R quais casos nós queremos substituir, e isso é incrivelmente fácil: assim como podemos usar log.y[1] para escolher o primeiro elemento do vetor log.y, podemos colocar entre colchetes uma condição, por exemplo os elementos do vetor log.y que são iguais a infinito. Já tratamos disso: log.y == -Inf faz esse trabalho. Assim, se quisermos substituir os elementos de log.y que são iguais a -Inf por NA, basta usar o seguinte trecho de código:

 log.y[log.y == -Inf] <- NA

Probits e Logits

Probits e logits também são úteis quando nossa variável dependente é uma dummy. Sempre podemos usar um modelo de probabilidade linear, e nesse caso o comando a ser utilizado é o lm. Mas em casos que queremos usar um probit ou logit, precisamos recorrer ao glm, que tem sintaxe muito parecida com o lm. Mas além de especificar as variáveis dependentes e independentes, também precisamos especificar o tipo de regressão, basicamente a distribuição da variável dependente. No caso de probits e logits, a variável dependente tem distribuição binomial. Depois, temos que especificar se a função de probabilidade da variável dependente é probit ou logit. Os comandos para estimar probits e logits são ilustrados abaixo:

mod.1 <- glm(y ~x.1 + x.2 + x.3, family = binomial(link = 'probit'))
mod.2 <- glm(y ~x.1 + x.2 + x.3, family = binomial(link = 'logit'))

E como de praxe, podemos usar o comando summary para obter os coeficientes, desvios padrões e estatísticas t.

Variável instrumental

Métodos de variável instrumental são muito úteis e populares, especialmente em casos de endogenidade. Existem várias implementações, mas para o mínimos quadrados de dois estágios usual, o pacote AER oferece um comando ivreg. A sintaxe é similar ao lm, mas com uma alteração na formula para inserir os instrumentos, que são separados dos regressores por |.

Por exemplo, suponha que temos a variável dependente y, as variáveis endogenas x.1 e x.2, a variável exogena x.3, e os instrumentos z.1 e z.2. Nesse caso, o ivreg seria usado:

modelo <- ivreg(y ~x.1 + x.2 + x.3|z.1+z.2+x.3)

E podemos usar o summary para ver o valor dos coeficientes, erros padrão e estatísticas t. Observe que o pacote não mostra o valor do teste F para o primeiro estágio nem de teste de sobre identificação.

Dados em painel

Em muitas aplicações usamos dados em painel - i.e., com dimensão temporal e cross section. Em geral, esse tipo de aplicação acaba envolvendo o uso de efeitos fixos. Existem duas maneiras de fazer: usando o pacote plm ou “na mão”, usando o lm. Não há nenhuma vantagem de usar o lm “na mão”, em geral, exceto em casos que temos mais de duas dimensões ou por algum motivo o plm não funciona. Exploraremos primeiro o uso do plm, que deve satisfazer a maioria dos usuários. A solução na mão vem depois e pode ser ignorada sem perda de continuidade.

Supondo que o pacote já foi instalado e carregado, precisamos (i) explicar para o pacote quais colunas são as colunas com efeitos fixos de tempo e unidade e (ii) rodar a regressão propriamente dita. Suponha, como de praxe, que temos um dataframe carregado no R com nome “dados”. Suponha que as colunas com as datas e um índice para unidades se chamam datas e unidades, respectivamente. Veja que essas colunas podem estar em formato de carácter, e que isso não deve nos preocupar no momento: poderia ser o caso de a unidade ser estados do país e o código da unidade ser o código do estado (DF,RJ,SP,…). O pacote plm disponibiliza o comando plm.data, que converte um data frame de forma que quando rodarmos a regressão, o R saiba quem são os efeitos fixos. Assim, vamos criar um dataframe chamado dd:

dd <- plm.data(dados,c('unidades','datas'))

O primeiro argumento da função é o data.frame a ser convertido, que contém as colunas para as quais criaremos efeitos fixos. O segundo argumento da função são as colunas com os efeitos fixos. Veja que a ordem é unidade e depois a variável temporal. Agora, a estimação do modelo pode ser feita usando o comando plm, que tem sintaxe muito similar ao lm: passamos uma formula com a variável dependente e as independentes, informamos a base de dados - que nesse caso é o objeto dd, não o objeto dados. Mas temos algumas novas opções: o método de estimação (em geral estamos interessados em within, o padrão), mas mais relevante é que efeitos fixos queremos colocar: só para indivíduo, só para tempo ou ambos. Abaixo, mostramos a sintaxe para cada um dos casos, respectivamente:

mod.1 <- plm(y ~x.1 + x.2,data = dd, effect = 'individual')
mod.2 <- plm(y ~x.1 + x.2,data = dd, effect = 'time')
mod.3 <- plm(y ~x.1 + x.2,data = dd, effect = 'twoways')

Como de praxe, podemos usar o comando summary para obter um sumário da regressão.

modelo <- lm(y ~x +tempo + unidade, data = dados)

Só podemos fazer análise de dados se tivermos… dados. Este capítulo ensina a colocar os dados no R e a manipular eles como para criar variáveis necessárias ou limpar os dados antes de iniciarmos a análise.

Esse capítulo não tenta ser enciclopédico nem detalhista: pelo contrário, ele omite muitas coisas. A omissão mais grave é, sem dúvida alguma, os pacotes do Tidyverse. A omissão se deve a ignorância do autor em usar estes pacotes.

read.csv2("C:/Usuário/Autor/Manual/dados.csv")

dados <- read.csv2("C:/Usuário/Autor/Manual/dados.csv")

dados <- read.csv2(file.choose())

dados <- read.dta("C:/Usuário/Autor/Manual/dados.dta")

library(BETS)

busca <- BETSsearch("IPCA")

## 
## BETS-package: Found 51 out of 18706 time series.

#View(busca)
IPCA <- BETSget(433, from = "2002-01-01")
plot(IPCA)

sal.min <- w < 678

O básico

Suponha que você, usando os métodos do capítulo 2, inseriu uma série temporal no R. O R não sabe, a priori, que os dados são uma série temporal. Você precisa contar isso a ele, e o comando que faz isso é o ts(). O ts recebe a série, a data de ínicio e a frequência. A frequência é como você dividiu o ano: 4 se o dado for trimestral, 12 se for mensal…

Por exemplo, podemos gerar uma série de variáveis aleatórias da normal (um ruído branco) e transformar em série temporal mensal começando em janeiro de 2000:

serie <- rnorm(1000)
serie <- ts(serie,start =c(2000,01), freq = 12)

Veja que o comando ts() é excessivamente engessado: os dados tem que ter uma frequência fixa, expressa como uma fração do ano. O pacote zoo extende bastante as capacidades do R em lidar com séries temporais, inclusive com séries irregulares.

Atenção: Séries Temporais e o lm Você não deve passar um objeto de série temporal para o lm(), já que o lm vai ignorar o formato de série temporal. Assim, estimar um modelo AR(1) usando lm(y ~ lag(y)) vai gerar uma regressão com coeficiente 1 para lag(y) e $R^2=1$. De fato, a regressão feita foi y em y - o que não é uma regressão muito emocionante.

##ARIMAs

Com uma série devidamente construída para ser um objeto ts - como nós fizemos acima- podemos tentar estimar algum modelo. O modelo base de séries temporais é o ARIMA. O R base já vem com muitas funções para lidar com isso, mas o pacote forecast extende bastante as capacidades do R em lidar com esse tipo de série. O primeiro passo para estimar um modelo Arima é obter a função de autocorrelação (FAC) e a função de autocorrelação parcial (FACP): elas são Acf e Pacf. Veja que existem versões na base do R que se chamam acf e pacf (notem que lá é com maiúscula e aqui com minúscula). A diferença fundamental entre os dois é que a Acf e a acf (e também a Pacf e a pacf) é que a primeira exclui a autocorrelação no momento 0 - que é trivialmente 1.

Uma vez conhecendo o formato da FAC e da FACP, podemos estimar o ARIMA. O comando para estimar um modelo ARIMA é Arima - e novamente, existe um arima com minúscula que é da base do R. O Arima basicamente recebe duas coisas, a série e a ordem do modelo(isso é, se o modelo é um ARMA(1,1), AR(1), MA(1) etc). Vamos gerar um exemplo de um AR(1) com dados simulados e obter a FAC e FACP e estimar o modelo Arma sugerido:

library(forecast)

u <- rnorm(500)
y <- rep(0,500) #nossa futura série

y[1] <- u[1]

for(i in 2:500){
  y[i] <- 0.6*y[i-1] + u[i] #um AR(1) com coeficiente 0.6
}

Acf(y)

Pacf(y)

modelo <- Arima(serie,order=c(1,0,0))
summary(modelo)

## Series: serie 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1     mean
##       -0.0312  -0.0157
## s.e.   0.0316   0.0315
## 
## sigma^2 estimated as 1.06:  log likelihood=-1447.14
## AIC=2900.28   AICc=2900.31   BIC=2915
## 
## Training set error measures:
##                        ME     RMSE       MAE    MPE     MAPE      MASE          ACF1
## Training set 1.500234e-05 1.028603 0.8091188 98.585 111.1065 0.7103052 -0.0002764195

Veja que o objeto modelo trás os coeficientes estimados, o erro padrão e alguns diagnósticos úteis como critérios de informação.

Em algumas situações pode ser muito difícil inferir o modelo certo a partir da FAC e da FACP. O comando auto.arima, do pacote forecast, seleciona um modelo a partir de algum critério de informação. Vamos ilustrar o ponto gerando uma série x que é um ARMA(1,2):

e <- rnorm(1000)

x <- rep(0,1000)

for(j in 1:998){
  x[j+2] <- 0.5*x[j+1] + e[j+2] - 0.3*e[j+1] + 0.4*e[j]
}

x <- x[500:1000]
x <- ts(x, start = c(1999,05), freq = 12)

Acf(x)

Pacf(x)

auto.arima(x,ic = "bic")

## Series: x 
## ARIMA(1,0,2)(1,0,0)[12] with zero mean 
## 
## Coefficients:
##          ar1      ma1     ma2     sar1
##       0.4648  -0.2850  0.4401  -0.0243
## s.e.  0.0766   0.0738  0.0418   0.0455
## 
## sigma^2 estimated as 0.8746:  log likelihood=-675.65
## AIC=1361.29   AICc=1361.41   BIC=1382.37

Nesse caso o auto.arima acertou, mas nem sempre isso ocorre.

VARs

Um VAR, teoricamente, é apenas uma generalização do AR. Ainda assim, do ponto de vista computacional, eles são distintos, e o VAR tem seu próprio conjunto de pacotes no R. O mais importante deles é o vars.

Começamos juntando todas as séries que queremos estimar o VAR em uma matriz (use o cbind() para isso). O passo seguinte é escolher a ordem do VAR - geralmente usando algum critério de informação. O comando VARselect faz isso e apresenta alguns critérios de informação e a quantidade de lags que minimizam cada um.

O comando que faz a estimativa per se é o VAR. Ele recebe a matriz com as séries e quantos lags você quer que sejam usados - ou o critério de informação a ser usado na hora de fazer a estimativa.

Por último, queremos recuperar a resposta dinâmica de cada uma das variáveis a um choque (não só um choque na própria variável, como o efeito cruzado de um choque em outra variável).

Vamos gerar um VAR(1) com duas variáveis apenas para ilustrar o uso do pacote:

library(vars)

T <- 500 #número de períodos
N <- 2 #número de variáveis

u <- matrix(rnorm(T*N),nrow = N,ncol = T)

x <- matrix(0,nrow = N, ncol = T)

A <- rbind(c(0.3,-0.2),c(0.6,0.2))

for(j in 2:T){
  x[,j] <- A%*%x[,j-1] + u[,j]
}

x <- x[,400:500]
x <- t(x)
colnames(x) <- c("x","y")

VARselect(x) 

## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      1      1      1      1 
## 
## $criteria
##                1         2         3         4         5        6         7         8         9        10
## AIC(n) 0.5073838 0.5664775 0.6364448 0.7148971 0.8004834 0.804686 0.8853781 0.9301667 0.9188877 0.9939262
## HQ(n)  0.5741733 0.6777934 0.7922872 0.9152658 1.0453785 1.094108 1.2193260 1.3086411 1.3418885 1.4614533
## SC(n)  0.6729349 0.8423961 1.0227309 1.2115506 1.4075043 1.522074 1.7131339 1.8682901 1.9673785 2.1527844
## FPE(n) 1.6610195 1.7624397 1.8909018 2.0466301 2.2319134 2.244812 2.4386605 2.5575410 2.5379366 2.7479537

modelo <- VAR(x,p = 1)
plot(irf(modelo, impulse = "x", response = c("x","y")))

Eu só pedi o plot do choque da primeira variável sobre as duas variáveis por que isso é uma ilustração. Fazer plot(irf(modelo)) em uma seção do R, ele vai plotar os choques de todas as variáveis sobre todas as variáveis.

Raiz unitária

O pacote urca nos trás testes de raiz unitária. O teste Dickey-Fuller, um dos mais populares, é chamado pelo ur.df(). Vamos gerar um passeio aleatório para mostrar:

 x <- cumsum(rnorm(1000))
summary(ur.df(x))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0769 -0.7404 -0.0143  0.6389  3.4537 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## z.lag.1    -2.254e-05  7.896e-04  -0.029    0.977
## z.diff.lag -1.212e-02  3.169e-02  -0.383    0.702
## 
## Residual standard error: 1.001 on 996 degrees of freedom
## Multiple R-squared:  0.0001483,  Adjusted R-squared:  -0.001859 
## F-statistic: 0.07388 on 2 and 996 DF,  p-value: 0.9288
## 
## 
## Value of test-statistic is: -0.0285 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Veja que o valor crítico para o teste Dick Fuller não é o valor usual da estatística t, mas sim o valor exibido na parte debaixo da tabela do sumário do teste. Nesse caso, a qualquer nível de significância, nós não rejeitamos a hipótese de raiz unitária. Vamos testar para um caso estacionário:

u <- rnorm(2000)
y <- rep(0,2000)

for(i in 2:2000){
  y[i] <- 0.5*y[i-1] + u[i]
}

y <- y[1000:2000]

summary(ur.df(y))

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7812 -0.6434  0.0092  0.6571  3.2324 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1    -0.501212   0.031742 -15.790   <2e-16 ***
## z.diff.lag -0.002658   0.031612  -0.084    0.933    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9926 on 997 degrees of freedom
## Multiple R-squared:  0.2511, Adjusted R-squared:  0.2496 
## F-statistic: 167.1 on 2 and 997 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -15.7901 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Desassonalizando

Dados de séries temporais, não raramente, apresentam sazonalidade. Por exemplo, gasto de energia elétrica tende a ser maior nos meses de dezembro a fevereiro, devido ao verão. Retirar sazonalidade é importante em muitas análises.

Um método padrão é colocar dummies para as unidades de tempo (uma para cada mês se o dado for mensal, uma para cada trimestre se for trimestral etc) e usar o resíduo dessa regressão somado a média da série (já que o resíduo tem média zero, por construção). Criar as dummies “no braço” pode ser tedioso, mas felizmente o pacote forecast trás o comando seasonaldummy() que cria as dummies automaticamente para a série.

energia <- BETSget(1406, from = "2002-01-01")
dum <- seasonaldummy(energia)
mod <- Arima(energia, xreg = dum)
des <- resid(mod) + mean(energia)
plot(energia)
lines(des, col = 2)
legend("topleft", legend = c("C/Sazonalidade", "Sem Sazonalidade"), lty = c(1,1), col = c(1,2))

Outra maneira comum de dessazonalizar é usando o X13, um programa do governo americano. O X13 pode ser acessado direto do R usando o pacote seasonal. O comando que acessa o X13 é o seas. O X13 são, na verdade, dois programas: um que é o X13 e o outro que é o SEATS. Ambos tem a mesma função: dessazonalizar. O X13 vem com todo tipo de método automático para detectar outliers, fazer transformações nas séries e uma infinidade de outras coisas. Nesse caso, nós vamos desligar todas essa opções:

library(seasonal)
modelo2 <- seas(energia, transform.function = "none", regression.aictest = NULL, outlier = NULL)

O comando final obtém a série dessazonalizada:

plot(energia)
lines(final(modelo2), col = 2)
legend("topleft", legend = c("C/Sazonalidade", "Sem Sazonalidade"), lty = c(1,1), col = c(1,2))

Filtro Hodrick-Prescott

O filtro Hodrick-Prescott (HP) é utilizado em séries não estacionárias quando queremos separar a tendência do componente ciclíco. Ele é polêmico, mas ainda é amplamente usado. No R, o pacote mFilter implementa ele e alguns outros. Vamos aplicar na série de energia:

library(mFilter)
filtrado <- mFilter(energia, "HP")
plot(filtrado)

Eu posso acessar a tendência e o componente cíclico usando o filtrado$trend e filtrado$cycle, no exemplo acima - que seria particularmente útil se eu quisesse utilizar os dados de cíclo para alguma estimação.