Afts Elec

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("Electric_Production.csv")
df.head()

df.info()

df['DATE']=pd.to_datetime(df['DATE'])
df = df.set_index('DATE')

df.head()

#H0 : Time Series is not stationary (it has time dependent structure and does not have constant variance over time)
#HA : Time Series is stationary

from statsmodels.tsa.stattools import adfuller, kpss

# Perform the ADF test
adf_result = adfuller(df['IPG2211A2N'])
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
print(f"Critical Values: {adf_result[4]}")
if adf_result[1] < 0.05:
    print("Reject H0: The series is stationary (p-value < 0.05)")
else:
    print("Fail to reject H0: The series is not stationary (p-value >= 0.05)\n")

# Perform the KPSS test
kpss_result = kpss(df['IPG2211A2N'], regression='c', nlags="auto")
print(f"KPSS Statistic: {kpss_result[0]:.4f}")
print(f"p-value: {kpss_result[1]:.4f}")
print(f"Critical Values: {kpss_result[3]}")
if kpss_result[1] < 0.05:
    print("Reject H0: The series is not stationary (p-value < 0.05)")
else:
    print("Fail to reject H0: The series is stationary (p-value >= 0.05)")

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['IPG2211A2N'], ax = plt.gca(), lags=30)
plt.show()

from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.figure(figsize=(20,5));
plt.grid();
plot_pacf(df['IPG2211A2N'], ax = plt.gca(), lags=20)
plt.show()

diff = df['IPG2211A2N'].diff().dropna()
diff.plot();

from statsmodels.tsa.stattools import adfuller, kpss

# Perform the ADF test
adf_result = adfuller(diff)
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
print(f"Critical Values: {adf_result[4]}")
if adf_result[1] < 0.05:
    print("Reject H0: The series is stationary (p-value < 0.05)\n")
else:
    print("Fail to reject H0: The series is not stationary (p-value >= 0.05)\n")

# Perform the KPSS test
kpss_result = kpss(diff, regression='c', nlags="auto")
print(f"KPSS Statistic: {kpss_result[0]:.4f}")
print(f"p-value: {kpss_result[1]:.4f}")
print(f"Critical Values: {kpss_result[3]}")
if kpss_result[1] < 0.05:
    print("Reject H0: The series is not stationary (p-value < 0.05)")
else:
    print("Fail to reject H0: The series is stationary (p-value >= 0.05)")

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=140)
plt.show()

from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
plt.figure(figsize=(20,5));
plt.grid();
plot_pacf(diff, ax = plt.gca(), lags=20)
plt.show()

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['IPG2211A2N'],order=(1,1,1),seasonal_order=(1,1,1,12))
model_1_fit=model_1.fit()
print(model_1_fit.summary())

#Note: here in code order (1,1,1) --> represents the parameter as d,p,q

#to get the best model of arima model --> use in-built auto_arima
from pmdarima.arima import auto_arima
model_auto_arima=auto_arima(df['IPG2211A2N'],seasonal=True,trace=True,m=12)

now using Best model: ARIMA(1,1,2)(1,0,1)[12] as suggested by auto_arima as it has the lowest AIC value
import statsmodels.api as sm
model_2=sm.tsa.statespace.SARIMAX(train_df['IPG2211A2N'],order=(1,1,2),seasonal_order=((1,0,1,12)))
model_2_fit=model_2.fit()
print(model_2_fit.summary())

interpretation of summary:
when we want to compare 2 ARIMA model compare their accuracy using AIC
here p-value for all component in low compared to above
test prediction --> we are doing this to understand whether our predicted model is accurate or Software predicted model
-using mean_absolute_error for first model

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

#### test prediction on software model 
#- using mean_absolute_error for second model

test_predictions2=model_2_fit.forecast(len(test_df))
forecast_further2=model_2_fit.forecast(len(test_df['IPG2211A2N'])+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['IPG2211A2N'],test_predictions1)
print('Mape for Test data:',mape_test1)

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

now results are subjective
if you are considering AIC to measure accuracy then Our predicted model performed well (mostly used)
if you are considering Mean_absolute_error to measure accuracy given by auto_arima model then in this case software predicted model performed better
if we want to compare 4 model Exponential Smoothing, ARIMA, SARIMA then better check Mean_absolute_error as result

#Visualizing SARIMA MOdel
# for custom model
plt.figure(figsize=(14,3))
plt.plot(train_df['IPG2211A2N'],label="Actual")
plt.plot(test_df['IPG2211A2N'],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['IPG2211A2N'],label="Actual")
plt.plot(test_df['IPG2211A2N'],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 study distribution of the residual (error) --> use diagnostic plot
it should be normal distribution
we want correlation to near zero to say model is working well --> errors are independent --> so, model is stationary

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

model_2_fit.plot_diagnostics(figsize=(15,12))
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,1,1)x(1,1,1,12)','SARIMAX(1,1,2)x(0,1,0,12)']
mape_test=[0.035345048590872966,0.033930556104136325]
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()