<jupyter_start><jupyter_text>Sistemas Híbridos Atividade referente ao assunto Sistemas Híbridos. O objetivo dessa atividade foi implementar um Sistema Híbirdo, conforme Zhang (2003), para realizar previsão do logaritmo da taxa de mortalidade da idade de 40 anos da população francesa. A metodologia consiste em realizar a modelagem da série $Z_t$ a partir da equação $Z_t=L_t + N_t$.  Para tanto, utilizou-se um Arima(p,d,q) para modelar $L_t$, enquanto que a estrutura não linear, $N_t$, foi modelada pelas Redes Neurais Artificiais, especificamente por uma Multilayer Perceptron (MLP).Metodologia:1. Treinar um ARIMA2. Calcular o erro do ARIMA nos dados de treinamento3. Treinar uma MLP no erro do ARIMA4. Realizar a Previsão do ARIMA: h = 205. Realizar a Previsão por MLP com 20 saídasO trabalho está dividido em: * Funções implementadas* Carregamento das bibliotecas* Importação dos dados* Análise gráfica e estatística da série temporal* ARIMA* MLP* Aplicação do Sistema Híbrido - ARIMA + MLP  Funções implementadas<jupyter_code>def normalizar_serie(serie):
    minimo = min(serie)
    maximo = max(serie)
    y = (serie - minimo) / (maximo - minimo)
    return y
def desnormalizar(serie_atual, serie_real):
    minimo = min(serie_real)
    maximo = max(serie_real)
    
    serie = (serie_atual * (maximo - minimo)) + minimo
    
    return pd.DataFrame(serie)
def gerar_janelas(tam_janela, serie):
    # serie: vetor do tipo numpy ou lista
    tam_serie = len(serie)
    tam_janela = tam_janela +1 # Adicionado mais um ponto para retornar o target na janela
    
    janela = list(serie[0:0+tam_janela]) #primeira janela p criar o objeto np
    janelas_np = np.array(np.transpose(janela))    
       
    for i in range(1, tam_serie-tam_janela):
        janela = list(serie[i:i+tam_janela])
        j_np = np.array(np.transpose(janela))        
        
        janelas_np = np.vstack((janelas_np, j_np))
            
    return janelas_np
def diferenciar_serie(serie):
    #serie: obj pandas
    return serie.diff().dropna()
def incrementar_serie(serie_real, serie_diff):
    return serie_real[0:-1] + serie_diff
def select_lag_acf(serie, max_lag):
    from statsmodels.tsa.stattools import acf
    x = serie[0: max_lag+1]
    
    acf_x, confint = acf(serie, nlags=max_lag, alpha=.05, fft=False,
                             unbiased=False)
       
    limiar_superior = confint[:, 1] - acf_x
    limiar_inferior = confint[:, 0] - acf_x
    
    lags_selecionados = []
    
    for i in range(1, max_lag+1):

        
        if acf_x[i] >= limiar_superior[i] or acf_x[i] <= limiar_inferior[i]:
            lags_selecionados.append(i-1)  #-1 por conta que o lag 1 em python é o 0
    
    #caso nenhum lag seja selecionado, essa atividade de seleção para o gridsearch encontrar a melhor combinação de lags
    if len(lags_selecionados)==0:


        print('NENHUM LAG POR ACF')
        lags_selecionados = [i for i in range(max_lag)]

    print('LAGS', lags_selecionados)

    #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    #inverte o valor dos lags para usar na lista de dados se os dados forem de ordem [t t+1 t+2 t+3]
    lags_selecionados = [max_lag - (i+1) for i in lags_selecionados]
    #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    return lags_selecionados
