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))

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')

 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')

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)}")