Análise Exploratória de Séries

PMC - Pesquisa Mensal de Comércio

A Pesquisa Mensal de Comércio disponibiliza indicadores que permitem acompanhar o comportamento do comércio varejista no País.

Fonte dos dados: https://www.ibge.gov.br/estatisticas/economicas/comercio/9227-pesquisa-mensal-de-comercio.html?=&t=o-que-e

# Bibliotecas utilizadas

library(e1071)      # estatísticas 
library(forecast)   # séries temporais

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(seasonal)   # determinar a sazonalidade

## Warning: package 'seasonal' was built under R version 4.1.1

library(readxl)     # para ler tabelas em xls
library(dplyr)      # para transformação de dados  

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyverse)  # cap 5. pre processamento de dados

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.4     v purrr   0.3.4
## v tibble  3.1.2     v stringr 1.4.0
## v tidyr   1.1.3     v forcats 0.5.1
## v readr   1.4.0

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x tibble::view()  masks seasonal::view()

library(tsibble)    # cap 5. pre processamento de dados

## Warning: package 'tsibble' was built under R version 4.1.1

## 
## Attaching package: 'tsibble'

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union

library(fpp2)

## Warning: package 'fpp2' was built under R version 4.1.1

## -- Attaching packages ---------------------------------------------- fpp2 2.4 --

## v fma       2.4     v expsmooth 2.3

## Warning: package 'fma' was built under R version 4.1.1

## Warning: package 'expsmooth' was built under R version 4.1.1

## 

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following object is masked from 'package:tsibble':
## 
##     interval

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(lmtest)

## Carregando pacotes exigidos: zoo

## 
## Attaching package: 'zoo'

## The following object is masked from 'package:tsibble':
## 
##     index

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

# ================================================================================


# Dados

# PMC - Pesquisa Mensal de Comércio - IBGE


dados <- read_xlsx("PMC.xlsx") 

# Transformando em série


# volume de vendas
volume <- ts(dados$volume,start = 2003, frequency = 12)
volume

##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2003  44.3  42.7  42.8  43.7  45.6  43.0  46.2  46.2  46.3  49.3  48.8  61.8
## 2004  47.2  44.8  50.9  49.0  51.4  50.3  52.3  51.8  50.8  52.9  53.3  68.4
## 2005  50.7  45.2  53.0  50.6  52.2  51.9  52.7  54.7  51.8  53.8  55.2  70.7
## 2006  52.7  46.9  55.1  52.0  56.8  52.8  56.1  59.0  57.1  59.6  60.8  74.8
## 2007  58.1  52.5  62.2  59.8  64.7  61.9  63.7  67.9  63.9  69.4  69.5  82.8
## 2008  66.4  62.2  69.7  69.2  72.3  70.6  74.0  72.6  74.0  71.8  66.6  83.8
## 2009  68.3  63.1  74.0  68.7  74.3  78.0  74.7  76.7  80.9  79.9  77.4  95.4
## 2010  75.3  71.8  90.5  77.0  81.4  80.5  84.0  87.4  85.7  88.7  90.6 109.8
## 2011  83.9  82.3  88.0  86.3  91.9  88.0  90.1  92.1  89.7  90.0  93.4 114.5
## 2012  90.8  84.9  97.1  88.8  96.5  99.0  99.3 106.5  91.5 103.1 100.1 120.3
## 2013  97.1  85.9 100.2  96.9 100.7  97.0 103.0 105.6  98.6 105.4 106.1 123.7
## 2014 101.7  92.9  94.6  96.9 101.4  91.2  97.9  98.5  97.4 102.9 103.6 121.0
## 2015  96.7  83.3  93.9  88.9  90.9  87.9  91.1  89.0  86.2  90.7  89.9 107.8
## 2016  83.0  78.7  86.4  80.7  81.6  80.9  81.4  82.2  78.9  81.6  85.2 100.5
## 2017  83.0  74.9  84.8  80.2  85.6  84.4  86.0  88.5  86.1  87.8  92.6 107.4
## 2018  88.3  78.9  92.3  87.2  87.4  87.6  88.5  94.6  88.1  93.3  98.1 109.3
## 2019  91.4  85.0  89.2  90.0  93.5  89.3  95.3  96.0  92.0  98.5 101.9 113.8
## 2020  94.5  87.5  83.5  65.3  78.6  87.3  96.9  99.6  98.8 104.5 106.1 117.0
## 2021  91.6  85.8  91.7  92.0  99.2  97.3