def split_serie_with_lags(serie, perc_train, perc_val = 0):
    
    #faz corte na serie com as janelas já formadas 
    
    x_date = serie[:, 0:-1]
    y_date = serie[:, -1]        
       
    train_size = np.fix(len(serie) *perc_train)
    train_size = train_size.astype(int)
    
    if perc_val > 0:        
        val_size = np.fix(len(serie) *perc_val).astype(int)
              
        
        x_train = x_date[0:train_size,:]
        y_train = y_date[0:train_size]
        print("Particao de Treinamento:", 0, train_size  )
        
        x_val = x_date[train_size:train_size+val_size,:]
        y_val = y_date[train_size:train_size+val_size]
        
        print("Particao de Validacao:",train_size, train_size+val_size)
        
        x_test = x_date[(train_size+val_size):-1,:]
        y_test = y_date[(train_size+val_size):-1]
        
        print("Particao de Teste:", train_size+val_size, len(y_date))
        
        return x_train, y_train, x_test, y_test, x_val, y_val
        
    else:
        
        x_train = x_date[0:train_size,:]
        y_train = y_date[0:train_size]

        x_test = x_date[train_size:-1,:]
        y_test = y_date[train_size:-1]

        return x_train, y_train, x_test, y_test<jupyter_output><empty_output><jupyter_text>Carregando biblioteca<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.api import SimpleExpSmoothing, Holt, ExponentialSmoothing
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.stattools import adfuller
from sklearn.metrics import mean_squared_error
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import mean_squared_error as MSE
from pmdarima import auto_arima
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] =15, 6
import os
os.chdir("/Users/filipeduarte/Documents/Mortality_data")
os.getcwd()<jupyter_output><empty_output><jupyter_text>Importação dos dados<jupyter_code># importing data
lnmx_series = pd.read_csv('lnmx_series.csv', sep = ',', index_col = 0)
lnmx_series.head()
import os
os.chdir("/Users/filipeduarte/Desktop/TimeSeriesCin-master")
os.getcwd()<jupyter_output><empty_output><jupyter_text>Análise gráfica e estatística<jupyter_code>idades = lnmx_series.columns
idades
# Começar com a idade 40 anos
lnmx = lnmx_series['40']
# gráfico 
lnmx.plot(label='log(m40)')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Verifica-se um comportamento muito distinto a partir do ano de 1950. Sendo assim, optou-se por utilizar a série a partir dessa data até 2010.<jupyter_code># Antes de aplicar o Arima, vamos dividir os dados em treinamento e teste. 
# Como vou aplicar uma janela com 2 lags passados e 20 anos para previsão, vamos ficar com os anos 1816:1990 para treinamento
# e 1991:2010 para teste
lnmx = lnmx[lnmx.index<2011]
treino, teste = lnmx[lnmx.index <= 1990], lnmx[lnmx.index > 1990]<jupyter_output><empty_output><jupyter_text>* Treinamento: 1816-1990* Teste: 1991-2010<jupyter_code># Visualizar o gráfico 
treino.plot(label = 'Treinamento')
teste.plot(label='Teste')
plt.legend()
plt.show()
# calcular as estatísticas
treino.describe()
teste.describe()
# Analisar a estacionariedade da série de treinamento
result = adfuller(treino)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('Critical Values:')
for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))<jupyter_output>ADF Statistic: 0.180194
p-value: 0.971129
Critical Values:
	1%: -3.469
	5%: -2.879
	10%: -2.576<jupyter_text>Resultados indicam que a série é não estacionária. Dessa maneira, temos que realizar a diferenciação pelo modelo Arima.Próximo passo será analisar o autocorrelograma/parcial da série de treinamento diferenciada.<jupyter_code>plot_acf(treino.diff().dropna(), lags=30)
plot_pacf(treino.diff().dropna(), lags = 30)
plt.show()<jupyter_output><empty_output><jupyter_text>Autocorrelação de ordem 2  ARIMAUtilizou-se o auto_arima da biblioteca pyramid.arima para modelar a série de treinamento. O auto_arima busca o modelo com o menor AIC, a partir de uma busca iniciando com os parâmetros pelos modelos mais parcimoniosos.<jupyter_code>stepwise_model = auto_arima(treino, start_p=1, start_q=1,
                           max_p=10, max_q=10, m=12,
                           start_P=0, seasonal=False,
                           d=1, D=1, trace=True,
                           error_action='ignore',  
                           suppress_warnings=True, 
                           stepwise=True)
