So I've been playing with some different forecasting methods on a data set that I have done some more basic analyses for in the past. Without going into to much detail, it's population data over time driven by a variety of unknown factors. I don't expect to be able to predict it very well, and, historically, I've seen MAPE values in the 25% range using more naive models like ARIMA (though even Prophet models haven't worked great). So today, for fun, I decided to try the LSTM code from https://machinelearningmastery.com/time-series-prediction-lstm-recurrent-neural-networks-python-keras/ to see how it'd do. Suffice it to say, it performed way too well. So well, I assume I've done something wrong. For fun, I decided to restrict the number of parameters to the smallest possible number. This model literally has 9 parameters:
I disabled the bias parameters on the LSTM and Dense layers, and I specified a linear activation function. I don't think it's possible to have LSTM+Dense with fewer parameters.
To make it harder on the model, I trained it for 1 step and set the lookback to 1. The resulting code looks like this:
# LSTM for international airline passengers problem with memory import numpy import matplotlib.pyplot as plt from pandas import read_csv import math from keras.models import Sequential from keras.layers import Dense from keras.layers import LSTM from sklearn.preprocessing import MinMaxScaler from sklearn.metrics import mean_squared_error # convert an array of values into a dataset matrix def create_dataset(ds, lb=1): dataX, dataY = ,  for i in range(len(ds)-lb-1): a = ds[i:(i+lb), 0] dataX.append(a) dataY.append(ds[i + lb, 0]) return numpy.array(dataX), numpy.array(dataY) # fix random seed for reproducibility numpy.random.seed(7) # load the dataset dataframe = df_forecasting dataset = dataframe.values dataset = dataset.astype('float32') # normalize the dataset scaler = MinMaxScaler(feature_range=(0, 1)) scaler.fit([, ]) # Maximum expected population is ~50 for problem-related reasons dataset = scaler.transform(dataset) # split into train and test sets train_size = list(df_forecasting.index > pd.to_datetime('20200101')).index(True) # I want to see about forecasting 2020 test_size = len(dataset) - train_size train, test = dataset[0:train_size,:], dataset[train_size:len(dataset),:] # reshape into X=t and Y=t+1 look_back = 1 trainX, trainY = create_dataset(train, look_back) testX, testY = create_dataset(test, look_back) # reshape input to be [samples, time steps, features] trainX = numpy.reshape(trainX, (trainX.shape, trainX.shape, 1)) testX = numpy.reshape(testX, (testX.shape, testX.shape, 1)) # create and fit the LSTM network batch_size = 1 model = Sequential() model.add(LSTM(1, batch_input_shape=(batch_size, look_back, 1), stateful=True, use_bias=False)) model.add(Dense(1, activation='linear', use_bias=False)) model.compile(loss='mean_squared_error', optimizer='adam') for i in range(1): model.fit(trainX, trainY, epochs=1, batch_size=batch_size, verbose=2, shuffle=False) model.reset_states() # make predictions trainPredict = model.predict(trainX, batch_size=batch_size) model.reset_states() testPredict = model.predict(testX, batch_size=batch_size) # invert predictions trainPredict = scaler.inverse_transform(trainPredict) trainY = scaler.inverse_transform([trainY]) testPredict = scaler.inverse_transform(testPredict) testY = scaler.inverse_transform([testY]) # calculate root mean squared error trainScore = math.sqrt(mean_squared_error(trainY, trainPredict[:,0])) print('Train Score: %.2f RMSE' % (trainScore)) testScore = math.sqrt(mean_squared_error(testY, testPredict[:,0])) print('Test Score: %.2f RMSE' % (testScore)) # shift train predictions for plotting trainPredictPlot = numpy.empty_like(dataset) trainPredictPlot[:, :] = numpy.nan trainPredictPlot[look_back:len(trainPredict)+look_back, :] = trainPredict # shift test predictions for plotting testPredictPlot = numpy.empty_like(dataset) testPredictPlot[:, :] = numpy.nan testPredictPlot[len(trainPredict)+(look_back*2)+1:len(dataset)-1, :] = testPredict # plot baseline and predictions plt.figure(figsize=(15, 10)) plt.plot(scaler.inverse_transform(dataset), c="g", label="Original Data") plt.plot(trainPredictPlot, c="b", label="Train Prediction") plt.plot(testPredictPlot, c="r", label="Test Prediction") plt.legend() plt.show() print(model.summary())
I'm stumped. I've traced through the code and I simply don't understand how it's working so well compared to other models I've tried. One theory I have is that this data is actually a moving average, but it's a causal average, so I don't see how that could leak information. Even if it could, 9 parameters can use that?
Would love to understand better what's going on here. The map is 17, which I find incredible compared to the other things I tried. When I un-break the model (101 params, 10 steps, 7day lookback) I get the MAPE down to 9. I think it would go substantially lower if I let it, but I want to believe first...