Séries Temporais

Série Mensal:

dataset_mensal é a base com os valores de fluxo mensais de 1992 a 2009.

# Frequency --> número de observações por unidade de tempo 
# define a unidade de tempo (e.g. 12: unidade de tempo = ano)
tsMensal <- ts(dataset_mensal, frequency=12, start=1992)

Série Diária:

dataset_diario é a base com valores de fluxos diários de 2002 a 2009.

tsDiario <- msts(dataset_diario, start=c(2002,1,1), seasonal.periods=c(7, 365.25))

Nota: é importante observar os anos bissextos, 2004 e 2008.

Inspeção dos Dados

Série para um dia específico da semana

Série para um dia e mês específico da semana

Observar a influência dos feriados (outliers). Exemplo, 12 de outubro foi uma quinta-feira em 2006.

Inspeção dos Dados - 2

Pode-se observar que as séries são:

• Tendenciosas: apresentam padrão de crescimento ao longo do tempo.
• Sazonais: apresentam padrão recorrente em um determinado período de tempo.

Fluxo Mensal

Fluxo Diário

Observações

1. Sazonalidade de 1 ano na série mensal.
2. Na série diária, há também sazonalidade semanal.

Inspeção dos Dados - 3

Fluxo Diário para um intervalo arbitrário de 8 semanas

Inspeção dos Dados - 4

Auto-correlação

acf(dataset_diario, col="red")

Parcial:

pacf(dataset_diario, col="red")

partial_autocorr <- pacf(dataset_diario, col="red", lag.max = 600, plot = FALSE)
which((abs(partial_autocorr$acf) > 0.08) == TRUE)

##  [1]   1   3   4   5   6   7   8  13  14  15  20  21  22  29  35  36  49
## [18]  56 365 366

Objetivos

Realizar prospecções de curto, médio e longo prazo.

Curto Prazo

• 30 Dias
• 45 Dias

Para curto prazo, será utilizada a série temporal com dados diários (tsDiario).

Médio Prazo

• 4 Meses
• 6 Meses

Longo Prazo

• 1 Ano
• 2 Anos

Para médio e longo prazo, será utilizada a série temporal com dados mensais (tsMensal).

Separação dos dados

Separação das séries temporais em conjuntos de treinamento e de testes.

Curto Prazo (tsDiario)

Conjunto de Treinamento

tsDiarioTrain <- window(tsDiario, end=c(2008,183)) # até 30/06/2008

Conjuntos de Teste

30 dias:

# De 01/07/2008 a 30/07/2008
tsDiarioTest30Days <- window(tsDiario, start=c(2008,184), end=c(2008,213)) 

45 dias:

# 45 dias a partir De 1/1/2008
tsDiarioTest45Days <- window(tsDiario, start=c(2008,184), end=c(2008,228)) 

Médio/Longo Prazo (tsMensal)

Conjunto de Treinamento

tsMensalTrain <- window(tsMensal, end=c(2007, 12)) # De Jan 1992 a Dez 2007

Conjuntos de Teste

4 meses:

tsMensalTest4mth <- window(tsMensal, start=2008, end=c(2008, 4)) # De Jan 2006 a Mar 2006

6 meses:

tsMensalTest6mth <- window(tsMensal, start=2008, end=c(2008, 6)) # De Jan 2006 a Jun 2006

1 ano:

tsMensalTest1yr <- window(tsMensal, start=2008, end=c(2008, 12)) # De Jan 2006 a Jan 2007

2 anos:

tsMensalTest2yr <- window(tsMensal, start=2008, end=c(2009, 12)) # De Jan 2006 a Jan 2008

Transformação Box-Cox

É possível observar um fator multiplicativo na componente sazonal. Por isso, a transformação Box-Cox será utilizada em alguns dos métodos apresentados na sequência.

• \(W_t = \begin{cases} \log\left(y_t\right) & \lambda = 0\\ \frac{(y_t^\lambda - 1)}{\lambda} & \lambda \neq 0\end{cases}\)

Série original

plot(tsMensalTrain, xlab="Anos", ylab="Fluxo Mensal")

Série transformada

Calcula-se o fator \(\lambda\) para a transformação:

lam <- BoxCox.lambda(tsMensalTrain)
lam

## [1] -0.06928

Série Resultante:

tsMensalBoxCox <- BoxCox(x = tsMensalTrain, lam)

plot(tsMensalBoxCox, col="red")

Bibliotecas Utilizada

Prospecções Estatísticas

ver (R. J. Hyndman 2014) e (Leek 2014):

library(forecast)

Redes Neurais:

library(nnet)

Predição a Médio Prazo

• ARIMA
• Exponential Smoothing

Modelo ARIMA

ARIMA (Autoregressive Integrated Moving Average)