# Plotando a série
plot(volume)

# Características gerais 

summary(volume)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   42.70   65.58   84.95   80.21   93.38  123.70

var(volume)

## [1] 359.3059

sd(volume)

## [1] 18.95537

kurtosis(volume)

## [1] -0.7775273

skewness(volume)

## [1] -0.3195193

hist(volume, main="Histograma do Volume ",xlab="Volume de Vendas",ylab="Frequência",col=rgb(0.2,0.4,0.6,0.4))

#plot(volume)
boxplot(volume,main="Boxplot do Volume ",col=rgb(0.2,0.4,0.6,0.4))

# Organizando meses

jan <- subset(dados, mes == 1)
jan <- ts(jan, start = 2003)
voljan <- jan[,3]



fev <- subset(dados, mes == 2)
fev <- ts(fev, start = 2003)
volfev <- jan[,3]

mar <- subset(dados, mes == 3)
mar <- ts(mar, start = 2003)
volmar <- mar[,3]

abr <- subset(dados, mes == 4)
abr <- ts(abr, start = 2003)
volabr <- abr[,3]

mai <- subset(dados, mes == 5)
mai <- ts(mai, start = 2003)
volmai <- mai[,3]

jun <- subset(dados, mes == 6)
jun <- ts(jun, start = 2003)
voljun <- jun[,3]

jul <- subset(dados, mes == 7)
jul <- ts(jul, start = 2003)
voljul <- jul[,3]

ago <- subset(dados, mes == 8)
ago <- ts(ago, start = 2003)
volago <- ago[,3]


set <- subset(dados, mes == 9)
set <- ts(set, start = 2003)
volset <- set[,3]


out <- subset(dados, mes == 10)
out <- ts(out, start = 2003)
volout <- out[,3]

nov <- subset(dados, mes == 11)
nov <- ts(nov, start = 2003)
volnov <- nov[,3]

dez <- subset(dados, mes == 12)
dez <- ts(dez, start = 2003)
voldez <- dez[,3]

# Plotando variação de volume em meses

plot(voljan, type = "b", frame = FALSE, pch = 8,ylim = c(0,120), xlim = c(2001,2020),
     col = "red", xlab = "tempo", ylab = "volume",lty = 22)

# Add proxiimo mes
lines(volfev, pch = 18, col = "blue", type = "b", lty = 2)
# Add proxiimo mes
lines(volmar, pch = 18, col = "green", type = "b", lty = 2)
# Add proxiimo mes
lines(volabr, pch = 18, col=c("#B0E2FF"), type = "b", lty = 2)
# Add proxiimo mes
lines(volmai, pch = 18, col=c("#FFFF00"), type = "b", lty = 2)
# Add proxiimo mes
lines(voljun, pch = 18, col=c("#FF1493"), type = "b", lty = 2)
# Add proxiimo mes
lines(voljul, pch = 18, col=c("#836FFF"), type = "b", lty = 2)
# Add proxiimo mes
lines(volago, pch = 18, col=c("#FFC1C1"), type = "b", lty = 2)
# Add proxiimo mes
lines(volset, pch = 18, col=c("#6B8E23"), type = "b", lty = 2)
# Add proxiimo mes
lines(volout, pch = 18, col=c("#FFD700"), type = "b", lty = 2)
# Add proxiimo mes
lines(volnov, pch = 18, col=c("#F0FFF0"), type = "b", lty = 2)
# Add proxiimo mes
lines(voldez, pch = 18, col=c("#2F4F4F"), type = "b", lty = 2)
# Add proxiimo mes

legend("bottomright", legend=c("Janeiro", "Fevereiro","Março","Abril","Maio", "Junho", "Julho","Agosto","Setembro", "Outubro","Novembro","Dezembro"),
       col=c("red", "blue","green","#B0E2FF","#FFFF00","#FF1493","#836FFF","#FFC1C1","#6B8E23","#FFD700","#F0FFF0","#2F4F4F"), lty = 3:2, cex=0.6, bty="n")

#volume
#ano <- volume[,1]

#library(ggplot2)

#ggplot(volume, aes(x=ano, y=value, col=variable)) + geom_line()

