时间序列预测——线性回归（上下界、异常检测），异常检测时候历史数据的输入选择是关键，使用过去历史值增加模型健壮性

data download:

https://github.com/nicolasmiller/pyculiarity/blob/master/tests/raw_data.csv

数据集样子：

                           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

做了shift处理前后：

                           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.816
                           y    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

代码：

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

# coding: utf-8
 
import 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)) * 100
 
 
def 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 starts
    test_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_test
 
 
def 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 = 20
        lower = 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.DataFrame
            dataframe with timeseries
        lag_start: int
            initial step back in time to slice target variable
            example - lag_start = 1 means that the model
                      will see yesterday's values to predict today
        lag_end: int
            final step back in time to slice target variable
            example - lag_end = 4 means that the model
                      will see up to 4 days back in time to predict today
        test_size: float
            size of the test dataset after train/test split as percentage of dataset
        target_encoding: boolean
            if True - add target averages to the dataset
 
    """
    # copy of the initial dataset
    data = pd.DataFrame(series.copy())
    data.columns = ["y"]
 
    # lags of series
    for 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.hour
    data["weekday"] = data.index.weekday
    data['is_weekend'] = data.weekday.isin([5, 6]) * 1
 
    if target_encoding:
        # calculate averages on train set only
        test_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 split
    y = data.dropna().y
    X = 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_test
 
 
def 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 XGBRegressor
    lr = 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)
 
    """
    IMPORTANT
    Generally 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()

可以看到相关系数！

重构代码，使其可以预测未来：

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

# coding: utf-8
 
import 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
 
 
def mean_absolute_percentage_error(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
 
 
def 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 starts
    test_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_predict
 
 
def predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler):
    # for predict
    # not OK for abnormal real value
    y_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_predict
 
 
def 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 - intervals
        upper = 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.DataFrame
            dataframe with timeseries
        lag_start: int
            initial step back in time to slice target variable
            example - lag_start = 1 means that the model
                      will see yesterday's values to predict today
        lag_end: int
            final step back in time to slice target variable
            example - lag_end = 4 means that the model
                      will see up to 4 days back in time to predict today
        test_size: float
            size of the test dataset after train/test split as percentage of dataset
        target_encoding: boolean
            if 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 the
        requested number of periods.
        """
        dates = pd.date_range(
            start=last_date,
            periods=periods + 1,  # An extra in case we include start
            freq=freq)
        dates = dates[dates > last_date]  # Drop start if equals last_date
        return dates[:periods]  # Return correct number of periods
 
    predict_points = days_to_predict * 1440 # 1 day = 60*24 minutes
 
    future_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 series
    for 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.hour
    data["weekday"] = data.index.weekday
    data['is_weekend'] = data.weekday.isin([5, 6]) * 1
 
    if target_encoding:
        # calculate averages on train set only
        test_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 split
    y = data.dropna().y
    X = 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_predict
 
 
def 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 = 30
    return  mae + scale * deviation
 
 
def plt_linear():
    data = pd.read_csv('raw_data.csv', usecols=['timestamp', 'count'])
    # input format
    data['timestamp'] = pd.to_datetime(data['timestamp'])
    data = data.sort_values('timestamp')
    data.rename(columns={'count': 'y'}, inplace=True)
 
    lag_start = 1
    lag_end = 100
    X_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 folds
    tscv = 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()

绘图：

尤其关键的是，lag_start, lag_end参数，如果lag_start=1，则表示使用前一个时刻输入，会导致模型过拟合，因此上面设置了lag_start=60表示使用一个小时前的数据来预测，防止过拟合。

来lag_start=1的效果：

可以看到，前一个时刻数据值的重要性！因此，最后做趋势预测的时候出现了重大失误，模型过拟合了。

bug修复：

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

# coding: utf-8
 
import 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 GradientBoostingRegressor
 
 
 
def mean_absolute_percentage_error(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
 
 
def predict_future(lr, X_predict, y_predict, lag_start, lag_end, scaler):
    # for predict
    y_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_predict
 
 
def 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 - intervals
        upper = 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.DataFrame
            dataframe with timeseries
        lag_start: int
            initial step back in time to slice target variable
            example - lag_start = 1 means that the model
                      will see yesterday's values to predict today
        lag_end: int
            final step back in time to slice target variable
            example - lag_end = 4 means that the model
                      will see up to 4 days back in time to predict today
        test_size: float
            size of the test dataset after train/test split as percentage of dataset
        target_encoding: boolean
            if 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 the
        requested number of periods.
        """
        dates = pd.date_range(
            start=last_date,
            periods=periods + 1,  # An extra in case we include start
            freq=freq)
        dates = dates[dates > last_date]  # Drop start if equals last_date
        return dates[:periods]  # Return correct number of periods
 
    predict_points = days_to_predict * 1440 # 1 day = 60*24 minutes
 
    future_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 series
    for 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.hour
    data["weekday"] = data.index.weekday
    data['is_weekend'] = data.weekday.isin([5, 6]) * 1
 
    test_index = int(len(series) * (1 - test_size))
    if target_encoding:
        # calculate averages on train set only
        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 split
    y = data.dropna().y
    X = data.dropna().drop(['y'], axis=1)
 
 
    total_size = len(X)
    # get the index after which test set starts
    X_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_predict
 
 
def 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 * std
 
 
def 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 * deviation
 
 
def 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 variable
            example - lag_start = 1 means that the model
                      will see yesterday's values to predict today
    :param lag_end:
            final step back in time to slice target variable
            example - lag_end = 4 means that the model
                      will 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.columns
    assert "y" in data_frame.columns
 
    X_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 folds
    tscv = 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, intervals
 
 
def 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 format
    data['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)　　