<jupyter_start><jupyter_text>Problem FormulationInput: time series data from a lahge retailer where each point is associated with a specific sales volumeOutput: prediction for unit sales volume 2 months (60 days) from nowRequired data mining function: determine patterns in fluctuating time data such as trend, seasonality, and stationary variancePotential challenges: noisy data that leads to ambiguous predictions, bias + overfitting as usualImpakt: more $$ in revenue for company, convert tough customers into happy onesIdeal solution: highly accurate projection from small dataset with minimal train/test time<jupyter_code>#regular schmegular importation

import numpy as np
from datetime import datetime
import pandas as pd

def rmsle(y_hat, y):
    """Compute Root Mean Squared Logarithmic Error"""
    metric = np.sqrt(sum((np.array(list(map(lambda x : np.log(x + 1), y_hat)))
                         - np.array(list(map(lambda x : np.log(x + 1), y))))**2)/len(y))

    return round(metric, 4)

from pandas import read_csv
#from pandas import datetime
from matplotlib import pyplot
# load data (monthly sales)

# use datetime (the Month column) as index
df = read_csv('train.csv', header=0)
df.head()

#change to unambiguous formatting
format = "%Y-%m-%d" 
df['ds']=df.ds.apply(lambda x: datetime.strptime(x, format).date())
df['ds']=pd.to_datetime(df['ds'])

df.info()
# plot the time series (monthly sale, the data was in month any way)
df.groupby(df["ds"].dt.month)["y"].mean().plot(
    kind='bar', rot=0, 
)
pyplot.show()
# plot the time series (yearly average sale)
df.groupby(df["ds"].dt.year)["y"].mean().plot(
    kind='bar', rot=0, 
)
pyplot.show()
# any other way to zoom in/zoom out on the time scale by resample?
# But we have to set Month (the datetime) column as index.
df_3m = df.set_index('ds')
# aggregation based on 3 months
df_3m.resample('3M').mean().head()

# aggregation based on 2 years
df_3m.resample('2Y').mean().head()
# Comparing multiple time series
import matplotlib.pyplot as plt
s1 = df.groupby(df["ds"].dt.year)["y"].mean()
s2 = df.groupby(df["ds"].dt.year)["y"].mean() + 1000
s1.plot()
s2.plot()
plt.legend(['s1','s2'])
plt.show()
# Rolling 90D average
rolling_df = df.set_index('ds')
rolling_df.y.rolling('90D').mean().plot()
plt.show()
df['y']
#attempt autocorrelation
from statsmodels.graphics.tsaplots import plot_acf
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (20,7)

acf_plot = plot_acf(df.y, lags = 70)
plt.xlabel('lag')
plt.ylabel('correlation')
plt.show()

#All lags having a value between 1.0 and 0.25 indicate that there is some monthly seasonality in this dataset
#attempt partial autocorrelation
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (20,7)

acf_plot = plot_pacf(df.y, lags = 10)
plt.xlabel('lag')
plt.ylabel('correlation')
plt.show()

#No statistically significant result after 9 lags<jupyter_output>C:\Users\Owner\AppData\Roaming\Python\Python39\site-packages\statsmodels\graphics\tsaplots.py:348: FutureWarning: The default method 'yw' can produce PACF values outside of the [-1,1] interval. After 0.13, the default will change tounadjusted Yule-Walker ('ywm'). You can use this method now by setting method='ywm'.
  warnings.warn(<jupyter_text>ARIMA<jupyter_code>from statsmodels.tsa.stattools import adfuller
from numpy import log

result = adfuller(df.y)
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])

#What effect will differencing have?

import numpy as np, pandas as pd
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt
plt.rcParams.update({'figure.figsize':(9,7), 'figure.dpi':120})

# Original Series
fig, axes = plt.subplots(4, 2, sharex=True)
axes[0, 0].plot(df.y); axes[0, 0].set_title('Original Series')
plot_acf(df.y, ax=axes[0, 1])

