y
timestamp
1980-09-25 14:01:00  182.478
1980-09-25 14:02:00  176.231
1980-09-25 14:03:00  183.917
1980-09-25 14:04:00  177.798
1980-09-25 14:05:00  165.469
1980-09-25 14:06:00  181.878
1980-09-25 14:07:00  184.502
1980-09-25 14:08:00  183.303
1980-09-25 14:09:00  177.578
1980-09-25 14:10:00  171.641
1980-09-25 14:11:00  191.014
1980-09-25 14:12:00  184.068
1980-09-25 14:13:00  188.457
1980-09-25 14:14:00  175.739
1980-09-25 14:15:00  175.524
1980-09-25 14:16:00  189.128
1980-09-25 14:17:00  176.885
1980-09-25 14:18:00  167.140
1980-09-25 14:19:00  173.723
1980-09-25 14:20:00  168.460
1980-09-25 14:21:00  177.623
1980-09-25 14:22:00  183.888

                           y
timestamp
1980-09-25 14:01:00  182.478
1980-09-25 14:02:00  176.231
1980-09-25 14:03:00  183.917
1980-09-25 14:04:00  177.798
1980-09-25 14:05:00  165.469
1980-09-25 14:06:00  181.878
1980-09-25 14:07:00  184.502
1980-09-25 14:08:00  183.303
1980-09-25 14:09:00  177.578
1980-09-25 14:10:00  171.641
1980-09-25 14:11:00  191.014
1980-09-25 14:12:00  184.068
1980-09-25 14:13:00  188.457
1980-09-25 14:14:00  175.739
1980-09-25 14:15:00  175.524
1980-09-25 14:16:00  189.128
1980-09-25 14:17:00  176.885
1980-09-25 14:18:00  167.140
1980-09-25 14:19:00  173.723
1980-09-25 14:20:00  168.460
1980-09-25 14:21:00  177.623
1980-09-25 14:22:00  183.888
1980-09-25 14:23:00  167.487
1980-09-25 14:24:00  165.572
1980-09-25 14:25:00  170.480
1980-09-25 14:26:00  172.474
1980-09-25 14:27:00  166.448
1980-09-25 14:28:00  163.098
1980-09-25 14:29:00  163.544
1980-09-25 14:30:00  163.816y    lag_6   ...      lag_23   lag_24
timestamp                               ...
1980-09-25 14:01:00  182.478      NaN   ...         NaN      NaN
1980-09-25 14:02:00  176.231      NaN   ...         NaN      NaN
1980-09-25 14:03:00  183.917      NaN   ...         NaN      NaN
1980-09-25 14:04:00  177.798      NaN   ...         NaN      NaN
1980-09-25 14:05:00  165.469      NaN   ...         NaN      NaN
1980-09-25 14:06:00  181.878      NaN   ...         NaN      NaN
1980-09-25 14:07:00  184.502  182.478   ...         NaN      NaN
1980-09-25 14:08:00  183.303  176.231   ...         NaN      NaN
1980-09-25 14:09:00  177.578  183.917   ...         NaN      NaN
1980-09-25 14:10:00  171.641  177.798   ...         NaN      NaN
1980-09-25 14:11:00  191.014  165.469   ...         NaN      NaN
1980-09-25 14:12:00  184.068  181.878   ...         NaN      NaN
1980-09-25 14:13:00  188.457  184.502   ...         NaN      NaN
1980-09-25 14:14:00  175.739  183.303   ...         NaN      NaN
1980-09-25 14:15:00  175.524  177.578   ...         NaN      NaN
1980-09-25 14:16:00  189.128  171.641   ...         NaN      NaN
1980-09-25 14:17:00  176.885  191.014   ...         NaN      NaN
1980-09-25 14:18:00  167.140  184.068   ...         NaN      NaN
1980-09-25 14:19:00  173.723  188.457   ...         NaN      NaN
1980-09-25 14:20:00  168.460  175.739   ...         NaN      NaN
1980-09-25 14:21:00  177.623  175.524   ...         NaN      NaN
1980-09-25 14:22:00  183.888  189.128   ...         NaN      NaN
1980-09-25 14:23:00  167.487  176.885   ...         NaN      NaN
1980-09-25 14:24:00  165.572  167.140   ...     182.478      NaN
1980-09-25 14:25:00  170.480  173.723   ...     176.231  182.478
1980-09-25 14:26:00  172.474  168.460   ...     183.917  176.231
1980-09-25 14:27:00  166.448  177.623   ...     177.798  183.917
1980-09-25 14:28:00  163.098  183.888   ...     165.469  177.798
1980-09-25 14:29:00  163.544  167.487   ...     181.878  165.469
1980-09-25 14:30:00  163.816  165.572   ...     184.502  181.878

