随着社会的现代化发展, 人类生活质量不断提高. 因此, 选择公共交通工具出行在日常通勤中的占比日益下降, 城市道路交通网络的运输压力亦会增加. 在居民选择出行路线以及有关部门对道路交通进行管理时, 如何提高城市道路交通网络的运输效率成为社会关注的重要问题之一.

由于城市道路交通符合网络数据集的典型结构特征, 本文将一般网络数据的研究方法, 运用于城市道路交通网络的分析中, 即将道路视为顶点, 路口的每个可行驶方向视为连接顶点的有向边, 构建道路交通网络的有向图结构.

为反映各条道路上的交通运输效率, 本文采用车辆平均行驶速度作为响应变量. 各道路在某一时段中的车辆平均行驶速度随时间的变化而变化, 为典型的时间序列的数据结构, 故可构建多元时间序列模型[1, 2, 3, 4, 5]以研究车辆平均行驶速度的变化规律. 多元时间序列模型的经典研究方法是分别对各分量上的时间序列进行建模, 分析与预测[6, 7, 8]. 然而, 在交通网络数据集中, 各道路在某一时段中的车辆平均行驶速度亦与前一时段中相邻道路的车辆平均行驶速度有关, 即在空间上亦具有变化规律. 故此经典研究方法不足以反映网络数据集中的响应变量在空间上的相关性. Bedrick and Tsai [9]与Lütkepohl [10]提出了向量自回归(VAR)模型以反映不同时间序列中响应变量间的相关性. 进一步, 结合网络数据的图结构, Zhu et al. [11]提出网络向量自回归(NAR)模型, 且Zhu and Pan [12]提出分组网络向量自回归(GNAR)模型, 来解释多元时间序列在空间上的联系. 然而, 在本文所处理的数据集中, 道路的数量, 亦即时间序列模型的维数远大于样本量, 因此若直接采用VAR, NAR抑或GNAR模型会面临因参数维数过大而产生的维数灾难带来的巨大挑战.

因此, 本文改进了NAR模型并提出有向图上的高维时间序列的线性自回归模型来研究交通运输效率在时间与空间上的变化规律. 在参数估计的方法上, 基于最小二乘法, 通过引入$ \ell_1 $–正则项以解决降低参数维数的问题. 进一步, 数值模拟的结果表明所提出的估计方法在有限样本下具有良好的表现. 本文通过处理成都市道路网络的城市链路速度数据集, 分别分析了位于成都市中心的交通运输压力较大的春熙路路段, 以及位于成都市市郊的交通运输压力较小的犀浦立交桥的各时段交通情况, 揭示了不同区域的交通运输效率在时间与空间上的变化规律, 并提供了疏导交通堵塞以及缓解交通压力的统计建议.

本文的章节安排如下. 在第二节中, 介绍了VAR模型及其平稳性条件并提出有向图上的高维时间序列模型. 在第三节中, 建立了有向图上的高维时间序列模型的参数估计方法, 并对调节参数的选择进行讨论. 在第四节中, 对本文提出的有向图上的高维时间序列模型的参数估计方法进行模拟实验, 并展示其结果. 在第五节中, 运用所提出的方法分析成都市道路网络的城市链路速度数据集. 第六节对全文进行总结与讨论.

本节首先介绍NAR模型及其参数的解释, 再由此引申出本文提出的有向图上的高维时间序列模型, 最后讨论其平稳性条件.

