### Bibliotecas utilizadas

library(tibble)
library(ggplot2)
library(seasonal)
library(dplyr)
library(TSA)
library(forecast)
library(fpp3)
library(stringr)
library(readr)

Banco de dados

O Brasil é o maior produtor mundial de café, mas ainda importa uma grande quantidade de café. O conjunto de dados nesta análise, inclui o valor mensal das importações de café no Brasil por 32,5 anos, de janeiro de 1990 a junho de 2022. Os dados foram retirados do site do Departamento de Agricultura dos Estados Unidos, Serviços de Agricultura Estrangeira. Para esse projeto, foi realizada uma análise de série temporais nos dados sobre Importação de Café no Brasil, disponível em : https://www.kaggle.com/datasets/srishti1209/coffee-imports-in-brazil

dados <- read_csv("Coffee_Time Series3.csv")
## Rows: 390 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Month Year
## dbl (1): Value of Import
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#head(dados)

Tratamento dos dados

Renomeando as variáveis

dados<- rename(dados, datacf = `Month Year`, vlimp= `Value of Import`)

Transformando a variável data

Criando uma sequência de datas seguindo os dados base e adicionando ao banco de dados A criação da variavel “dtcf” deve-se ao fato da dificuldade de trabalhar com os dados da varíavel data que veio no banco de dados,sendo assim criamos uma nova variavel com o mesmo tempo decorrido da database

dados$dtcf<- seq(as.Date("1990-01-01"),as.Date("2022-06-01"),by="month")
head(dados)
## # A tibble: 6 × 3
##   datacf vlimp dtcf      
##   <chr>  <dbl> <date>    
## 1 01/90  26173 1990-01-01
## 2 02/90  19159 1990-02-01
## 3 03/90  27451 1990-03-01
## 4 04/90  26171 1990-04-01
## 5 05/90  28544 1990-05-01
## 6 06/90  17922 1990-06-01
dados<-dados[,-c(1)]
head(dados)
## # A tibble: 6 × 2
##   vlimp dtcf      
##   <dbl> <date>    
## 1 26173 1990-01-01
## 2 19159 1990-02-01
## 3 27451 1990-03-01
## 4 26171 1990-04-01
## 5 28544 1990-05-01
## 6 17922 1990-06-01

Análise exploratória

str(dados)
## tibble [390 × 2] (S3: tbl_df/tbl/data.frame)
##  $ vlimp: num [1:390] 26173 19159 27451 26171 28544 ...
##  $ dtcf : Date[1:390], format: "1990-01-01" "1990-02-01" ...

Análise de dados faltantes

colSums(is.na(dados))
## vlimp  dtcf 
##     0     0

Podemos verificar que não temos dados faltantes nesses dados

Análise Descritiva dos dados, valor mínimo, 1º quartil, mediana, média, 3º quartil,valor máximo

summary(dados)
##      vlimp             dtcf           
##  Min.   :  9066   Min.   :1990-01-01  
##  1st Qu.: 28588   1st Qu.:1998-02-08  
##  Median : 54886   Median :2006-03-16  
##  Mean   : 65607   Mean   :2006-03-17  
##  3rd Qu.: 94362   3rd Qu.:2014-04-23  
##  Max.   :268247   Max.   :2022-06-01

Podemos observar então que a quantidade mínima de importação de café foi 9066, a quantidade média foi de 65607 e a quantidade máxima foi de 268247

Transformando os dados em tstibble para indicar a série temporal

Transformando os dados para o formato tstibble

serie <- dados %>%
  mutate(dtcf = yearmonth(dtcf)) %>%
  as_tsibble(index = dtcf)
serie
## # A tsibble: 390 x 2 [1M]
##    vlimp     dtcf
##    <dbl>    <mth>
##  1 26173 1990 jan
##  2 19159 1990 fev
##  3 27451 1990 mar
##  4 26171 1990 abr
##  5 28544 1990 mai
##  6 17922 1990 jun
##  7 26481 1990 jul
##  8 23562 1990 ago
##  9 31341 1990 set
## 10 43549 1990 out
## # … with 380 more rows
tail(serie)
## # A tsibble: 6 x 2 [1M]
##    vlimp     dtcf
##    <dbl>    <mth>
## 1 146751 2022 jan
## 2 177513 2022 fev
## 3 185460 2022 mar
## 4 197431 2022 abr
## 5 119075 2022 mai
## 6 181707 2022 jun

Visualizando a Série temporal

Vamos visualizar graficamente o comportamento da série temporal da quantidade importação de café no Brasil de 1990 a 2022