# coding: utf-8import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit  # you have everything done for you
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, RidgeCV# for time-series cross-validation set 5 folds
tscv = TimeSeriesSplit(n_splits=5)def mean_absolute_percentage_error(y_true, y_pred):return np.mean(np.abs((y_true - y_pred) / y_true)) * 100def timeseries_train_test_split(X, y, test_size):"""Perform train-test split with respect to time series structure"""# get the index after which test set startstest_index = int(len(X) * (1 - test_size))X_train = X.iloc[:test_index]y_train = y.iloc[:test_index]X_test = X.iloc[test_index:]y_test = y.iloc[test_index:]return X_train, X_test, y_train, y_testdef plotModelResults(model, X_train, X_test, y_train, y_test, plot_intervals=False, plot_anomalies=False):"""Plots modelled vs fact values, prediction intervals and anomalies"""prediction = model.predict(X_test)plt.figure(figsize=(15, 7))plt.plot(prediction, "g", label="prediction", linewidth=2.0)plt.plot(y_test.values, label="actual", linewidth=2.0)if plot_intervals:cv = cross_val_score(model, X_train, y_train,cv=tscv,scoring="neg_mean_absolute_error")mae = cv.mean() * (-1)deviation = cv.std()scale = 20lower = prediction - (mae + scale * deviation)upper = prediction + (mae + scale * deviation)plt.plot(lower, "r--", label="upper bond / lower bond", alpha=0.5)plt.plot(upper, "r--", alpha=0.5)if plot_anomalies:anomalies = np.array([np.NaN] * len(y_test))anomalies[y_test < lower] = y_test[y_test < lower]anomalies[y_test > upper] = y_test[y_test > upper]plt.plot(anomalies, "o", markersize=10, label="Anomalies")error = mean_absolute_percentage_error(prediction, y_test)plt.title("Mean absolute percentage error {0:.2f}%".format(error))plt.legend(loc="best")plt.tight_layout()plt.grid(True);plt.savefig("linear.png")def plotCoefficients(model, X_train):"""Plots sorted coefficient values of the model"""coefs = pd.DataFrame(model.coef_, X_train.columns)coefs.columns = ["coef"]coefs["abs"] = coefs.coef.apply(np.abs)coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)plt.figure(figsize=(20, 12))coefs.coef.plot(kind='bar')plt.grid(True, axis='y')plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')plt.savefig("linear-cov.png")def code_mean(data, cat_feature, real_feature):"""Returns a dictionary where keys are unique categories of the cat_feature,and values are means over real_feature"""return dict(data.groupby(cat_feature)[real_feature].mean())def prepareData(series, lag_start, lag_end, test_size, target_encoding=False):"""series: pd.DataFramedataframe with timeserieslag_start: intinitial step back in time to slice target variableexample - lag_start = 1 means that the modelwill see yesterday's values to predict todaylag_end: intfinal step back in time to slice target variableexample - lag_end = 4 means that the modelwill see up to 4 days back in time to predict todaytest_size: floatsize of the test dataset after train/test split as percentage of datasettarget_encoding: booleanif True - add target averages to the dataset"""# copy of the initial datasetdata = pd.DataFrame(series.copy())data.columns = ["y"]# lags of seriesfor i in range(lag_start, lag_end):data["lag_{}".format(i)] = data.y.shift(i)# datetime features# data.index = data.index.to_datetime()data["hour"] = data.index.hourdata["weekday"] = data.index.weekdaydata['is_weekend'] = data.weekday.isin([5, 6]) * 1if target_encoding:# calculate averages on train set onlytest_index = int(len(data.dropna()) * (1 - test_size))data['weekday_average'] = list(map(code_mean(data[:test_index], 'weekday', "y").get, data.weekday))data["hour_average"] = list(map(code_mean(data[:test_index], 'hour', "y").get, data.hour))# drop encoded variables# data.drop(["hour", "weekday"], axis=1, inplace=True)# train-test splity = data.dropna().yX = data.dropna().drop(['y'], axis=1)X_train, X_test, y_train, y_test = \timeseries_train_test_split(X, y, test_size=test_size)return X_train, X_test, y_train, y_testdef plt_linear():data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])data['timestamp'] = pd.to_datetime(data['timestamp'])data.set_index("timestamp", drop=True, inplace=True)data.rename(columns={'count': 'y'}, inplace=True)X_train, X_test, y_train, y_test = \prepareData(data, lag_start=6, lag_end=25, test_size=0.3, target_encoding=True)scaler = StandardScaler()X_train_scaled = scaler.fit_transform(X_train)X_test_scaled = scaler.transform(X_test)# lr = LinearRegression()lr = LassoCV(cv=tscv)# lr = RidgeCV(cv=tscv)# lr.fit(X_train_scaled, y_train)"""from xgboost import XGBRegressorlr = XGBRegressor()# lr = xgb.XGBRegressor(max_depth=5, learning_rate=0.1, n_estimators=160, silent=False, objective='reg:gamma')"""lr.fit(X_train_scaled, y_train)"""IMPORTANTGenerally tree-based models poorly handle trends in data, compared to linear models, so you have to detrend your series first or use some tricks to make the magic happen. Ideally make the series stationary and then use XGBoost, for example, you can forecast trend separately with a linear model and then add predictions from xgboost to get final forecast."""plotModelResults(lr, X_train=X_train_scaled, X_test=X_test_scaled,  y_train=y_train, y_test=y_test, plot_intervals=True, plot_anomalies=True)plotCoefficients(lr, X_train=X_train)plt_linear()