记网络的有向图结构中顶点的数量为$ d $, 称之为网络大小. 记$ \mathbb{Y}=\{Y_{i, t}:i=1, \ldots, d;t=1, \ldots, T\} $为有向图上的时间序列, 其中$ Y_{i, t}\in \mathbb{R}^{1} $为$ t $时刻于顶点$ i $处的连续型响应变量. 记$ {{\bf{Y}}}_{\cdot, t}=(Y_{1, t}, \ldots, Y_{d, t})^\top $, $ {{\bf{Y}}}=({{\bf{Y}}}_{\cdot, 1}, \ldots, {{\bf{Y}}}_{\cdot, T}) $. $ {\bf{A}}=\left ( a_{ij} \right )_{d\times d} $为有向图的邻接矩阵. Zhu et al. [11]在提出的NAR模型假设中, 考虑了网络的图结构框架, 即响应变量在时间与空间上具有相关性, 亦即在时间上$ Y_{i, t} $会与该顶点在前一时段的响应变量值$ Y_{i, t-1 } $具有相关性, 在空间上$ Y_{i, t} $会与相邻顶点在前一时段的响应变量值$ Y_{j, t-1}, j\in\left\lbrace j:a_{ij}=1\right\rbrace $具有相关性. 因此, 响应变量之间有如下的自回归关系:

$ \begin{equation} Y_{i, t}=\beta _{0}+\beta _{1}Y_{i, t-1 }+\beta_2n_i^{-1}\sum\limits_{j=1}^{d}a_{ij}Y_{j, t-1 }+\varepsilon _{i, t}, \end{equation} $

(2.1)

其中, $ \beta _{0}, \beta _{1}, \beta_2 $为待估参数, $ n_i=\sum_{j=1}^d a_{ij} $为顶点$ i $的出度, $ \{\varepsilon _{i, t}:i=1, \ldots, d;t=1, \ldots, T\} $为一系列独立的随机噪声项. 然而, 在Zhu et al. [11]的NAR模型假设中, 默认某顶点处某时段的响应变量值在空间上与相邻顶点处在前一时段的相关性是相同的. 但是在本文考虑的成都市道路网络的城市链路速度数据集中, 此假设不足以反映各条道路因地理位置的不同而导致交通网络提供的运输效率上的差异. 因此, 考虑空间上的差异性, 对(2.1)式进行改进, 本文提出有向图上的时间序列模型：

$ \begin{equation} Y_{i, t}=\beta _{0}+\beta _{1}Y_{i, t-1 }+\sum\limits_{j=1}^{d}a_{ij}\eta _{j}Y_{j, t-1 }+\varepsilon _{i, t}, \end{equation} $

(2.2)

其中, 与(2.1)不同的是, (2.2)中相邻顶点处响应变量对应的回归系数$ \boldsymbol{\eta}=(\eta _1, \ldots, \eta_d)^\top $为$ d $维待估参数向量.

在模型的有向图中, 若顶点个数$ d $固定, 则$ \mathbb{Y} $是多元时间序列. 若其满足平稳性, 则该序列的统计性质不随时间的推移发生改变. 为方便表述, 记$ \Lambda_{\max}({\bf{M}}) $为矩阵$ {\bf{M}} $的特征值模的最大值, $ {\bf{M}}_{l, \cdot} $与$ {\bf{M}}_{\cdot, m} $分别为矩阵$ {\bf{M}} $的第$ l $行与第$ m $列, $ {\bf{M}}_{-l, \cdot} $与$ {\bf{M}}_{\cdot, -m} $分别为矩阵$ {\bf{M}} $删去第$ l $行与删去第$ m $列得到的子矩阵. 记$ {\rm diag}({\bf{v}}) $为向量$ {\bf{v}} $中的所有分量形成的对角阵, $ {\rm vec}(\cdot) $为矩阵的按列拉直算子, $ {\bf{I}}_d $为$ d $阶单位矩阵, $ {\bf{J}}_{l\times m} $为$ l\times m $阶分量全为$ 1 $的矩阵. 记$ {\bf{H}}={\rm diag}(\boldsymbol{\eta}) $, $ {\bf{W}}=\beta_1{\bf{I}}_d+{\bf{A}}{\bf{H}} $.

有向图上的时间序列模型(2.2)存在唯一严格平稳解的充分必要条件为

$ \begin{equation} \Lambda_{\max}( {\bf{W}} )=\Lambda_{\max}( \beta _{1}{\bf{I}}_d+{\bf{A}}{\bf{H}} )< 1, \end{equation} $

(2.3)

其解的形式为:

