# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: env
#     language: python
#     name: env
# ---

# # P03_01_ARIMA_AR

section_fig_main_title = 'P03_01_ARIMA_AR'

import tensorflow as tf
tf.__version__

import warnings
warnings.filterwarnings('ignore')
warnings.simplefilter(action='ignore', category=FutureWarning)

# %config InlineBackend.figure_format = 'retina'


import gif

import sys
# print(sys.executable)

# +
# # !/Users/amin/Desktop/proj/env/bin/python -m pip install matplotlib
# -

import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras

import pandas as pd
import numpy as np
import seaborn as sns

from statsmodels.tsa.ar_model import AutoReg
from sklearn.metrics import mean_squared_error
# import statsmodels.api as sm

# ## Read data

csv_path = "../Data/jena_climate_2009_2016.csv"
df = pd.read_csv(csv_path, index_col=0, parse_dates=True)

df.head()

df.columns

df.info()

# +
# df.corr().T
# -

# I choose the first 20000 timestamps as sub-data and used 720 as lag. The model checks out the 720 timestamps and uses them for future predictions. (this selection is based on keras website example for time series modeling). later I apply the model for 72 timestamps prediction ahead (12 hours ahead).

# +
column_name = 'T (degC)'
start = 0
end = 20000

training_df = df[column_name][start:end]
# -

# ## TO BE stationary OR NOT TO BE stationary
#
# Let's use Augmented Dickey Fuller test (ADF Test) to see if data is stationary or not.

from statsmodels.tsa.stattools import adfuller

adfuller_test = adfuller(df[column_name][start: end], autolag='AIC')

test_details = ['ADF', 'P-Vale', 'NumofLags', 'NumOfObser_ADF_regr', 'critical_values']
for i,name in enumerate(test_details):
    print(name, '  ', adfuller_test[i])

# Base on ADF test, the Data is stationary (p-value < 0.05) so null hypothesis (data is not stationary) is rejected. also as was observed in EDA parts, no trend is observed. 

plt.figure(figsize=(20,5))
df['T (degC)'].plot(style='k.')
plt.show()

# plt.figure(figsize=(20,5))
training_df.plot(style='k.')
plt.show()

# However I see a seanality in data! Let's split the data and check out the mean and variance and compare if there are any changes.

# +
# source https://machinelearningmastery.com/time-series-data-stationary-python/
X = df['T (degC)']

split = int(len(X) / 2)
X1, X2 = X[0:split], X[split:]
mean1, mean2 = X1.mean(), X2.mean()
var1, var2 = X1.var(), X2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))
print('variance1=%f, variance2=%f' % (var1, var2))
# -

# The mean and variance are so close! let's plot the histogram for data and check out the data distribution.

df.hist(column = ['T (degC)'], figsize=(10,5), bins=50, legend=True)

# the data distribution is normal. I think we can accept that the data is stationary (taking the ADF test, mean/var, and data distribution into account).



# # ARIMA
#
# I choose the ARIMA model which combines three methods including autoregressive, moving average, and integration.
#
# ### Searching for order of ARIMA model

# ARIMA model is composed of three parts including autoregressive, moving average, and integrated.
# auto regressive means the model considers the relationship between an obervation and some lagged observations (history). The integrated part use difference (subtraction) method to make data stationary. Moving average part, first applies a moving average on lagged observationas, second it calculates the difference between the moving average model and the residuals and, finally uses a relationship between an obervation and residuals (predicted by moving average model).
#
# The ARIMA model has three hyper-parameters (p,d,q) which are repersentative for number of lagged observation to be used in autoregressive (so called laged order)(p), the number of times the observations are differenced (d), and the size of moving average window (q).
#
# If the order of p,q are 5, that means the model used last 5 observation to predict the future time-stamps.
#
#
# ARMA model is used for stationary data, while if the data is not stationary, the ARIMA comes to play.

# +
# # !/Users/amin/Desktop/proj/env/bin/python -m pip install pmdarima
# -

from pmdarima import auto_arima

arima_model_fit = auto_arima(training_df, trace=True, suppress_warnings=True)

# It is seen that the best ARIMA model has (2,1,3) orders. The summary is provided at the following.

arima_model_fit.summary()

# The ARIMA model coeficinet are:
#
# ```ar.L1 = 1.9857, ar.L2 = -0.9876, 
# ma.L1 = -1.4953, ma.L2 = 0.1089, ma.L3 = 0.3896```
#

# ## Prediction
#
# ### ARIMA prediction for 72 timestamps ahead

