时间序列预测算法——MQR(C)NN

论文传送门

概述：

类似于DeepAR，MQR$($C$)$NN也是一种可以对批量时间序列统一建模和预测的算法，采用的也是seq2seq模型框架，即encoder-decoder结构。MQR$($C$)$NN代表MQRNN和MQCNN两个算法，两个算法唯一的不同就是encode部分，MQRNN的encode部分用的是RNN，而MQCNN用的是CNN，更确切的说应该是WaveNet，其他部分完全一样，decode部分采用MLP。另外，与DeepAR不同的是，MQR$($C$)$NN利用了fork的机制，且在预测时候并不是递归式的预测，而是一次性的预测decode_len长度的值。

算法原理及实现过程

上图是模型的网络结构。

原理

MQR$($C$)$NN是基于历史encoder_len长度的时间序列值及特征，预测未来decoder_len长度的序列值。用公式表示如下：

其中，$y_{\cdot ,i}$表示第$i$条时间序列；$x_{:t,i}^{(h)}$表示历史的时间特征，作为encoder部分的输入；$x_{:t,i}^{(f)}$表示未来的时间特征，作为decoder部分的输入；$x_{:t,i}^{(s)}$表示一些静态的不变的特征，比如每条时间序列的category编码，或者周期性平稳性特征等。

下面解析模型的网络结构：

关于Encoder，首先该部分可以使用LSTM（模型就叫MQRNN），也可以使用WaveNet（模型就叫MQCNN）。Encoder在每个时间点的输入有两个部分：一是当前时刻的特征，二是序列值，将这些历史输入通过LSTM或WaveNet编码得到最终的hidden_state$h_t$，hidden_state需要传到Deceder部分。另外，上面网络结构图中有一个很容易忽略的点就是在encoder部分的每一个时间点后都接着一个decoder，这是论文中提到一个重要机制：fork机制，代码实现时要重点注意，如果不在每一个时间点decoder，模型训练的loss，mse，rmse等指标会较差。

Decoder部分采用两个MLP，即global_mlp和local_mlp。global_mlp用来summarize encoder部分的输出$h_t$和未来的特征输入到2个contexts：$c_{t+K}$和$c_a$，k=decoder_len，$c_{t+1}$…$c_{t+K}$代表未来每个时间点specific horizon，并分别作为decoder每一个时间点的输入，而ca是计算捕获到的公共的信息horizon-agnostic，并作为decoder所有时间点的输入；local_mlp组合“未来的特征输入”和“来自global_mlp的两个上下文输出”作为输入，最终得到未来每个时间点上在各个分位数上的预测值，用公式表达：

其中$q_{(\cdot )}$代表分位数列表里的每一个分位数，模型训练时，是作为一个参数传入。

关于算法的loss function，结合了分位数回归的思想，MQR$($C$)$NN采用分位数损失函数，即：

模型通过最小化total分位数损失${\Sigma }_t{\Sigma }_q{\Sigma }_k{L}_q(y_{t+k},{y^{\wedge }}_{t+k}^q)$训练模型。

实现过程

原理部分已经对模型的输入，网络如何传递，输出，还有损失函数等做了大致的讲解，其实差不多就是算法的整个实现过程了，下面是encoder,decoder,loss_function的代码实现，代码是从完整的代码摘取出来的，显得有些冗余，但基本上涵盖了该算法的实现逻辑，在下一章(MQRNN代码实现)给大家示例用完整代码训练一个模型。

encoder

# MQRNN的encoder部分，MQCNN可以直接用WaveNet作为encoder
def lstm_cell(layer_num,keep_prob,hidden_size):
    stacked_rnn = []
    for i in range(layer_num):
        lstm = tf.contrib.rnn.LSTMCell(hidden_size, forget_bias=1.0)
        drop = tf.nn.rnn_cell.DropoutWrapper(lstm, output_keep_prob=keep_prob)
        stacked_rnn.append(drop)
        lstm_multi = tf.contrib.rnn.MultiRNNCell(stacked_rnn)
    return lstm_multi

def encoder(encoder_input, encoder_seq_len, keep_prob):
    cells = lstm_cell(keep_prob)
    with tf.variable_scope('encoder', reuse=tf.AUTO_REUSE):
        output, _ = tf.nn.dynamic_rnn(cell=cells, inputs=encoder_input, sequence_length=encoder_seq_len, dtype=tf.float32)
        return output

decoder

def add_layer(inputs, n_hidden, keep_prob):
    '''
    自定义神经层添加函数
    '''
    outputs = tf.contrib.layers.fully_connected(inputs, n_hidden, activation_fn=tf.nn.relu, biases_initializer=tf.zeros_initializer())
    #　outputs_dropout = tf.nn.dropout(outputs, keep_prob)  
    return outputs

def mlp(inputs, mlp_layers, n_hidden, keep_prob):
    for i in range(mlp_layers):
        outputs_dropout = add_layer(inputs, n_hidden, keep_prob)
        inputs = outputs_dropout
    outputs_dropout = inputs
    return outputs_dropout