$ \begin{equation} {\bf{Y}}_{\cdot, t}=\beta_0\left ( {\bf{I}}_d-{\bf{W}} \right )^{-1}{\bf{J}}_{d\times 1}+\sum\limits_{j=0}^{\infty }{\bf{W}}^{j}\boldsymbol{\varepsilon}_{t-j}, \end{equation} $

(2.4)

其中, $ \boldsymbol{\varepsilon}_t=(\varepsilon_{1, t}, \ldots, \varepsilon_{d, t})^\top $. 特别地, 若$ \varepsilon _{i, t} $独立同分布于标准正态分布, 即$ \varepsilon _{i, t}\sim N\left ( 0, \sigma ^{2} \right ) $, 严格平稳解服从正态分布, 均值$ \boldsymbol{\mu} $和协方差矩阵$ \boldsymbol{\Sigma} $分别为:

$ \begin{equation*} \boldsymbol{\mu} =\beta_0\left ( {\bf{I}}_d-{\bf{W}} \right )^{-1}{\bf{J}}_{d\times 1}=\beta_0\left ( {\bf{I}}_d-\beta_{1}{\bf{I}}_d-{\bf{A}}\bf{H} \right )^{-1}{\bf{J}}_{d\times 1}, \end{equation*} $

$ \begin{equation*} {\rm vec}\left ( \boldsymbol{\Sigma} \right )=\sigma ^{2}\left ( {\bf{I}}_{d\times d}-{\bf{W}}\otimes {\bf{W}} \right )^{-1}{\rm vec}\left ( {\bf{I}}_d \right ). \end{equation*} $

进一步, 若

$ \begin{equation} \left \vert \beta _{1} \right \vert+\sum\limits_{j=1}^{d}\left \vert \eta _{j} \right \vert<1, \end{equation} $

(2.5)

则有向图上的时间序列模型(2.2)存在唯一严格平稳解, 称(2.5)为(2.2)存在唯一严格平稳解的充分条件.

此外, 本文将此性质推广至有向图上的高维时间序列的平稳性条件, 即有向图中, 顶点个数$ d $远大于样本量$ T $, 且二者均趋于无穷的情形. 首先定义有向图上的高维时间序列的平稳性：令$ \left \{ {\bf{Y}}_{\cdot, t}\in \mathbb{R}^{d} \right \} $是$ d\rightarrow \infty $时的$ d $维时间序列.

定义$ \boldsymbol{\Omega}=\left \{ \boldsymbol{\Omega} \in \mathbb{R}^{\infty }:\sum_{i=1}^\infty \left \vert \omega _{i} \right \vert< \infty \right \} $, 其中$ \boldsymbol{\Omega}=\left ( \omega _{i}\in \mathbb{R}^{1}:1\leq i\leq \infty \right )^{\top}\in \mathbb{R}^{\infty } $. 对于任意$ \boldsymbol{\Omega} \in \boldsymbol{\Omega} $, 令$ \boldsymbol{\Omega}_{d}=\left ( \omega _{1}, \cdots , \omega _{d} \right )^{\top}\in \mathbb{R}^{d} $是截断的$ d $维向量. 若$ {\bf{Y}}_{\cdot, t} $满足下列条件: $ \forall \boldsymbol{\Omega} \in \boldsymbol{\Omega} $,

(ⅰ) $ Y_{t}^{\omega }=\lim_{d\rightarrow \infty }\boldsymbol{\Omega}_{d}^{\top}{\bf{Y}}_{\cdot, t} $几乎必然存在,

(ⅱ) $ \left \{ Y_{t}^{\omega } \right \} $是严格平稳的,

(ⅲ) $ {\bf{Y}}_{\cdot, t} $具有有限的一阶矩, 即$ \underset{1\leq i\leq \infty }{\max}E\left \vert Y_{it} \right \vert< \infty $,

