library(seastests)
library(urca)
library(randtests)
library(TSA)
library(fpp3)
library(tibble)
library(readr)
library(reshape2)  
library(ggplot2)
library(forecast)
library(fpp2)
library(seasonal)
library(dplyr)

Banco de dados

Nesse projeto vamos utilizar o banco de dados da Divida mobiliária que está disponivel no site: https://dadosabertos.bcb.gov.br/dataset/4176-divida-mobiliaria---participacao-por-indexador---posicao-em-carteira---igp-di

Conceito: A dívida mobiliária federal é constituída dos títulos emitidos pelo Tesouro Nacional e pelo Banco Central do Brasil e contempla os valores que estão em poder do mercado ou em carteira no Banco Central.

Possui 268 informações em dados mensais sobre o valor da dívida na data referente

dados <- read.csv2("/cloud/project/bcdata.sgs.4176.csv")
head(dados)
##         data valor
## 1 01/01/2000  2.04
## 2 01/02/2000  2.05
## 3 01/03/2000  1.96
## 4 01/04/2000  4.47
## 5 01/05/2000  4.43
## 6 01/06/2000  4.39

Indicando ao R a variável data e a formatação dela

dados$data<- as.Date(dados$data,"%d/%m/%Y")

Análise exploratória dos dados

str(dados)
## 'data.frame':    268 obs. of  2 variables:
##  $ data : Date, format: "2000-01-01" "2000-02-01" ...
##  $ valor: num  2.04 2.05 1.96 4.47 4.43 4.39 4.49 4.51 4.39 4.38 ...
summary(dados)
##       data                valor       
##  Min.   :2000-01-01   Min.   :0.0300  
##  1st Qu.:2005-07-24   1st Qu.:0.0500  
##  Median :2011-02-15   Median :0.3000  
##  Mean   :2011-02-15   Mean   :0.9173  
##  3rd Qu.:2016-09-08   3rd Qu.:1.2975  
##  Max.   :2022-04-01   Max.   :4.5100

Análise de dados faltantes

colSums(is.na(dados))
##  data valor 
##     0     0

Indicando ao R que nosso banco de dados é uma série temporal

dadost<- ts(dados,start=c(2000,1,1),frequency=12,end =c(2022,4,1))

Séries Temporais

Caracterização, Modelagem e Previsão de uma série temporal

Uma série temporal (ST) é um conjunto de observações, geralmente equiespaçadas, obtidas a partir da medição/observação de uma variável ao longo do tempo;

A frequência de observação/medição de uma ST pode variar dependendo do fenômeno observado: minuto, diária, semanal, mensal, anual etc;

A partir da análise de ST, é possível obter subsídios para a escolha de um modelo adequado para modelar a série, escolhido dentro de uma classe de modelos pré-existentes;

Uma vez construído, um modelo de ST pode ser utilizado para efetuar previsões probabilísticas sobre o futuro da série e/ou analisar a sensibilidade de algumas variáveis à variável objetivo;

A capacidade de realizar previsões é fundamental no processo de tomada de decisões em diversos contextos e em diversos lugares como orgãos públicos e empresas.

dados %>% 
  ggplot(aes(x = data, y = valor)) +
  geom_line() +
  # adicionar curva de tendencia
  geom_smooth(se = FALSE) +
  theme_bw() +
  # quebrar eixo x em 1 mes
  scale_x_date(date_labels = "%Y/%m",
               minor_breaks = NULL) +
  # inverter eixos
  theme(axis.text.x = element_text(angle = 90))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Analisando esse gráfico inicial podemos supor algumas características da série temporal: Tendência: Está muito visivel a tendência de decrescimento dos dados, mais ou menos em 2004 começou a “cair” e a série não voltou mais a “subir” Sazonalidade: Não é possivel identificar na série temporal padrões de comportamento, ou seja oscilções de subida e descida em períodos especificos, ainda mais após a queda de 2004 Estacionaridade:Uma série temporal estacionária é aquela cujas propriedades estatísticas, como a média, a variância e a auto correlação, são constantes ao longo do tempo, mas na série estacionária não possui sazonalidade ou tendência, como podemos observar claramente a série tmeporal para esses dados possuim uma tendência de descrescimo por isso não é estacionária

Correlação

Correlação é um indicador estatístico que mede o nível de dependência linear entre duas variáveis. Está definida no intervalo [−1,+1]. Se a correlação é negativa, indica que as variáveis são inversamente proporcinais: quando uma aumenta, a outra diminui. Se é positiva, indica que as variáveis são diretamente proporcionais. Mas para esse banco de dados não precisamos calcular a correlação já que temos apenas a variável valor relacionada a data.

Autocorrelação

Uma correlação nos diz o quão relacionado é um valor com outro valor e uma autocorrelação (ACF) é quando comparamos o valor do presente com valores do passado da mesma série. A diferença entre a autocorrelação e a autocorrelação parcial (PACF) é quase um detalhe: em uma ACF temos a correlação direta e indireta e em uma PACF apenas a correlação direta

##Funcao de autocorrelacao e autocorrelacao parcial
par(mfrow=c(1,2))
acf(dados$valor, main = "FAC")
pacf(dados$valor, main = "FACP")

Para as séries temporais esperamos que cada autocorrelação seja próxima de zero, eles não serão exatamente iguais a zero, pois há alguma variação aleatória. Podemos observar claramente a tendência da série pelo gráfico ACF. A linha tracejada azul no gráfico ACF parcial indica onde é significativamente diferente de zero, praticamente todos os valores ACF estão dentro do limite da linha tracejada azul, ou seja, autocorrelação igual a zero, indicando que a série é aleatória, conforme o esperado.