• Combinação dos modelos:
  • auto-regressivo (AR)
  • médias-móveis (MA)

• Método ARMA, acrescido da etapa de integração (I).

• Caso a série temporal seja estacionária, o método se reduz ao ARMA.

Equação de Predição

\[y(t) = \sum_{i=1}^{p}\varphi_iy(t-i) + \sum_{j=1}^{q}\theta_j\epsilon(t-j)\]

Notação

\(ARIMA(p,d,q)\), onde:

• p = número de valores passados da variável (AR)

• d = número de diferenças não-sazonais (I)

• q = número de erros de predição passados (MA)

Função auto.arima no R

• A função auto.arima ajusta o melhor modelo ARIMA à série temporal, com base no menor parâmetro AIC.

• Protótipo com os principais parâmetros:

auto.arima(x, max.p=\(5\), max.d=\(2\), max.q=\(5\), stationary=\(FALSE\), seasonal=\(TRUE\), approximation=\(TRUE\), xreg=\(NULL\))

• Com todos os parâmetros automáticos, a função auto.arima forneceu predições ruins.

Regressores Adicionais

• Aumentar o número de termos na equação de regressão do modelo ARIMA.

• Adicionar ao modelo informações que possuem alta correlação com a série temporal em questão.

Regressores incorporados ao modelo ARIMA

Indicação dos meses dos conjuntos mensais de treinamento/teste.

Regressor que pode assumir \(12\) categorias possíveis: janeiro, fevereiro, …, dezembro.

Quantas segundas, terças, …, domingos há em cada mês.

Base dataset_diario_1992-2009.csv, com as datas diárias de Jan de 1992 a Dez de 2009.

Treinamento do Modelo ARIMA em prospecções de Médio Prazo

Regressores externos utilizados no treinamento

# datas diárias do conjunto de treinamento, de Jan 1992 a Dez 2007
datasDiarioTrain_Mensal <- window(index(dataset_diario_full), end=5844)

# Série contendo as datas do conjunto de treinamento, de Jan 1992 a Dez 2007
datasMensalTrain <- dataset_mensal[1:192]

# Indica quais são os meses (janeiro, ..., dezembro) do conjunto de treinamento
meses <- getMonths( format(index(datasMensalTrain),"%m") )

# conta quantas segundas,...,domingos há em cada mês do conjunto de treinamento
diasSemanasMes_treino <- countWeekdaysPerMonth(datasDiarioTrain_Mensal)

Modelo ARIMA utilizando os regressores: Indicação dos meses e Quantidade dos dias da semana de cada mês

arima_model_mensal <- auto.arima(tsMensalTrain, approximation=FALSE, lambda=BoxCox.lambda(tsMensalTrain), xreg=cbind(meses,diasSemanasMes_treino))

Método ARIMA em prospecções de Médio Prazo

Para 4 meses

mesesAPrever <- 4

# Série contendo as datas do conjunto de testes, de Jan 2008 a Abr 2008
datasMensal4mth <- dataset_mensal[193:196]

# Indica quais são os meses (janeiro, ..., dezembro) do conjunto de testes
mesesFuturo4 <- getMonths( format(index(datasMensal4mth),"%m") )


# datas diárias do conjunto de teste (4 meses)
datasDiarioTeste4mth_Mensal <- window(index(dataset_diario), start=2192, end=2312)

# conta quantas segundas,...,domingos há em cada mês do conjunto de teste (4 meses)
diasSemanasMes_futuro_4mth <- countWeekdaysPerMonth(datasDiarioTeste4mth_Mensal)

Predição:

fcast_arima_4mth <- forecast(arima_model_mensal, xreg=cbind(mesesFuturo4,diasSemanasMes_futuro_4mth), h=mesesAPrever)

plot(fcast_arima_4mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest4mth, col="red")

Acurácia:

accuracy(fcast_arima_4mth, tsMensalTest4mth)

##                    ME RMSE  MAE       MPE  MAPE   MASE    ACF1 Theil's U
## Training set   -15.24 8727 6389 -0.005435 1.256 0.1598 -0.0472        NA
## Test set     -7837.87 9728 8071 -0.980798 1.012 0.2019 -0.5471     0.162

Método ARIMA em prospecções de Médio Prazo

Para 6 meses

mesesAPrever <- 6

Predição:

datasMensal6mth <- window(index(dataset_mensal), start=193, end=198)
mesesFuturo6 <- getMonths( format(datasMensal6mth,"%m") )

# quantas segundas, tercas, ..., domingos há em cada mês do conjunto de teste (6 meses)
datasDiarioTeste6mth_Mensal <- window(index(dataset_diario), start=2192, end=2373)

# conta quantas segundas,...,domingos há em cada mês do conjunto de teste (6 meses)
diasSemanasMes_futuro_6mth <- countWeekdaysPerMonth(datasDiarioTeste6mth_Mensal)

