基于SARIMA、XGBoost和CNN-LSTM的时间序列预测对比

 from statsmodels.tsa.stattools import adfuller
 
 result = adfuller(plant2_dcpower.values)
 
 print('ADF Statistic: %f' % result[0])
 print('p-value: %f' % result[1])
 print('Critical Values:')
 for key, value in result[4].items():
    print('\t%s: %.3f' % (key, value))
1.
2.
3.
4.
5.
6.
7.
8.
9.

from time import time
 from sklearn.metrics import mean_squared_error
 from statsmodels.tsa.statespace.sarimax import SARIMAX
 
 configg = [(2, 1, 4), (2, 1, 6, 24), 'n']
 
 def train_test_split(data, test_len=48):
    """
    Split data into training and testing.
    """
    train, test = data[:-test_len], data[-test_len:]
    return train, test
 
 def sarima_model(data, cfg, test_len, i):
    """
    SARIMA model which outputs prediction and model.
    """
    order, s_order, t = cfg[0], cfg[1], cfg[2]
    model = SARIMAX(data, order=order, seasonal_order=s_order, trend=t,
                    enforce_stationarity=False, enfore_invertibility=False)
    model_fit = model.fit(disp=False)
    yhat = model_fit.predict(len(data))
     
    if i + 1 == test_len:
        return yhat, model_fit
    else:
        return yhat
 
 def walk_forward_val(data, cfg):
    """
    A walk forward validation technique used for time series data. Takes current value of x_test and predicts
    value. x_test is then fed back into history for the next prediction.
    """
    train, test = train_test_split(data)
    pred = []
    history = [i for i in train]
    test_len = len(test)
     
    for i in range(test_len):
        if i + 1 == test_len:
            yhat, s_model = sarima_model(history, cfg, test_len, i)
            pred.append(yhat)
            mse = mean_squared_error(test, pred)
            return pred, mse, s_model
        else:
            yhat = sarima_model(history, cfg, test_len, i)
            pred.append(yhat)
            history.append(test[i])
    pass
 
 if __name__ == '__main__':
    start_time = time()
    sarima_pred_plant2, sarima_mse, s_model = walk_forward_val(plant2_dcpower, configg)
    time_len = time() - start_time
 
    print(f'SARIMA runtime: {round(time_len/60,2)} mins')
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.

 import numpy as np
 import pandas as pd
 import xgboost as xgb
 from sklearn.preprocessing import MinMaxScaler
 from time import time
 
 def train_test_split(df, test_len=48):
    """
    split data into training and testing.
    """
    train, test = df[:-test_len], df[-test_len:]
    return train, test
 
 def data_to_supervised(df, shift_by=1, target_var='DC_POWER'):
    """
    Convert data into a supervised learning problem.
    """
    target = df[target_var][shift_by:].values
    dep = df.drop(target_var, axis=1).shift(-shift_by).dropna().values
    data = np.column_stack((dep, target))
    return data
 
 
 def xgb_forecast(train, x_test):
    """
    XGBOOST model which outputs prediction and model.
    """
    x_train, y_train = train[:,:-1], train[:,-1]
    xgb_model = xgb.XGBRegressor(learning_rate=0.01, n_estimators=1500, subsample=0.8,
                                  colsample_bytree=1, colsample_bylevel=1,
                                  min_child_weight=20, max_depth=14, objective='reg:squarederror')
    xgb_model.fit(x_train, y_train)
    yhat = xgb_model.predict([x_test])
    return yhat[0], xgb_model
 
 def walk_forward_validation(df):
    """
    A walk forward validation approach by scaling the data and changing into a supervised learning problem.
    """
    preds = []
    train, test = train_test_split(df)
     
    scaler = MinMaxScaler(feature_range=(0,1))
    train_scaled = scaler.fit_transform(train)
    test_scaled = scaler.transform(test)
     
    train_scaled_df = pd.DataFrame(train_scaled, columns = train.columns, index=train.index)
    test_scaled_df = pd.DataFrame(test_scaled, columns = test.columns, index=test.index)
     
    train_scaled_sup, test_scaled_sup = data_to_supervised(train_scaled_df), data_to_supervised(test_scaled_df)
    history = np.array([x for x in train_scaled_sup])
     
    for i in range(len(test_scaled_sup)):
        test_x, test_y = test_scaled_sup[i][:-1], test_scaled_sup[i][-1]
        yhat, xgb_model = xgb_forecast(history, test_x)
        preds.append(yhat)
        np.append(history,[test_scaled_sup[i]], axis=0)
         
    pred_array = test_scaled_df.drop("DC_POWER", axis=1).to_numpy()
    pred_num = np.array([pred])
    pred_array = np.concatenate((pred_array, pred_num.T), axis=1)
    result = scaler.inverse_transform(pred_array)
     
    return result, test, xgb_model
 
 if __name__ == '__main__':
    start_time = time()
    xgb_pred, actual, xgb_model = walk_forward_validation(dropped_df_cat)
    time_len = time() - start_time
 
    print(f'XGBOOST runtime: {round(time_len/60,2)} mins')
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.
57.
58.
59.
60.
61.
62.
63.
64.
65.
66.
67.
68.
69.
70.
71.