#p1 <- ggplot(volume) +
 # geom_line(aes(x = fev, y = jan, colour = "Valor máximo")) 

#p1

#plot(volume$Jan)

# Sazonalidade

seasonplot(volume,col = rainbow(18), year.labels =TRUE, type = "o",pch = 16)

Decomposicção da série | Cap. 3

# Decomposição da série

# Método 11 - Baseia-se na decomposição clássica, mas inclui muitos passos e características extras para superar as desvantagens da decomposição clássica que foram discutidas na seção anterior.

fit <- seas(volume, x11 = "")
# varejoms %>% seas(x11='') -> fit # uso a serie 'varejoms', aplico o 'seas' x11
# e gero 'fit'
library(fpp2)


#decomposição multiplicativa
decomposicao.mult <- decompose(volume, type = "multiplicative")  # multiplicativa
plot(decomposicao.mult)

#decomposição aditiva
decomposeAditiva <- decompose(   x = volume,  type = "additive")
plot(  decomposeAditiva)

# Sazonalidade Ajustada

autoplot(volume, series = "Data") + autolayer(trendcycle(fit), series = "Trend") + 
  autolayer(seasadj(fit), series = "Seasonally Adjusted") + xlab("Ano") + ylab("Index") + 
  ggtitle("Índice de volume de vendas no varejo Total Brasil") + 
  scale_colour_manual(values = c("gray", "blue", "red"), breaks = c("Data", "Seasonally Adjusted", 
                                                                    "Trend"))

#volume %>% diff() %>% ggtsdisplay(main = "Série varejo do Brasil em primeira diferença")

# Checagem de Resíduos
checkresiduals(fit)

## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

Tentativa de Previsão | Cap. 5

# Método Média

media <- meanf(volume, h = 5)
plot(meanf(volume, h = 5), main="Método Média")

# Método Ingênuo 

ingenuo <- naive(volume, h = 5)
plot(naive(volume, h = 5),  main="Método Ingênuo")

# Método Ingênuo Sazonal | Cap 5
ingenuo_saz <- snaive(volume, h = 5)
plot(snaive(volume, h = 5),  main="Método Ingênuo Sazonal")

# Método Drift
drift_m <- rwf(volume, h = 5)
plot(rwf(volume, h = 8, drift = T),  main="Método Drift")

# Tentativa de comparar os 3 considerando h=10

autoplot(volume) +
  autolayer(meanf(volume, h=10),
    series="Mean", PI=FALSE) +
  autolayer(naive(volume, h= 10),
    series="Naïve", PI=FALSE) +
  autolayer(snaive(volume, h=10),
    series="Seasonal naïve", PI=FALSE) +
  ggtitle("Tentativa de previsão com os 3 métodos") +
  xlab("Year") + ylab("volume") +
  guides(colour=guide_legend(title="Métodos"))

# Checagem de Resíduos

checkresiduals(media)

## 
##  Ljung-Box test
## 
## data:  Residuals from Mean
## Q* = 2866.8, df = 23, p-value < 2.2e-16
## 
## Model df: 1.   Total lags used: 24

checkresiduals(ingenuo)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 491.37, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

checkresiduals(ingenuo_saz)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 675.95, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

checkresiduals(drift_m)

## 
##  Ljung-Box test
## 
## data:  Residuals from Random walk
## Q* = 491.37, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

SUAVIZAÇÃO EXPONENCIAL SIMPLES| Cap. 8

modelo <- ets(volume,"ANN")

modelo

## ETS(A,N,N) 
## 
## Call:
##  ets(y = volume, model = "ANN") 
## 
##   Smoothing parameters:
##     alpha = 0.26 
## 
##   Initial states:
##     l = 44.4041 
## 
##   sigma:  7.7862
## 
##      AIC     AICc      BIC 
## 2114.629 2114.739 2124.837

plot(forecast(modelo, h=5))

modelo2 <- ses(volume, alpha = 0.2, h =3)

modelo3 <- ses(volume, alpha = 0.6, h =3)

modelo4 <- ses(volume,alpha = 0.9, h =3)

plot(modelo2, plot.conf=FALSE, ylab="Volume", xlab="Tempo", main="", fcol="white", type="o")

## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" não é um parâmetro
## gráfico

## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" não é
## um parâmetro gráfico

## Warning in axis(1, ...): "plot.conf" não é um parâmetro gráfico

