import math
import numpy as np
import pandas as pd

# Date & Time
import datetime
import matplotlib.dates as mdates

# TSA from Statsmodels
import statsmodels.api as sm
import statsmodels.tsa.api as smt
import statsmodels.formula.api as smf

#Ploting
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline

def test_stationarity(timeseries, maxlag=None, regression=None, autolag=None,
                      window=None, plot=False, verbose=False):
    
    if regression is None:
        regression = 'c'
    
    if verbose:
        print('Running Augmented Dickey-Fuller test with paramters:')
        print('maxlag: {}'.format(maxlag))
        print('regression: {}'.format(regression))
        print('autolag: {}'.format(autolag))
    
    if plot:
        if window is None:
            window = 4
        rolmean = timeseries.rolling(window=window, center=False).mean()
        rolstd = timeseries.rolling(window=window, center=False).std()
        
        orig = plt.plot(timeseries, color='blue', label='Original')
        mean = plt.plot(rolmean, color='red', label='Rolling Mean ({})'.format(window))
        std = plt.plot(rolstd, color='black', label='Rolling Std ({})'.format(window))
        plt.legend(loc='best')
        plt.title('Rolling Mean & Standard Deviation')
        plt.show(block=False)
    
    dftest = smt.adfuller(timeseries, maxlag=maxlag, regression=regression, autolag=autolag)
    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
    if verbose:
        print('Results of Augmented Dickey-Fuller Test:')
        print(dfoutput)
    return dfoutput

"""
Explanation: Airlines Ticket Forecasting 1949-1960
Time Series Analysis Tutorial
<br>
End of explanation
"""


air_data = pd.read_csv("data/AirPassengers.csv", header = 0, index_col=0)
air_data.head(13)

"""
Explanation: Reading The Data
End of explanation
"""


n = len(air_data.index)
temp_time = pd.DataFrame({'year': [int(math.floor(air_data['time'][i])) for i in range (1, n+1)],
                          'month': [int(round((air_data['time'][i]-int(math.floor(air_data['time'][i])))*12, 0))+1 
                                    for i in range (1, n+1)],
                          'day': [1]*n}, index=air_data.index)

air_data['time'] = pd.to_datetime(temp_time)
air_data.head(13)

"""
Explanation: Data Preprocessing:
the formating of the need to be transformed into dates
End of explanation
"""


x_axis = pd.DataFrame([mdates.date2num(air_data['time'][i]) for i in range(1,n+1)])
# [mdates.date2num(air_data['time'][i]) for i in range(1,n)]

fig, ax = plt.subplots(figsize=(9,5));
ax.set_title('International airline passengers, 1949-1960');
sns.tsplot(data=air_data['AirPassengers'], time=x_axis, ax=ax)
ax.xaxis.set_major_locator(mdates.AutoDateLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
fig.autofmt_xdate()
plt.show()

"""
Explanation: Exploring the Time Series:
To determine its componenets: Error, Trend and Seasonality (ETS)
End of explanation
"""


moving_avg = air_data['AirPassengers'].rolling(window=12,center=False).mean()

fig, ax = plt.subplots(figsize=(8,4));
ax.set_title('Series Trend Component');
moving_avg.plot()
plt.show()

"""
Explanation: ETS Approch:
1- Check the components of the series (Trend and Seasonality).<br>
2- Stationarize the series.<br>
3- Apply the appropriate ETS model.<br>
End of explanation
"""


air_data['lnAirPassengers'] = np.log(air_data['AirPassengers'])
air_data['firstDiff'] = air_data['lnAirPassengers'].diff()
air_data['secondDiff'] = air_data['firstDiff'].diff()

fig = plt.figure(figsize=(15, 5))
plt.subplot(121)
plt.title("First Difference - Standard Dev. = {}".format(round(np.std(air_data['firstDiff']),2)));
pd.Series(air_data['firstDiff']).plot()
plt.subplot(122)
plt.title("Second Difference - Standard Dev. = {}".format(round(np.std(air_data['secondDiff']),2)));
pd.Series(air_data['secondDiff']).plot()
plt.show()

"""
Explanation: From the 2 figures above, one can conclude the the series has an additive trend with a multiplicative seasonality. Hence, Seasonal Holts winter method can be utilized.
Box-Jenkins Approach:
1- Stationarize the Series:


There are several diffrenceing methods can be used to stationarize the series such as: Deflation by CPI, Logarithmic Function, Difference, Seasonal Difference and Seasonal Adjustment.


Usually Dickey-Fuller test is used to assess the stationarity of the series and compare the previous methods (the test statistic must be greater than the critical value).


However I used the first and second Seasonal Difference and choosed the one with the minimum variance.
End of explanation
"""


print test_stationarity(air_data['firstDiff'].dropna())
print "------------------------------------------"
print test_stationarity(air_data['secondDiff'].dropna())

"""
Explanation: Testing Stationarity:
End of explanation
"""


fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(121)
fig = smt.graphics.plot_acf(air_data['secondDiff'].iloc[2:], lags=40, ax=ax1)
ax2 = fig.add_subplot(122)
fig = smt.graphics.plot_pacf(air_data['secondDiff'].iloc[2:], lags=40, ax=ax2)

"""
Explanation: Exploring the ACF and PACF:
To determine the SARIMA components
End of explanation
"""


# air_data = air_data.set_index('time')
mod = sm.tsa.statespace.SARIMAX(air_data['AirPassengers'][1:120], order=(0,0,0), seasonal_order=(3,2,0,12))
results = mod.fit()
print results.summary()

results.plot_diagnostics(figsize=(16, 12));

"""
Explanation: Building & Selecting the SARIMA Model
Note That, I kept the last two years out of the trainnig process for the testing phase
End of explanation
"""


pred_begin = air_data.index[-35]
pred_end = air_data.index[-1]
air_data['forecast'] = results.get_prediction(start=pred_begin.strftime('%Y-%m-%d'),
                                              end=pred_end.strftime('%Y-%m-%d')).predicted_mean


air_data[['AirPassengers', 'forecast']].plot(figsize=(10, 5)).set_title(
          'Forecasting International airline passengers, 1949-1960');

"""
Explanation: Model Selection
Future Prediction
End of explanation
"""