import pandas as pd
 import numpy as np
 from sklearn.metrics import mean_squared_error
 from sklearn.preprocessing import MinMaxScaler
 import keras
 from keras.models import Sequential
 from keras.layers.convolutional import Conv1D, MaxPooling1D
 from keras.layers import LSTM, TimeDistributed, RepeatVector, Dense, Flatten
 from keras.optimizers import Adam
 
 n_steps = 1
 subseq = 1
 
 def train_test_split(df, test_len=48):
    """
    Split data in training and testing. Use 48 hours as testing.
    """
    train, test = df[:-test_len], df[-test_len:]
    return train, test
 
 def split_data(sequences, n_steps):
    """
    Preprocess data returning two arrays.
    """
    x, y = [], []
    for i in range(len(sequences)):
        end_x = i + n_steps
         
        if end_x > len(sequences):
            break
        x.append(sequences[i:end_x, :-1])
        y.append(sequences[end_x-1, -1])
         
    return np.array(x), np.array(y)
   
  def CNN_LSTM(x, y, x_val, y_val):
    """
    CNN-LSTM model.
    """
    model = Sequential()
    model.add(TimeDistributed(Conv1D(filters=14, kernel_size=1, activation="sigmoid",
                                      input_shape=(None, x.shape[2], x.shape[3]))))
    model.add(TimeDistributed(MaxPooling1D(pool_size=1)))
    model.add(TimeDistributed(Flatten()))
    model.add(LSTM(21, activation="tanh", return_sequences=True))
    model.add(LSTM(14, activation="tanh", return_sequences=True))
    model.add(LSTM(7, activation="tanh"))
    model.add(Dense(3, activation="sigmoid"))
    model.add(Dense(1))
     
    model.compile(optimizer=Adam(learning_rate=0.001), loss="mse", metrics=['mse'])
    history = model.fit(x, y, epochs=250, batch_size=36,
                        verbose=0, validation_data=(x_val, y_val))
     
    return model, history
   
 # split and resahpe data
 train, test = train_test_split(dropped_df_cat)  
 
 train_x = train.drop(columns="DC_POWER", axis=1).to_numpy()
 train_y = train["DC_POWER"].to_numpy().reshape(len(train), 1)
 
 test_x = test.drop(columns="DC_POWER", axis=1).to_numpy()
 test_y = test["DC_POWER"].to_numpy().reshape(len(test), 1)
 
 #scale data  
 scaler_x = MinMaxScaler(feature_range=(-1,1))
 scaler_y = MinMaxScaler(feature_range=(-1,1))
 
 train_x = scaler_x.fit_transform(train_x)
 train_y = scaler_y.fit_transform(train_y)
 
 test_x = scaler_x.transform(test_x)
 test_y = scaler_y.transform(test_y)
 
 
 # shape data into CNN-LSTM format [samples, subsequences, timesteps, features] ORIGINAL
 train_data_np = np.hstack((train_x, train_y))
 x, y = split_data(train_data_np, n_steps)
 x_subseq = x.reshape(x.shape[0], subseq, x.shape[1], x.shape[2])
 
 # create validation set
 x_val, y_val = x_subseq[-24:], y[-24:]
 x_train, y_train = x_subseq[:-24], y[:-24]
 
 n_features = x.shape[2]
 actual = scaler_y.inverse_transform(test_y)
 
 # run CNN-LSTM model
 if __name__ == '__main__':
    start_time = time()
 
    model, history = CNN_LSTM(x_train, y_train, x_val, y_val)
    prediction = []
     
    for i in range(len(test_x)):
        test_input = test_x[i].reshape(1, subseq, n_steps, n_features)
        yhat = model.predict(test_input, verbose=0)
        yhat_IT = scaler_y.inverse_transform(yhat)
        prediction.append(yhat_IT[0][0])
 
    time_len = time() - start_time
    mse = mean_squared_error(actual.flatten(), prediction)
     
    print(f'CNN-LSTM runtime: {round(time_len/60,2)} mins')
    print(f"CNN-LSTM MSE: {round(mse,2)}")
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.
57.
58.
59.
60.
61.
62.
63.
64.
65.
66.
67.
68.
69.
70.
71.
72.
73.
74.
75.
76.
77.
78.
79.
80.
81.
82.
83.
84.
85.
86.
87.
88.
89.
90.
91.
92.
93.
94.
95.
96.
97.
98.
99.
100.
101.
102.
103.
104.
105.
106.