Explicando o comportamento de depósitos totais (Depósito à vista, Depósito a Prazo e Poupança) do Brasil nos ultimos cinco anos uitlizando o PIB, Selic, Inflação e Desemprego

Origem dos dados utilizados

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

Definindo o timeframe da série de dados

# 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. 

Consultando tabelas de variáveis do pacote Sidrar

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.

Obtendo e tratando os dados da SELIC

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

Baixando os Dados da Inflação

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 dados para o nível de Desemprego

# 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

Baixando dados do PIB direto da tabela do IPEAData via Webscrapping

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

Juntando os dados em um único dataframe

#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

Extraindo rendimentos da poupança

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

Baixando dados para depósitos totais

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

Inspeção visual dos dados!

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")

Inspecionando a correlação entre os dados

# 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

Avaliando estacionariedade dos dados pelo teste ADF (Augmented Dickey-Fuller)

# 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.

Para Depósitos Totais

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.