AFTS Year

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

df = pd.read_csv("Yearly_revenue.csv")
df.head()

years = range(1975, 1975 + len(df))
df['Date']=pd.to_datetime([f'{year}-01-01' for year in years])
df.head()

df['Date'] = pd.to_datetime(df['Date'], format='%Y-%m-%d')
df.set_index('Date', inplace=True)

df.head()

result=seasonal_decompose(df[['Yearly_revenue']],model='multiplicative')
result.plot()
plt.show()

yearly_result = seasonal_decompose(df['Yearly_revenue'], model='additive')
yearly_result.plot()
plt.show()

import pymannkendall as mk
mk.original_test(df['Yearly_revenue'])

train_df=df[:int(df.shape[0]*0.7)]
test_df=df[:int(df.shape[0]*0.7):]

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(df, 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()

from statsmodels.tsa.api import Holt
model_double = Holt(train_df)
model_double_fit = model_double.fit()

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

plt.plot(df, 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()

from statsmodels.tsa.api import ExponentialSmoothing
model_triple = ExponentialSmoothing(train_df, seasonal_periods=12, trend='add', seasonal='add')
model_triple_fit =model_triple.fit()

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

plt.plot(df, 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()

from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error
mape_train = mean_absolute_percentage_error(train_df['Yearly_revenue'],model_triple_fit.fittedvalues)
print('Mape Train',mape_train)

mape_test = mean_absolute_percentage_error(test_df['Yearly_revenue'],model_triple_fit.fittedvalues)
print('Mape Test',mape_test)

# ADF Test
from statsmodels.tsa.stattools import adfuller
print("Result of Dickeu Fuller Test")
dftest=adfuller(df['Yearly_revenue'])
dfoutput=pd.Series(dftest[0:4],index=["Test Statistic","p-value",'#Lags used','Number of Observation'])
for key,value in dftest[4].items():
    dfoutput['Critical Value(%s)' %key]= value
print(dfoutput)

#KPSS Test
from statsmodels.tsa.stattools import kpss
kp=kpss(df['Yearly_revenue'])
p=kp[1]
print("p-value for Kps test(untransformed)=",p)

# acf and pacf 
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(df['Yearly_revenue'], ax=plt.gca(),lags=30)
plt.show()

plt.figure(figsize=(20,5))
plt.grid()
plot_pacf(df['Yearly_revenue'], ax=plt.gca(),lags=20)
plt.show()

# First order differencing
diff=df['Yearly_revenue'].diff().dropna()
plt.figure(figsize=(14,3))
plt.grid()
plt.plot(diff)
plt.show()

#Now check the stationarity and for that run adf, kpss, and acf plot again
from statsmodels.tsa.stattools import adfuller
print("Result of Dickeu Fuller Test")
dftest=adfuller(diff)
dfoutput=pd.Series(dftest[0:4],index=["Test Statistic","p-value",'#Lags used','Number of Observation'])
for key,value in dftest[4].items():
    dfoutput['Critical Value(%s)' %key]= value
print(dfoutput)

from statsmodels.tsa.stattools import kpss
kp=kpss(diff)
p=kp[1]
print("p-value for Kps test(untransformed)=",p)

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(diff, ax=plt.gca(),lags=30)
plt.show()

# second order diffrencing 
diff2=df['Yearly_revenue'].diff().diff().dropna()
plt.figure(figsize=(14,3))
plt.grid()
plt.plot(diff)
plt.show()

from statsmodels.tsa.stattools import adfuller
print("Result of Dickeu Fuller Test")
dftest=adfuller(diff2)
dfoutput=pd.Series(dftest[0:4],index=["Test Statistic","p-value",'#Lags used','Number of Observation'])
for key,value in dftest[4].items():
    dfoutput['Critical Value(%s)' %key]= value
print(dfoutput)

from statsmodels.tsa.stattools import kpss
kp=kpss(diff2)
p=kp[1]
print("p-value for Kps test(untransformed)=",p)

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(diff2, ax=plt.gca(),lags=30)
plt.show()

We still didn't get the data as stationary.
Since it's a yearly data and there's no seasonality present we can't go for Seasonal Decomposition hence we will continue with the non seasonal decomposition and do the 3rd order differencing.
If it still does not lead to stationarity then we will stop here.

# 3rd order non seasonal - differencing
diff3=df['Yearly_revenue'].diff().diff().diff().dropna()
plt.figure(figsize=(14,3))
plt.grid()
plt.plot(diff)
plt.show()

from statsmodels.tsa.stattools import adfuller
print("Result of Dickeu Fuller Test")
dftest=adfuller(diff3)
dfoutput=pd.Series(dftest[0:4],index=["Test Statistic","p-value",'#Lags used','Number of Observation'])
for key,value in dftest[4].items():
    dfoutput['Critical Value(%s)' %key]= value
print(dfoutput)

from statsmodels.tsa.stattools import kpss
kp=kpss(diff3)
p=kp[1]
print("p-value for Kps test(untransformed)=",p)

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(diff3, ax=plt.gca(),lags=30)
plt.show()

from statsmodels.tsa.stattools import adfuller
print("Results of Dickey-Fuller Test: ")
dftest = adfuller(df['Yearly_revenue'])
dfoutput = pd.Series

#Results of Dickey-Fuller Test: 
#Even though the adf and kpss test says that it is stationary, the graph still has 1 point very outside the range. But that can be neglected.
#We give more priority to statistical outputs rather than the graphical visualisations
#Hence since the adf and kpss test says that it is Stationary we will ASSUME that yes the data has finally become stationary on the basis of the statistical tests
#ARIMA -> When there is Seasonal component we will be doing SARIMA
#But since there is no seasonal component here we will do ARIMA

q from ACF plot (the lag that crosses the significant level the most) =1 
p from PACF plot (the lag that crosses the significant level the most) =1 
d which order of differencing gave stationarity = 3

# order = p,d,q
import statsmodels.api as sm
from sklearn.model_selection import train_test_split
train_df,test_df=train_test_split(df,test_size=0.3,shuffle=False)
model_1 = sm.tsa.statespace.SARIMAX(train_df['Yearly_revenue'],order=(1,3,1),seasonal_order=(0,0,0,0)) # no seasonality therefore 0
model_1_fit = model_1.fit()
print(model_1_fit.summary())

#AUTO ARIMA
from pmdarima.arima import auto_arima
model_auto_arima=auto_arima(df['Yearly_revenue'], seasonal=False, trace=True, m=1)
model_auto_arima_fit = model_auto_arima.fit(df['Yearly_revenue'])

# we get the best model as 0,1,0 so now we'll re reun the model1 with the updated parameters and change it's name to model 2
model_2 = sm.tsa.statespace.SARIMAX(train_df['Yearly_revenue'],order=(0,1,0),seasonal_order=(0,0,0,0)) # no seasonality therefore 0
model_2_fit = model_2.fit()
print(model_2_fit.summary())

#1st model AIC = 489.039
#2nd model AIC = 528.186
#This suggests that thw 1st model is better since it has less AIC. The reason why Auto arima didn't give a better result is because it doesn't consider AIC. It gives the result on the basis of the MAPE value
#Hence now we'll find the mape for that we need to find the test predictions

test_predictions1=model_1_fit.forecast(len(test_df))
forecast_further1=model_1_fit.forecast(len(test_df['Yearly_revenue'])+90)
forecast_further1=forecast_further1[test_df.index.max():]
forecast_further1

test_predictions2=model_2_fit.forecast(len(test_df))
forecast_further2=model_2_fit.forecast(len(test_df['Yearly_revenue'])+90)
forecast_further2=forecast_further2[test_df.index.max():]
forecast_further2

### now comparing using mean_absolute error of both models
from sklearn.metrics import mean_absolute_percentage_error
mape_test1=mean_absolute_percentage_error(test_df['Yearly_revenue'],test_predictions1)
print('Mape for Test data 1:',mape_test1)

from sklearn.metrics import mean_absolute_percentage_error
mape_test2=mean_absolute_percentage_error(test_df['Yearly_revenue'],test_predictions2)
print('Mape for Test data 2:',mape_test2)

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

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

#Visualizing SARIMA Model
# for custom model
plt.figure(figsize=(14,3))
plt.plot(train_df['Yearly_revenue'],label="Actual")
plt.plot(test_df['Yearly_revenue'],label="test")
plt.plot(test_predictions1,label="Prediction")
plt.plot(forecast_further1,label="future")
plt.title('SARIMAX model')
plt.xlabel('date')
plt.ylabel('ep')
plt.legend()
plt.show()

## for auto_arima model
plt.figure(figsize=(14,3))
plt.plot(train_df['Yearly_revenue'],label="Actual")
plt.plot(test_df['Yearly_revenue'],label="test")
plt.plot(test_predictions2,label="Prediction")
plt.plot(forecast_further2,label="future")
plt.title('SARIMAX model')
plt.xlabel('date')
plt.ylabel('ep')
plt.legend()
plt.show()


#to display plot for Mean_absolute_Error
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
models=['SARIMAX(1,3,1)x(0,0,0,0)','SARIMAX(0,1,0)x(0,0,0,0)']
mape_test=[0.54849386356058,0.19009645478424367]
table=pd.DataFrame({'Model':models,'MAPE TEST':mape_test})
plt.figure(figsize=(10,6))
sns.barplot(x='Model',y='MAPE TEST',data=table)
plt.title('MAPE OF DIFF MODEL')
plt.xlabel('MODEL')
plt.ylabel('MAPE')
plt.show()