Banco de dados

library(tidyverse)
library(fpp)
library(keras)
library(forecast)
library(tseries)
library(gridExtra)
library(astsa)
brazil<-read_rds("/opt/datasets/brazil.rds")
serie<-brazil %>%  filter(Subject == "Currency exchange rates, monthly average", TIME>="2000-01-01")

Análise descritiva

Antes de apicarmos qualquer tipo de modelo pra previsão da variável resposta, uma análise descritiva precisa ser feita, afim de estudar o máximo possível a respeito da variabilidade dos dados.

Observando o gráfico da série original, podemos dizer que os dados não apresentam característica de estacionariedade, apresentando variações ao longo do tempo. Como para série, precisa de estabelecer uma base de treinamento e teste, então no segundo gráfico representa justamente a parte da séries que vai ser usada para treinar o modelo, com a cor azul, e por final, a parte com cor vermelha vai ser usada para verificar o ajuste dos modelos impostos, sendo denominada de base de teste.

N<-round(.7*nrow(serie))
g1<- serie %>%dplyr::select(Time = TIME, Value) %>%
  ggplot() + geom_point(aes(x = Time, y = Value), alpha = .5) +
  geom_line(aes(x = Time, y = Value), color = "red", alpha = .5)  
  

g2<-serie %>%dplyr::select(Time = TIME, Value) %>%
  mutate(Base = if_else(row_number() >= N, "test", "train")) %>% 
  ggplot() + geom_point(aes(x = Time, y = Value, color=Base), alpha = .5) +
  geom_line(aes(x = Time, y = Value), alpha = .5)

gridExtra::grid.arrange(g1,g2,ncol=1,nrow=2)

Continuando a análise, o próximo passo é de aplicar um teste de não-estacionaridade, e dependendo do resultado, saber quantas diferenças precisam ser feitas afim de tornar essa série em estacionária.

myts<- ts(serie%>%dplyr::select(Value), start=c(2000,1),end=c(2019,8), frequency=12)
#print(myts)

#Número de diferenças necessárias para a série se tornar estacionária:#
ndiffs(myts)
## [1] 1

Utilizando o teste de dickey-fuller, afim de estudar os aspcectos de não normalidade dos dados, temos que para um nível \(\alpha = 5\%\), temos a evidência de não rejeitar a hipótese ue os dados são não normais.