则称$ {\bf{Y}}_{\cdot, t} $具有严格平稳性. 在此定义下, 平稳性条件可以由如下命题给出: 若$ \Lambda_{\max}( {\bf{W}} )=\Lambda_{\max}( \beta _{1}{\bf{I}}_d+{\bf{A}}{\bf{H}} )< 1 $, 则(2.4)为一阶矩有限的有向图上的高维时间序列模型(2.2)的唯一严格平稳解.

本节介绍有向图上的高维时间序列模型(2.2)的参数估计方法. 由于在成都市道路网络的城市链路速度数据集中, 顶点个数$ d $远大于观测时间段数$ T $, 即模型参数的维数远大于样本量, 因此选用正则化方法进行参数估计.

在高维时间序列的框架下, 本文通过对最小二乘损失函数加$ \ell_1 $-惩罚项, 得到正则化估计. 具体如下,

$ \begin{equation} (\widehat{\boldsymbol{\beta}}, \widehat{\boldsymbol{\eta}})\in\arg\min\limits_{\boldsymbol{\beta}\in\mathbb{R}^2, \boldsymbol{\eta}\in\mathbb{R}^d} Q(\boldsymbol{\beta}, \boldsymbol{\eta}), \end{equation} $

(3.1)

其中

$ \begin{equation} Q(\boldsymbol{\beta}, \boldsymbol{\eta})= \sum\limits_{i=1}^{d}\sum\limits_{t=2}^{T}\left ( Y_{i, t}-\beta _{0}-\beta _{1}Y_{i, t-1 }-\sum\limits_{j=1}^{d}a_{ij}\eta _{j}Y_{j, t-1 } \right )^{2}+\lambda \left \vert \boldsymbol{\eta} \right \vert_{1}, \end{equation} $

(3.2)

$ \boldsymbol{\beta} =\left ( \beta _{0}, \beta _{1} \right )^\top\in \mathbb{R}^{2} $, 且$ \lambda>0 $为调节参数, $ \vert\cdot\vert_r $表示向量的$ \ell_r $-范数, $ 1\leq r\leq \infty $. 根据优化问题(3.1)的Karush–Kuh–Tucker条件, 可以得到如下等式:

$ \begin{eqnarray} \frac{\partial Q(\widehat{\boldsymbol{\beta}}, \widehat{\boldsymbol{\eta}})}{\partial \boldsymbol{\beta}}&=&2\sum\limits_{t=2}^{T}{\bf{X}}_{\cdot, t-1}^{\top}{\bf{X}}_{\cdot, t-1} \hat{\boldsymbol{\beta} }-2\sum\limits_{t=2}^{T}{\bf{X}}_{\cdot, t-1}^{\top} \left ({\bf{Y}}_{\cdot, t}-{\bf{A}} \widehat{{\bf{H}}}{\bf{Y}}_{\cdot, t-1} \right )={\bf 0}, \end{eqnarray} $

(3.3)

$ \begin{eqnarray} \frac{\partial Q(\widehat{\boldsymbol{\beta}}, \widehat{\boldsymbol{\eta}})}{\partial \boldsymbol{\eta}}&=&-2\sum\limits_{t=2}^{T}{\rm diag}({{\bf{Y}}}_{\cdot, t}){\bf{A}}^\top({{\bf{Y}}}_{\cdot, t-1}-{{\bf{X}}}_{\cdot, t}\widehat{\boldsymbol{\beta}}-{\bf{A}}\widehat{{\bf{H}}}{{\bf{Y}}}_{\cdot, t-1}) +\lambda{\rm sign}(\boldsymbol{\widehat\eta})={\bf 0}, \end{eqnarray} $

(3.4)

其中$ {{\bf{X}}}_{i, t }=\left(1, Y_{i, t} \right)^{T} $, $ {\bf{X}}_{\cdot, t}=\left({{\bf{X}}}_{1, t}, \cdots , {{\bf{X}}}_{d, t} \right)^{\top} $, 且$ \widehat{{\bf{H}}}={\rm diag}(\widehat{\boldsymbol{\eta}}) $. 由于(3.1)中的目标函数是凸函数, 通过(3.3)和(3.4)中的关系, 本文运用交替迭代的方式结合坐标下降算法以求优化问题的最优解. 记$ \Pi_{\theta, m}({\bf{v}})=(v_1, \ldots, v_{m-1}, 0, v_{m+1}, \ldots, v_d)^{T} $为作用在向量$ {\bf{v}}=(v_1, \ldots, v_d)^{T} $上的投影算子. 所提算法的迭代过程如下:

