import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from pmdarima import auto_arima
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from sklearn.metrics import mean_squared_error, mean_absolute_error

import warnings
warnings.filterwarnings('ignore')

df = pd.read_csv('./Data/EnergyProduction.csv', index_col='DATE', parse_dates=True)

df.plot();

len(df)

n_obs = 24

train_data = df.iloc[:-n_obs]
test_data = df.iloc[-n_obs:]

print(train_data.shape, test_data.shape)

plt.plot(train_data, label='train')
plt.plot(test_data, label='test')
plt.legend()
plt.show()

def exp_smoothing(train_data, trend, seasonal, test_data = test_data):
    fitted_model = ExponentialSmoothing(train_data, trend=trend, seasonal=seasonal,
                                       seasonal_periods=12).fit()
    test_predictions = fitted_model.forecast(24)
    plt.figure(figsize=(12,5))
    plt.plot(test_data, label='test')
    plt.plot(test_predictions, label='prediction')
    
    mae = mean_absolute_error(test_data, test_predictions)
    mse = mean_squared_error(test_data, test_predictions)
    rmse = np.sqrt(mse)
    
    print('mae :', mae)
    print('mse :', mse)
    print('rmse :', rmse)

exp_smoothing(train_data, 'mul', 'mul')

exp_smoothing(train_data, 'add','add')

final_model = ExponentialSmoothing(df, trend='add',seasonal='add',
                                  seasonal_periods=12).fit()

forecast_prediction = final_model.forecast(36)

plt.figure(figsize=(12,8))
plt.plot(df)
plt.plot(forecast_prediction)
plt.show()

from statsmodels.tsa.stattools import adfuller
from pmdarima import auto_arima

from statsmodels.tsa.arima_model import ARIMA, ARIMAResults
from statsmodels.tsa.statespace.sarimax import SARIMAX, SARIMAXResults
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.statespace.tools import diff
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tools.eval_measures import rmse

result = seasonal_decompose(df, model='add')

result.plot();

adfuller(df)[1]

df_diff = df.diff().dropna()

adfuller(df_diff)[1]

stepwise_fit = auto_arima(df, seasonal=True, trace=True, m=12)

stepwise_fit.summary()

model = SARIMAX(train_data, order=(3,0,1), seasonal_order=(0,1,1,12))

results = model.fit()

results.summary()

start = len(train_data)
end = len(train_data) + len(test_data) - 1

pred = results.predict(start, end, typ='level')

plt.plot(pred)
plt.plot(test_data)
plt.show()

error = rmse(test_data, pred)
error

model = SARIMAX(df, order=(3,0,1), seasonal_order=(0,1,1,12))

results = model.fit()

results.summary()

forecast = results.predict(len(df)-1, len(df)+11, typ='levels')

plt.plot(df['1983-01-01':])
plt.plot(forecast)

from tensorflow.keras.preprocessing.sequence import TimeseriesGenerator
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Dense, LSTM, Input

nobs = 12

train = df.iloc[:-nobs]
test = df.iloc[-nobs:]

mm_scaler = MinMaxScaler()
st_scaler = StandardScaler()

mm_scaler.fit(train)
st_scaler.fit(train)

mm_scaled_train = mm_scaler.transform(train)
mm_scaled_test = mm_scaler.transform(test)
st_scaled_train = st_scaler.transform(train)
st_scaled_test = st_scaler.transform(test)

n_input = 12
n_features = 1

mm_gen = TimeseriesGenerator(mm_scaled_train, mm_scaled_train, length=n_input,
                         batch_size=1)
st_gen = TimeseriesGenerator(st_scaled_train, st_scaled_train, length=n_input,
                         batch_size=1)

inputs = Input(shape=(n_input, n_features))
x = LSTM(100, activation='relu')(inputs)
x = Dense(1)(x)

mm_model = Model(inputs, x)
st_model = Model(inputs, x)

mm_model.compile(optimizer='adam', loss='mse')
st_model.compile(optimizer='adam', loss='mse')

mm_hist = mm_model.fit_generator(mm_gen, epochs=50, verbose=0)
st_hist = st_model.fit_generator(st_gen, epochs=50, verbose=0)

plt.plot(mm_hist.history['loss'], label='MinMax')
plt.plot(st_hist.history['loss'], label='Standard')
plt.legend()
plt.show()

first_eval_batch = mm_scaled_train[-12:]
first_eval_batch = first_eval_batch.reshape((1, n_input, n_features))

mm_model.predict(first_eval_batch)

test_predictions = []
first_eval_batch = mm_scaled_train[-n_input:]
current_batch = first_eval_batch.reshape((1, n_input, n_features))

for i in range(len(test)):
    current_pred = model.predict(current_batch)[0]
    test_predictions.append(current_pred)
    current_batch = np.append(current_batch[:, 1:, :], [[current_pred]], axis=1)

test_predictions

true_predictions = mm_scaler.inverse_transform(test_predictions)

test['predictions'] = true_predictions

test.plot(figsize=(12,9))

print(rmse(test['EnergyIndex'], test['predictions']))