adf.test(x = myts, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  myts
## Dickey-Fuller = -1.5054, Lag order = 6, p-value = 0.7838
## alternative hypothesis: stationary

Com essa função abaixo, podemos fazer o gráfico de autocorrelação e de autocorrelação parcial. Ao obseervar o primeiro gráfico, temos que os dados realmente não são estacionários.

acf2(myts)

##         ACF  PACF
##  [1,]  0.98  0.98
##  [2,]  0.95 -0.10
##  [3,]  0.92 -0.07
##  [4,]  0.88 -0.10
##  [5,]  0.85 -0.02
##  [6,]  0.81 -0.03
##  [7,]  0.78  0.05
##  [8,]  0.75  0.06
##  [9,]  0.72  0.00
## [10,]  0.69  0.01
## [11,]  0.67  0.01
## [12,]  0.64 -0.12
## [13,]  0.61  0.04
## [14,]  0.59  0.03
## [15,]  0.56 -0.02
## [16,]  0.54 -0.02
## [17,]  0.52  0.07
## [18,]  0.50  0.03
## [19,]  0.48 -0.01
## [20,]  0.47 -0.02
## [21,]  0.45 -0.06
## [22,]  0.43 -0.03
## [23,]  0.41  0.00
## [24,]  0.39  0.03
## [25,]  0.37  0.00
## [26,]  0.35 -0.01
## [27,]  0.33 -0.02
## [28,]  0.32 -0.03
## [29,]  0.30  0.03
## [30,]  0.28  0.01
## [31,]  0.27  0.03
## [32,]  0.26  0.00
## [33,]  0.25 -0.03
## [34,]  0.23  0.01
## [35,]  0.22 -0.04
## [36,]  0.21 -0.02
## [37,]  0.19 -0.02
## [38,]  0.18  0.05
## [39,]  0.17 -0.09
## [40,]  0.15 -0.14
## [41,]  0.12 -0.01
## [42,]  0.10 -0.04
## [43,]  0.07 -0.12
## [44,]  0.04 -0.01
## [45,]  0.01  0.03
## [46,] -0.02 -0.03
## [47,] -0.04 -0.01
## [48,] -0.07 -0.03

Após isso, podemos vizualizar como a série se comporta após a primeira diferença.

serie_diff <- serie %>% dplyr::select(Time = TIME, Value) %>% 
  mutate(Diff = Value - lag(Value)) %>%
  dplyr::filter(row_number() >= 2) %>%
  mutate( Base = if_else(row_number() >= N, "test", "train")) %>% 
  dplyr::select(Time, Value, Diff, Base)

g3<-serie_diff  %>% 
  ggplot() + geom_point(aes(x = Time, y = Diff, color=Base), alpha = .5) +
  geom_line(aes(x = Time, y = Diff, color=Base), alpha = .5)

g3

Aplicando um teste de para verificar a estacionariedade da nova série, apenas com a base de treinamento, temos que de fato, bastou apenas uma diferença para tornar a série em estacionária.

newts<- ts(serie_diff  %>% filter(Base=="train") %>% dplyr::select(Diff))

#nNúmero de diferenças necessárias para a série se tornar estacionária:
ndiffs(newts)
## [1] 0
#Teste de Dickey-Fuller para não estacionariedade da série
adf.test(x = newts, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  newts
## Dickey-Fuller = -5.4662, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Agora, basta normalizar a série e dividi-la em base de treino com base de teste.

scales <- serie_diff %>%
  filter(Base == "train") %>%
  summarise(min = min(Diff), max = max(Diff), Range = max - min)

serie_norm <- serie_diff %>%
  mutate(DiffPadrao = 2*((Diff - scales$min)/scales$Range - .5)) 

g4<-serie_norm %>% ggplot(aes(x = Time, y = DiffPadrao, color = Base)) + geom_point(alpha = .5) +
  geom_line(alpha = .5)
g4

newts_norm <- serie_norm %>% 
  filter(Base == "train") %>% 
  dplyr::select(DiffPadrao) %>% 
  ts(start=c(2000,1),end=c(2013,9), frequency=12)

Agora podemos análisar os gráfico da função de autocorrelação (ACF) e autocorrelação parcial (PACF) da base de treino normalizada, já que a gente só vai precisar dela para obter as nossas previsões, no qual, em prática a gente nunca pode utilizar a base de teste pois geralmente não se tem acesso a ela.

No gráfico do ACF, apresenta que o primero lag (\(h>0\)) foi significativo, o segundo lag pode ser considerado na margem de erro, não apresentando um lag significativo. No

acf2(newts_norm)

##         ACF  PACF
##  [1,]  0.33  0.33
##  [2,]  0.20  0.10
##  [3,]  0.08 -0.02
##  [4,]  0.10  0.07
##  [5,]  0.04 -0.02
##  [6,] -0.17 -0.23
##  [7,] -0.25 -0.17
##  [8,] -0.18 -0.02
##  [9,] -0.08  0.04
## [10,] -0.08 -0.01
## [11,] -0.11 -0.03
## [12,]  0.02  0.09
## [13,]  0.05 -0.02
## [14,]  0.12  0.04
## [15,]  0.12  0.06
## [16,]  0.12  0.04
## [17,]  0.10 -0.01
## [18,]  0.02 -0.07
## [19,]  0.12  0.14
## [20,]  0.08  0.06
## [21,] -0.01 -0.07
## [22,]  0.07  0.16
## [23,]  0.00  0.00
## [24,] -0.01 -0.08
## [25,]  0.03  0.09
## [26,] -0.12 -0.13
## [27,] -0.04  0.03
## [28,] -0.04  0.02
## [29,]  0.05  0.08
## [30,] -0.07 -0.07
## [31,] -0.11 -0.13
## [32,] -0.07 -0.02
## [33,]  0.03  0.08
## [34,]  0.08  0.01
## [35,] -0.05 -0.09
## [36,] -0.05 -0.02
## [37,] -0.05 -0.10
## [38,]  0.00 -0.06
## [39,]  0.02  0.03
## [40,] -0.11 -0.03
## [41,] -0.07 -0.07
## [42,]  0.04  0.08
## [43,]  0.04  0.02
## [44,]  0.02 -0.09
## [45,]  0.00  0.04
## [46,] -0.08 -0.10
## [47,] -0.02 -0.04
## [48,] -0.03  0.00

Modelo normal

Como estamos trabalhando na classe dos modelos arima, então pode-se trabalhar com a série original, selecionando apenas a base de treinamento.

serie_models_train <- serie %>% dplyr::select(Time = TIME, Value) %>%
  mutate(Base = if_else(row_number() >= N, "test", "train")) %>% 
  filter(Base =="train") %>% 
  dplyr::select(Value) %>%
  ts(start=c(2000,1),end=c(2013,8), frequency=12) 

serie_models_test<- serie %>% dplyr::select(Time = TIME, Value) %>%
  mutate(Base = if_else(row_number() >= N, "test", "train")) %>% 
  filter(Base =="test") %>% dplyr::select(Value)

Como foi visto anteriormente os possíveis parâmetros para o modelo arima são (1,1,0).

m1<-Arima(serie_models_train,order=c(1,1,0))
checkresiduals(m1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 28.174, df = 23, p-value = 0.2093
## 
## Model df: 1.   Total lags used: 24
autoplot(forecast(m1,71), xlim=c(2000,2019))

a<-forecast(m1,71) %>% as.data.frame()
pred_arima<-a$`Point Forecast`


EQM_arima<-sum((pred_arima - serie_models_test)^2)
#Soma dos quadrados dos resíduos para o modelo arima
#EQM_arima

Holt-Winters

No caso dos modelos de suavização exponencial, podemos trabaçhar com a série original

m2<-HoltWinters(serie_models_train)
checkresiduals(m2)

autoplot(forecast(m2,71), xlim=c(2000,2019))

a<-forecast(m2,71) %>% as.data.frame()
pred_holts<-a$`Point Forecast`

EQM_holts<-sum((pred_holts - serie_models_test)^2)
#Soma dos quadrados dos resíduos para o modelo de suavização exponenciall
#EQM_holts

Modelo de rede neural

Primeiramente, precisa estabelecer as bases que serão utilizadas no modelo. Todos essas quantidades foram antes aplicadas uma transformação de uma diferença e posteriormente normalizada em relação a base de trinamento.

serie_norm2 <- serie_norm %>%
    transmute(Time, Y = DiffPadrao, X1 = lag(Y), Base) %>% 
    filter(row_number() >= 2)

x_train <- serie_norm2 %>%
  filter(Base == "train") %>%
  dplyr::select(X1) %>%
  as.matrix()

y_train <- serie_norm2 %>%
  filter(Base == "train") %>%
  dplyr::select(Y) %>%
  as.matrix()

x_test <- serie_norm2 %>%
  filter(Base == "test") %>%
  dplyr::select(X1) %>%
  as.matrix()

y_test <- serie_norm2 %>%
  filter(Base == "test") %>%
  dplyr::select(Y) %>%
  as.matrix()
# Reshape the input to 3-dim
dim(x_train) <- c(length(x_train), 1, 1)

# specify required arguments
X_shape2 = dim(x_train)[2]
X_shape3 = dim(x_train)[3]
batch_size = 1 # must be a common factor of both the train and test samples
#units = 10 # can adjust this, in model tuninig phase

No modelo, foi utilizado apenas duas camadas sem contar com acamda de saída. A primeira foi utilizado 10 neurônios e na segunda apenas 5 neurônios.

model <- keras_model_sequential()
model%>%
layer_lstm(units=10, batch_input_shape = c(batch_size, X_shape2, X_shape3), stateful= TRUE) %>%
layer_dense(units=5) %>% 
layer_dense(units=3) %>% 
layer_dense(units = 1)


model %>% compile(
loss = 'mean_squared_error',
optimizer = optimizer_adam(lr= 0.02, decay = 1e-6 ),
metrics = c('mse')
)

Obtendo as previsões da base de teste.

Epochs = 50
for(i in 1:Epochs ){
  model %>% fit(x_train, y_train, epochs = 1, batch_size = batch_size,
verbose = 1, shuffle = FALSE)
model %>% reset_states()
}

Colocando as previsões obtidas anteriormente na escala original.

L = length(x_test)
predictions = numeric(L)
X = x_train[length(x_train)]
for(i in 1:L){
  #cat(" Valor X inicial:",X)
  dim(X) = c(1,1,1)
  y_sdt = model %>% predict(X, batch_size=batch_size)
  # invert scaling
  yhat = (y_sdt/2 + .5)*scales$Range + scales$min
  #cat("\n")
  #cat("Valor yhat previsto:", yhat)
  X<-y_sdt
  # store
  predictions[i] <- yhat
}

Obtendo um data frame, nele tem o valor real da base de teste, a previsão da base de teste, e o tempo correspondente desses dois valores. Com esse data frame pode ser obter o erro quadrático médio.

predict_df <- serie_norm %>% filter(row_number() >= N) %>%
  mutate(yhat = predictions) %>% 
  mutate(Vhat = yhat + lag(Value)) %>% 
  select(Time, Value, Vhat) %>%
  pivot_longer(cols = -Time,
               names_to = "Serie",
               values_to = "Valor")

g_comparativo<-predict_df  %>%
  ggplot(aes(x = Time, y = Valor, color = Serie)) +
  geom_point() +
  geom_line()

g_comparativo

Cálculo do EQM, foi criado um vetores coluna com a base de predição feita anteriormente, apenas pra fazer o comparativo com a previsão do modelo.

df_value_test<-predict_df %>%  filter(Serie=="Value") %>% select(Time, Valor, -Serie)
df_predict_test<-predict_df %>%  filter(Serie=="Vhat") %>% select(Time, Valor,-Serie)

#Primeiro dado da predição é NA, não sei direito o erro.

EQM_LSTM<-sum((df_predict_test$Valor[2:71] -df_value_test$Valor)^2)

#Soma dos quadrados dos resíduos para o modelo de Rede Neural
#EQM_LSTM

Observando a previsão na fase de teste, vizualmente ficou um bom ajuste, no qual as estimativas da base de teste ficaram próximas.

g_final<-serie %>% dplyr::select(Time = TIME, Value) %>%
  mutate(Base = if_else(row_number() >= N, "test", "train")) %>% 
  ggplot() + geom_point(aes(x = Time, y = Value, color=Base), alpha = .5) +
  geom_line(aes(x = Time, y = Value,color=Base), alpha = .5) +
  geom_point(aes(x = Time, y = Valor), color="#E7B800", alpha = .5, data=df_predict_test) +
  geom_line(aes(x= Time, y=Valor), color="#E7B800", alpha=.5, data=df_predict_test, size = 0.6)

g_final

Comparativo

Nesta seção será feita um comparativo de todos os modelos discutidos anteriormente. O critério de seleção será feito pelo erro quadrático médio. Com isso, podemos dizer que o modelo de série temporal por rede neural foi extramamente superior em relação aos demais.

##   EQM-ARIMA EQM-HOLTWINTERS EQM-LSTM
## 1  74.33742        12.36104 3.361468