Neste relatório, os dados utilizados foram obtidos a partir das bibliotecas sidrar , ipeaDataR, GetBCBData e rvest
Caso não tenha instalado a biblioteca, será necessário descomentar comando comando install.packages(“sidrar”) e install.packages(“ipeadatar”) no código R abaixo:
#install.packages("sidrar")
#install.packages("ipeadatar")
#install.packages("GetBCBData")
#install.packages("lubridate")
#install.packages("ts")
#install.packages("XML")
#install.packages("rvest")
#install.packages("bcb")
#install.packages("MASS")
library(readxl)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(urca)
library(vars)
## Carregando pacotes exigidos: MASS
## Carregando pacotes exigidos: strucchange
## Carregando pacotes exigidos: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Carregando pacotes exigidos: sandwich
## Carregando pacotes exigidos: lmtest
library(lmtest)
library(tsDyn)
#
#rm(list = ls()) # código limpeza do ambiente
# Definir as datas inicial e hoje dinâmicamente sempre que o código for atualizado
rm(list=ls());
DataInicial <- Sys.Date() - (365*8) # calculada a data de anos atrás
DataHoje <- Sys.Date() # Recebe do sistema do dia atual
DataFiltro <- "2022-12-01" # necessário porque o PIB esta apenas disponível até o ano de 2022 com dados mensais completos.
Para utilização do pacote Sidrar é necessário passar como parâmetro uma tabela de referência.
Esta pode ser obtida pelo comando search_sidrar(c(Busca)). O método utilizado permite ao pesquisador produzir uma lista de tabelas para extrair os dados que comporão as séries históricas.
De modo análogo, __search_sidra c((“Busca”)) pode ser utilizado para o pacote sidrar.
library(ipeadatar)
#ls('package:ipeadatar')
# Procurando tabela com dados de inflação.
search_series(c('Taxa de juros - Over / Selic - acumulada no mês'))
## Warning: `filter_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `filter()` instead.
## ℹ See vignette('programming') for more help
## ℹ The deprecated feature was likely used in the ipeadatar package.
## Please report the issue at
## <https://github.com/gomesleduardo/ipeadatar/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## # A tibble: 1 × 7
## code name theme source freq lastupdate status
## <chr> <chr> <fct> <fct> <fct> <date> <fct>
## 1 BM12_TJOVER12 Taxa de juros - Over / Sel… Macr… Bacen… Mont… 2023-05-19 Active
# selecionando os dados do ipeadata
SELIC <- ipeadata("BM12_TJOVER12")
# aplicando o timeframe desejado
SELIC <- subset(SELIC, date >= DataInicial & date <= DataFiltro)
# Adicionar uma coluna com a marcação do trimestre
SELIC$Tri <- quarters(SELIC$date)
# Adicionar uma coluna com a marcação do ano
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
SELIC$Ano <- year(as.Date(SELIC$date))
# Utilizando a mediana das médias mensais, agregando e ordenando
SELIC <- aggregate(value ~ Ano + Tri, data = SELIC, FUN = mean)
SELIC <- SELIC[order(SELIC$Ano, SELIC$Tri), ]
data_SELIC_in <- as.Date("2015-04-01")
data_SELIC_fim <- as.Date("2022-10-01")
# Adicionando a sequência trimestral ao dataframe.
seq <- seq(data_SELIC_in, data_SELIC_fim, by = "3 months")
SELIC$date <- c(seq)
# Formantando SELIC para tamanho de dados desejados.
SELIC <- subset(SELIC, date >= "2016-07-01" & date <= DataFiltro)
head(SELIC, 1)
## Ano Tri value date
## 17 2016 Q3 1.146667 2016-07-01
tail(SELIC,1)
## Ano Tri value date
## 31 2022 Q4 1.053333 2022-10-01
summary(SELIC)
## Ano Tri value date
## Min. :2016 Length:26 Min. :0.1567 Min. :2016-07-01
## 1st Qu.:2018 Class :character 1st Qu.:0.4100 1st Qu.:2018-01-23
## Median :2019 Mode :character Median :0.5233 Median :2019-08-16
## Mean :2019 Mean :0.6027 Mean :2019-08-16
## 3rd Qu.:2021 3rd Qu.:0.8342 3rd Qu.:2021-03-09
## Max. :2022 Max. :1.1467 Max. :2022-10-01
nrow(SELIC)
## [1] 26
library(ipeadatar)
# Procurando tabela com dados de inflação.
# search_series(c('IPCA'))
# Dentre as possíveis tabelas, escolhi BM12_IPCA20N12
# selecionando os dados do ipeadata
IPCA <- ipeadata("BM12_IPCA20N12")
# aplicando o timeframe desejado
IPCA <- subset(IPCA, date >= DataInicial & date <= DataFiltro)
# Adicionar uma coluna com a marcação do trimestre
IPCA$Tri <- quarters(IPCA$date)
# Adicionar uma coluna com a marcação do ano
library(lubridate)
IPCA$Ano <- year(as.Date(IPCA$date))
# Utilizando a mediana das médias mensais, agregando e ordenando
IPCA <- aggregate(value ~ Ano + Tri, data = IPCA, FUN = mean)
IPCA <- IPCA[order(IPCA$Ano, IPCA$Tri), ]
data_IPCA_in <- as.Date("2015-04-01")
data_IPCA_fim <- as.Date("2022-10-01")
# Adicionando a sequência trimestral ao dataframe.
seq <- seq(data_IPCA_in, data_IPCA_fim, by = "3 months")
IPCA$date <- c(seq)
# Formantando IPCA para tamanho de dados desejados.
IPCA <- subset(IPCA, date >= "2016-07-01" & date <= DataFiltro)
head(IPCA, 1)
## Ano Tri value date
## 17 2016 Q3 0.33 2016-07-01
tail(IPCA, 1)
## Ano Tri value date
## 31 2022 Q4 0.4166667 2022-10-01
summary(IPCA)
## Ano Tri value date
## Min. :2016 Length:26 Min. :0.04667 Min. :2016-07-01
## 1st Qu.:2018 Class :character 1st Qu.:0.22167 1st Qu.:2018-01-23
## Median :2019 Mode :character Median :0.30000 Median :2019-08-16
## Mean :2019 Mean :0.35987 Mean :2019-08-16
## 3rd Qu.:2021 3rd Qu.:0.46000 3rd Qu.:2021-03-09
## Max. :2022 Max. :0.86000 Max. :2022-10-01
nrow(IPCA)
## [1] 26
# Baixando os dados com início em 2016
library(sidrar)
# Caminho da tabela definida em 'sidra.ibge.gov.br/tabela/6379'
cod_sidraDes <- '/t/6397/n1/all/v/4099/p/last%2024/c58/95253/d/v4099%201'
# Baixando as infos da tabela pela API
IBGE <- sidrar::get_sidra(api = cod_sidraDes)
## All others arguments are desconsidered when 'api' is informed
tail(IBGE,1)
## Nível Territorial (Código) Nível Territorial Unidade de Medida (Código)
## 25 1 Brasil 2
## Unidade de Medida Valor Brasil (Código) Brasil Variável (Código)
## 25 % 8.8 1 Brasil 4099
## Variável
## 25 Taxa de desocupação, na semana de referência, das pessoas de 14 anos ou mais de idade
## Trimestre (Código) Trimestre Grupo de idade (Código) Grupo de idade
## 25 202301 1º trimestre 2023 95253 Total
# convertendo a coluna Trimestre para formato data válido
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
#Obtendo os anos
ano <- substr(IBGE$Trimestre, nchar(IBGE$Trimestre) - 3, nchar(IBGE$Trimestre))
# Obtendo os trimestres
trimestre <- substr(IBGE$Trimestre, 1, 1)
# Substituindo strings 1, 2, 3 4, respectivamente, pelos meses 01, 04, 07 e 10
mes <- ifelse(trimestre == "1", "01",
ifelse(trimestre == "2", "04",
ifelse(trimestre == "3", "07",
ifelse(trimestre == "4", "10", NA))))
# Criando datas no padrão trimestral
dataDES <- paste(ano, mes, "01", sep = "-")
# Adicionando datas válidas ao dataframe
IBGE$Trimestre <- as.Date(dataDES)
tail(IBGE$Trimestre)
## [1] "2019-10-01" "2020-01-01" "2022-04-01" "2022-07-01" "2022-10-01"
## [6] "2023-01-01"
# Criando dataframe somente com info necessária
DES <- data.frame(Date = IBGE$Trimestre, Des = IBGE$Valor)
# Formatando série trimestral para mesmo tamanho de dados
DES <- subset(DES, Date >= '2016-07-01' & Date <= DataFiltro)
head(DES,1)
## Date Des
## 6 2016-07-01 11.9
tail(DES, 1)
## Date Des
## 23 2022-10-01 7.9
nrow(DES)
## [1] 18
# Aparantemente os daods trimestrais de 2020 e 2021 estão faltantes na base do IBGE.
# Assim foram manualmente substituidos pelos dados disponívels eim: https://agenciabrasil.ebc.com.br/economia/noticia/2021-03/desemprego-registrou-taxa-media-de-135-em-2020#:~:text=No%20%C3%BAltimo%20trimestre%20de%202020,em%20apenas%20em%20cinco%20estados.
#https://agenciabrasil.ebc.com.br/economia/noticia/2020-08/pnad-desemprego-vai-133-no-segundo-trimestre
# Preenchendo os dados faltantes de 2020
DES <- rbind(DES,c('2020-04-01', 13.3))
DES <- rbind(DES,c('2020-07-01', 14.6))
DES <- rbind(DES,c('2020-10-01', 13.9))
# Preenchendo os dados faltantes de 2021
DES <- rbind(DES,c('2021-01-01',14.9))
DES <- rbind(DES,c('2021-04-01',12.85))
DES <- rbind(DES,c('2021-07-01',12.6))
DES <- rbind(DES,c('2021-10-01', 11.1))
# Preenchendo primeiro trimestre de 2020
DES <- rbind(DES,c('2022-01-01', 11.1))
# Ficou estável ao tri do ano anterior.
# ordendando novamente a série temporal
DES <- DES[order(DES$Date), ]
nrow(DES)
## [1] 26
tail(DES,5)
## Date Des
## 25 2021-10-01 11.1
## 26 2022-01-01 11.1
## 21 2022-04-01 9.3
## 22 2022-07-01 8.7
## 23 2022-10-01 7.9
library(rvest)
# Local dos quais dados serão baixados
url <- "http://www.ipeadata.gov.br/ExibeSerie.aspx?serid=38415"
# Realizar o scraping dos dados
page <- read_html(url)
dados_html <- page %>%
html_nodes(xpath = '//table[@class="dxgvTable"]') %>%
html_table(fill = TRUE)
# Converter os dados em um dataframe
PIB <- data.frame(dados_html[[1]])
# Alterando nomes das colunas
PIB <- data.frame(Data = PIB$X1, PIB = PIB$X2)
# Limpando dataframe para somente valores numéricos
PIB <- PIB[4:nrow(PIB), ]
# Convertendo String do tipo XXXX - XT para dado trimestral do tipo data.
# Exemplo de dado na coluna
# Extrair o ano e o trimestre do dado
anos <- substring(PIB$Data, 1, 4)
trimestre <- substring(PIB$Data, 6, 8)
# Definir o mês baseado no trimestre
meses <- ifelse(trimestre == "T1", "01",
ifelse(trimestre == "T2", "04",
ifelse(trimestre == "T3", "07",
ifelse(trimestre == "T4", "10", NA))))
# Criar as datas no formato "yyyy-mm-dd"
datas <- paste(anos, meses, "01", sep = "-")
# Converter coluna de datas para o formato de data trimestral
PIB$Data <- as.Date(datas)
# Aplicando o filtro para os ultimos 5 anos
PIB <- subset(PIB, Data >= '2016-07-01' & Data <= DataFiltro)
# Por fim, convertendo os valores String para strings capazes de serem posteiormente interpretados como numéricos
PIB$PIB <- as.numeric(gsub("\\.", "", PIB$PIB))
tail(PIB, 5)
## Data PIB
## 107 2021-10-01 2309564
## 108 2022-01-01 2315709
## 109 2022-04-01 2471837
## 110 2022-07-01 2543645
## 111 2022-10-01 2584126
nrow(PIB)
## [1] 26
#adicionando as variáveis disponíveis para o dataframe dado
dados <- data.frame(Date = as.Date(DES$Date, format = "%Y-%m-%d"),
PIBmi = as.numeric(PIB$PIB),
IPCA = as.numeric(IPCA$value),
DES = as.numeric(DES$Des),
SELIC = as.numeric(SELIC$value))
tail(dados, 10)
## Date PIBmi IPCA DES SELIC
## 17 2020-07-01 1929703 0.2266667 14.60 0.1700000
## 18 2020-10-01 2048979 0.5366667 13.90 0.1566667
## 19 2021-01-01 2152622 0.4200000 14.90 0.1600000
## 20 2021-04-01 2182049 0.4733333 12.85 0.2633333
## 21 2021-07-01 2254492 0.6600000 12.60 0.4100000
## 22 2021-10-01 2309564 0.7333333 11.10 0.6166667
## 23 2022-01-01 2315709 0.7800000 11.10 0.8066667
## 24 2022-04-01 2471837 0.8600000 9.30 0.9600000
## 25 2022-07-01 2543645 0.4800000 8.70 1.0900000
## 26 2022-10-01 2584126 0.4166667 7.90 1.0533333
library(GetBCBData)
cod_poup <- "7828"
# tabela para a rentabilidade da poupança
rentpoup <- gbcbd_get_series(cod_poup)
##
## Fetching id = 7828 [7828] from BCB-SGS with cache
## Found 121 observations
# Filtrando rendimentos da poupança
rentpoup <- subset(rentpoup, ref.date >= '2016-07-01' & ref.date <= DataFiltro)
rentpoup <- data.frame(data = rentpoup$ref.date, valor = rentpoup$value)
# Extraindo o trimestre e o ano das datas
library(lubridate)
# Converter a coluna de data para o formato de data
rentpoup$data <- as.Date(rentpoup$data)
# Adicionar a coluna de trimestre
rentpoup$trimestre <- quarter(rentpoup$data)
# Adicionar a coluna de ano
rentpoup$ano <- year(rentpoup$data)
head(rentpoup, 5)
## data valor trimestre ano
## 1 2016-07-01 0.6629 3 2016
## 2 2016-08-01 0.7558 3 2016
## 3 2016-09-01 0.6583 3 2016
## 4 2016-10-01 0.6609 4 2016
## 5 2016-11-01 0.6435 4 2016
anos <- rentpoup$ano
trimestre <- rentpoup$trimestre
# Definir o mês baseado no trimestre
meses <- ifelse(trimestre == "1", "01",
ifelse(trimestre == "2", "04",
ifelse(trimestre == "3", "07",
ifelse(trimestre == "4", "10", NA))))
# Criar as datas no formato "yyyy-mm-dd"
datastri <- paste(anos, meses, "01", sep = "-")
# Converter coluna de datas para o formato de data trimestral
rentpoup$DatasTRI <- c(datastri)
# Agrupando pela média
rentpoup <- aggregate(valor ~ datastri, data = rentpoup, FUN = mean)
# juntando rentabilidade da poupança ao dataframe de dados
dados$SavI = c(rentpoup$valor)
tail(dados, 5)
## Date PIBmi IPCA DES SELIC SavI
## 22 2021-10-01 2309564 0.7333333 11.1 0.6166667 0.5163333
## 23 2022-01-01 2315709 0.7800000 11.1 0.8066667 0.5528000
## 24 2022-04-01 2471837 0.8600000 9.3 0.9600000 0.6240000
## 25 2022-07-01 2543645 0.4800000 8.7 1.0900000 0.6958000
## 26 2022-10-01 2584126 0.4166667 7.9 1.0533333 0.6699333
library(GetBCBData)
cod_prazo <- "27805" # Depósitos a prazo
#Meios de pagamento amplos - Depósitos a prazo (saldo em final de período) - Novo u.m.c. (mil)
dep_prazo <- gbcbd_get_series(cod_prazo)
##
## Fetching id = 27805 [27805] from BCB-SGS with cache
## Found 119 observations
head(dep_prazo, 10)
## ref.date value id.num series.name
## 1 2013-05-01 654729170 27805 id = 27805
## 2 2013-06-01 661146880 27805 id = 27805
## 3 2013-07-01 648627183 27805 id = 27805
## 4 2013-08-01 640591070 27805 id = 27805
## 5 2013-09-01 636467388 27805 id = 27805
## 6 2013-10-01 632614042 27805 id = 27805
## 7 2013-11-01 629290367 27805 id = 27805
## 8 2013-12-01 631550861 27805 id = 27805
## 9 2014-01-01 627715639 27805 id = 27805
## 10 2014-02-01 623862192 27805 id = 27805
cod_vista <- "27787"
# Meios de pagamento - Depósitos à vista (saldo em final de período) - Novo u.m.c. (mil)
dep_vista <- gbcbd_get_series(cod_vista)
##
## Fetching id = 27787 [27787] from BCB-SGS with cache
## Found 119 observations
head(dep_vista,10)
## ref.date value id.num series.name
## 1 2013-05-01 164917324 27787 id = 27787
## 2 2013-06-01 170053202 27787 id = 27787
## 3 2013-07-01 172438147 27787 id = 27787
## 4 2013-08-01 167912251 27787 id = 27787
## 5 2013-09-01 169761564 27787 id = 27787
## 6 2013-10-01 168983698 27787 id = 27787
## 7 2013-11-01 171308235 27787 id = 27787
## 8 2013-12-01 187720429 27787 id = 27787
## 9 2014-01-01 178135098 27787 id = 27787
## 10 2014-02-01 170502813 27787 id = 27787
cod_savings <- "1835"
# 1835 Meios de pagamento amplos - Depósito de Poupança (saldo em final de #período) u.m.c. (mil)
savings <- gbcbd_get_series(cod_savings)
##
## Fetching id = 1835 [1835] from BCB-SGS with cache
## Found 119 observations
head(savings,10)
## ref.date value id.num series.name
## 1 2013-05-01 527860458 1835 id = 1835
## 2 2013-06-01 539315239 1835 id = 1835
## 3 2013-07-01 551158753 1835 id = 1835
## 4 2013-08-01 558449217 1835 id = 1835
## 5 2013-09-01 567882037 1835 id = 1835
## 6 2013-10-01 575369015 1835 id = 1835
## 7 2013-11-01 584781093 1835 id = 1835
## 8 2013-12-01 599825987 1835 id = 1835
## 9 2014-01-01 604824604 1835 id = 1835
## 10 2014-02-01 609877474 1835 id = 1835
# aplicando o filtro desejado
savings <- subset(savings, ref.date >= '2016-07-01' & ref.date <= DataFiltro)
dep_vista <- subset(dep_vista, ref.date >= '2016-07-01' & ref.date <= DataFiltro)
dep_prazo <- subset(dep_prazo, ref.date >= '2016-07-01' & ref.date <= DataFiltro)
# verificando a quantidade de linhas nas trÊs tabelas filtradas
nrow(dep_prazo)
## [1] 78
nrow(savings)
## [1] 78
nrow(dep_vista)
## [1] 78
# Formatando para série trimestral
library(lubridate)
anos <- year(savings$ref.date)
trimestre <- quarters(savings$ref.date)
# Definir o mês baseado no trimestre
meses <- ifelse(trimestre == "Q1", "01",
ifelse(trimestre == "Q2", "04",
ifelse(trimestre == "Q3", "07",
ifelse(trimestre == "Q4", "10", NA))))
# Criar as datas no formato "yyyy-mm-dd"
datastri <- paste(anos, meses, "01", sep = "-")
# Converter coluna de datas para o formato de data trimestral
savings$DatasTRI <- c(datastri)
dep_prazo$DatasTRI <- c(datastri)
dep_vista$DatasTRI <- c(datastri)
savings <- aggregate(value ~ DatasTRI, data = savings, FUN = mean)
dep_prazo <- aggregate(value ~ DatasTRI, data = dep_prazo, FUN = mean)
dep_vista <- aggregate(value ~ DatasTRI, data = dep_vista, FUN = mean)
dados <- cbind(dados, Savings = savings$value, DepVista = dep_vista$value, DepPrazo = dep_prazo$value)
dados$DepTotais <- rowSums(data.frame(dados$DepPrazo, dados$DepVista, dados$Savings))
head(dados, 10)
## Date PIBmi IPCA DES SELIC SavI Savings DepVista
## 1 2016-07-01 1577172 0.3300000 11.9 1.1466667 0.6923333 644299807 148085053
## 2 2016-10-01 1632766 0.3566667 12.2 1.0700000 0.6634000 656256520 156864085
## 3 2017-01-01 1585585 0.3066667 13.9 1.0033333 0.6179667 663075626 153017312
## 4 2017-04-01 1630638 0.2133333 13.1 0.8433333 0.5435667 670974684 153672040
## 5 2017-07-01 1648631 0.1800000 12.5 0.7466667 0.5379333 690841486 155032100
## 6 2017-10-01 1720626 0.2966667 11.9 0.5833333 0.5000000 710506319 163379596
## 7 2018-01-01 1682460 0.1733333 13.2 0.5266667 0.5000000 729469370 163855257
## 8 2018-04-01 1734454 0.2133333 12.6 0.5200000 0.5000000 745119657 168568061
## 9 2018-07-01 1767868 0.2433333 12.0 0.5266667 0.5000000 768784127 173730991
## 10 2018-10-01 1819359 0.2266667 11.7 0.5066667 0.5000000 787989308 177216453
## DepPrazo DepTotais
## 1 659675119 1452059979
## 2 684942790 1498063395
## 3 711482590 1527575528
## 4 789986433 1614633157
## 5 826600944 1672474529
## 6 838968408 1712854323
## 7 843486825 1736811452
## 8 895058350 1808746067
## 9 942237342 1884752460
## 10 988249872 1953455632
par(mfrow=c(2,3))
plot(dados$Date, dados$PIBmi, col = 'black', type = 'l',
ylim = c(min(dados$PIBmi), max(dados$PIBmi)),
main = "PIB", xlab = "Data", ylab = "PIB")
plot(dados$Date, dados$SELIC, col = 'black', type = 'l',
ylim = c(min(dados$SELIC), max(dados$SELIC)),
main = "Taxa de juros - Selic", xlab = "Data", ylab = "SELIC")
plot(dados$Date, dados$DES, col = 'black', type = 'l',
ylim = c(min(dados$DES), max(dados$DES)),
main = "Desemprego", xlab = "Data", ylab = "Desemprego")
plot(dados$Date, dados$IPCA, col = 'black', type = 'l',
ylim = c(min(dados$IPCA), max(dados$IPCA)),
main = "IPCA", xlab = "Data", ylab = "IPCA")
plot(dados$Date, dados$DepTotais, col = 'black', type = 'l',
ylim = c(min(dados$DepTotais), max(dados$DepTotais)),
main = "Depósitos totais", xlab = "Data", ylab = "Depósitos totais")
# Criando a matriz de correlação para as taxas
matriz_taxas <- data.frame(dados$PIB, dados$DepTotais, dados$SELIC, dados$IPCA, dados$DES)
matriz_correl <- cor(matriz_taxas)
# exibindo a relação gráfica
matriz_correl
## dados.PIB dados.DepTotais dados.SELIC dados.IPCA dados.DES
## dados.PIB 1.0000000 0.9462268 0.0680610 0.7440503 -0.6004160
## dados.DepTotais 0.9462268 1.0000000 -0.1476187 0.6890658 -0.3694546
## dados.SELIC 0.0680610 -0.1476187 1.0000000 0.2234004 -0.6140568
## dados.IPCA 0.7440503 0.6890658 0.2234004 1.0000000 -0.3934822
## dados.DES -0.6004160 -0.3694546 -0.6140568 -0.3934822 1.0000000
# PIB e depósitos totais é 0,94, PIB e IPCA é 0,74. Correlação desemprego e PIB é negativa.
# Inflação e depósitos totais é positiva (expansão monetária)
# Taxa de juros e de semprego é negativa - 0,61
# Inflação e desemprego é negativa, embora seja mais fraca
# Visualização gráfica das correlações
pairs(matriz_correl)
## Preparando séire para modelos de VEC/VAR
# Conversão para log natural para diminuir o impacto de depósitos totais e PIB serem séires nominais, enquanto as demais serem taxas. O Ideal seria ter todos os dados nominais ou converter PIB e Des para um índice, mas isto funcionaria como uma transformação. Espera-se que o log suavize esta questão.
# Excluir a coluna de data
logdados <- subset(dados, Date >= "2017-01-01")
logdados <- subset(dados, select = -1)
head(logdados, 1)
## PIBmi IPCA DES SELIC SavI Savings DepVista DepPrazo DepTotais
## 1 1577172 0.33 11.9 1.146667 0.6923333 644299807 148085053 659675119 1452059979
# converter para log
logdados <- log(logdados)
sum(is.na(logdados))
## [1] 0
dadosts<-ts(as.data.frame(cbind(logdados$DepTotais,logdados$PIB,logdados$DES, logdados$IPCA,logdados$SELIC),start=c(1997,1), frequency = 4))
head(dadosts, 5) # visualizar os 5 primeiros dados da série temporal
## Time Series:
## Start = 1
## End = 5
## Frequency = 1
## V1 V2 V3 V4 V5
## 1 21.09625 14.27114 2.476538 -1.108663 0.13685918
## 2 21.12744 14.30579 2.501436 -1.030954 0.06765865
## 3 21.14695 14.27646 2.631889 -1.181994 0.00332779
## 4 21.20237 14.30448 2.572612 -1.544899 -0.17039299
## 5 21.23757 14.31546 2.525729 -1.714798 -0.29213642
# Regras de decisão para os testes ADF:
# Se Tau3 calculado < q ue tau3 crítico com α = 5% , unicaudal à esquerda, não há raiz unitária, caso contrário testar M3 X M2
# Se phi3 calculado maior que phi3 crítico com α = 5%, M3 melhor que M2, caso contrário M2 melhor que M3.
# Se phi2 calculado maior que phi2 crítico com α = 5%, modelo M3 é melhor que M1.
# Se M3 melhor que (M2 & M1) -> M3 vence!
# Se M2 melhor que M2 melhor que M3 & M3 melhor que M1 -> M2 vence! Fazer ADF2 com trend para confirmar a relação de M1 e M2.
# Se tau2 calculado < tau2 crítico com α = 5%, tem raiz unitária
# Se phi1 calculado > phi1 crítico com α = 5%, M2 melhor que M1. Caso contrário, M1 melhor que M2.
library(urca)
# Tem ou não RU & M3 X M2 & M3 X M1, para M3 com tendência, M2 com drift e# # M1 = None
DepTotais_adf3 <-ur.df(logdados$DepTotais, type = "trend", selectlags = c("BIC"))
summary(DepTotais_adf3)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.040130 -0.019747 -0.004929 0.012561 0.107438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.35943 2.58445 2.461 0.0231 *
## z.lag.1 -0.30180 0.12312 -2.451 0.0236 *
## tt 0.01195 0.00500 2.389 0.0268 *
## z.diff.lag 0.50597 0.19279 2.624 0.0162 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03176 on 20 degrees of freedom
## Multiple R-squared: 0.3302, Adjusted R-squared: 0.2298
## F-statistic: 3.287 on 3 and 20 DF, p-value: 0.04194
##
##
## Value of test-statistic is: -2.4513 4.0649 3.0092
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.15 -3.50 -3.18
## phi2 7.02 5.13 4.31
## phi3 9.31 6.73 5.61
# tau3 -2.4513 > -3.50 -> Tem RU
# phi3 3.0092 < 6.73 -> vence M2, modelo com trend
# phi2 4.0649 < 5.13 -> Vence M1, modelo com none
# Conclucão: Testar M2 contra M1 por ADF2
# Tem ou não RU & M2 X M1 com drift para M2 = Drift e M1 = None
DepTotais_adf2<-ur.df(logdados$DepTotais, type = "drift", selectlags = c("BIC"))
summary(DepTotais_adf2)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression drift
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.041455 -0.016962 -0.002084 0.006697 0.135675
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29903 0.54945 0.544 0.5920
## z.lag.1 -0.01282 0.02552 -0.502 0.6207
## z.diff.lag 0.36777 0.20349 1.807 0.0851 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03514 on 21 degrees of freedom
## Multiple R-squared: 0.139, Adjusted R-squared: 0.05702
## F-statistic: 1.695 on 2 and 21 DF, p-value: 0.2077
##
##
## Value of test-statistic is: -0.5022 2.6486
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau2 -3.58 -2.93 -2.60
## phi1 7.06 4.86 3.94
# tau2 -0.5022 > -2.93 -> Tem RU
# phi1 2.64 < 4.86 -> M1 vence. Modelo none!
# Conclusão: vence modelo com none
# Tem ou não RU!
DepTotais_adf1 <-ur.df(logdados$DepTotais, type = "none", selectlags = c("BIC"))
summary(DepTotais_adf1)
##
## ###############################################
## # 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
## -0.045240 -0.012672 -0.001222 0.004993 0.136225
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 0.0010707 0.0004711 2.273 0.0331 *
## z.diff.lag 0.3600333 0.1997219 1.803 0.0852 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03457 on 22 degrees of freedom
## Multiple R-squared: 0.5736, Adjusted R-squared: 0.5348
## F-statistic: 14.8 on 2 and 22 DF, p-value: 8.471e-05
##
##
## Value of test-statistic is: 2.273
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.62 -1.95 -1.61
# conclusão: Tau1 2.273 > -1.95; Tem raiz unitária. None
# Repetir ADF3 para variávei na primeira diferença.
Diff_DepTotais_adf3<-ur.df(diff(logdados$DepTotais), type = "trend", selectlags = c("BIC"))
summary(Diff_DepTotais_adf3)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.037467 -0.013380 -0.002297 0.006319 0.135662
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0307432 0.0190632 1.613 0.12330
## z.lag.1 -0.7359481 0.2558418 -2.877 0.00966 **
## tt -0.0002695 0.0011491 -0.235 0.81710
## z.diff.lag 0.1456260 0.2261724 0.644 0.52736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03655 on 19 degrees of freedom
## Multiple R-squared: 0.3397, Adjusted R-squared: 0.2354
## F-statistic: 3.258 on 3 and 19 DF, p-value: 0.04436
##
##
## Value of test-statistic is: -2.8766 2.7867 4.1791
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.38 -3.60 -3.24
## phi2 8.21 5.68 4.67
## phi3 10.61 7.24 5.91
# Continua tendo raiz unitária.
# tau3: 2.8766 > -3.60 -> Tem RU
# phi3: 4.1791 < 7.24 -> M2 melhor M3
# phi2: 2.7867 < 5.68 -> M1 melhor M3
# Conclusão: Testar segunda diferença e reavalir modelos. Como sabemos que None é o melhor, verifica-se a presença de raiz unitária pelo ADF, pior cenário.
Diff2_DepTotais_adf3<-ur.df(diff(logdados$DepTotais,2), type = "none", selectlags = c("BIC"))
summary(Diff2_DepTotais_adf3)
##
## ###############################################
## # 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
## -0.08213 -0.01555 0.01716 0.03518 0.14732
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z.lag.1 -0.2116 0.1117 -1.895 0.0727 .
## z.diff.lag 0.3935 0.2055 1.914 0.0700 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04806 on 20 degrees of freedom
## Multiple R-squared: 0.2198, Adjusted R-squared: 0.1417
## F-statistic: 2.817 on 2 and 20 DF, p-value: 0.0836
##
##
## Value of test-statistic is: -1.8947
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau1 -2.66 -1.95 -1.6
# tau1 -1.8947 < -1.95. Nova diferenciação necessária
Diff3_DepTotais_adf3<-ur.df(diff(logdados$DepTotais,3), type = "trend", selectlags = c("AIC"))
summary(Diff3_DepTotais_adf3)
##
## ###############################################
## # Augmented Dickey-Fuller Test Unit Root Test #
## ###############################################
##
## Test regression trend
##
##
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.097909 -0.018331 0.001742 0.013973 0.131515
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0440990 0.0290548 1.518 0.1474
## z.lag.1 -0.4169144 0.1476159 -2.824 0.0117 *
## tt 0.0002188 0.0018052 0.121 0.9049
## z.diff.lag 0.5619940 0.2042036 2.752 0.0136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04998 on 17 degrees of freedom
## Multiple R-squared: 0.4004, Adjusted R-squared: 0.2945
## F-statistic: 3.783 on 3 and 17 DF, p-value: 0.03019
##
##
## Value of test-statistic is: -2.8243 2.6606 3.9906
##
## Critical values for test statistics:
## 1pct 5pct 10pct
## tau3 -4.38 -3.60 -3.24
## phi2 8.21 5.68 4.67
## phi3 10.61 7.24 5.91
# A partir daqui o trabalho começou a não fazer sentido. Vou tentar uma transformação Box-Cox para resolver a estacionariedade.