Forecasting.py

# -*- coding: utf-8 -*-
"""Aula 2 - Forecasting

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1qMk3Ro_I_OAdWWf6FZaRT95JYlQgFvS0

Quanto maior no futuro que querer ver algo, menor será a certeza, ou seja, a minha capacidade de inferir sobre algo

Sempre que eu queira prever algo, algum erro terá

Para isso, há técnicas e métodos de estatística mais avançada
"""

# Commented out IPython magic to ensure Python compatibility.
#importando as bibliotecas
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import mean_squared_error

sns.set_style()

#Deixa o gráfico melhor, no formato svg
# %config InlineBackend.figure_format = 'svg'

#URL
dataset_path = "https://raw.githubusercontent.com/carlosfab/escola-data-science/master/datasets/electricity_consumption/Electric_Production.csv"

#Lendo o URL
df = pd.read_csv(dataset_path)

#Importar o csv para um dataframe
df.DATE = pd.to_datetime(df.DATE, format="%m-%d-%Y")
df.set_index('DATE', inplace=True)

df.head()

#Dividir um grupo de treino e validação (as datas foram arbitrárias)

train = df[df.index <= '2012-8-1']
valid = df[df.index > '2012-8-1']

#Aqui iremos criar um dataframe para armazenar as previsões
y_hat = valid.copy()

#Utilizando períodos amplos para que a aplicação dos métodos seja perceptível visualmente
y_hat['Naive'] = train.iloc[-1].values[0]

#Plotar train e valid
fig, ax = plt.subplots(figsize=(8,4))
train.plot(ax=ax)
valid.plot(ax=ax)

y_hat['Naive'].plot(ax=ax)

'''O gráfico pegou todo o conjunto de dados e separou até 2014. 
Em laranja mostra como se fosse os dias atuais
O método Naive pegou o último valor e replicou o valor para todas
as previsões'''

# Erro do método Naive
print('Erro do método Naive')

# (valor real, valor que eu previ, tira o quadrado)
mean_squared_error(y_hat.Value, y_hat.Naive, squared=True)

'''Média móvel
Ao contrário do método Naive que pega somente
o último valor, a média móvel permite usar uma janela de intervalo.

Elas são boas para suavisar as curvas, siminuindo dispersões/ruídos
para criar novas variáveis - feature engineering'''

# Calcular a média dos últimos 7 valores disponíveis
# Conjunto de treino

# A função .rolling - do pandas - permite pegar a média do período entre ()
y_hat['m7'] = train.Value.rolling(7).mean().iloc[-1]

# Nas 7 primeiras entradas aparece NaN, pois não tem nada para comparar

# Plot
fig, ax = plt.subplots(figsize=(8,4))
train.plot(ax=ax)
valid.plot(ax=ax)
y_hat['m7'].plot(ax=ax)
plt.show()

''' O gráfico não pegou mais o último valor, agora ele pegou os 7
períodos anteriores'''

# Cálculo do erro
mean_squared_error(y_hat.Value, y_hat.m7, squared=True)

''' Agora trabalharemos com uma técnica que também lida com as tendências.
O nome é Holt's Linear TRend Model'''

# Importando as bibliotecas necessárias
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.holtwinters import Holt

# Salvando os dados de treino em uma variável

result = seasonal_decompose(train)

# Plotando as componentes
fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(8,10))
result.observed.plot(ax=ax1)
result.trend.plot(ax=ax2)
result.seasonal.plot(ax=ax3)
result.resid.plot(ax=ax4)

# Serve para ajustar os gráficos e não deixá-los desalinhados
plt.tight_layout()

'''Salvando os valores para criar uma reta que se encaixe na inclinação do gráfico 1
para não ficar simplesmente uma reta horizontal'''

# Smoothing level trabalha um ajuste de altura enquanto o slope é a inclinação da reta
y_hat['holt'] = Holt(train.Value).fit(smoothing_level=0.1,
                                      smoothing_slope=0.1).forecast(len(valid))

# Plot train e valid
fig, ax = plt.subplots(figsize=(8,4))
train.plot(ax=ax)
valid.plot(ax=ax)
y_hat['holt'].plot(ax=ax)
ax.legend(['Train', 'Valid', 'Holt'])
plt.show()

# Calcular o erro para ver se a altura e slope estão mais adequados
print('Erro do mean squared error: {}'.format(mean_squared_error(y_hat.Value, y_hat.holt, squared=True)))

''' Aparentemente o erro reduziu... quando se altera o smoothing level
e smoothing slope o erro aumenta ou diminui'''

"""Para saber se uma Série Temporal é uma série estacionária, basta verificar se as propriedades são constantes. Conceitos básicos:

media e variância são constantes em relação ao tempo
covariância entre os termos Ti e Ti+m também são constantes

Para aplicar um modelo preditivo, devemos saber se a série é estacionária. Para verificar isso, usamos o teste de Dickey-Fuller

Vamos entender um pouco sobre o que está sendo verificado.

A Hipótese Nula ($H_0$) do teste é que a TS não é estacionária. Ou seja, possui algum tipo de dependência em relação ao tempo.
A Hipótese Alternativa ($H_1$) rejeita a hipótese nula, ou seja, que a TS é estacionária.
Interpretamos o resultado do teste usando o valor-p com um threshold que traga indícios suficientes para rejeitarmos a hipótese nula. Caso o valor do teste fique acima desse threshold, falharemos em rejeitar a hipótese nula e manteremos a premissa de que ela não é estacionária.

Os valores do threshold que irei usar são:

$\text{Valor-p } \leq 0.05$: Rejeitamos $H_0$ e a TS é estacionária; e
$\text{Valor-p } >= 0.05$: Falhamos em rejeitar $H_0$ e a TS é não-estacionária.
"""

from statsmodels.tsa.stattools import adfuller

# Importando os dados como sendo uma coluna única
df_est = pd.read_csv(dataset_path, index_col=0, squeeze=True)

df_est.head()

# Extrair somente os valores
x = df_est.values
print(x)

# Aplicar o ADP e imprimir
# Se p-Valor > 5% rejeita-se a hipótese de que a série era estacionária
result = adfuller(x)
print('Dickey-Fuller Aumentado')
print('Teste Estatístico {:.4f}'.format(result[0]))
print('Valor-p: {:.4f}'.format(result[1]))
print('Valores Críticos: ')

for key, value in result[4].items():
  print("\t{}: {:.4f}".format(key, value))