# melhor AIC - modelo (1,1,1)
stepwise_model.aic()<jupyter_output><empty_output><jupyter_text>O modelo selecionado será o ARIMA(p=1,q=1,d=1), ou seja, diferenciado em 1ª ordem e ARMA(1).<jupyter_code># treinamento do arima(1,1,1)
stepwise_model.fit(treino)
# extração do fit e plotando o gráfico
arima_fit_treino = stepwise_model.predict_in_sample()
#arima_fit_treino = incrementar_serie(serie_diff=arima_fit_treino, serie_real=treino)
arima_fit_treino = pd.Series(arima_fit_treino, index = treino[1:].index)
treino.plot(label='Treinamento - Real')
arima_fit_treino.plot(label='Treinamento - Fitted arima(0,1,1)')
plt.legend()
plt.show()
# extração dos resíduos e plotando o gráfico
resid_arima = stepwise_model.resid()
resid_arima = pd.Series(resid_arima, index = treino[1:].index)
resid_arima.plot(label='resíduos arima')
plt.legend()
plt.show()
# Análise do autocorrelograma dos resíduos
plot_acf(resid_arima, lags=30)
plot_pacf(resid_arima, lags=30)
plt.show()
from pandas import DataFrame
from pandas import concat
 
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
    """
    Frame a time series as a supervised learning dataset.
    Arguments:
        data: Sequence of observations as a list or NumPy array.
        n_in: Number of lag observations as input (X).
        n_out: Number of observations as output (y).
        dropnan: Boolean whether or not to drop rows with NaN values.
    Returns:
        Pandas DataFrame of series framed for supervised learning.
    """
    n_vars = 1 if type(data) is list else data.shape[1]
    df = DataFrame(data)
    cols, names = list(), list()
    # input sequence (t-n, ... t-1)
    for i in range(n_in, 0, -1):
        cols.append(df.shift(i))
        names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
    # forecast sequence (t, t+1, ... t+n)
    for i in range(0, n_out):
        cols.append(df.shift(-i))
        if i == 0:
            names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
        else:
            names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
    # put it all together
    agg = concat(cols, axis=1)
    agg.columns = names
    # drop rows with NaN values
    if dropnan:
        agg.dropna(inplace=True)
    return agg
# Criar as janelas com 20 saídas para y e 2 lags
values = list(resid_arima)
data = series_to_supervised(values, 2, 20)
data.head()<jupyter_output><empty_output><jupyter_text>MLP<jupyter_code>def treinar_mlp(x_train, y_train, x_val, y_val, num_exec):
    
    
    neuronios =  [1, 2, 3, 5, 10]    #[1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 150, 170, 200]
    func_activation =  ['tanh', 'relu']   #['identity', 'tanh', 'relu']
    alg_treinamento = ['lbfgs', 'sgd', 'adam']#, ['lbfgs','sgd', 'adam']
    max_iteracoes = [10000] #[100, 1000, 10000]
    learning_rate = ['adaptive']  #['constant', 'invscaling', 'adaptive']
    qtd_lags_sel = len(x_train[0])
    best_result = np.Inf
    for i in range(0,len(neuronios)):
        for j in range(0,len(func_activation)):
            for l in range(0,len(alg_treinamento)):
                for m in range(0,len(max_iteracoes)):
                    for n in range(0,len(learning_rate)):
                        for qtd_lag in range(1, len(x_train[0]+1)): #variar a qtd de pontos utilizados na janela 
                            
                            print('QTD de Lags:', qtd_lag, 'Qtd de Neuronios' ,neuronios[i], 'Func. Act', func_activation[j])


                            for e in range(0,num_exec):
                                mlp = MLPRegressor(hidden_layer_sizes=neuronios[i], activation=func_activation[j], solver=alg_treinamento[l], max_iter = max_iteracoes[m], learning_rate= learning_rate[n])
                                
                                
                                mlp.fit(x_train[:,-qtd_lag:], y_train)
                                predict_validation = mlp.predict(x_val[:,-qtd_lag:])
                                rmse = np.sqrt(MSE(y_val, predict_validation))

                                if rmse < best_result:
                                    best_result = rmse
                                    print('Melhor RMSE:', best_result)
                                    select_model = mlp
                                    qtd_lags_sel = qtd_lag
    
    
    return select_model, qtd_lags_sel