Predição:

fcast_arima_6mth <- forecast(arima_model_mensal, xreg=cbind(mesesFuturo6,diasSemanasMes_futuro_6mth), h=mesesAPrever)

plot(fcast_arima_6mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest6mth, col="red")

Acurácia:

accuracy(fcast_arima_6mth, tsMensalTest6mth)

##                    ME RMSE  MAE       MPE   MAPE   MASE    ACF1 Theil's U
## Training set   -15.24 8727 6389 -0.005435 1.2559 0.1598 -0.0472        NA
## Test set     -7345.11 9042 7501 -0.908452 0.9291 0.1876 -0.5392     0.184

Exponential Smoothing

• Médias móveis simples pondera igualmente as observações passadas em uma série temporal.
• Exponential smoothing atribui pesos que decrescem exponencialmente ao longo do tempo, levando em consideração:
  • Erro
  • Tendência
  • Sazonalidade
• Estes parâmetros podem ser escolhidos como (e.g. tendência \(T_h\), assumindo um termo de nível \(l\) e termo de crescimento \(b\), para \(h\) períodos temporais e \(0 < \phi < 1\)):
  • Nulo: \(T_h = l\)
  • Aditivo: \(T_h = l + bh\)
  • Aditivo amortecido: \(T_h = l + \left( \phi + \phi^2 + \cdots + \phi^h \right)b\)
  • Multiplicativo: \(T_h = l b^h\)
  • Multiplicativo amortecido: \(T_h = l b^{\left( \phi + \phi^2 + \cdots + \phi^h \right)}\)

Exponential Smoothing - 2

Ver (R. J. Hyndman and Athanasopoulos 2013). O autor possui um livro exclusivamente sobre o assunto em (R. Hyndman et al. 2008).

Para cada um dos 15 métodos, há dois modelos possíveis: com erro aditivo e com erro multiplicativo.

Há, portanto, 30 métodos distintos.

Exponential Smoothing no R

• A função ets(ts) sem argumentos além da série temporal ts determina automaticamente o método apropriado.
• O algoritmo avalia os 30 métodos e busca aquele que tiver o menor AIC.
  • AIC - “Akaike Information Criterion”: uma estimativa da informação perdida quando um determinado modelo é usado para representar o verdadeiro processo que gera a série temporal.
• Ao final, o objeto retornado por ets apresenta o modelo escolhido e o AIC resultante.
• A função ets não trabalha com séries de sazonalidade superior a 24 unidades temporais. Portanto, será utilizada somente para prospecções realizadas com dados mensais (sazonalidade = 12).
• O modelo ets parametrizado automaticamente deve fornecer resultados muito precisos para prospecções de poucos pontos (por exemplo, 4 pontos).

Método Exponential Smoothing em Médio Prazo

Para 6 meses

mesesAPrever <- 6

Modelo:

etsMensal <- ets(tsMensalTrain)

Prospecção:

fcastMensal6mth <- forecast(etsMensal, h=mesesAPrever)

Comparar predição com valores reais do conjunto de teste:

