# AFTS Prac

```python
AIR PASSENGERS
# Import packages
import pandas as pd
import numpy as np
import os 
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose 
import seaborn as sns

time_series_data = pd.read_csv("AirPassengers.csv")
time_series_data.head()

time_series_data.shape

# Convert column year in datetime format
time_series_data['Month'] = pd.to_datetime(time_series_data['Month'])
time_series_data = time_series_data.set_index('Month')

time_series_data.head()

sns.lineplot(time_series_data)
plt.ylabel("#passenger")

plt.figure(figsize=(10, 6))
plt.plot(time_series_data.index , time_series_data['#Passengers'])
plt.xlabel('Month')
plt.ylabel('#passenger')
plt.show()

#We can see that the number of passengers appears to be increasing over time, with some seasonality as well.
# Decomposition of time series - Multiplicative model
result = seasonal_decompose(time_series_data[['#Passengers']], model='multiplicative')
result.plot()
plt.show()

result = seasonal_decompose(time_series_data[['#Passengers']], model='additive', period = 10)
result.plot()
plt.show()

import pymannkendall as mk
# Perform Mann-Kendall test
# H0: There is no monotonic trend in the series
mk.original_test(time_series_data['#Passengers'])

#pvalue <0.5 => reject Ho. There is a trend present
train_df=time_series_data[:int(time_series_data.shape[0]*0.7)]
test_df=time_series_data[int(time_series_data.shape[0]*0.7):]

train_df.head(5)

train_df.shape

test_df.shape

#SINGLE EXPONENTIAL SMOOTHNING
from statsmodels.tsa.api import SimpleExpSmoothing
model = SimpleExpSmoothing(train_df)
model_single_fit = model.fit()

forecast_single = model_single_fit.forecast(6)
print(forecast_single)

forecast_single = model_single_fit.forecast(40)

model_single_fit.params

plt.plot(time_series_data, label='Original Data')
plt.plot(model_single_fit.fittedvalues, label='Fitted Values')
plt.plot(forecast_single, label='Forecast')
plt.xlabel('Year')
plt.ylabel('Number of Passengers')
plt.title('Single Exponential Smoothing')
plt.legend()
plt.show()

#DOUBLE EXPONENTIAL SMOOTHNING
from statsmodels.tsa.api import Holt
model_double = Holt(train_df)
model_double_fit = model_double.fit()

forecast_double = model_double_fit.forecast(6)
print(forecast_double)

forecast_double = model_double_fit.forecast(40)

model_double_fit.params

plt.plot(time_series_data, label='Original Data')
plt.plot(model_double_fit.fittedvalues, label='Fitted Values')
plt.plot(forecast_double, label='Forecast')
plt.xlabel('Year')
plt.ylabel('Number of Passengers')
plt.title('Double Exponential Smoothing')
plt.legend()
plt.show()

#TRIPLE EXPONENTIAL SMOOTHNING

from statsmodels.tsa.api import ExponentialSmoothing
# Create an instance of ExponentialSmoothing class
model_triple = ExponentialSmoothing(train_df, seasonal_periods=12, trend='add', seasonal='add')
# Fit the model to the data
model_triple_fit = model_triple.fit()

forecast_triple = model_triple_fit.forecast(6)
print(forecast_triple)

forecast_triple = model_triple_fit.forecast(40)

model_triple_fit.params

plt.plot(time_series_data, label='Original Data')
plt.plot(model_triple_fit.fittedvalues, label='Fitted Values')
plt.plot(forecast_triple, label='Forecast')
plt.xlabel('Year')
plt.ylabel('Number of Passengers')
plt.title('Triple Exponential Smoothing')
plt.legend()
plt.show()

residuals = model_triple_fit.resid
residuals.head()

test_predictions = model_triple_fit.forecast(len(test_df)).rename('HW Test Forecast')

test_predictions.tail()

plt.figure(figsize=(16,8))
plt.plot(train_df['#Passengers'],label='Train')
plt.plot(test_df['#Passengers'], label = 'Test')
plt.plot(model_triple_fit.fittedvalues , label='fitted')
plt.plot(test_predictions , label ='forecast')
plt.xlabel('Year')
plt.ylabel('Number of passengers')
plt.legend()
plt.show()

from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error
mape_train = mean_absolute_percentage_error(train_df['#Passengers'], model_triple_fit.fittedvalues)
print('MAPE for Train Data:', mape_train)

mape_test = mean_absolute_percentage_error(test_df['#Passengers'],test_predictions)
print('MAPE for Test Data:', mape_test)

#ADF TEST
#Null Hypothesis (HO): Series is non-stationary, or series has a unit root.
#Alternate Hypothesis(HA): Series is stationary, or series has no unit root.
#p-value < 0.05 – Reject Null Hypothesis(HO)
from statsmodels.tsa.stattools import adfuller
print("Results of Dickey-Fuller Test:")
dftest = adfuller(time_series_data['#Passengers'])
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)
#Here the condition is not satisfied so we fail to reject the null hypothesis i.e Series is NOT stationary

#KPSS TEST
#Null Hypothesis (HO): Series is trend stationary or series has no unit root.
#Alternate Hypothesis(HA): Series is non-stationary, or series has a unit root.
#p-value < 0.05 –  Reject Null Hypothesis(HO)
from statsmodels.tsa.stattools import kpss
kp = kpss(time_series_data['#Passengers'])
p = kp[1]
print("p-value for KPSS test (untransformed) = ", p)
#Here p-value is less than 0.05 so reject Null Hypothesis i.e data is NOT trend stationary

from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.figure(figsize=(20, 5))
plt.grid()
plot_acf(time_series_data['#Passengers'] , ax=plt.gca(), lags =30 )
plt.show()
#ACF decays slowly, a sinusoidal or linear decay, or remains high for many lags --> NOT STATIONARY.
#Values close to 1 or -1 at a lag suggest strong positive or negative correlation.
#Values near 0 indicate no correlation.

plt.figure(figsize=(20, 5))
plt.grid()
plot_pacf(time_series_data['#Passengers'] , ax=plt.gca(), lags = 30)
plt.show()

#PACF may show significant correlations across multiple lags
#With the help of ACF and PACF we can see that values of lag are falling outside the significant region so we can conclude that our data is non-stationary

#NON-SEASONAL DIFFERENCING
diff = time_series_data['#Passengers'].diff().dropna()     # TRYING TO MAKE DATA TREND STATIONARY
plt.figure(figsize=(14, 3))
plt.grid()
plt.plot(diff)
plt.show()

#SEASONAL DIFFERENCING
sdiff = time_series_data['#Passengers'].diff(periods = 12).dropna()          # TRYING TO REMOVE SEASONALITY
plt.figure(figsize=(14, 3))
plt.grid()
plt.plot(sdiff)
plt.show()

#SEASONAL AND NON-SEASONAL DIFFERENCING
sdiff = time_series_data['#Passengers'].diff(periods = 12).dropna()  # seasonal differencing
sddiff = sdiff.diff().dropna()                          # non-seasonal differencing  # 1ST ORDER DIFFERENCING
plt.figure(figsize=(14, 3))
plt.grid()
plt.plot(sddiff)
plt.show()

#ADF
#Null Hypothesis (HO): Series is non-stationary, or series has a unit root.
#Alternate Hypothesis(HA): Series is stationary, or series has no unit root.
#p-value < 0.05 – Reject Null Hypothesis(HO)

print("Results of Dickey-Fuller Test:")
dftest = adfuller(sddiff)
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)
## SERIES IS STATIONARY

#KPSS TEST
#Null Hypothesis (HO): Series is trend stationary or series has no unit root.
#Alternate Hypothesis(HA): Series is non-stationary, or series has a unit root.
#p-value < 0.05 –  Reject Null Hypothesis(HO)
from statsmodels.tsa.stattools import kpss
kp = kpss(sddiff)
p = kp[1]
print("p-value for KPSS test (untransformed) = ", p)
## Data is trend stationary

plt.figure(figsize=(20, 5))
plt.grid()
plot_acf(sddiff , ax=plt.gca(), lags = 30)
plt.show()

plt.figure(figsize=(20, 5))
plt.grid()
plot_pacf(sddiff , ax=plt.gca(), lags = 30)
plt.show()
#The ACF and PACF graphs also give the same result that the data is stationory after first order differencing

from pmdarima.arima import auto_arima
model_auto_arima = auto_arima(time_series_data['#Passengers'], seasonal=True, trace = True , m=12)
model_auto_arima_fit = model_auto_arima.fit(time_series_data['#Passengers'])
print(model_auto_arima_fit.summary())

#INTERPRETATION ON ACF AND PACF
import statsmodels.api as sm
model_autosarima_dom1 = sm.tsa.statespace.SARIMAX(train_df['#Passengers'], order=(1,1,1), seasonal_order = (1,1,1,12))
model_autosarima_dom1_fit = model_autosarima_dom1.fit()
print(model_autosarima_dom1_fit.summary())

model_autosarima_dom1_fit.plot_diagnostics(figsize=(15,12))
plt.show()

test_predictions1 = model_autosarima_dom1_fit.forecast(len(test_df))
forecast_further1=model_autosarima_dom1_fit.forecast(len(test_df['#Passengers'])+90)
forecast_further1 = forecast_further1[test_df.index.max():]
forecast_further1

plt.figure(figsize=(12, 8))
plt.plot(train_df['#Passengers'], label='Actual')
plt.plot(test_df['#Passengers'], label='Test')
plt.plot(test_predictions1, label='Prediction')
plt.plot(forecast_further1, label='Future')
plt.title('SARIMAX Model Forecast')
plt.xlabel('Year')
plt.ylabel('#Passengers')
plt.legend()
plt.show()

mape_train = mean_absolute_percentage_error(train_df['#Passengers'], model_autosarima_dom1_fit.fittedvalues)
print('MAPE for Train Data:', mape_train)

mape_test = mean_absolute_percentage_error(test_df['#Passengers'],test_predictions1)
print('MAPE for TEST Data:', mape_test)

#BASED ON AUTO ARIMA BEST
model_autosarima_dom = sm.tsa.statespace.SARIMAX(time_series_data['#Passengers'], order=(2,1,1), seasonal_order = (0,1,0,12))
model_autosarima_dom_fit = model_autosarima_dom.fit()
print(model_autosarima_dom_fit.summary())

model_autosarima_dom_fit.plot_diagnostics(figsize=(15,12))
plt.show()

test_predictions = model_autosarima_dom_fit.forecast(len(test_df))
forecast_further=model_autosarima_dom_fit.forecast(len(test_df['#Passengers'])+90)
forecast_further = forecast_further[test_df.index.max():]
forecast_further

plt.figure(figsize=(12, 8))
plt.plot(train_df['#Passengers'], label='Actual')
plt.plot(test_df['#Passengers'], label='Test')
plt.plot(test_predictions, label='Prediction')
plt.plot(forecast_further, label='Future')
plt.title('SARIMAX Model Forecast')
plt.xlabel('Year')
plt.ylabel('#Passengers')
plt.legend()
plt.show()

from sklearn.metrics import mean_absolute_percentage_error
mape_test = mean_absolute_percentage_error(test_df['#Passengers'],test_predictions)
print('MAPE for TEST Data:', mape_test)

# prompt: plot the mape of all models make a table
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Create some data for the table
models = ['Exponential Smoothing', 'SARIMAX (1,1,1)x(1,1,1,12)', 'SARIMAX (2,1,1)x(0,1,0,12)']
mape_test = [0.05397259855132907, 0.05292896562888002, 0.3820631602301439]
# Create the table
table = pd.DataFrame({'Model': models, 'MAPE Test': mape_test})
# Plot the table
plt.figure(figsize=(10, 6))
sns.barplot(x='Model', y='MAPE Test', data=table)
plt.xlabel('Model')
plt.ylabel('MAPE')
plt.title('MAPE of Different Models')
plt.show()
```