Traduzido de: Complete guide to create a Time Series Forecast (with Codes in Python)

Autor: Aarshay Jain

Introdução

Séries Temporais (referido a partir de agora como TS) é considerada uma das técnicas menos conhecidas no espaço de análise (eu mesmo tinha poucas pistas sobre isso até alguns dias atrás). Mas como você sabe, o nosso Mini Hackathon inaugural foi baseado nela, então iniciei uma jornada para aprender os passos básicos para a resolução de um problema de TS e aqui vou compartilha-lo com você. Estas informações irão ajudá-lo a obter um modelo decente no nosso hackaton de hoje.

Antes de começar o artigo, recomendo a leitura do tutorial completo sobre Time Series Modelagem em R, que é como um pré-requisito para este artigo. O tutorial concentra-se em conceitos fundamentais baseados no R nos quais vou me basear para resolver um problema juntamente com códigos em Python. Existem muitos recursos para TS em R mas muito poucos estão lá para Python, então eu vou usar o Python neste artigo.

Nossa jornada irá passar pelas seguintes etapas:

1. O que torna TS algo tão especial?
2. Carregando e manuseando TS com Pandas
3. Como verificar a Estacionaridade de uma TS?
4. Como fazer uma TS estacionária?
5. Projetando uma série temporal

1. O que torna TS algo tão especial?

Como o nome sugere, TS é um conjunto de pontos de dados coletados em intervalos de tempo constante. Estes são analisados para determinar a tendência de longo prazo, de modo a prever o futuro ou executar alguma outra forma de análise. Mas o que faz um TS diferente, por exemplo, de um problema de regressão comum? São 2 coisas:

1. É dependente do tempo. Assim, o pressuposto básico de um modelo de regressão linear, que as observações são independentes, não se sustenta neste caso.
2. Junto com uma tendência crescente ou decrescente, a maioria dos TS têm alguma forma de tendência de sazonalidade, ou seja, variações específicas a um determinado período de tempo. Por exemplo, se você olhar as vendas de um casaco de lã ao longo do tempo, você vai sempre encontrar picos de vendas nas estações de inverno.

2. Carregando e Manipulando TS com Pandas

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6

data = pd.read_csv('AirPassengers.csv') 
print data.head() 
print '\n Data Types:' 
print data.dtypes

dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m') 
data = pd.read_csv('AirPassengers.csv', parse_dates='Month', index_col='Month',date_parser=dateparse) 
print data.head()

Vamos entender a argumentos um por um:

1. parse_dates: especifica a coluna que contém as informações de data e hora. Como dizemos acima, o nome da coluna é “Mês”.
2. index_col: uma ideia central por trás do uso do Pandas para dados de TS é que o índice tem de ser a informação que descreve a data em tempo variável. Portanto, este argumento diz para o pandas usar a coluna “Mês” como índice.
3.  date_parser: especifica uma função que converte uma cadeia de entrada na variável de data e hora. O padrão Pandas lê os dados no formato AAAA-MM-DD HH: MM: SS ‘. Se os dados não estiverem neste formato, o formato tem de ser definido manualmente. Algo semelhante à função dataparse definida aqui pode ser utilizado para este fim.

Agora que podemos ver os dados objeto de tempo como o índice e #Passengers como a coluna, podemos confrontar o tipo de dados do índice com o seguinte comando:

data.index

Observe o dtype = ‘datetime [ns] “, que confirma que é um objeto de data e hora. Como uma preferência pessoal, eu converteria a coluna em um objeto Series para evitar referir-se a nomes de colunas cada vez que usar os TS. Sinta-se livre para usar como dataframe, se assim funcionar melhor para você.