## Warning in axis(2, ...): "plot.conf" não é um parâmetro gráfico

## Warning in box(...): "plot.conf" não é um parâmetro gráfico

lines(fitted(modelo2), col="blue", type="o")

lines(fitted(modelo3), col="red", type="o")

lines(fitted(modelo4), col="green", type="o")

lines(modelo2$mean, col="blue", type="o")

lines(modelo3$mean, col="red", type="o")

lines(modelo4$mean, col="green", type="o")

legend("topleft",lty=1, col=c(1,"blue","red","green"), c("Série Original", expression(alpha == 0.2), expression(alpha == 0.6), expression(alpha == 0.9)),pch=1)

Linear de Holt

Holt criou uma extensão ao método de suavização simples que permite prever dados com tendência que possui dois parâmetros α e β. Matematicamente, temos:

modelo5 <- holt(volume, alpha = 0.6, beta = 0.4)
modelo6 <- holt(volume)
modelo5$modelo6

## NULL

plot(modelo5)
lines(fitted(modelo5), col = "blue")
lines(fitted(modelo6), col = "red")

Holt-winter aditivo e multiplicativo

Uma evolução do modelo linear de Holt foi criado por Holt e Winter para possibilitar a modelagem de séries temporais por suavização exponencial que também possuam um componente sazonal. O método de Holt-Winters possui três equações para calcular os componentes lt de nível, Tt de tendência e st de sazonalidade, com os parâmetros α , β e γ. .

# testando o metodo aditivo quanto o multiplicativo para a serie de exemplo
ajuste_ad <- hw(volume, seasonal = "additive")
ajuste_mult <- hw(volume, seasonal = "multiplicative")


plot(volume)
lines(fitted(ajuste_ad), col = "blue")
lines(fitted(ajuste_mult), col = "red")

Tentativa de ARIMA | Cap. 9

# Dividindo a série (testagem e validação)

volume_test <- ts(volume[1:170], frequency = 12)
volume_validacao <- ts(volume[171:222], frequency = 12)


fit <- Arima(y = volume_test, order = c(1,1,0), seasonal = c(1,1,0), lambda = TRUE)

# order = c(1,1,0) | ordem dentro dos componentes não sazonais,1 componente regressiva, 1 diferenciação e 0 componentes de medias moveis
# seasonal " só que para componentes sazonais

summary(fit)

## Series: volume_test 
## ARIMA(1,1,0)(1,1,0)[12] 
## Box Cox transformation: lambda= TRUE 
## 
## Coefficients:
##           ar1     sar1
##       -0.4676  -0.3083
## s.e.   0.0706   0.0752
## 
## sigma^2 estimated as 11.61:  log likelihood=-414.95
## AIC=835.9   AICc=836.06   BIC=845.07
## 
## Training set error measures:
##                       ME     RMSE     MAE        MPE     MAPE      MASE
## Training set -0.07426952 3.253366 2.32795 -0.1261029 2.922217 0.3978111
##                    ACF1
## Training set -0.1314016

coeftest(fit)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.467599   0.070593 -6.6239 3.499e-11 ***
## sar1 -0.308349   0.075212 -4.0998 4.136e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#comparar a série original com os valores ajustados

plot(volume_test)

lines(fit$fitted, col = "blue")

# perdeu algumas observações no início pela diferenciação sazonal

# avaliando métricas

accuracy(volume_test, fit$fitted)

##                  ME     RMSE     MAE         MPE     MAPE       ACF1 Theil's U
## Test set 0.07426952 3.253366 2.32795 -0.03075505 2.910793 -0.1314016 0.3744947

# comparando com a base de validação nas 51 observações restantes


predi <- forecast(fit, h = 51)

plot(predi)

# plotando só as previsões contra a base de validação para averiguar se fez previu bem

plot(as.numeric(volume_validacao), type = 'l')
lines(as.numeric(predi$mean), col = 'blue')

accuracy(as.numeric(volume_validacao), as.numeric(predi$mean))

##                 ME     RMSE      MAE       MPE     MAPE
## Test set -18.06565 20.34778 18.12062 -25.24913 25.33154

# MAPE = 2.910793 2 %? (ver com a professora)


# Diagnóstico de Resíduos

tsdiag(fit)

qqnorm(fit$residuals)
qqline(fit$residuals)