def global_mlp(encoder_state,decoder_input,time_series_dim,decoder_len,hidden_size,mlp_layers, n_hidden,　keep_prob):
    with tf.variable_scope('global_mlp', reuse=tf.AUTO_REUSE):
        encoder_state = tf.tile(encoder_state, [1, decoder_len, 1])
        global_mlp_input = tf.concat([encoder_state, decoder_input], -1)
        global_mlp_input = tf.reshape(global_mlp_input,　[-1, (hidden_size + time_series_dim) * decoder_len])
        n_hidden = 64 * (decoder_len + 1)
        global_mlp_output = mlp(global_mlp_input, mlp_layers, n_hidden,　keep_prob)
        global_mlp_output = tf.reshape(global_mlp_output,　[-1, (decoder_len + 1), 64])
        ck = tf.slice(global_mlp_output, [0, 0, 0], [-1, decoder_len, -1])
        ca = tf.slice(global_mlp_output, [0, decoder_len, 0], [-1, 1, -1])
    return ck, ca

def local_mlp(inputs,mlp_layers,Q_list_len,keep_prob):
    '''
    Q_list_len:分位数列表的长度，也就是分位数的个数
    '''
    with tf.variable_scope('local_mlp', reuse=tf.AUTO_REUSE):
        quantile_pred = mlp(inputs, mlp_layers, Q_list_len, keep_prob)
    return quantile_pred

def decoder(encoder_state,decoder_input,time_series_dim,decoder_len):
    time_series_dim
    with tf.variable_scope('decoder', reuse=tf.AUTO_REUSE):
        global_mlp_output, ca = global_mlp(encoder_state, decoder_input, time_series_dim)
        ca_tile = tf.tile(ca, [1, decoder_len, 1])
        local_mlp_input = tf.concat([global_mlp_output, ca_tile, decoder_input],-1)
        local_mlp_output = local_mlp(local_mlp_input)
    return local_mlp_output

def fork_transformer(h_tensor,y_tensor,x_tensor,encoder_len,decoder_len,is_training):
    decoder_y_ta = tf.TensorArray(dtype=tf.float32, size=encoder_len)
    decoder_output_ta = tf.TensorArray(dtype=tf.float32, size=encoder_len)

    def cond_fn(time, decoder_y_ta, decoder_output_ta):
        return time <= encoder_len

    def loop_fn(time, decoder_y_ta, decoder_output_ta):
        decoder_y = tf.slice(y_tensor, [0, time], [-1, decoder_len])
        decoder_x = tf.slice(x_tensor, [0, time, 0], [-1, decoder_len, -1])
        decoder_h = tf.slice(h_tensor, [0, time - 1, 0], [-1, 1, -1]) 
        decoder_output = decoder(decoder_h, decoder_x)
        decoder_y_ta = decoder_y_ta.write(time - 1, decoder_y)
        decoder_output_ta = decoder_output_ta.write(time - 1, decoder_output)
        return time + 1, decoder_y_ta, decoder_output_ta
    if is_training:
        # 训练的时候使用fork机制,每个时间点进行decoder
        loop_init = [tf.constant(1, dtype=tf.int32),
                     decoder_y_ta,
                     decoder_output_ta]
    else:
        # 预测的时候不使用fork机制，直接将encoder最后state及x_f输入decoder，计算得到decoder_output
        loop_init = [tf.constant(encoder_len, dtype=tf.int32),
                     decoder_y_ta,
                     decoder_output_ta]
    _, decoder_y_loop, decoder_output_loop = tf.while_loop(cond_fn, loop_fn, loop_init)
    decoder_y_tensor = tf.transpose(decoder_y_loop.stack(), perm=[1, 0, 2])
    decoder_output_tensor = tf.transpose(decoder_output_loop.stack(), perm=[1, 0, 2, 3])
    return decoder_y_tensor, decoder_output_tensor

loss function

def quantile_loss(y_true,y_hat,Q_list,Q_list_len,encoder_len,decoder_len):
    '''
    Q_list：分位数列表
    '''
    with tf.name_scope("loss"):
        y_true_expand = tf.tile(tf.expand_dims(y_true, -1), [1, 1, 1, Q_list_len]) 
        batch_size = tf.shape(y_true_expand)[0]
        q_tensor = tf.tile(tf.expand_dims(Q_list, 0), [decoder_len, 1])
        q_tensor = tf.tile(tf.expand_dims(q_tensor, 0), [encoder_len, 1, 1])
        q_tensor = tf.tile(tf.expand_dims(q_tensor, 0), [batch_size, 1, 1, 1])
        loss_op = tf.math.reduce_mean(tf.multiply(q_tensor, tf.math.maximum(0., y_true_expand - y_hat)) + tf.multiply((1 - q_tensor), tf.math.maximum(0., y_hat - y_true_expand)))
    return loss_op