# coding: utf-8import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, RidgeCVdef mean_absolute_percentage_error(y_true, y_pred):return np.mean(np.abs((y_true - y_pred) / y_true)) * 100def timeseries_train_test_split(X, y, test_size, predict_size):"""Perform train-test split with respect to time series structure"""total_size = len(X)# get the index after which test set startstest_index = int(total_size * (1 - test_size))X_train = X.iloc[:test_index]y_train = y.iloc[:test_index]X_test = X.iloc[test_index:total_size-predict_size]y_test = y.iloc[test_index:total_size-predict_size]X_predict = X.iloc[-predict_size:]y_predict = y.iloc[-predict_size:]return X_train, X_test, y_train, y_test, X_predict, y_predictdef predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler):# for predict# not OK for abnormal real valuey_predict[0:lag_start] = lr.predict(scaler.transform(X_predict.iloc[0:lag_start]))# last_line = X_test.iloc[-1]for i in range(lag_start, len(X_predict)):# for i in range(0, len(X_predict)):last_line = X_predict.iloc[i-1]index = X_predict.index[i]for j in range(lag_end-1, lag_start):X_predict.at[index, "lag_{}".format(j)] = last_line["lag_{}".format(j-1)]X_predict.at[index, "lag_{}".format(lag_start)] = y_predict[i-1]y_predict[i] = lr.predict(scaler.transform([X_predict.iloc[i]]))[0]return y_predictdef plot_results(y_predict, y, intervals, img_filename, plot_intervals=False, plot_anomalies=False, extra_plot=None):"""Plots modelled vs fact values, prediction intervals and anomalies"""assert len(y_predict) == len(y)plt.figure(figsize=(15, 7))# plt.plot(y.index, y_predict, "g", label="prediction", linewidth=3.0)# plt.plot(y.index, y.values, label="actual", linewidth=1.0)plt.plot(y.index, y_predict, ls='-', c='#0072B2', label='predicted y')plt.plot(y.index, y.values, 'k.', label='y')if extra_plot is not None:# plt.plot(extra_plot.index, extra_plot.values, "y", label="future predict", linewidth=3.0)plt.plot(extra_plot.index, extra_plot.values, 'y', label='predicted y')if plot_intervals:lower = y_predict - intervalsupper = y_predict + intervals# plt.plot(y.index, lower, "r--", label="upper bond / lower bond", alpha=0.5)# plt.plot(y.index, upper, "r--", alpha=0.5)plt.fill_between(y.index, lower, upper, color='#0072B2', alpha=0.2, label='predicted upper/lower y')if extra_plot is not None:# plt.plot(extra_plot.index, extra_plot.values-intervals, "r--", label="upper bond / lower bond", alpha=0.5)# plt.plot(extra_plot.index, extra_plot.values+intervals, "r--", alpha=0.5)plt.fill_between(extra_plot.index, extra_plot.values-intervals, extra_plot.values+intervals,color='#0072B2', alpha=0.2, label='predicted upper/lower y')if plot_anomalies:anomalies_lower = y[y < lower]anomalies_upper = y[y > upper]# plt.plot(anomalies_lower.index, anomalies_lower.values, "ro", markersize=10, label="Anomalies(+)")# plt.plot(anomalies_upper.index, anomalies_upper.values, "ro", markersize=10, label="Anomalies(-)")plt.plot(anomalies_lower.index, anomalies_lower.values, "rX", label='abnormal points')plt.plot(anomalies_upper.index, anomalies_upper.values, "rX")error = mean_absolute_percentage_error(y_predict, y)plt.title("Mean absolute percentage error {0:.2f}%".format(error))plt.legend(loc="best")plt.tight_layout()plt.grid(True);plt.savefig(img_filename)def plot_arg_importance(model, X_train):"""Plots sorted coefficient values of the model"""coefs = pd.DataFrame(model.coef_, X_train.columns)coefs.columns = ["coef"]coefs["abs"] = coefs.coef.apply(np.abs)coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)plt.figure(figsize=(20, 12))coefs.coef.plot(kind='bar')plt.grid(True, axis='y')plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')plt.savefig("linear-cov.png")def code_mean(data, cat_feature, real_feature):"""Returns a dictionary where keys are unique categories of the cat_feature,and values are means over real_feature"""return dict(data.groupby(cat_feature)[real_feature].mean())def prepare_data(series, lag_start, lag_end, test_size, target_encoding=False, days_to_predict=1):"""series: pd.DataFramedataframe with timeserieslag_start: intinitial step back in time to slice target variableexample - lag_start = 1 means that the modelwill see yesterday's values to predict todaylag_end: intfinal step back in time to slice target variableexample - lag_end = 4 means that the modelwill see up to 4 days back in time to predict todaytest_size: floatsize of the test dataset after train/test split as percentage of datasettarget_encoding: booleanif True - add target averages to the dataset"""last_date = series["timestamp"].max()def make_future_date(periods, freq='D'):"""Simulate the trend using the extrapolated generative model.Parameters----------periods: Int number of periods to forecast forward.freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.Returns-------pd.Dataframe that extends forward from the end of self.history for therequested number of periods."""dates = pd.date_range(start=last_date,periods=periods + 1,  # An extra in case we include startfreq=freq)dates = dates[dates > last_date]  # Drop start if equals last_datereturn dates[:periods]  # Return correct number of periodspredict_points = days_to_predict * 1440 # 1 day = 60*24 minutesfuture_dates = make_future_date(periods=predict_points, freq='T')df_future = pd.DataFrame({"timestamp": future_dates, "y": np.zeros(len(future_dates))})data = pd.concat([series, df_future])data.set_index("timestamp", drop=True, inplace=True)# data = pd.DataFrame(series.copy())# data.columns = ["y"]print(data[:30])# lags of seriesfor i in range(lag_start, lag_end):data["lag_{}".format(i)] = data.y.shift(i)print(data[:30])# datetime features# data.index = data.index.to_datetime()data["hour"] = data.index.hourdata["weekday"] = data.index.weekdaydata['is_weekend'] = data.weekday.isin([5, 6]) * 1if target_encoding:# calculate averages on train set onlytest_index = int(len(data.dropna()) * (1 - test_size))data['weekday_average'] = list(map(code_mean(data[:test_index], 'weekday', "y").get, data.weekday))data["hour_average"] = list(map(code_mean(data[:test_index], 'hour', "y").get, data.hour))# drop encoded variables# data.drop(["hour", "weekday"], axis=1, inplace=True)# train-test splity = data.dropna().yX = data.dropna().drop(['y'], axis=1)X_train, X_test, y_train, y_test, X_predict, y_predict = \timeseries_train_test_split(X, y, test_size=test_size, predict_size=predict_points)return X_train, X_test, y_train, y_test, X_predict, y_predictdef calculate_intevals(lr, X_train, y_train, tscv):cv = cross_val_score(lr, X_train, y_train,cv=tscv,scoring="neg_mean_absolute_error")mae = cv.mean() * (-1)deviation = cv.std()scale = 30return  mae + scale * deviationdef plt_linear():data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])# input formatdata['timestamp'] = pd.to_datetime(data['timestamp'])data = data.sort_values('timestamp')data.rename(columns={'count': 'y'}, inplace=True)lag_start = 1lag_end = 100X_train, X_test, y_train, y_test, X_predict, y_predict = \prepare_data(data, lag_start=lag_start, lag_end=lag_end, test_size=0.3, target_encoding=True)scaler = StandardScaler()X_train_scaled = scaler.fit_transform(X_train)X_test_scaled = scaler.transform(X_test)# for time-series cross-validation set 5 foldstscv = TimeSeriesSplit(n_splits=5)# lr = LinearRegression()lr = LassoCV(cv=tscv)# lr = RidgeCV(cv=tscv)lr.fit(X_train_scaled, y_train)intervals = calculate_intevals(lr, X_train, y_train, tscv)"""y_test_predict = lr.predict(X_test_scaled)plot_results(y_predict=y_test_predict, y=y_test, intervals=intervals, img_filename="linear-test-result.png", plot_intervals=True, plot_anomalies=True)"""y2 = lr.predict(np.concatenate((X_train_scaled, X_test_scaled)))y = pd.concat([y_train, y_test])# plot_results(y_predict=y2, y=y, intervals=intervals, img_filename="linear-all-result.png", plot_intervals=True, plot_anomalies=True)y_future = predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler)plot_results(y_predict=y2, y=y, intervals=intervals, img_filename="linear-all.png", plot_intervals=True, plot_anomalies=True, extra_plot=y_future)plot_arg_importance(lr, X_train=X_train)plt_linear()