plot_serie= serie %>%
  autoplot(vlimp)+
  labs(title=" Quantidade de café importado ",
       y=" Quantidade ", x="Tempo")+
  geom_line()+
  geom_smooth(se = FALSE) +
  theme_bw()
plot_serie
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Ao analisar inicialmente a série podemos supor algumas características: * Tendência: A série mostra tendência de crescimento após 2000 * Sazonalidade: A série não mostra uma sazonalidade explicita, mas como a série é muito longa (1990 a 2022), pode ser que dificulte a visualização dessas sazonalidades * Estacionariedade: A série não é estacionária, até 2005 poderiamos supor isso, que os dados estariam em torno de uma média mas após esse ano há claramente uma mudança na quantidade de importação, saindo dessa média de 2005 a anos anteriores

Análise de autocorrelação

acf = serie %>%
  ACF(vlimp) %>% 
  autoplot() + labs(title="ACF")
#acf

pacf = serie%>%
  PACF(vlimp) %>% 
  autoplot() + labs(title="PACF") 
#pacf

gridExtra::grid.arrange(plot_serie, acf, pacf,layout_matrix=rbind(c(1,1),c(2,3) ))
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

### Decomposição da série

serie %>%
  model(
    classical_decomposition(vlimp, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Decomposição classica aditiva")
## Warning: Removed 6 rows containing missing values (`geom_line()`).

serie %>%
  model(
    classical_decomposition(vlimp, type = "multiplicative")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Decomposição classica multiplicativa")
## Warning: Removed 6 rows containing missing values (`geom_line()`).

Estudando a sazonalidade dos dados

serie %>%
  gg_season(vlimp, labels = "both") +
  labs(y = "Quantidade de importação",
       title = "Sazonalidade:Importação de Café",
       x="Tempo")

Teste de estacionariedade

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

serie %>%
  features(vlimp, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1      4.44        0.01

O valor do teste foi menor que 0,05 portanto há indicios que esta série não é estacionária

Aplicando o operador diferença

plot_serieD1 = serie %>%
  autoplot(difference(vlimp) ) +
  labs(title="Série diferenciada",
       y="% of GDP")
#plot_serieD1
acfd = serie %>%
  ACF(difference(vlimp)) %>% 
  autoplot() + labs(title="ACF")
#acf

pacfd = serie%>%
  PACF(difference(vlimp)) %>% 
  autoplot() + labs(title="PACF") 
#pacf

gridExtra::grid.arrange(plot_serieD1, acfd, pacfd,layout_matrix=rbind(c(1,1),c(2,3) ))
## Warning: Removed 1 row containing missing values (`geom_line()`).

Ao analisar o gráfico ACF e PACF podemos supor os modelos MA(4), AR(1),AR(4) No ACF também podemos indicar uma sazonalidade nos dados a cada 6 meses, podendo ser modelo SARIMA

teste para saber se a série com 1 diferenciação se tornou estacionária

serie %>%
  features(difference(vlimp), unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0377         0.1

O valor do teste foi maior que 0,05 portanto há indicios que esta série é estacionária

Definindo base de treino e teste para criação dos modelos

train <- serie %>%
  filter_index("1990" ~ "2020")
## Warning: `yearmonth()` may yield unexpected results.
## ℹ Please use arg `format` to supply formats.
## `yearmonth()` may yield unexpected results.
## ℹ Please use arg `format` to supply formats.
tail(train)
## # A tsibble: 6 x 2 [1M]
##    vlimp     dtcf
##    <dbl>    <mth>
## 1  95533 2019 ago
## 2  96513 2019 set
## 3  92604 2019 out
## 4  99221 2019 nov
## 5 107555 2019 dez
## 6 101334 2020 jan
test <- serie %>%
  filter_index("2021" ~ . );test
## Warning: `yearmonth()` may yield unexpected results.
## ℹ Please use arg `format` to supply formats.
## # A tsibble: 18 x 2 [1M]
##     vlimp     dtcf
##     <dbl>    <mth>
##  1 109252 2021 jan
##  2 112305 2021 fev
##  3 115099 2021 mar
##  4  94671 2021 abr
##  5 118536 2021 mai
##  6 108001 2021 jun
##  7  81897 2021 jul
##  8  98620 2021 ago
##  9  85452 2021 set
## 10 131631 2021 out
## 11 120896 2021 nov
## 12 149556 2021 dez
## 13 146751 2022 jan
## 14 177513 2022 fev
## 15 185460 2022 mar
## 16 197431 2022 abr
## 17 119075 2022 mai
## 18 181707 2022 jun

Analisando os dados de treino

Comportamento da série, acf,pacf

plot_seriet = train %>%
  autoplot(vlimp) +
  labs(title="Importação de Café no Brasil",
       y="Quantidade")
acft = train %>%
  ACF(vlimp) %>% 
  autoplot() + labs(title="ACF")


pacft = train %>%
  PACF(vlimp) %>% 
  autoplot() + labs(title="PACF") 


gridExtra::grid.arrange(plot_seriet, acft, pacft,layout_matrix=rbind(c(1,1),c(2,3) ))

Aplicando o teste de estacionariedade nos dados de treino

train %>%
  features(vlimp, unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1      3.94        0.01

O valor do teste foi menor que 0,05 portanto há indicios que esta série não é estacionária

Aplicando o operador diferença nos dados de treino

plot_trainD1 = train %>%
  autoplot(difference(vlimp) ) +
  labs(title="Série diferenciada",
       y="% of GDP")
#plot_serieD1
acfdt = train %>%
  ACF(difference(vlimp)) %>% 
  autoplot() + labs(title="ACF")
#acf

pacfdt = serie%>%
  PACF(difference(vlimp)) %>% 
  autoplot() + labs(title="PACF") 
#pacf

gridExtra::grid.arrange(plot_trainD1, acfdt, pacfdt,layout_matrix=rbind(c(1,1),c(2,3) ))
## Warning: Removed 1 row containing missing values (`geom_line()`).

* Ao analisar o gráfico ACF podemos supor o modelo MA(4) * Ao analisar o gráfico PACF podemos supor os modelo AR(1) 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).

Aplicando o teste de estacionariedade nos dados de treino com 1 diferenciação

train %>%
  features(difference(vlimp), unitroot_kpss)
## # A tibble: 1 × 2
##   kpss_stat kpss_pvalue
##       <dbl>       <dbl>
## 1    0.0167         0.1

Estimando modelos

fit <- train %>% model(
  ETS1 = ETS( vlimp ~ error("A") + trend("N") + season("N")),
  ETS2 = ETS( vlimp ~ error("A") + trend("A") + season("N")),
  ETS3 = ETS( vlimp~ error("A") + trend("M") + season("N")),
  ETS4 = ETS( vlimp ~ error("A") + trend( "N" ) + season("A")),
  ETS5 = ETS( vlimp ~ error("A") + trend( "N" ) + season("M")),
  ARIMA014 = ARIMA(vlimp ~ pdq(0,1,4)), #MA(4)
  ARIMA110 = ARIMA(vlimp ~ pdq(1,1,0)), #AR(1)
  #ARIMA114 = ARIMA(vlimp ~ pdq(1,1,4)), #AR(1)MA(4)
  SARIMA014012 = ARIMA(vlimp ~ pdq(0,1,4)+ PDQ(0,1,2)), #MA(4)MA(2)
  SARIMA110110 = ARIMA(vlimp ~ pdq(1,1,0)+PDQ(1,1,0)), #AR(1)#AR(1)
  #SARIMA114112 = ARIMA(vlimp ~ pdq(1,1,4)+PDQ(1,1,2)), #AR(1)MA(4) MA(2)AR(1)
  SARIMA014110 = ARIMA(vlimp ~ pdq(0,1,4)+ PDQ(1,1,0)), #MA(4)AR(1)
  SARIMA110012 = ARIMA(vlimp ~ pdq(1,1,0)+PDQ(0,1,2)), #AR(1)#MA(2)
  auto = ARIMA(vlimp),
  auto1 = ETS(vlimp))
## Warning in sqrt(diag(best$var.coef)): NaNs produzidos
#fit

Escolha do melhor modelo

glance(fit)%>% arrange(AICc) %>% select(.model:BIC)
## # A tibble: 13 × 6
##    .model        sigma2 log_lik   AIC  AICc   BIC
##    <chr>          <dbl>   <dbl> <dbl> <dbl> <dbl>
##  1 SARIMA014012 2.13e+8  -3837. 7688. 7688. 7715.
##  2 SARIMA110012 2.18e+8  -3842. 7692. 7692. 7708.
##  3 SARIMA014110 2.95e+8  -3886. 7785. 7785. 7808.
##  4 SARIMA110110 3.02e+8  -3892. 7790. 7790. 7801.
##  5 ARIMA014     2.23e+8  -3968. 7950. 7951. 7977.
##  6 ARIMA110     2.29e+8  -3974. 7957. 7957. 7972.
##  7 auto         2.30e+8  -3975. 7961. 7961. 7984.
##  8 auto1        5.29e-2  -4418. 8867. 8868. 8925.
##  9 ETS5         1.85e+8  -4492. 9013. 9015. 9072.
## 10 ETS4         2.10e+8  -4514. 9059. 9060. 9117.
## 11 ETS1         2.42e+8  -4546. 9099. 9099. 9110.
## 12 ETS3         2.42e+8  -4545. 9101. 9101. 9120.
## 13 ETS2         2.45e+8  -4547. 9105. 9105. 9124.

Para esses dados o melhor modelo foi o SARIMA014012

Análise de resíduos do modelo SARIMA014012

fit %>%
  select(SARIMA014012) %>%
  gg_tsresiduals()

#### Forecast a 40 passos

fc <- fit %>% forecast(h = 40)
#fc
# Grafico da serie com as previsoes 
fc %>%
  autoplot(train, level = NULL) +
  autolayer(
    test,
    colour = "black"
  ) +
  labs(
    y = "Quantidade ",
    x="Tempo",
    title = "Importação de Café no Brasil"
  ) +
  guides(colour = guide_legend(title = "Forecast"))
## Plot variable not specified, automatically selected `.vars = vlimp`

Analisando o desempenho dos modelos

accuracy(fc, test)%>% arrange(ME) 
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing. 
## 22 observations are missing between 2020 fev and 2023 mai
## # A tibble: 13 × 10
##    .model       .type     ME   RMSE    MAE   MPE  MAPE  MASE RMSSE  ACF1
##    <chr>        <chr>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 ETS2         Test  12476. 34325. 25007.  3.89  17.8   NaN   NaN 0.491
##  2 ETS1         Test  27514. 44195. 32840. 15.8   22.0   NaN   NaN 0.543
##  3 auto         Test  29445. 45246. 33779. 17.4   22.5   NaN   NaN 0.539
##  4 ARIMA014     Test  31665. 46698. 34860. 19.3   23.1   NaN   NaN 0.527
##  5 ARIMA110     Test  31741. 46758. 34898. 19.3   23.1   NaN   NaN 0.525
##  6 ETS4         Test  36397. 51269. 38811. 23.1   25.9   NaN   NaN 0.515
##  7 ETS5         Test  39433. 52521. 39850. 26.0   26.5   NaN   NaN 0.409
##  8 SARIMA014012 Test  39485. 53663. 40691. 25.8   27.2   NaN   NaN 0.448
##  9 auto1        Test  40677. 56155. 43031. 26.3   29.1   NaN   NaN 0.541
## 10 SARIMA110012 Test  41008. 54812. 42057. 27.1   28.3   NaN   NaN 0.450
## 11 SARIMA014110 Test  47156. 61102. 47220. 31.7   31.8   NaN   NaN 0.428
## 12 ETS3         Test  51468. 64436. 51468. 34.8   34.8   NaN   NaN 0.606
## 13 SARIMA110110 Test  53952. 67190. 53952. 37.1   37.1   NaN   NaN 0.447

Analisando a acurácia, temos que que o melhor modelo para esses dados foi o ETS2(A,A.N) que obteve o menor valor de Erro médio (ME) e assim seguindo para outras medidas de erro, apresentou ser o menor valor entre todos os modelos estimados. No modelo ETS(A,A,N), a primeira letra se refere ao componente de erro A(aditivo), a segunda ao componente de tendência A(aditivo), a terceira ao componente de sazonalidade N( não possui sazonalidade)

Análise de Resíduos para o modelo ETS2(A,A.N)

fit %>%
  select(ETS2) %>%
  gg_tsresiduals()

Criando um fit com os 2 melhores modelos

fitf <- serie %>%
  model(
    SARIMA014012  = ARIMA(vlimp ~ pdq(0,1,4)+ PDQ(0,1,2)), #MA(4)MA(2),
    ETS2 = ETS( vlimp ~ error("A") + trend( "A" ) + season("M")) 
  )

Previsão de um ano da série com os modelos finais

frf <- fitf %>% forecast(h = 12)
#frf

# Grafico da serie com as previsoes 
frf %>%
  autoplot(serie, level = NULL) +
  autolayer(
    test,
    colour = "black"
  ) +
  labs(
    y = "Quantidade ",
    x="Tempo",
    title = "Previsão de um ano da importação de café no Brasil"
  ) +
  guides(colour = guide_legend(title = "Forecast"))
## Plot variable not specified, automatically selected `.vars = vlimp`