Algorithm 1优化问题(3.1)的迭代流程图

Require: $ k=0 $, $ \widehat{\boldsymbol{\eta}}^{\left ( 0 \right )}=\left ( \widehat{\eta_{1}}^{\left (0 \right )}, \cdots \widehat{\eta_{d}}^{\left (0 \right )} \right )^{\top} $, $ {\rm tol}=10^{-2} $.   $ \left( i\right) $由(3.3)式更新参数$ \widehat{\boldsymbol{\beta} }^{(k)}=(\widehat{\beta}^{(k)}_0, \widehat{\beta}^{(k)}_1)^\top $的值:      $ \widehat{\boldsymbol{\beta} }^{(k)}=\left ( \sum\limits_{t=2}^{T}{\bf{X}}_{\cdot, t-1}^{\top}{\bf{X}}_{\cdot, t-1} \right )^{-1}\sum\limits_{t=2}^{T}{\bf{X}}_{\cdot, t-1}^{\top} \left ({\bf{Y}}_{\cdot, t}-{\bf{A}} {\rm diag}\left(\widehat{\boldsymbol{\eta}}^{\left ( k \right )}\right){\bf{Y}}_{\cdot, t-1} \right ); $ $ \left( ii\right) $计算矩阵 $\widetilde{{\bf{G}}}^{(k)}= {\bf{A}}_{\cdot, j} {\bf{Y}}_{j, -T}\odot \left( {\bf{Y}}_{\cdot, -T}-\widehat{\beta}^{(k)}_0{\bf{J}}_{N\times \left( T-1\right)}-\widehat{\beta}^{(k)}_1{\bf{Y}}_{\cdot, -T} -{\bf{A}}{\rm diag}\left\{\Pi_{0, j}\left(\widehat{\boldsymbol{\eta}}^{\left ( k \right )}\right)\right\}{\bf{Y}}_{\cdot, -T} \right); $ 其中$ \odot $为矩阵的Hadamard乘积; $ \left( iii\right) $计算数值        $ G^{(k)}= {\bf{J}}_{1\times d}\widetilde{{\bf{G}}}^{(k)}{\bf{J}}_{(T-1)\times1}; $ $ \left( iv\right) $ for $ j=1, \ldots , p $ do   由(3.4)式更新参数$ \widehat{\boldsymbol{\eta}}^{\left ( k+1 \right )}=\left ( \widehat{\eta}_{1}^{\left (k+1 \right )}, \cdots, \widehat{\eta}_{d}^{\left (k+1 \right )} \right )^{\top} $的值:    if $ \left \vert G^{(k)} \right \vert<\lambda/2, $ then     $ \widehat{\eta}_{j}^{(k+1)}=0 $;    else      $ \widehat{\eta}_{j}^{(k+1)}=\left\{G^{(k)}-\lambda{\rm sign}\left(G^{(k)}\right)/2\right\} /\{\vert{{\bf{Y}}}_{j, -T}\vert^2_2\vert{\bf{A}}_{\cdot, j}\vert^2_2\} $;   end if end for $ \left( v\right) $ if $ \left\vert\widehat{\boldsymbol{\eta}}^{(k+1)}-\widehat{\boldsymbol{\eta}}^{(k)}\right\vert_\infty<{\rm tol}, $ then   stop; else    $ k\leftarrow k+1 $,   return to $ \left( i\right) $; end if Output $ \widehat{\boldsymbol{\eta}}=\widehat{\boldsymbol{\eta}}^{(k+1)} $且$ \widehat{\boldsymbol{\beta}}=\widehat{\boldsymbol{\beta}}^{(k)} $; Return $ \widehat{\boldsymbol{\eta}} $与$ \widehat{\boldsymbol{\beta}} $.