# +
# from statsmodels.tsa.arima.model import ARIMA

# +
# # fit model
# model = ARIMA(df['T (degC)'][0:20000], order=(720,1,0))
# model_fit = model.fit()
# # summary of fit model
# print(model_fit.summary())

# +
# start = len(training_df)
# ahead = 72 # 12 hours ahead
# end = start + ahead


# expected = df['T (degC)'][start:end]
# predicted = arima_model_fit.predict(n_periods=ahead)

# +
start = len(training_df)
ahead = 72 # 12 hours ahead
end = start + ahead


expected = df['T (degC)'][start:end]
predicted = arima_model_fit.predict(n_periods=ahead)

diff = predicted - expected

predicted = pd.DataFrame(predicted)
predicted.index = df[column_name][start:end].index
diff = pd.DataFrame(diff)


plt.figure(figsize=(20,8))
plt.plot(predicted, label='ARIMA prediction')
plt.plot(expected, 'r^', label='Data - Expected')


plt.xlabel('Date')
plt.ylabel('Temperature c')
plt.title("Temperature prediction")

plt.ylim(0,)

plt.legend()
plt.savefig("../Figures/" + section_fig_main_title + "ARIMA_Temp_esti_72ahead_1.png")

# +
start = len(training_df)
ahead = 72
end = start + ahead

past = 100


expected = df['T (degC)'][start:end]
predicted = arima_model_fit.predict(n_periods=ahead)

diff = predicted - expected

predicted = pd.DataFrame(predicted)
predicted.index = df[column_name][start:end].index
diff = pd.DataFrame(diff)


plt.figure(figsize=(20,8))
plt.plot(predicted, label='ARIMA prediction')
plt.plot(expected, 'r^', label='Data - Expected')



plt.plot(df[column_name][start-past:end], c='r')


plt.xlabel('Date')
plt.ylabel('Temperature c')
plt.title("Temperature prediction")


star = df['T (degC)'][start:start+1].index
plt.axvline(star,color='k', linestyle='--', label='test_line')
# plt.ylim(0,)

plt.legend()

# plt.show()

plt.savefig("../Figures/" + section_fig_main_title + "ARIMA_Temp_esti_72ahead.png")
# -

# ### ARIMA prediction for 720 timestamps ahead

# +
start = len(training_df)
ahead = 720
end = start + ahead

past = 100


expected = df['T (degC)'][start:end]
predicted = arima_model_fit.predict(n_periods=ahead)

diff = predicted - expected

predicted = pd.DataFrame(predicted)
predicted.index = df[column_name][start:end].index
diff = pd.DataFrame(diff)


plt.figure(figsize=(20,8))
plt.plot(predicted, label='ARIMA prediction')
plt.plot(expected, 'r^', label='Data - Expected')



plt.plot(df[column_name][start-past:end], c='r')


plt.xlabel('Date')
plt.ylabel('Temperature c')
plt.title("Temperature prediction")


# star = df['T (degC)'][start].index
star = df['T (degC)'][start:start+1].index
plt.axvline(star,color='k', linestyle='--', label='test_line')

# plt.ylim(0,)

plt.legend()
plt.savefig("../Figures/" + section_fig_main_title + "ARIMA_Temp_esti_720ahead.png")
# -

star = df['T (degC)'][start:start+1].index
star


# as it is seen from above plot, the model error increases by time. 
#
# The ARIMA model successfully estimates the temperature for future 72 time steamps. 
# Applying the ARIMA model for higher time steamps decreases the model estimation accuracy. 

# ### GIF visualization ARIMA model

# let's write a function for plot and use it make an animation visualization.

# +
@gif.frame
def plott(ahead, past):
    
    start = len(training_df)
    ahead = ahead
    end = start + ahead

    past = past


    expected = df['T (degC)'][start:end]
    predicted = arima_model_fit.predict(n_periods=ahead)

    diff = predicted - expected

    predicted = pd.DataFrame(predicted)
    predicted.index = df[column_name][start:end].index
    diff = pd.DataFrame(diff)


    plt.figure(figsize=(20,8))
    plt.plot(predicted, label='ARIMA prediction')
    plt.plot(expected, 'r^', label='Data - Expected')


    plt.plot(df[column_name][start-past:end], c='r')

    plt.xlabel('Date')
    plt.ylabel('Temperature c')
    plt.title("Temperature prediction")


    star = df['T (degC)'][start:start+1].index
    plt.axvline(star,color='k', linestyle='--', label='test_line')
    # plt.ylim(0,)

    plt.legend()

    # plt.show()

