本文主要參考了 Jason Brownlee 的博文 Time Series Prediction With Deep Learning in Keras
原文使用 python 實現模型，這里是用 R
基於 Keras 用深度學習預測時間序列
時間序列預測一直以來是機器學習中的一個難題。
在本篇文章中，將介紹如何在 R 中使用 keras 深度學習包構建神經網絡模型實現時間序列預測。
文章的主要內容：
- 如何將時間序列預測問題表示成為一個回歸問題，並建立對應的神經網絡模型。
- 如何使用滯后時間的數據實現時間序列預測，並建立對應的神經網絡模型。
問題描述
“航班旅客數據”是一個常用的時間序列數據集，該數據包含了 1949 至 1960 年 12 年間的月度旅客數據，共有 144 個觀測值。
下載鏈接：international-airline-passengers.csv
多層感知機回歸
時間序列預測中最簡單的思路之一便是尋找當前和過去數據（\(X_t, X_{t-1}, \dots\)）與未來數據（$ X_{t+1}$）之間的關系，這種關系通常會表示成為一個回歸問題。
下面着手將時間序列預測問題表示成一個回歸問題，並建立神經網絡模型用於預測。
首先，加載相關 R 包。
library(keras)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(lubridate)
神經網絡模型在訓練時存在一定的隨機性，所以要為計算統一隨機數環境
set.seed(7)
畫出整體數據的曲線圖，對問題有一個直觀的認識。
dataframe <- read.csv(
'international-airline-passengers.csv')
dataframe$Month <- paste0(dataframe$Month,'-01') %>%
ymd()
ggplot(
data = dataframe,
mapping = aes(
x = Month,
y = passengers)) +
geom_line() +
geom_point() +
theme_economist() +
scale_color_economist()
圖1
很顯然，數據體現出“季節性”，同時存在線性增長和波動水平增大的趨勢。
將數據集分成兩部分：訓練集和測試集，比例分別占數據集的 2/3 和 1/3。
dataset <- dataframe$passengers
train_size <- as.integer(length(dataset) * 0.67)
test_size <- length(dataset) - train_size
train <- dataset[1:train_size]
test <- dataset[(train_size + 1):length(dataset)]
cat(length(train), length(test))
96 48
為訓練神經網絡對數據做預處理，用數據構造出兩個矩陣，分別是“歷史數據”（作為預測因子）和“未來數據”（作為預測目標）。這里用最近一個月的歷史數據做預測。
create_dataset <- function(dataset,
look_back = 1)
{
l <- length(dataset)
dataX <- matrix(nrow = l - look_back, ncol = look_back)
for (i in 1:ncol(dataX))
{
dataX[, i] <- dataset[i:(l - look_back + i - 1)]
}
dataY <- matrix(
data = dataset[(look_back + 1):l],
ncol = 1)
return(
list(
dataX = dataX,
dataY = dataY))
}
look_back <- 1
trainXY <- create_dataset(train, look_back)
testXY <- create_dataset(test, look_back)
下面構造神經網絡的框架結構並用處理過的訓練數據訓練。
model <- keras_model_sequential()
model %>%
layer_dense(
units = 8,
input_shape = c(look_back),
activation = 'relu') %>%
layer_dense(units = 1) %>%
compile(
loss = 'mean_squared_error',
optimizer = 'adam') %>%
fit(
trainXY$dataX,
trainXY$dataY,
epochs = 200,
batch_size = 2,
verbose = 2)
訓練結果如下。
trainScore <- model %>%
evaluate(
trainXY$dataX,
trainXY$dataY,
verbose = 0)
testScore <- model %>%
evaluate(
testXY$dataX,
testXY$dataY,
verbose = 0)
sprintf(
'Train Score: %.2f MSE (%.2f RMSE)',
trainScore,
sqrt(trainScore))
sprintf(
'Test Score: %.2f MSE (%.2f RMSE)',
testScore,
sqrt(testScore))
[1] "Train Score: 538.50 MSE (23.21 RMSE)"
[1] "Test Score: 2342.33 MSE (48.40 RMSE)"
把訓練數據的擬合值、測試數據的預測值和原始數據畫在一起。
trainPredict <- model %>%
predict(trainXY$dataX)
testPredict <- model %>%
predict(testXY$dataX)
df <- data.frame(
index = 1:length(dataset),
value = dataset,
type = 'raw') %>%
rbind(
data.frame(
index = 1:length(trainPredict) + look_back,
value = trainPredict,
type = 'train')) %>%
rbind(
data.frame(
index = 1:length(testPredict) + look_back + length(train),
value = testPredict,
type = 'test'))
ggplot(data = df) +
geom_line(
mapping = aes(
x = index,
y = value,
color = type)) +
geom_point(
mapping = aes(
x = index,
y = value,
color = type)) +
geom_vline(
xintercept = length(train) + 0.5) +
theme_economist() +
scale_color_economist()
圖2
黑線左邊是訓練部分，右邊是測試部分。
從圖中可以看出，神經網絡模型抓住了數據線性增長和波動率逐漸增加的兩大趨勢，在不做數據轉換的前提下，這是經典的時間序列分析模型不容易做到的；但是很可能沒有識別出“季節性”的結構特點，因為訓練和預測結果和原始數據之間存在“平移錯位”。
多層感知機回歸結合“窗口法”
前面的例子可以看出，如果僅使用\(X_{t-1}\)來預測\(X_t\)，很難讓神經網絡模型識別出“季節性”的結構特征，因此有必要嘗試增加“窗口”寬度，使用更多的歷史數據（包含一個完整的周期）訓練模型。
下面將數 create_dataset 中的參數 look_back 設置為 12，用來包含過去 1 年的歷史數據，重新訓練模型。
look_back <- 12
trainXY <- create_dataset(train, look_back)
testXY <- create_dataset(test, look_back)
model <- keras_model_sequential()
model %>%
layer_dense(
units = 8,
input_shape = c(look_back),
activation = 'relu') %>%
layer_dense(units = 1) %>%
compile(
loss = 'mean_squared_error',
optimizer = 'adam') %>%
fit(
trainXY$dataX,
trainXY$dataY,
epochs = 200,
batch_size = 2,
verbose = 2)
trainScore <- model %>%
evaluate(
trainXY$dataX,
trainXY$dataY,
verbose = 0)
testScore <- model %>%
evaluate(
testXY$dataX,
testXY$dataY,
verbose = 0)
sprintf(
'Train Score: %.2f MSE (%.2f RMSE)',
trainScore,
sqrt(trainScore))
sprintf(
'Test Score: %.2f MSE (%.2f RMSE)',
testScore,
sqrt(testScore))
trainPredict <- model %>%
predict(trainXY$dataX)
testPredict <- model %>%
predict(testXY$dataX)
df <- data.frame(
index = 1:length(dataset),
value = dataset,
type = 'raw') %>%
rbind(
data.frame(
index = 1:length(trainPredict) + look_back,
value = trainPredict,
type = 'train')) %>%
rbind(
data.frame(
index = 1:length(testPredict) + look_back + length(train),
value = testPredict,
type = 'test'))
ggplot(data = df) +
geom_line(
mapping = aes(
x = index,
y = value,
color = type)) +
geom_point(
mapping = aes(
x = index,
y = value,
color = type)) +
geom_vline(
xintercept = length(train) + 0.5) +
theme_economist() +
scale_color_economist()
[1] "Train Score: 157.17 MSE (12.54 RMSE)"
[1] "Test Score: 690.69 MSE (26.28 RMSE)"
圖3
新的模型基本上克服了“平移錯位”的現象，同時依然能夠識別出線性增長和波動率逐漸增加的兩大趨勢。
改進方向
- 目前對“季節性”的識別是靠增加歷史數據實現的，能否從神經網絡結構的方向入手。
- 目前的模型中幾乎沒有用到“特征工程”，如何用特征工程表示數據中存在的主要趨勢和結構化特征。
- DNN + ARIMA：一方作為另外一方的“特征工程”手段。