Forecasting

# Import dependencies

import pandas as pd
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from pmdarima import auto_arima
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
from fbprophet import Prophet
from sklearn.metrics import mean_absolute_error

warnings.filterwarnings("ignore")

metric_df = pd.read_pickle("../data/raw/ts.pkl")
### 30 mins because it gets very difficult to model otherwise
ts = metric_df["value"].astype(float).resample("30min").mean()
sns.set()
ts.plot(figsize=(15, 10))
plt.title("Visualize time series")
plt.ylabel("Node memory active bytes")
plt.show()

# Divide into train and test set
train = ts[:"2021-02-07"]
test = ts["2021-02-08":]

# Define helper dates
first_date = train.index[0]
sample_date = pd.to_datetime("2021-02-06", utc=True)
start_date = pd.to_datetime("2021-02-08", utc=True)
end_date = pd.to_datetime("2021-02-09", utc=True)

# Plot train and test sets
train.plot(figsize=(15, 10), color="green")
test.plot(figsize=(15, 10), color="blue")

def plot_predicts_aam(pred_train, pred):
    """
    Plots train, test, prediction of training set, and prediction of test set
    for auto arima and fbprophet
    """
    train.plot(figsize=(15, 10), color="green", label="Train actual")

    pd.Series(pred_train, index=train[:start_date].index).plot(
        figsize=(15, 10), color="red", label="Train prediction"
    )

    test.plot(figsize=(15, 10), color="blue", label="Test actual")

    pd.Series(pred, index=test[start_date:].index).plot(
        figsize=(15, 10), color="orange", label="Test prediction"
    )
    plt.legend()
    plt.show()


def plot_predicts_man(pred_train, pred):
    """
    Plots train, test, prediction of training set, and prediction of test set
    for manual ARIMA
    """
    train.plot(figsize=(15, 10), color="green", label="Train actual")

    pred_train.plot(figsize=(15, 10), color="red", label="Train prediction")

    test.plot(figsize=(15, 10), color="blue", label="Test actual")

    pred.plot(figsize=(15, 10), color="orange", label="Test prediction")
    plt.legend()
    plt.show()

aam_default = auto_arima(train)
pred_aam_default = aam_default.predict(n_periods=len(test))
pred_train_aam_default = aam_default.predict(n_periods=len(train[:start_date]))

aam_default.summary()

plot_predicts_aam(pred_train_aam_default, pred_aam_default)

mae_train_aam_default = mean_absolute_error(pred_train_aam_default, train)
mae_aam_default = mean_absolute_error(pred_aam_default, test)
print(f"train error: {mae_train_aam_default}, test error: {mae_aam_default}")

aam_tuned = auto_arima(
    train,
    #            exogenous=,
    m=4,  # SARIMAX s
    max_order=None,
    max_p=5,  # Search till p=6
    max_q=5,  # Search till q=6
    max_d=2,  # Search till d=2
    max_P=4,  # Search till P=2
    max_Q=4,  # Search till Q=2
    max_D=2,  # Search till D=2
    maxiter=30,  # Increase if you see no convergence
    njobs=7,  # Number of parallel processes
    trend="ct",  ##ctt for quadratic; accounts for trend in data
    information_criterion="aic",  # out of bag aic, aicc, bic, hqic
    # out_of_sample_size=int(len(ts) * 0.2),  ## Validation set of 20% for oob
)

aam_tuned.summary()

pred_aam_tuned = aam_tuned.predict(n_periods=len(test))
pred_train_aam_tuned = aam_tuned.predict(n_periods=len(train[:start_date]))
plot_predicts_aam(pred_train_aam_tuned, pred_aam_tuned)

mae_train_aam_tuned = mean_absolute_error(pred_train_aam_tuned, train)
mae_aam_tuned = mean_absolute_error(pred_aam_tuned, test)
print(f"train error: {mae_train_aam_tuned}, test error: {mae_aam_tuned}")

first_date = train.index[1]
end_date = test.index[-1]

arima_312 = ARIMA(train, order=(3, 1, 2)).fit()
pred_arima_312 = arima_312.predict(start=start_date, end=end_date)
pred_train_arima_312 = arima_312.predict(start=first_date, end=start_date)
arima_312.summary()

plot_predicts_man(pred_train_arima_312, pred_arima_312)

mae_train_arima_312 = mean_absolute_error(pred_train_arima_312, train)
mae_arima_312 = mean_absolute_error(pred_arima_312, test)
print(f"train error: {mae_train_arima_312}, test error: {mae_arima_312}")