#     plt.savefig("../Figures/" + section_fig_main_title + "Temp_esti_72ahead.png")

# +
# # !/Users/amin/Desktop/proj/env/bin/python -m pip install gif

# +
# import gif

# +
frames = []    
for i in range(1,720):

    
    frame = plott(ahead = i, past=100)
    frames.append(frame)
    
# gif.save(frames, "../Figures/GIF.gif", duration=2)
# -

gif.save(frames, "../Figures/GIF.gif", duration=2)

# ## Autoregressive (AR)
#
# I apply AR and see how close I can get compare to ARIMA. similarly, i use 20000 sub-data with lags=720 for training and later I apply the model for 72 timestamps ahead (12 hours).

# train autoregression
model = AutoReg(training_df, lags=720)
ar_model_fit = model.fit()
# print('Coefficients: %s' % ar_model_fit.params)

# ### AR prediction for 72 timestamps ahead

start = len(training_df)
ahead = 72
end = start + ahead
print('start, ahead, end\n', start, ahead, end)
test = df['T (degC)'][start:end]

# +
# make predictions
predictions = ar_model_fit.predict(start=len(training_df), 
                                end=len(training_df)+len(test)-1, 
                                dynamic=False)

# apply same index for prediction df
predictions.index = test.index

# +
plt.figure(figsize=(20,8))
plt.plot(predictions, label='AR prediction')
plt.plot(test, 'r^', label='Data - Expected')

past = 100

plt.xlabel('Date')
plt.ylabel('Temperature c')
plt.title("Temperature prediction")

star = df['T (degC)'][start:start+1].index
plt.axvline(star,color='k', linestyle='--', label='test_line')
plt.plot(df[column_name][start-past:end], c='r')

plt.ylim(0,)

plt.legend()
plt.savefig("../Figures/" + section_fig_main_title + "AR_Temp_esti_72ahead_1.png")
# -

# ### AR prediction for 720 timestamps ahead

start = len(training_df)
ahead = 720
end = start + ahead
print('start, ahead, end\n', start, ahead, end)
test = df['T (degC)'][start:end]

# +
# make predictions
predictions = ar_model_fit.predict(start=len(training_df), 
                                end=len(training_df)+len(test)-1, 
                                dynamic=False)

# apply same index for prediction df
predictions.index = test.index

# +
plt.figure(figsize=(20,8))
plt.plot(predictions, label='AR prediction')
plt.plot(test, 'r^', label='Data - Expected')


plt.xlabel('Date')
plt.ylabel('Temperature c')
plt.title("Temperature prediction")


plt.plot(df[column_name][start-past:end], c='r')
star = df['T (degC)'][start:start+1].index
plt.axvline(star,color='k', linestyle='--', label='test_line')
plt.ylim(0,)




plt.legend()
plt.savefig("../Figures/" + section_fig_main_title + "AR_Temp_esti_720ahead_1.png")
# -



@gif.frame
def plott(ahead, past):
    
    start = len(training_df)
    ahead = ahead
    end = start + ahead

    past = past


    expected = df['T (degC)'][start:end]
#     predicted = arima_model_fit.predict(n_periods=ahead)
    predicted =  ar_model_fit.predict(start=len(training_df), 
                                end=end-1, 
                                dynamic=False)

    diff = predicted - expected

    predicted = pd.DataFrame(predicted)
    predicted.index = df[column_name][start:end].index
    diff = pd.DataFrame(diff)


    plt.figure(figsize=(20,8))
    plt.plot(predicted, label='AR prediction')
    plt.plot(expected, 'r^', label='Data - Expected')


    plt.plot(df[column_name][start-past:end], c='r')

    plt.xlabel('Date')
    plt.ylabel('Temperature c')
    plt.title("Temperature prediction")


    star = df['T (degC)'][start:start+1].index
    plt.axvline(star,color='k', linestyle='--', label='test_line')

    plt.legend()

# +
# plott(2, 100)
# -

frames = []    
for i in range(1,720):

    
    frame = plott(ahead = i, past=100)
    frames.append(frame)

gif.save(frames, "../Figures/AR_GIF.gif", duration=2)

# Note that the ARIMA model uses 2 lags (d=2) for modeling while in the AR model I applied 72 and 720 lags which results in better estimation. The ARIMA model applied d=1 to make sure data is stationary (suggested by the autoarima library). The data as was observed in the P03_01 and previous EDA part is almost stationary therefore there is no need for applying differencing method.





StarcoderdataJupyterScriptsDedupFiltered