ts = data[‘#Passengers’] ts.head(10)

Antes de prosseguir, eu vou discutir algumas técnicas de indexação para os dados de TS. Vamos começar por selecionar um determinado valor no objeto Series. Isto pode ser feito de 2 formas:

#1. Specific the index as a string constant: 
ts['1949-01-01'] 

#2. Import the datetime library and use 'datetime' function: 
from datetime import datetime ts[datetime(1949,1,1)]

Ambos irão retornar o valor “112”, que também pode ser confirmado a partir da saída anterior. Suponha que queremos todos os dados até maio de 1949. Isso pode ser feito de 2 maneiras:

#1. Specify the entire range:
ts['1949-01-01':'1949-05-01']

#2. Use ':' if one of the indices is at ends:
ts[:'1949-05-01']

Ambos renderiam o seguinte resultado:

Deve-se notar duas coisas:

1. Ao contrário de indexação numérica, o índice final é incluído aqui. Por exemplo, se indexarmos uma lista como um [: 5], em seguida ele iria retornar os valores em índices – [0,1,2,3,4]. Mas aqui o índice ‘1949/05/01’ foi incluído na saída.
2. Os índices têm de ser ordenados para funcionar. Se você embaralhar aleatoriamente o índice, não vai funcionar.
   

Considere um outro exemplo onde você precisa de todos os valores do ano 1949. Isto pode ser feito com:

ts['1949']

A parte do mês foi omitida. Da mesma forma, se você tiver todos os dias de um determinado mês, a parte do dia pode ser omitida.

Agora, vamos passar para a análise do TS.

3. Como verificar se uma TS é estacionária?

A TS é dita estacionária se as suas propriedades estatísticas, tais como a média e variância permanecem constantes ao longo do tempo. Mas por que isso é importante? A maioria dos modelos TS trabalha com o pressuposto de que a TS é estacionária. Intuitivamente, podemos entender que se uma TS tem um comportamento particular ao longo de um tempo, há uma probabilidade muito elevada de que seguirá o mesmo comportamento no futuro. Além disso, as teorias relacionadas com séries estacionárias são mais maduras e mais fáceis de implementar, em comparação com séries não-estacionárias.

1. Uma média constante
2. Uma variância constante
3. Uma autocovariância que não dependa do tempo

Eu vou pular os detalhes, uma vez que eles são muito claramente definido neste artigo. Vamos passar para os modos de teste de Estacionaridade.

plt.plot(ts)

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
    
    #Determing rolling statistics
    rolmean = pd.rolling_mean(timeseries, window=12)
    rolstd = pd.rolling_std(timeseries, window=12)

    #Plot rolling statistics:
    orig = plt.plot(timeseries, color='blue',label='Original')
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    #Perform Dickey-Fuller test:
    print 'Results of Dickey-Fuller Test:'
    dftest = adfuller(timeseries, autolag='AIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical Value (%s)'%key] = value
    print dfoutput

O código é bastante simples, vamos executá-lo para a nossa série de entrada:

test_stationarity(ts)

Embora a variação de desvio padrão seja pequena, a média está claramente aumentando com o tempo e isto não é uma característica de série estacionária. Além disso, a estatística de teste é muito mais do que os valores críticos. Note que os valores assinalados devem ser comparados e não os valores absolutos.

A seguir, vamos discutir as técnicas que podem ser usadas para levar este TS à direção estacionária.

4. Como fazer TS estacionária?

Embora o pressuposto de estacionaridade seja considerado em muitos modelos TS, quase nenhuma série temporal na prática é estacionária. Assim, os estatísticos descobriram maneiras de fazer séries estacionárias que vamos discutir agora. Na verdade, é quase impossível fazer uma série perfeitamente estacionária, mas tentaremos chegar o mais próximo possível.

Existem 2 principais razões por trás da não estacionaridade de uma TS:

Estimando & Eliminando Tendência

ts_log = np.log(ts)
plt.plot(ts_log)

Neste caso mais simples, é fácil ver para tendência futura nos dados. Mas não é muito intuitivo no caso de haver ruídos. Assim, podemos usar algumas técnicas para estimar ou modelar a tendência e, em seguida, removê-la da série. Podem haver muitas maneiras de fazer isso e algumas das mais comumente utilizadas são:

1- Agregação – usando a média para um período de tempo como médias semanais ou mensais
2- Suavização – usando médias moveis
3- Ajuste polinomial – encaixando um modelo de regressão

Vou usar suavização aqui, mas você deve tentar outras técnicas que podem funcionar para outros problemas. Suavização refere-se às estimativas contínuas, ou seja, considerando os últimos casos. Pode ser feito de várias maneiras, aqui vou discutir duas delas.

Média Móvel

Nessa abordagem, consideramos a média de ‘K’ valores consecutivos, dependendo da frequência da série temporal. Aqui nós podemos tomar a média do último ano, ou seja os últimos 12 valores.

O Pandas tem funções específicas definidas para determinar as estatísticas de rolagem.

moving_avg = pd.rolling_mean(ts_log,12)
plt.plot(ts_log)
plt.plot(moving_avg, color='red')

A linha vermelha mostra a média móvel. Iremos subtrair isso da série original. Note-se que uma vez que estamos tomando média dos últimos 12 valores, a média móvel não está definida para 11 primeiros valores. Isto pode ser observado como:

ts_log_moving_avg_diff = ts_log - moving_avg
ts_log_moving_avg_diff.head(12)

Observe os primeiro 11 como sendo NaN. Vamos usar esses valores NaN para testar a estacionaridade.

ts_log_moving_avg_diff.dropna(inplace=True)
test_stationarity(ts_log_moving_avg_diff)

Parece uma série muito melhor. Os valores parecem ter variado um pouco, mas não há nenhuma tendência específica. Além disso, a estatística do teste é menor do que os valores críticos de 5% para que possamos dizer com 95% de confiança que esta é uma série estacionária.

expwighted_avg = pd.ewma(ts_log, halflife=12)
plt.plot(ts_log)
plt.plot(expwighted_avg, color='red')

ts_log_ewma_diff = ts_log - expwighted_avg 
test_stationarity(ts_log_ewma_diff)

Eliminando Tendência e Sazonalidade

As técnicas simples de redução de tendência discutidos anteriormente não funcionam em todos os casos, particularmente naqueles com alta sazonalidade. Vamos discutir duas maneiras de remoção de tendência e de sazonalidade:

1- Diferenciação – tomando a diferença com um determinado intervalo de tempo

2- Decomposição – modelando tanto a tendência quanto a sazonalidade e removendo-as a partir do modelo.

Diferenciação

Um dos métodos mais comuns de lidar tanto com tendência quanto sazonalidade é a diferenciação. Nessa técnica, tomamos a diferença de observação em um instante particular com a do instante anterior. Isso funciona particularmente bem para estacionaridade crescente. Diferenciação de primeira ordem pode ser feito em pandas com:

ts_log_diff = ts_log - ts_log.shift()
plt.plot(ts_log_diff)

Parece ter reduzido a tendência consideravelmente. Vamos verificar usando nossos gráficos:

ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)