在高维模型中, 调节参数的大小控制了模型参数的稀疏度, 因而如何选择是一个非常重要的环节. 本节对于有向图上高维时间序列的参数估计中涉及到的调节参数$ \lambda $的不同选择方法进行比较与讨论. 一般而言, $ K $-折交叉验证法和$ hold $-$ out $交叉验证法是两种最常见的选择调节参数的方法. 然而, 对于本文提出的有向图上的高维时间序列模型而言, 之后发生的事件与之前的事件有关, 因此须保证测试集数据在时间顺序上是位于训练数据之后的, 故在此情况下$ K $-折交叉验证不能用于本文的调节参数选择. 至于$ hold $-$ out $交叉验证法, 它是将数据集根据时间顺序分割为训练集和测试集两个子集. 然而, 经典的$ hold $-$ out $交叉验证法对于测试集的选择具有不确定性, 会导致在测试集上预测能力不足. 为了克服这些不足, 本文提出日向前链交叉验证法, 即按时间进行分割, 将该天的数据作为测试集, 并将此前的所有时段的数据分配到训练集中. 其与$ hold $-$ out $交叉验证的区别如下图所示:

若数据集有$ T $天, 则将生成$ T-1 $个不同的训练集和测试集分割, 如图(b)所示. 这样既可以保证测试集的时间顺序在训练子集后面, 又避免了分割测试集时的不确定性. 因此, 在本文的模拟中采用日向前链交叉验证法对调节参数进行选择.

将时间序列数据集按图(b)方法进行$ T-1 $次分割, 结合式子$ \left( 3.2\right) $, 记$ \widehat{\boldsymbol{\beta} }^{\left (m \right )}\left ( \lambda \right ), \widehat{\boldsymbol{\eta}}^{\left (m \right )}\left ( \lambda \right ) $为第$ m $次分割$ \left( m=1, \ldots, T-1\right) $求解的参数估计值. 给出交叉验证得分的定义如下:

$ {\rm CV}\left ( \lambda \right )=\frac{1}{T-1}\sum\limits_{m=1}^{T-1}L^{\left( m\right) }\left (\widehat{\boldsymbol{\beta} }^{\left (m \right )}\left ( \lambda \right ), \widehat{\boldsymbol{\eta}}^{\left (m \right )}\left ( \lambda \right ) \right ), $

其中$ L^{\left ( m \right )}\left (\cdot\right ) $为第$ m $次分割的损失函数

$ L^{\left( m\right) }\left (\widehat{\boldsymbol{\beta} }^{\left (m \right )}\left ( \lambda \right ), \widehat{\boldsymbol{\eta}}^{\left (m \right )}\left ( \lambda \right ) \right )=\sum\limits_{i=1}^{d}\sum\limits_{t=m+1} ( Y_{i, t}-\widehat{\beta }_{0}^{\left (m \right )}\left ( \lambda \right )-\widehat{\beta }_{1}^{\left (m \right )}\left ( \lambda \right )Y_{i, t-1 }-\sum\limits_{j=1}^{d}a_{ij}\widehat{\eta }_{j}^{\left (m \right )}\left ( \lambda \right )Y_{j, t-1 } )^{2}. $

本文通过求取$ {\rm CV}\left ( \lambda \right ) $的最小值点以得到调节参数$ \lambda $的选取值.

本节对有向图上的高维时间序列模型进行模拟实验, 用于评估所提出方法的实用性. 利用2.2节中提到的充分必要条件(2.3)和充分条件(2.5), 可控制生成时间序列数据的平稳性. 分别考虑(2.3)和(2.5)成立时, $ d=50 $, 100和200以及$ T=30 $, 100, 150和200时参数的估计, 模拟次数设定为500次. 结果的展示中省略了$ \eta_{j}=0 $的估计值$ \hat{\eta}_{j} $, 缺失项NA表示该$ \eta_{j} $无法估计. 采用日向前链时间序列交叉验证法选取$ \lambda $为0.01.