arima_313 = ARIMA(train, order=(3, 1, 3)).fit()
pred_arima_313 = arima_313.predict(start=start_date, end=end_date)
pred_train_arima_313 = arima_313.predict(start=first_date, end=start_date)
arima_313.summary()

plot_predicts_man(pred_train_arima_313, pred_arima_313)

mae_train_arima_313 = mean_absolute_error(pred_train_arima_313, train)
mae_arima_313 = mean_absolute_error(pred_arima_313, test)
print(f"train error: {mae_train_arima_313}, test error: {mae_arima_313}")

### SARIMAX

sarima_3133134 = SARIMAX(
    train, order=(3, 1, 3), seasonal_order=(3, 1, 3, 4)
).fit()
pred_sarima_3133134 = sarima_3133134.predict(start=start_date, end=end_date)
pred_train_sarima_3133134 = sarima_3133134.predict(
    start=first_date, end=start_date
)
sarima_3133134.summary()

plot_predicts_man(pred_train_sarima_3133134, pred_sarima_3133134)

mae_train_sarima_3133134 = mean_absolute_error(
    pred_train_sarima_3133134, train
)
mae_sarima_3133134 = mean_absolute_error(pred_sarima_3133134, test)
print(
    f"train error: {mae_train_sarima_3133134}, test error: {mae_sarima_3133134}"
)

## Convert data to prophet type
train_fb = (
    train.tz_localize(None)
    .reset_index()
    .rename(columns={"timestamp": "ds", "value": "y"})
)
test_fb = (
    test.tz_localize(None)
    .reset_index()
    .rename(columns={"timestamp": "ds", "value": "y"})
)
train_fb

proph_model = Prophet()
pm = proph_model.fit(train_fb)

future = pd.concat([train_fb[["ds"]], test_fb[["ds"]]])
forecast = proph_model.predict(future)
forecast[["ds", "yhat", "yhat_lower", "yhat_upper"]].tail()

plot_predicts_aam(forecast["yhat"][:278].values, forecast["yhat"][278:].values)

## Prophet finds components as well and forecasts the trend
fig2 = proph_model.plot_components(forecast)

mae_train_prophet = mean_absolute_error(forecast.yhat.iloc[:278], train)
mae_prophet = mean_absolute_error(forecast.yhat.iloc[278:], test)
print(f"train error: {mae_train_prophet}, test error: {mae_prophet}")

sns.set()
plt.rc("figure", figsize=(15, 10))
sd_add = seasonal_decompose(train, model="additive", freq=2 * 24)
sd_add.plot()
plt.show()

train_adj = train - sd_add.seasonal
train_adj_fb = (
    train_adj.tz_localize(None)
    .reset_index()
    .rename(columns={"timestamp": "ds", 0: "y"})
)
test_fb = (
    test.tz_localize(None)
    .reset_index()
    .rename(columns={"timestamp": "ds", 0: "y"})
)

# Here we take one period from the train decomposed set
# with similar start and end point as the test set
test_set_period_adj = sd_add.seasonal["2021-02-03":"2021-02-04 4:30"]

## Prophet model
proph_model = Prophet()
pm = proph_model.fit(train_adj_fb)
future = pd.concat([train_adj_fb[["ds"]], test_fb[["ds"]]])
forecast = proph_model.predict(future)

train_predict_fb_adj = forecast["yhat"][:278].values + sd_add.seasonal.values
predict_fb_adj = forecast["yhat"][278:].values + test_set_period_adj.values
plot_predicts_aam(train_predict_fb_adj, predict_fb_adj)

mae_train_fb_adj = mean_absolute_error(train_predict_fb_adj, train)
mae_fb_adj = mean_absolute_error(predict_fb_adj, test)
print(f"train error: {mae_train_fb_adj}, test error: {mae_fb_adj}")

pd.DataFrame(
    data=[
        ["aam_default", mae_train_aam_default, mae_aam_default],
        ["aam_tuned", mae_train_aam_tuned, mae_aam_tuned],
        ["arima_312", mae_train_arima_312, mae_arima_312],
        ["ARIMA_313", mae_train_arima_313, mae_arima_313],
        ["sarima_3133134", mae_train_sarima_3133134, mae_sarima_3133134],
        ["prophet", mae_train_prophet, mae_prophet],
        ["decompose+ARIMA", mae_train_fb_adj, mae_fb_adj],
    ],
    columns=["Method", "Train error", "Test error"],
)