Podemos ver que a média e desvio padrão têm pequenas variações com o tempo. Além disso, a estatística de teste Dickey-Fuller é menor do que o valor crítico de 10%, assim a TS é estacionária com 90% de confiança. Nós também podemos ter diferenças de segunda ou terceira ordem que poderiam obter resultados ainda melhores em determinadas aplicações. Deixo isso para você experimentar.

Decomposição

Nesta abordagem, tanto tendência quanto sazonalidade são modeladas separadamente e a parte restante da série é devolvida. Vou pular as estatísticas e vou direto para os resultados:

from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log) 
trend = decomposition.trend 
seasonal = decomposition.seasonal 
residual = decomposition.resid 
plt.subplot(411) 
plt.plot(ts_log, label='Original') 
plt.legend(loc='best') 
plt.subplot(412) 
plt.plot(trend, label='Trend') 
plt.legend(loc='best') 
plt.subplot(413) 
plt.plot(seasonal,label='Seasonality') 
plt.legend(loc='best') 
plt.subplot(414) 
plt.plot(residual, label='Residuals') plt.legend(loc='best') 
plt.tight_layout()

Aqui podemos ver que a tendência e a sazonalidade são separadas  e podemos modelar os resíduos. Vamos verificar a estacionaridade dos resíduos:

ts_log_decompose = residual 
ts_log_decompose.dropna(inplace=True) 
test_stationarity(ts_log_decompose)

#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf

lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')

#Plot ACF: 
plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')

#Plot PACF: 
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray') 
plt.title('Partial Autocorrelation Function') 
plt.tight_layout()

from statsmodels.tsa.arima_model import ARIMA

Os valores p, d, q podem ser especificados usando o argumento de ordem de ARIMA que tomar uma tupla (p, d, q). Vamos modelar os 3 casos:

 Modelo AR

model = ARIMA(ts_log, order=(2, 1, 0))  
results_AR = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))

 Modelo MA

model = ARIMA(ts_log, order=(0, 1, 2))  
results_MA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))

Modelo Combinado

model = ARIMA(ts_log, order=(2, 1, 2)) 
results_ARIMA = model.fit(disp=-1) 
plt.plot(ts_log_diff) 
plt.plot(results_ARIMA.fittedvalues, color='red') 
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))

Aqui podemos ver que os modelos AR e MA têm quase a mesma RSS, mas combinado é significativamente melhor. Agora, ficamos com um último passo, ou seja, trazer esses valores de volta para a escala original.

Trazendo de volta à escala original

Uma vez que o modelo combinado deu melhor resultado, vamos traze-lo de volta para os valores originais e ver o quão bem ele executa lá. O primeiro passo seria para armazenar os resultados previstos como uma série separada e observá-la.

predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True) 
print predictions_ARIMA_diff.head()

predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum() 
print predictions_ARIMA_diff_cumsum.head()

Você pode rapidamente fazer algumas contas de cabeça usando a saída anterior para verificar se estão corretos. Em seguida temos que adicioná-los ao número base. Para isso, vamos criar uma série com todos os valores como número base e adicionar as diferenças a ela. Isto pode ser feito como:

predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index)
predictions_ARIMA_log=predictions_ARIMA_log.add
(predictions_ARIMA_diff_cumsum,fill_value=0) predictions_ARIMA_log.head()

Aqui, o primeiro elemento é o número base em si e a partir dele os valores cumulativamente somados. O último passo é tomar o expoente e comparar com a série original.

predictions_ARIMA = np.exp(predictions_ARIMA_log) 
plt.plot(ts) plt.plot(predictions_ARIMA) 
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))

Finalmente temos uma previsão na escala original. Não é uma projeção muito boa, eu diria, mas você entendeu o conceito? Agora, eu deixo pra você aperfeiçoar a metodologia daqui pra frente e fazer uma solução melhor!

Notas finais

Neste artigo eu tentei te dar uma abordagem padrão para resolver o problema de séries temporais. Cobrimos conceitos de estacionaridade, como tornar uma série de tempo mais perto da estacionaridade e, finalmente, prever os resíduos. Foi uma longa viagem e eu pulei alguns detalhes estatísticos que eu te encorajo a se aprofundar usando o material sugerido. Se você não quiser copiar e colar, você pode baixar o notebook ipython com todos os códigos de meu repositório GitHub.

Espero que este artigo o ajude a alcançar uma boa primeira solução.

Gostou do artigo? Algo incomodando você que deseja discutir mais? Por favor, sinta-se livre para deixar um comentário.

Veja também:

• Guia de Acesso Rápido Python – Numpy
• Guia de Acesso rápido ao Python
• Participantes do Mercado de Pagamentos