# Transformando em array
X, y = np.array(data.iloc[:,:2]), np.array(data.iloc[:,2:])
# Criando os conjuntos de treinamento, validação e teste
def divisao_dados_temporais(X,y, perc_treino, perc_val = 0):
    tam_treino = int(perc_treino * len(y))
    
    if perc_val > 0:        
        tam_val = int(len(y)*perc_val)
              
        X_treino = X[0:tam_treino,:]
        y_treino = y[0:tam_treino,:]
        
        print("Particao de Treinamento:", 0, tam_treino)
        
        X_val = X[tam_treino:tam_treino+tam_val,:]
        y_val = y[tam_treino:tam_treino+tam_val,:]
        
        print("Particao de Validacao:",tam_treino,tam_treino+tam_val)
        
        X_teste = X[(tam_treino+tam_val):-1,:]
        y_teste = y[(tam_treino+tam_val):-1,:]
        
        print("Particao de Teste:", tam_treino+tam_val, len(y))
        
        return X_treino, y_treino, X_teste, y_teste, X_val, y_val
        
    else:
        
        X_treino = X[0:tam_treino,:]
        y_treino = y[0:tam_treino,:]

        X_teste = X[tam_treino:-1,:]
        y_teste = y[tam_treino:-1,:]

        return X_treino, y_treino, X_teste, y_teste 
X_treino, y_treino, X_val, y_val, X_teste, y_teste = divisao_dados_temporais(X, y, perc_treino = 0.6, perc_val = 0.28)
# treinar MLP
mlp_model, lag_sel = treinar_mlp(X_treino, y_treino, X_val, y_val, num_exec = 5)
mlp_model
lag_sel<jupyter_output><empty_output><jupyter_text>Aplicação do Sistema Híbrido - ARIMA + MLP<jupyter_code># Realizando as previsões
previsao_mlp = mlp_model.predict(resid_arima.values[-1]).reshape(20,)

# Transformando em Série temporal 
previsao_mlp = pd.Series(previsao_mlp, index = teste.index)

# Prevendo com o Arima
arima_teste = stepwise_model.predict(n_periods=20)
arima_teste = pd.Series(arima_teste, index = teste.index)
# Criando o sistema híbrido z = arima + mlp
z = arima_teste+previsao_mlp
# plotando
plt.plot(teste, label = 'Série', color = 'blue')
plt.plot(arima_teste,label = 'ARIMA', color = 'coral')
plt.plot(z,label = 'Z = ARIMA + MLP', color = 'green')

plt.legend()
plt.show()
print('RMSE - ARIMA = ',np.sqrt(MSE(teste, arima_teste)))
print('RMSE - Z     = ', np.sqrt(MSE(teste, z)))<jupyter_output>RMSE - ARIMA =  0.136612432187195
RMSE - Z     =  0.14597662518655155
{'fraction_non_alphanumeric': 0.120716777238219, 'fraction_numerical': 0.02277545995485945, 'mean_word_length': 3.1932893604817894, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '185076', 'license': 'MIT', 'n_tokens_mistral': 5746, 'n_tokens_neox': 5428, 'n_words': 1259, 'path': 'sistemas_hibridos_sem_teste.ipynb', 'repo_name': 'ggcarvalho/hybridSystems', 'size': '14621'}