# coding: utf-8import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import TimeSeriesSplit
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV, RidgeCV
from sklearn.ensemble import GradientBoostingRegressordef mean_absolute_percentage_error(y_true, y_pred):return np.mean(np.abs((y_true - y_pred) / y_true)) * 100def predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler):# for predicty_predict[0:lag_start] = lr.predict(scaler.transform(X_predict.iloc[0:lag_start]))for i in range(lag_start, len(X_predict)):last_line = X_predict.iloc[i-1]index = X_predict.index[i]for j in range(lag_end-1, lag_start, -1):X_predict.at[index, "lag_{}".format(j)] = last_line["lag_{}".format(j-1)]X_predict.at[index, "lag_{}".format(lag_start)] = y_predict[i-1]y_predict[i] = lr.predict(scaler.transform([X_predict.iloc[i]]))[0]return y_predictdef plot_results(y_predict, y, intervals, img_filename, plot_intervals=False, plot_anomalies=False, extra_plot=None):"""Plots modelled vs fact values, prediction intervals and anomalies"""assert len(y_predict) == len(y)plt.figure(figsize=(15, 7))# plt.plot(y.index, y_predict, "g", label="prediction", linewidth=3.0)# plt.plot(y.index, y.values, label="actual", linewidth=1.0)plt.plot(y.index, y_predict, ls='-', c='#0072B2', label='predicted y')plt.plot(y.index, y.values, 'k.', label='y')if extra_plot is not None:# plt.plot(extra_plot.index, extra_plot.values, "y", label="future predict", linewidth=3.0)plt.plot(extra_plot.index, extra_plot.values, 'y', label='predicted y')if plot_intervals:lower = y_predict - intervalsupper = y_predict + intervals# plt.plot(y.index, lower, "r--", label="upper bond / lower bond", alpha=0.5)# plt.plot(y.index, upper, "r--", alpha=0.5)plt.fill_between(y.index, lower, upper, color='#0072B2', alpha=0.2, label='predicted upper/lower y')if extra_plot is not None:# plt.plot(extra_plot.index, extra_plot.values-intervals, "r--", label="upper bond / lower bond", alpha=0.5)# plt.plot(extra_plot.index, extra_plot.values+intervals, "r--", alpha=0.5)plt.fill_between(extra_plot.index, extra_plot.values-intervals, extra_plot.values+intervals,color='#0072B2', alpha=0.2, label='predicted upper/lower y')if plot_anomalies:anomalies_lower = y[y < lower]anomalies_upper = y[y > upper]# plt.plot(anomalies_lower.index, anomalies_lower.values, "ro", markersize=10, label="Anomalies(+)")# plt.plot(anomalies_upper.index, anomalies_upper.values, "ro", markersize=10, label="Anomalies(-)")plt.plot(anomalies_lower.index, anomalies_lower.values, "rX", label='abnormal points')plt.plot(anomalies_upper.index, anomalies_upper.values, "rX")error = mean_absolute_percentage_error(y_predict, y)plt.title("Mean absolute percentage error {0:.2f}%".format(error))plt.legend(loc="best")plt.tight_layout()plt.grid(True)plt.savefig(img_filename)def plot_arg_importance(model, X_train, img_filename="linear-cov.png"):"""Plots sorted coefficient values of the model"""coefs = pd.DataFrame(model.coef_, X_train.columns)coefs.columns = ["coef"]coefs["abs"] = coefs.coef.apply(np.abs)coefs = coefs.sort_values(by="abs", ascending=False).drop(["abs"], axis=1)plt.figure(figsize=(20, 12))coefs.coef.plot(kind='bar')plt.grid(True, axis='y')plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles='dashed')plt.savefig(img_filename)def code_mean(data, cat_feature, real_feature):"""Returns a dictionary where keys are unique categories of the cat_feature,and values are means over real_feature"""return dict(data.groupby(cat_feature)[real_feature].mean())def prepare_data(series, lag_start, lag_end, test_size, target_encoding=False, days_to_predict=2):"""series: pd.DataFramedataframe with timeserieslag_start: intinitial step back in time to slice target variableexample - lag_start = 1 means that the modelwill see yesterday's values to predict todaylag_end: intfinal step back in time to slice target variableexample - lag_end = 4 means that the modelwill see up to 4 days back in time to predict todaytest_size: floatsize of the test dataset after train/test split as percentage of datasettarget_encoding: booleanif True - add target averages to the dataset"""last_date = series["timestamp"].max()def make_future_date(periods, freq='D'):"""Simulate the trend using the extrapolated generative model.Parameters----------periods: Int number of periods to forecast forward.freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.Returns-------pd.Dataframe that extends forward from the end of self.history for therequested number of periods."""dates = pd.date_range(start=last_date,periods=periods + 1,  # An extra in case we include startfreq=freq)dates = dates[dates > last_date]  # Drop start if equals last_datereturn dates[:periods]  # Return correct number of periodspredict_points = days_to_predict * 1440 # 1 day = 60*24 minutesfuture_dates = make_future_date(periods=predict_points, freq='T')df_future = pd.DataFrame({"timestamp": future_dates, "y": np.zeros(len(future_dates))})data = pd.concat([series, df_future])data.set_index("timestamp", drop=True, inplace=True)# data = pd.DataFrame(series.copy())# data.columns = ["y"]# print(data[:30])# lags of seriesfor i in range(lag_start, lag_end):data["lag_{}".format(i)] = data.y.shift(i)# print(data[:30])# datetime features# data.index = data.index.to_datetime()data["hour"] = data.index.hourdata["weekday"] = data.index.weekdaydata['is_weekend'] = data.weekday.isin([5, 6]) * 1test_index = int(len(series) * (1 - test_size))if target_encoding:# calculate averages on train set onlydata['weekday_average'] = list(map(code_mean(data[:test_index], 'weekday', "y").get, data.weekday))data["hour_average"] = list(map(code_mean(data[:test_index], 'hour', "y").get, data.hour))# drop encoded variables# data.drop(["hour", "weekday"], axis=1, inplace=True)# train-test splity = data.dropna().yX = data.dropna().drop(['y'], axis=1)total_size = len(X)# get the index after which test set startsX_train = X.iloc[:test_index]y_train = y.iloc[:test_index]X_test = X.iloc[test_index:total_size-predict_points]y_test = y.iloc[test_index:total_size-predict_points]X_predict = X.iloc[-predict_points:]y_predict = y.iloc[-predict_points:]return X_train, X_test, y_train, y_test, X_predict, y_predictdef mean_absolute_error(y_true, y_pred):abs_err = np.abs(y_true-y_pred)return np.mean(abs_err), np.std(abs_err)def calculate_intervals2(y_true, y_pred, scale):mae, std = mean_absolute_error(y_true, y_pred)return mae + scale * stddef calculate_intervals(lr, X_train, y_train, tscv, scale):cv = cross_val_score(lr, X_train, y_train,cv=tscv,scoring="neg_mean_squared_error")mae = cv.mean() * (-1)deviation = cv.std()return mae + scale * deviationdef linear_predict(data_frame, interval_scale=30, lag_start=60, lag_end=100, days_to_predict=2):"""predict time series data using linear model.:param data_frame:  input data frame. Must have timestamp and y columns. Also set index with timestamp.:param interval_scale:  interval range scale.:param lag_start:initial step back in time to slice target variableexample - lag_start = 1 means that the modelwill see yesterday's values to predict today:param lag_end:final step back in time to slice target variableexample - lag_end = 4 means that the modelwill see up to 4 days back in time to predict today:param days_to_predict: predicted days in future/:return:history_predict: history predicted value (pandas.Series)y_future: predicted value for future (pandas.Series)intervals: upper and lower range (float)."""assert "timestamp" in data_frame.columnsassert "y" in data_frame.columnsX_train, X_test, y_train, y_test, X_predict, y_predict = \prepare_data(data_frame, lag_start=lag_start, lag_end=lag_end, test_size=0.3, target_encoding=True, days_to_predict=days_to_predict)scaler = StandardScaler()X_train_scaled = scaler.fit_transform(X_train)X_test_scaled = scaler.transform(X_test)# for time-series cross-validation set 5 foldstscv = TimeSeriesSplit(n_splits=5)# lr = LinearRegression()lr = GradientBoostingRegressor(n_estimators=100, learning_rate=0.1, max_depth=1, random_state=666)# lr = LassoCV(cv=tscv)# lr = RidgeCV(cv=tscv)lr.fit(X_train_scaled, y_train)# plot_arg_importance(lr, X_train=X_train, img_filename="linear-cov.png")y_future = predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler)y = lr.predict(np.concatenate((X_train_scaled, X_test_scaled)))y_history = pd.concat([y_train, y_test])# intervals = calculate_intervals(lr, X_train, y_train, tscv, scale=interval_scale)intervals = calculate_intervals2(y_history, y, interval_scale)# anomalies_lower = y_history[y_history<y-intervals]# anomalies_upper = y_history[y_history>y+intervals]assert len(y_history.index) == len(y)return pd.Series(data=y, index=y_history.index, name="history_predict"), y_future, intervalsdef get_anoms(y_real, y_predict, intervals, lower_ratio=1.0, upper_ratio=1.0):"""calculat anomal point using predicted value and interval:param y_real:  real value (pandas.Series):param y_predict: predicted value (pandas.Series):param intervals:  upper and lower bound range (float):param lower_ratio: lower ratio you want to scale:param upper_ratio:  upper ratio you want scale:return:anoms_lower, anoms_upper (pandas.Series)"""anomalies_lower_index,anomalies_lower_val = [], []anomalies_upper_index,anomalies_upper_val = [], []for timestamp, expect_val in zip(y_predict.index, y_predict.values):real_val = y_real.loc[timestamp]if (expect_val-intervals) * lower_ratio > real_val:anomalies_lower_index.append(timestamp)anomalies_lower_val.append(real_val)if (expect_val+intervals) * upper_ratio < real_val:anomalies_upper_index.append(timestamp)anomalies_upper_val.append(real_val)return pd.Series(data=anomalies_lower_val, index=anomalies_lower_index),\pd.Series(data=anomalies_upper_val, index=anomalies_upper_index)def plot_history_and_future(y_predict, y_real, intervals, anomalies_lower, anomalies_upper, predicted_future, img_filename, need_lower=True):"""Plots modelled vs fact values, prediction intervals and anomalies"""assert len(y_predict) < len(y_real)plt.figure(figsize=(15, 7))plt.plot(y_predict.index, y_predict.values, ls='-', c='#0072B2', label='predicted y')plt.plot(y_real.index, y_real.values, 'k.', label='y')if need_lower:plt.fill_between(y_predict.index, y_predict.values - intervals, y_predict.values + intervals, color='#0072B2', alpha=0.2, label='predicted upper/lower y')else:plt.fill_between(y_predict.index, 0, y_predict.values + intervals, color='#0072B2', alpha=0.2, label='predicted upper/lower y')plt.plot(predicted_future.index, predicted_future.values, 'y', label='predicted y')if need_lower:plt.fill_between(predicted_future.index, predicted_future.values-intervals, predicted_future.values+intervals,color='#0072B2', alpha=0.2)else:plt.fill_between(predicted_future.index, 0, predicted_future.values+intervals,color='#0072B2', alpha=0.2)if need_lower:plt.plot(anomalies_lower.index, anomalies_lower.values, "rX", label='abnormal points')plt.plot(anomalies_upper.index, anomalies_upper.values, "rX")else:plt.plot(anomalies_upper.index, anomalies_upper.values, "rX", label="abnormal points")error = mean_absolute_percentage_error(y_predict, y_real)plt.title("Mean absolute percentage error {0:.2f}%".format(error))plt.legend(loc="best")plt.tight_layout()plt.grid(True)plt.savefig(img_filename)if __name__ == "__main__":data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])# input formatdata['timestamp'] = pd.to_datetime(data['timestamp'])data = data.sort_values('timestamp')data.rename(columns={'count': 'y'}, inplace=True)data.set_index("timestamp", drop=False, inplace=True)y_predict, y_future, intervals = linear_predict(data, interval_scale=5, lag_start=60, lag_end=100, days_to_predict=3)anomalies_lower, anomalies_upper = get_anoms(data['y'], y_predict, intervals)plot_history_and_future(y_predict=y_predict, y_real=data['y'], intervals=intervals, anomalies_lower=anomalies_lower,anomalies_upper=anomalies_upper, predicted_future=y_future, img_filename="linear.png", need_lower=True)