Decomposição de série temporal

a função decompose() a qual realiza a decomposição sazonal clássica por médias móveis lidando com componentes sazonais aditivos ou multiplicativos.Neste gráfico pode-se verificar os dados observados, a tendência, a sazonalidade e a componente aleatória extraídos

Decomposição aditiva

decomp <- decompose(dadost[,2],type = "additive")
plot(decomp)

Pelo gráfico trend podemos reafirmr que há uma tendência de decrescimo a longo prazo nos dados. Já o componente sazonal reflete a influência de um determinado fator externo, que ocorre sempre no mesmo período mas para esses dados não temos sazonalidade. Pelo gráfico random podemos observar que os valores estão próximos de 0.

Decomposição multiplicativa

decomp <- decompose(dadost[,2],type = "multiplicative")
plot(decomp)

Tendência

A tendência de uma série temporal é definida como um padrão de crescimento/descrecimento da variável em um certo período de tempo.

Existem testes específicos para a identificação da tendência, como o Teste de Wald e o de Cox-Stuart.

Teste de tendência de Cox-stuart

\(H_0\): Não possui tendência \(H_1\): Possui tendência

cox.stuart.test(dadost) 
## 
##  Cox Stuart test
## 
## data:  dadost
## statistic = 0, n = 268, p-value < 2.2e-16
## alternative hypothesis: non randomness

Pelo teste de Cox Stuart foi obtido o p-valor < 0,05, logo ao nível de 5% de significância temos evidências que a série possui tendencia.

Uma técnica muito utilizada é o ajuste de uma Regressão Linear Simples para a identificação da inclinação da reta de tendência. Vale lembrar que o ajuste da Regressão Linear, neste caso, pode levar a resultados enviesados e, para evitar este problema, estimadores robustos à autocorrelação podem ser utilizados Para remover a tendência pode-se utilizar uma regressão linear. Usando a regressão linear para modelar os dados da série temporal com índices lineares. Os resíduos do modelo resultante são uma representação da série temporal desprovida de tendência.

Ajuste de tendência

# Decompor a série 
otdecomp<- stl(dadost[,2], s.window="periodic")

# pegar o componentes (sazonalidade)  
ajuste.sazonal <- dadost[,2]-otdecomp$time.series[,"seasonal"]  
ajuste.tendencia <- lm(ajuste.sazonal~c(1:length(ajuste.sazonal))) 
plot(resid(ajuste.tendencia), type="l")

Teste de estacionariedade

\(H_0\): Não é estacionária \(H_1\): É estacionária

ur.kpss(dadost)
## 
## ####################################### 
## # KPSS Unit Root / Cointegration Test # 
## ####################################### 
## 
## The value of the test statistic is: 5.4309

O valor do teste foi maior que 0,05 portanto há indicios que esta série não é estacionária Mas avaliando o gráfico da série novamente

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Podemos questionar se após a queda brusca, a partir de 2015 a série se tornaria estacionária, para testar vamos realizar um recorte da série

est<- ts(dados,start=c(2015,1,1),frequency=12,end =c(2022,4,1))
ur.kpss(est)
## 
## ####################################### 
## # KPSS Unit Root / Cointegration Test # 
## ####################################### 
## 
## The value of the test statistic is: 2.8424

Após realizar o teste para após o período de 2015, apesar do gráfico, podemos afirmar que a série é não estacionária mesmo com o recorte temporal

Teste de sazonalidade

Caso tenhamos dúvidas sobre a presença de sazonalidade após a inspeção gráfica, podemos realizar um teste de sazonalidade com o pacote seastests

isSeasonal(dadost[,2], freq = 12)
## [1] FALSE

Ajustando um modelo

modelo.ets <- ets(dadost[,2])
summary(modelo.ets)
## ETS(M,Ad,N) 
## 
## Call:
##  ets(y = dadost[, 2]) 
## 
##   Smoothing parameters:
##     alpha = 0.5189 
##     beta  = 0.1428 
##     phi   = 0.8662 
## 
##   Initial states:
##     l = 0.6399 
##     b = 1.4498 
## 
##   sigma:  0.0692
## 
##       AIC      AICc       BIC 
## -620.0601 -619.7383 -598.5142 
## 
## Training set error measures:
##                       ME      RMSE        MAE       MPE     MAPE      MASE
## Training set -0.02583864 0.1610694 0.04836649 -1.808264 3.855734 0.2005153
##                   ACF1
## Training set 0.1367532

um modelo ETS(M, Ad, N). A primeira letra se refere ao componente de erro M (multiplicativo), a segunda ao componente de tendência é Ad (aditivo amortecido), e a terceira ao componente de sazonalidade N ( não possui sazonalidade). O ajuste do modelo, graficamente, é:

plot(dadost[,2])
lines(fitted(modelo.ets), col = "red")

Forecast

Previsão pelo modelo forecast ets para os próximos 3 anos

fit <- ets(dadost[,2])
fit %>% forecast(h=36, level= 95) %>%
  autoplot() +
  ylab ("Valor")

##### Acuracia do modelo

accuracy(fit)
##                       ME      RMSE        MAE       MPE     MAPE      MASE
## Training set -0.02583864 0.1610694 0.04836649 -1.808264 3.855734 0.2005153
##                   ACF1
## Training set 0.1367532

As medidas calculadas são:

ME: Erro médio

RMSE: Erro quadrático médio da raiz

MAE: Erro Absoluto Médio

MPE: Erro de Porcentagem Média

MAPE: Erro de Porcentagem Absoluta Média

MASE: Erro Escalado Absoluto Médio

ACF1: Autocorrelação de erros na defasagem 1.