plot(fcastMensal6mth, xlab="Anos", ylab="Fluxo Mensal", include = 36);
lines(tsMensalTest6mth, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

Acurácia:

accuracy(fcastMensal6mth, tsMensalTest6mth)

##                  ME  RMSE  MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set  856.7  9849 7308  0.13394 1.442 0.1828 -0.02704        NA
## Test set     -755.7 10937 8653 -0.06425 1.104 0.2164 -0.63165    0.2848

Predição a Longo Prazo

• Prospecção com decomposição sazonal de Loess (STL)

• Rede Neural

Prospecção com decomposição sazonal de Loess - STLF

1. Aplica-se a decomposição STL
2. Remove-se a sazonalidade da série temporal
3. Um modelo é treinado na série resultante
4. Faz-se a prospecção
5. Adiciona-se novamente o último período da sazonalidade estimada nos resultados

mesesAPrever <- 24

Dados com sazonalidade removida:

tsMensal.stl <- stl(tsMensalTrain[,1], s.window=12)
# Seasonally adjusted data constructed by removing the seasonal component.
plot(seasadj(tsMensal.stl))

stlf_model <- stlf(tsMensalTrain[,1], lambda = lam)
stlf_fcast <- forecast(stlf_model, h=mesesAPrever)
plot(stlf_fcast, xlab="Anos", ylab="Fluxo Mensal")
lines(tsMensalTest2yr, col="red")

Acurácia:

accuracy(stlf_fcast, tsMensalTest1yr)

##                   ME  RMSE   MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set  -25.42  8186  6175 -0.01584 1.238 0.1544 -0.02348        NA
## Test set     8132.24 17138 13696  0.92760 1.619 0.3425 -0.11417    0.4831

Rede Neural

Predição a curto prazo

• Rede Neural
• TBATS

Rede Neural

Pacotes

require(quantmod)
require(nnet)

Escolha das entradas:

partial_autocorr <- pacf(dataset_diario, col="red", lag.max = 600, plot = FALSE)
which((abs(partial_autocorr$acf) > 0.1) == TRUE)

##  [1]   1   3   4   5   6   7   8  13  14  15  21  22  29  35  36 366

dat <- data.frame( y=dataset, 
                   Lag(dataset,1),                                           
                   Lag(dataset,3),
                   Lag(dataset,4),                   
                   Lag(dataset,7),
                   Lag(dataset,14),
                   Lag(dataset,15),                      
                   Lag(dataset,35),
                   Lag(dataset,364),
                   month=as.numeric(format(index(dataset), "%m")), 
                   year=(as.numeric(format(index(dataset), "%y")) - 1), 
                   holiday=isHoliday(dataset))

Modelo

model  <- nnet(y ~ ., data = datTrain, repeats=50, size=10, decay=0.127, linout=TRUE, trace=FALSE, 
               maxit=10000)
model

## a 11-10-1 network with 131 weights
## inputs: Lag.1 Lag.3 Lag.4 Lag.7 Lag.14 Lag.15 Lag.35 Lag.364 month year holidayTRUE 
## output(s): y 
## options were - linear output units  decay=0.127

Predição 30 dias

Acurácia

mape(datTest$y[1:nPredict], fcast)*100

## [1] 0.963

Predição 45 dias

Acurácia

mape(datTest$y[1:nPredict], fcast)*100

## [1] 1.786

Modelo TBATS

• Técnica recente, publicada em 2011 (De Livera, Hyndman, and Snyder 2011).

• A sigla TBATS é composta das iniciais das técnicas utilizadas: Trigonométrica, transformação Box-cox, correção de erro pela técnica ARMA, e componentes Trend e Sazonal.

• Desenvolvido para séries temporais com características especiais de sazonalidade.

Justificativa

• Por inspeção visual do fluxo diário, supõe-se que a série temporal apresenta sazonalidades semanal e anual.

Modelo TBATS, a partir da série multi-sazonal do conjunto de treinamento:

tbats_model <- tbats(tsDiarioTrain,                      
                     use.box.cox = TRUE, 
                     use.trend = TRUE, 
                     use.damped.trend = FALSE,
                     seasonal.periods = c(7, 365.25),
                     use.arma.errors = TRUE)                    

Decomposição do fluxo diário em termos trend+season pelo método TBATS

O método TBATS é indicado para a previsão do fluxo diário, pois:

1. Sazonalidades múltiplas: semanal e anual

2. Sazonalidade semanal é de alta frequência

Método TBATS em prospecções de Curto Prazo

Para 30 dias

diasAPrever <- 30

Predição:

tbats_fc30Days <- forecast(tbats_model, h=diasAPrever)

Acurácia:

##                    ME  RMSE   MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set    7.641 750.8 426.2 -0.07899 2.041 0.2136 -0.01679        NA
## Test set     -198.853 502.6 358.1 -0.75395 1.311 0.1795  0.71247     0.258

Método TBATS em prospecções de Curto Prazo

Para 45 dias

diasAPrever <- 45

Predição:

tbats_fc45Days <- forecast(tbats_model, h=diasAPrever)

Acurácia:

##                   ME  RMSE   MAE      MPE  MAPE   MASE     ACF1 Theil's U
## Training set   7.641 750.8 426.2 -0.07899 2.041 0.2136 -0.01679        NA
## Test set     -99.500 499.1 371.6 -0.41877 1.331 0.1862  0.69601    0.2385

Benchmark dos Métodos Utilizados

Curto prazo

##           Training.set Test..30d. Test..45d.
## Mean            14.656     22.506     23.334
## Näive            6.214      4.674      5.004
## N-Sazonal        8.720      9.066      9.304
## Drift            6.216      4.788      4.993
## Regressao        5.868      7.383      7.777
## ARIMA            5.330      5.415      6.111
## TBATS            2.036      1.339      1.327

Médio e Longo prazo

##           Training.set Test..4mths. Test..6mths. Test..1yr. Test..2yr.
## Mean            25.582       35.295       36.430     39.867     41.232
## Näive            4.440        5.639        3.982      5.708      6.801
## N-Sazonal        7.680        6.051        5.902      6.898      9.046
## Drift            4.470        6.542        5.016      4.722      5.210
## Regressao        4.741        6.894        7.387     10.023     10.398
## ETS              1.442        1.447        1.104      2.603      2.964
## Arima            1.625        1.933        1.539      1.841      2.390

Referências