表 1, 表 2, 表 3分别是在(2.3)成立时, $ d=50 $, 100和200时参数$ \boldsymbol{\beta} $和$ \boldsymbol{\eta} $的估计值. 在网络大小$ d $不变的条件下, 随着时间序列$ T $的增大, 估计值的标准差越来越小. 在时间序列$ T $相同的条件下, 随着$ d $的增大, 参数$ \boldsymbol{\eta} $的维数也随之增加, 尽管如此, 所提的估计方法依然表现良好.

表 4, 表 5, 表 6分别是在(2.5)成立时, $ d=50 $, 100和200时参数$ \boldsymbol{\beta} $和$ \boldsymbol{\eta} $的估计值, 从中可以得到类似结论, 不做赘述.

本节将本文所提出的方法应用于分析网络的城市链路速度数据集(https://www.nature.com/articles/s41597-019-0060-3). 该数据集包含成都市交通网络的$ d $=5943条道路, 所观测为2015年6月1日08:00–10:00, 12:00–14:00, 17:00–19:00, 21:00–23:00四个不同时间段内每条道路的汽车平均行驶速度, 每个时间段中以两分钟为一个时间间隔, 即$ T=60 $. 将本文提出的有向图上的高维时间序列模型分别应用于成都市交通网络的四个不同时间段, 以研究成都市交通网络在时间以及空间上的变化规律.

我们首先选取成都市的繁华区域春熙路附近进行分析, 如图 2所示. 以道路4为例, 它在四个时间段的交通平均速度的估计函数分别为:

$ \hat{Y}_{4t}=12.588+0.594\hat{Y}_{4\left ( t-1 \right )}-0.004\hat{Y}_{5\left ( t-1 \right )}-0.086\hat{Y}_{7\left ( t-1 \right )}-0.090\hat{Y}_{10\left ( t-1 \right )}, t\in \left ( 08:00, 10:00 \right ), $

$ \hat{Y}_{4t}=14.552+0.595\hat{Y}_{4\left ( t-1 \right )}-0.075\hat{Y}_{5\left ( t-1 \right )}-0.116\hat{Y}_{7\left ( t-1 \right )}+0.058\hat{Y}_{10\left ( t-1 \right )}, t\in \left ( 12:00, 14:00 \right ), $

$ \hat{Y}_{4t}=11.427+0.607\hat{Y}_{4\left ( t-1 \right )}-0.017\hat{Y}_{5\left ( t-1 \right )}-0.069\hat{Y}_{7\left ( t-1 \right )}-0.078\hat{Y}_{10\left ( t-1 \right )}, t\in \left ( 17:00, 19:00 \right ), $

$ \hat{Y}_{4t}=13.641+0.596\hat{Y}_{4\left ( t-1 \right )}-0.024\hat{Y}_{5\left ( t-1 \right )}-0.081\hat{Y}_{7\left ( t-1 \right )}-0.013\hat{Y}_{10\left ( t-1 \right )}, t\in \left ( 21:00, 23:00 \right ). $

可以看出, 道路4在$ t $时刻的平均速度受道路4在$ t-1 $时刻的平均速度, 以及其相邻道路5, 7, 10在$ t-1 $时刻的平均速度的影响. 不同时间段内, 道路4在$ t-1 $时刻的平均速度对道路4在$ t $时刻的平均速度影响都是正相关, 即道路4在$ t-1 $时刻的平均速度越快, 道路4在$ t $时刻的平均速度也越快. 不同时间段内, 道路5, 7在$ t-1 $时刻的平均速度负向影响道路4在$ t $时刻的平均速度, 这可能是因为道路5, 7周围是春熙路附近的商业街, 每个时间段内车流量都相对较大. 道路5, 7在$ t-1 $时刻车辆较少时, 车辆几乎都拥堵在道路4上, 于是道路5, 7在$ t-1 $时刻平均速度越快, 道路4在$ t $时刻的平均速度相对较慢. 但不同时间段内, 道路10在$ t-1 $时刻的平均速度对道路4在$ t $时刻的平均速度的影响不同. 这可能与道路10旁侧商店较少且位于居民楼附近有一定关系, 故中午12:00–14:00车流量较小, 其余时间段附近车流量大. 因此中午12:00–14:00道路10在$ t-1 $时刻速度越快, 道路4在$ t $时刻速度也相对较快. 但在其余时间段, 道路10在$ t-1 $时刻车辆较少时, 车辆均拥堵在道路4上, 于是道路10在$ t-1 $时刻速度越快, 道路4在$ t $时刻的平均速度也变慢.

进一步, 我们选取了成都市的郊区犀浦立交桥附近路段进行分析, 如图 3所示. 以道路16为例, 它在四个时间段的交通平均速度的估计函数分别为:

$ \hat{Y}_{16t}=12.588+0.594\hat{Y}_{16\left ( t-1 \right )}-0.044\hat{Y}_{12\left ( t-1 \right )}+0.099\hat{Y}_{13\left ( t-1 \right )}, t\in \left ( 08:00, 10:00 \right ), $

$ \hat{Y}_{16t}=14.552+0.595\hat{Y}_{16\left ( t-1 \right )}+0.070\hat{Y}_{12\left ( t-1 \right )}-0.002\hat{Y}_{13\left ( t-1 \right )}, t\in \left ( 12:00, 14:00 \right ), $

$ \hat{Y}_{16t}=11.427+0.607\hat{Y}_{16\left ( t-1 \right )}+0.017\hat{Y}_{12\left ( t-1 \right )}-0.009\hat{Y}_{13\left ( t-1 \right )}, t\in \left ( 17:00, 19:00 \right ), $

$ \hat{Y}_{16t}=13.641+0.596\hat{Y}_{16\left ( t-1 \right )}-0.031\hat{Y}_{12\left ( t-1 \right )}+0.246\hat{Y}_{13\left ( t-1 \right )}, t\in \left ( 21:00, 23:00 \right ). $

同样可以看出, 道路16在$ t $时刻的平均速度受道路16在$ t-1 $时刻的平均速度, 以及其相邻道路12, 13在$ t-1 $时刻的平均速度的影响. 不同时间段内, 道路16在$ t-1 $时刻的平均速度对道路16在$ t $时刻的平均速度影响都是正相关, 即道路16在$ t-1 $时刻的平均速度越快, 道路16在$ t $时刻的平均速度也越快. 但不同时间段内, 道路12, 13在$ t-1 $时刻的平均速度对道路16在$ t $时刻的平均速度的影响不同. 在早上08:00–10:00和晚上21:00–23:00, 道路12在$ t-1 $时刻的平均速度负向影响道路16在$ t $时刻的平均速度. 道路12前往成都电子科技大学, 可能与教职工在早上08:00–10:00上班和晚上21:00–23:00下班有一定关系, 故道路12在这两个时间段车流量大, 道路12在$ t-1 $时刻车辆较少时, 车辆均拥堵在道路16上. 因此道路12在$ t-1 $时刻平均速度越快, 道路16在$ t $时刻的平均速度相对较慢. 而道路13前往博雅体育俱乐部, 人们运动和学习的时间往往是相反的, 于是在这两个时间段, 道路13在$ t-1 $时刻的平均速度正向影响道路16在$ t $时刻的平均速度. 其余两个时间范围则恰恰相反.

值得注意的是, 由于正则化方法只给出了参数的估计值, 而对于其假设检验问题未曾涉及. 显然, 该问题具有更大的挑战性. 鉴于此, 这些统计结论均为描述性的, 其显著性如何尚未考究.

本文提出了有向图上的高维时间序列模型, 建立其正则化统计方法. 数值的结果表明所提出的方法具有良好的表现, 在成都市交通网络数据的应用上, 也得到一些可解释的统计结论. 进一步, 若考虑网络的组群效应, 可将group lasso [13]作为正则项引入, 相应的估计方法与算法实现均需进一步研究.