# 1st Differencing
axes[1, 0].plot(df.y.diff()); axes[1, 0].set_title('1st Order Differencing')
plot_acf(df.y.diff().dropna(), ax=axes[1, 1])

# 2nd Differencing
axes[2, 0].plot(df.y.diff().diff()); axes[2, 0].set_title('2nd Order Differencing')
plot_acf(df.y.diff().diff().dropna(), ax=axes[2, 1])

# 3rd Differencing
axes[3, 0].plot(df.y.diff().diff().diff()); axes[3, 0].set_title('3rd Order Differencing')
plot_acf(df.y.diff().diff().diff().dropna(), ax=axes[3, 1])

plt.show()

#Approximate stationarity is reached after 3 differences - VERY GOOD
plt.rcParams["figure.figsize"] = (20,7)

acf_plot = plot_pacf(df.y, lags = 30) #will increasing the number of lags have an effect?
plt.xlabel('lag')
plt.ylabel('correlation')
plt.show()

#The lack of statistically significant effects after the 15th lag suggests that seasonality is limited to approx. every 2 weeks
#next, investigate whether shifting the moving average order makes anyhting interesting happen

fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].plot(df.y.diff()); axes[0].set_title('1st Differencing')
axes[1].set()
plot_acf(df.y.diff().dropna(), ax=axes[1])

plt.show()
#INITIATE OPERATION CROSS_VALIDATE
import statsmodels.api as sm

df_ar = df.copy()

train = df_ar.y[:1200]
test = df_ar.y[1200:]

# order = (p, d, q) 
model = sm.tsa.ARIMA(train, (1, 2, 0))
model_fit = model.fit()
model_fit.summary()
# Plot residual errors
import pandas as pd

residuals = pd.DataFrame(model_fit.resid)
fig, ax = plt.subplots(1,2)
residuals.plot(title="Residuals", ax=ax[0])
residuals.plot(kind='kde', title='Density', ax=ax[1])
plt.show()
model_fit.plot_predict(dynamic=False)
plt.show()
# Forecast
fc, se, conf = model_fit.forecast(12, alpha=0.9)  # 95% conf

# Make as pandas series
fc_series = pd.Series(fc, index=test.index)
lower_series = pd.Series(conf[:, 0], index=test.index)
upper_series = pd.Series(conf[:, 1], index=test.index)

# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(train, label='training')
plt.plot(test, label='actual')
plt.plot(fc_series, label='forecast')
plt.fill_between(lower_series.index, lower_series, upper_series, 
                 color='k', alpha=.15)
plt.title('Forecast vs Actuals')
plt.legend(loc='upper left', fontsize=8)
plt.show()<jupyter_output><empty_output><jupyter_text>PROPHET<jupyter_code>from pandas import read_csv
from pandas import to_datetime
from pandas import DataFrame
from fbprophet import Prophet
from matplotlib import pyplot


# prophet requires fixed name of columns. 
# ds -> datetime
# y -> values

df_pp = df.copy()
df_pp.columns = ['ds', 'y']
df_pp['ds'] = to_datetime(df_pp['ds'])
df_pp.head()

train_pp = df_pp[:1300]
test_pp = df_pp[1300:]
# define the model
model = Prophet()
# fit the model
model.fit(df_pp)

# get the prediction dates
future = test_pp
forecast = model.predict(future)
# summarize the forecast
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].head())
# plot forecast
model.plot(forecast)
pyplot.show()
#DEKOMPOSE

model.plot_components(forecast)
# plot change points in the trend.
from fbprophet.plot import add_changepoints_to_plot

fig = model.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), model, forecast)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14424703891708968, 'fraction_numerical': 0.02143260011280316, 'mean_word_length': 5.543357933579336, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '307880', 'license': 'MIT', 'n_tokens_mistral': 2663, 'n_tokens_neox': 2478, 'n_words': 681, 'path': 'A3_SUBMIT.ipynb', 'repo_name': 'phiriv/CISC--371-Intro-to-Data-Sci-', 'size': '7092'}