Note：(Grouped) Network Autoregression Model (1)

发表于 2025-07-22 更新于 2025-07-25 分类于读笔 Waline：

传统的空间计量模型依赖于空间之间的相关性（例如计算Moran's I），其核心依旧在于外生冲击X对Y的影响（增加了空间扩散效应），也就是说，在传统的面板数据的，增加了或者这样的传导路径，然而，对于纯粹的空间溢出效应（或者说网络效应，不那么像地理），更多采用VAR模型来刻画（包括TVP-VAR等方法），但是一方面存在纬度诅咒，另一方面，我们可以获得一些直接可观测的网络结构，这篇blog介绍一些新的方法，核心参考文献：

Zhu, X., Pan, R., Li, G., Liu, Y., & Wang, H. (2017). Network Vector Autoregression. The Annals of Statistics, 45(3), 1096-1123.
Zhu, X., & Pan, R. (2020). Grouped network vector autoregression. Statistica Sinica, 30(3), 1437-1462.
Zhu, X., Xu, G., & Fan, J. (2025). Simultaneous estimation and group identification for network vector autoregressive model with heterogeneous nodes. Journal of Econometrics, 249, 105564.

1 模型设定

首先介绍 NAR 模型

其中： - ：节点在时间的观测值 - ：节点的外生协变量向量 - ：网络邻接矩阵元素，表示节点和之间的连接权重 - ：节点的度数（连接数） - ：独立同分布的随机误差项

因此，NAR 模型可以看作是对传统面板数据模型的扩展，引入了网络结构信息，从而更好地捕捉节点之间的相互影响。

此外，在传统的动态面板模型（个体很大，时间有限）中，引入节点异质性会导致 OLS 估计量的不一致（Nickell Bias）。对于 NAR 模型，Zhu et al.（2017）证明了在节点数量和观测时间同时趋于无穷时，即使节点性质控制变量簇和误差项弱相关，也可以获得无偏一致估计量。

Nickell（1981）给出了固定效应估计（Fixed Effects, FE）在包含滞后因变量的动态模型中存在的偏误（Nickell Bias）大小，即：当时，Nickell Bias 几乎可以忽略。

接下来介绍 GNAR 模型

假设第个节点携带一个潜在变量，其中当属于第个分组时，否则。因此，GNAR 模型可以写为：

其中为独立噪声项，服从标准正态分布。

进一步拓展，如果把节点的分组潜变量作为待估参数，同时考虑节点间的网络效应不仅由接受节点，也由发送节点的状态影响，则 GNAR 模型可以表示为：

考虑一个包含个节点的网络，节点编号为，其关系由邻接矩阵记录，其中表示第个节点关注第个节点，否则为。约定。对于第个节点，观测到一个连续变量的时间序列，以及一组节点特定的协变量。特别地，的第一个分量始终为，对应截距项。

为刻画网络异质性，假设网络节点可分为个组，每组内回归效应同质，记表示第个节点的分组。GNAR 模型可表示为：

其中，表示节点的出度，为均值为、方差为的独立同分布随机噪声。所有模型参数及节点分组均需估计。

本篇笔记主要介绍公式（3）的模型估计方法

2 模型估计

2.1 损失函数设定

首先设定损失函数

def loss_function(self, beta=None, group=None):
    """ Calculate the loss function based on the parameters beta
    beta: (G+p+1)*G matrix of coefficients, if None, use self.beta
    group: N vector of group labels, if None, use self.group
    return:
    loss: float, the loss value
    """
    if beta is None:
        beta = self.beta
    if group is None:
        group = self.group
    #按照N个节点对应的组，构建N*N的系数矩阵，其中i，j节点的系数为beta[group[j], group[i]]
    self.coef = np.zeros((self.N, self.N))
    for i in range(self.N):
        for j in range(self.N):
            self.coef[i, j] = beta[group[j], group[i]]
    # 按照N个节点对应的组，把G个动量系数v构建长度为N的动量系数
    self.v_coef = np.zeros(self.N)
    for i in range(self.N):
        self.v_coef[i] = beta[self.G, group[i]]
    # 按照N个节点对应的组，把G*p个CV系数gamma构建N*p的CV系数
    self.gamma_coef = np.zeros((self.N, self.p))
    for i in range(self.N):
        self.gamma_coef[i, :] = beta[self.G+1:self.G+self.p+1, group[i]]
    # 将N*N*T的network矩阵与系数矩阵相乘，得到N*N*T的网络系数矩阵
    self.network_coef = np.zeros((self.N, self.N, self.T))
    for t in range(self.T):
        self.network_coef[:, :, t] = self.network[:, :, t] * self.coef
    #计算每一天的 Yt = v@Yt-1 + network_coef@Yt-1 + gamma@CV + eps，注意是矩阵乘法
    Y_est = np.zeros((self.N, self.T-1))
    for t in range(1, self.T):
        # 计算上一天的Yt-1
        Y_prev = self.Y[:, t-1]
        # 计算动量部分
        momentum = self.v_coef * Y_prev
        # 计算网络部分
        network_part = self.network_coef[:, :, t-1] @ Y_prev
        # 计算CV部分（先乘再求和）
        CV_part = np.sum(self.gamma_coef * self.CV[:, :, t-1], axis=1)
        # 计算当天的Yt
        Y_est[:, t-1] = momentum + network_part + CV_part
    loss = np.sum((Y_est - self.Y[:, 1:]) ** 2) / (self.N * (self.T-1))
    return loss

其中， - ，其中 -

对于给定的分组，的最优解为其中为所有按堆叠而成的设计矩阵，为对应的响应变量。

def cul_para_OLS(self):
    """ Calculate the parameters using self.group
    return:
    beta_all: (G+p+1)*G matrix of coefficients
    """
    X = np.zeros((self.N, self.G + self.p + 1, self.T - 1))
    Y = np.zeros((self.N, self.T - 1))
    beta_all = np.zeros((self.G + self.p + 1, self.G))
    for i in range(self.N):
        # 取Network 第i列
        X_i = np.zeros((self.G+self.p+1, self.T-1))
        Y_i = self.Y[i, 1:]
        for g in range(self.G):
            WY = self.network[i, :, :-1] * self.Y[:, :-1] * np.repeat((self.group == g).reshape(-1, 1), self.T-1, axis=1)  # N*T-1
            X_i[g, :] = np.sum(WY, axis=0)
        X_i[self.G, :] = self.Y[i, :-1]
        X_i[self.G+1:self.G+self.p+1, :] = self.CV[i, :, :-1]
        X[i, :, :] = X_i
        Y[i, :] = Y_i
    for g in range(self.G):
        if np.sum(self.group == g) == 0:
            # 初始化系数的时候不会出现某个组为0的情况
            beta_all[:, g] = self.beta[:, g]
            continue
        Xg = np.concatenate(X[self.group == g],axis=1)  # N_obs*(G+p+1)*(T-1)
        yg = np.concatenate(Y[self.group == g])  # N_obs*(T-1)
        # 计算OLS系数
        beta = np.linalg.inv(Xg @ Xg.T) @ Xg @ yg
        beta_all[:, g] = beta
    return beta_all

如何理解上述操作？

当分组确定的时候，注意到对任何给定的，方程3中的所有待估计都只和这一组内的所有节点有关，因此可以单独提取出这一组节点，然后进行最小二乘估计。

2.2 分组初始化

为获得初始分组，可以使用 k-means 算法对节点特征进行聚类。具体步骤如下：

定义中心化变量：并对应地令 , 。

基于GNAR模型，有：其中，。

初始参数估计可通过对每个节点分别做最小二乘回归获得：令，其中。则有其中为回归参数。

def ridgeOLS(self, X, y):
    """ Ridge Ordinary Least Squares estimation """
    # 增加ridge tuning parameter
    XTX = X @ X.T
    # 计算lambda值并添加到对角线
    lambda_values = []
    for j in range(X.shape[0]):
        x_it_norm_squared = np.sum(X[j, :] ** 2)
        lambda_j = 0.01 * x_it_norm_squared / (X.shape[1] + 1) + 1e-6
        lambda_values.append(lambda_j)
    # 在XTX对角线上加入lambda值
    XTX += np.diag(lambda_values)
    beta = np.linalg.inv(XTX) @ X @ y
    return beta

def initialize_base_OLS(self):
    """ Initialize for kmean
    return:
    beta: N*N matrix of coefficients
    v: N vector of momentum coefficients
    fi: N vector of intercepts
    """
    # 删去第0天并对每个个体去均值
    Y_tire = self.Y[:, 1:] - np.mean(self.Y[:, 1:], axis=1, keepdims=True)
    # 删去最后一天并对每个个体去均值
    Y_tire_lag = self.Y[:, :-1] - np.mean(self.Y[:, :-1], axis=1, keepdims=True)
    # 初始化beta N+1* N
    beta = np.zeros((self.N + 1, self.N))
    fi = np.zeros(self.N)
    for i in range(self.N):
        Y_tire_lag_i = Y_tire_lag[i, :].reshape(1, -1) # 1*T-1
        Y_tire_i = Y_tire[i, :] # 1*T-1
        # 取Network
        WY_tire_lag_i = self.network[i, :, :-1] * Y_tire_lag   # N*T-1
        # 拼接
        X = np.concatenate([WY_tire_lag_i,Y_tire_lag_i], axis=0) # N+1*T-1
        # 计算beta
        beta[:, i] = self.ridgeOLS(X, Y_tire_i)  # T-1*1
        # 计算所有时间的网络均值
        net_mean = np.mean(self.network[i, :, :-1], axis=1)  # N*1
        fi[i] = np.mean(self.Y[i,1:]) - np.sum(beta[:-1, i] * net_mean * np.mean(self.Y[:, :-1], axis=1))- beta[-1, i] * np.mean(self.Y[i, :-1])
    return beta[:-1,:],beta[-1,:], fi

如何理解上述操作？

通过去均值消除不重要的参数，然后假设所有节点存在异质性，进行完全的最小二乘估计，有个网络参数。

在实际操作中，初始分组可通过以下三种 k-means 聚类方案获得：

基于各节点动量参数估计进行 k-means 聚类。
基于各节点固定效应估计进行 k-means 聚类。
基于各节点网络效应参数（）进行如下两步聚类：
1. 首先对所有进行 k-means 聚类，聚类数为，聚类标签记为。
2. 对每个节点，定义向量，其中，。
3. 令，对所有节点的进行 k-means 聚类，聚类数为，得到初始分组向量。

def k_means_clustering(self, method='networkeffect'):
    """ Perform k-means clustering based on the specified method
    method = networkeffect: 'momentum', 'fixedeffect', or 'networkeffect'
    ! networkeffect is recommended for better results
    """
    beta, v , fi = self.initialize_base_OLS()
    if method == 'momentum':
        # 使用动量系数进行k-means聚类
        Kres = KMeans(n_clusters=self.G, random_state=self.seed).fit(v.reshape(-1, 1))
        self.group = Kres.labels_
    elif method == 'fixedeffect':
        # 使用固定效应进行k-means聚类
        Kres = KMeans(n_clusters=self.G, random_state=self.seed).fit(fi.reshape(-1, 1))
        self.group = Kres.labels_
    elif method == 'networkeffect':
        # 使用网络效应进行k-means聚类, Two steps
        # 先对beta进行kmeans G^2个组
        Kres1 = KMeans(n_clusters=self.G**2, random_state=self.seed).fit(beta.reshape(-1, 1))
        c = Kres1.labels_.reshape(self.N, self.N)
        beta_net = np.zeros((self.N,1+self.G**2))
        for l in range(self.G**2):
            beta_ = beta.copy()
            # 只保留c=1的beta，其他设为0
            beta_[c != l] = 0
            # 对每个组的beta求均值
            beta_net[:,1+l] = np.mean(beta_, axis=0) / np.sum(c == l, axis=0)
        beta_net[:, 0] = v
        # fillna
        beta_net = np.nan_to_num(beta_net, nan=0.0)
        Kres2 = KMeans(n_clusters=self.G, random_state=self.seed).fit(beta_net)
        self.group = Kres2.labels_
    return self.group

2.3参数迭代估计方法

提出如下迭代算法以联合优化：

初始化分组：使用 k-means 算法获得初始分组估计，并据此用 (2.3) 得到初始参数和
更新分组成员：在第次迭代中，依次更新每个节点的分组估计。具体地，对每个节点，其分组成员通过如下方式更新：其中表示将第个节点分配到分组后，其余节点分组保持不变。对所有重复上述步骤，直到分组不再发生变化。具体地，即将第个节点分配到分组，其余节点分组保持当前迭代状态。

如何理解上述操作？ 对任何一组给定的系数（理论上分为个组就有组系数），依次按序调整每一个节点，保证对每一个节点，其最优分组是当前分组的最优解

更新参数估计：固定分组，利用 (2.3) 对每个分组分别进行最小二乘估计，得到更新后的参数和。
重复迭代：重复步骤 (b)-(c)，直到损失函数收敛为止。

def update_para_group(self):
    """ Update the group labels based on the loss function
    return:
    update_count: int, the number of nodes whose group labels are updated
    """
    update_count = 0
    if self.group is None:
        raise ValueError("Group is not initialized, Run k_means_clustering() first")
    # update beta
    self.beta = self.cul_para_OLS()
    for i in tqdm(range(self.N)):
        i_group_loss = []
        for g in range(self.G):
            group = self.group.copy()
            group[i] = g
            loss = self.loss_function(group=group)
            i_group_loss.append(loss)
        # 找到最小损失对应的组
        if self.group[i] != np.argmin(i_group_loss):
            update_count += 1
            self.group[i] = np.argmin(i_group_loss)
    return update_count

def update_all_para(self, max_iter=100):
    """ Update all parameters until convergence or max_iter reached
    return:
    group: N vector of group labels
    beta: (G+p+1)*G matrix of coefficients
    """
    for epoch in range(max_iter):
        update_count = self.update_para_group()
        print(f"Epoch {epoch}: {update_count} nodes updated")
        if update_count == 0:
            print("Convergence reached")
            break
    return self.group, self.beta

2.4 收敛性分析

原文对于迭代算法和参数的收敛性进行了详细分析，这里仅摘取（我认为）比较重要的几个条件：

Condition 2 真实参数

假设 (a) ；(b) 存在常数，使得

Condition 3 网络结构

对任意，定义比例，。

假设存在极限和，使得，，且存在常数，使得

这两个条件的核心点在于

约束了自回归过程的“温和和稳定性”，这一点在实际应用中不难满足，因为观测到的值通常不会是发散的
约束了系数的异质性，保证了不同组之间的系数差异是显著的，这一点比较难以满足，可以参考下一章节的模拟分析

3 模拟分析

这一部分主要是我自己的一些测试和不成熟想法，可能存在错误

3.1 测试

首先，参考Zhu et al.（2025）中的模拟分析，设定：

sim = simulator.GNAR_simulator(N=100, T=300, G=3)
beta = np.array([[0.15,0.2,-0.1],[0.1,0.3,-0.2],[0.15,0.1,0.3]]).T
v = np.array([0.2,0.4,0.6])
beta, v, gamma, group = sim.generate_para(beta=beta, v=v)
#beta, v, gamma, group = sim.generate_para()
Y, CV, network = sim.generate_data()
est = estimator.GNAR_estimator(Y, CV, network,G=3)
group_est_0 = est.k_means_clustering(method='networkeffect')
group_est, beta_est = est.update_all_para()
display(beta_est)
'''
100%|██████████| 100/100 [00:03<00:00, 26.65it/s]
Epoch 0: 1 nodes updated
100%|██████████| 100/100 [00:03<00:00, 26.93it/s]
Epoch 1: 0 nodes updated
Convergence reached
[[ 0.29871302  0.1506184   0.10031757]
 [-0.096042    0.14965717  0.19838028]
 [-0.20008931  0.10065114  0.30056093]
 [ 0.60005426  0.19921955  0.39979723]
 [ 0.61187271 -0.12828938 -0.56726863]
 [-0.15554652  0.06205307 -0.10291666]
 [-0.17400902  0.97779666  0.04877714]
 [ 0.31758294 -0.34173447 -0.8287787 ]]

拟合结果较好

然而，当选择随机系数的时候，在的条件下，模型的表现依旧良好，仍然能收敛到真实数值，然而，G=3的时候，模型的表现急剧下降，错误的分组往往导致模型的动量效应被错误的估计到接近1，而网络效应接近于0。

#beta, v, gamma, group = sim.generate_para()
Y, CV, network = sim.generate_data()
est = estimator.GNAR_estimator(Y, CV, network,G=3)
group_est_0 = est.k_means_clustering(method='networkeffect')
group_est, beta_est = est.update_all_para()
display(beta_est)
'''
100%|██████████| 100/100 [00:03<00:00, 26.39it/s]
Epoch 0: 45 nodes updated
100%|██████████| 100/100 [00:03<00:00, 27.15it/s]
Epoch 1: 26 nodes updated
100%|██████████| 100/100 [00:03<00:00, 26.38it/s]
Epoch 2: 9 nodes updated
100%|██████████| 100/100 [00:03<00:00, 27.12it/s]
Epoch 3: 0 nodes updated
Convergence reached
[[0.14955406 0.09261598 0.13903634]
 [0.23811479 0.03663893 0.24985683]
 [0.19167095 0.17679469 0.20424485]]
 [0.57343237 0.61152089 0.6912202 ]
[[ 2.05792403e-01  5.43748256e-03  2.27624382e-01]
 [ 1.79468545e-01 -1.70604967e-02  2.62645702e-01]
 [ 1.75468029e-01  1.09056927e-02 -1.67185548e+00]
 [ 6.91405397e-01  9.99097215e-01  5.02156734e-01]

注意到，上述迭代方程实际上只更新了一次节点，然而，在初始系数比较接近的时候，K-means的分组效果很差，会大量更新节点，在的时候，更新节点（大概率）能更新到正确的位置，但在的时候，更新节点的错误反而会导致网络项系数的错误估计，最终造成结果不收敛

3.2 Time Variant GNAR

这一部分完全基于我自己的一些数值模拟，尝试性的调整，无严格推证

我们进一步假设网络结构是时间变化的，控制变量也是时间变化的，模型变为：

在新的设定下，首先，原先的K-means分组方法不再适用，因为网络结构和控制变量都是时间变化的，因此无法通过简单的去均值来消除固定效应，需要引入改进的K-means初始化

首先，在初始参数估计方面，同样可通过对每个节点分别做最小二乘回归获得，但此时不再使用去均值的方式，而是直接使用原始数据进行回归，同时也需要纳入控制变量的完全系数！！令，其中。则有其中为回归参数。

如何理解上述操作？

新设定下，原先的所谓固定效应被细化成了控制变量的回归系数，因此，kmeans的聚类也从的空间转变为控制变量的回归系数空间。

新设定的好处如下： - 由于控制变量是时间变化的，在回归中，其实际有意义的样本数从增加到了，这使得模型能够更好地捕捉到时间变化带来的影响 - 由于网络效应同样随时间变化，模型能够更好地捕捉到网络结构的变化对节点行为的影响，对的估计也更为准确

使用随机系数进行数值模拟，结果如下：

sim = simulator.GNAR_simulator(N=100, T=300, G=3, network_time_varying=True, CV_time_varying=True)
beta = np.array([[0.15,0.2,-0.1],[0.1,0.3,-0.2],[0.15,0.1,0.3]]).T
v = np.array([0.2,0.4,0.6])
#beta, v, gamma, group = sim.generate_para(beta=beta, v=v)
beta, v, gamma, group = sim.generate_para()
Y, CV, network = sim.generate_data()
est = estimator.GNAR_estimator(Y, CV, network,G=3)
group_est_0 = est.k_means_clustering(method='networkeffect',time_varying=True)
group_est, beta_est = est.update_all_para()
print(beta)
print(v)
print(beta_est)
'''
100%|██████████| 100/100 [00:03<00:00, 26.68it/s]
Epoch 0: 0 nodes updated
Convergence reached
[[0.23457801 0.22121888 0.04764032]
 [0.33318105 0.02456799 0.37209593]
 [0.04396837 0.04075084 0.276216  ]]
[0.53957165 0.7045408  0.42471186]
[[ 0.27675469  0.37250482  0.04731419]
 [ 0.04031485  0.02641971  0.22076745]
 [ 0.04350052  0.33138305  0.23508159]
 [ 0.42439437  0.70449122  0.53957243]
 [-0.37458463  0.7889177  -0.28970793]
 [-0.49728413 -0.64000801  0.31258883]
 [ 0.24437562 -0.10160471  0.3148274 ]
 [ 0.02428095 -0.07410229 -0.07048587]]

此时，无论初始通过网络效应还是控制变量系数（固定效应）进行k-means聚类，模型都能收敛到真实参数，且分组效果非常良好。

上述测试所用的完整代码附于：https://github.com/Ardentem/GNAR_python，已经封装成了一个完整的包，可以直接进行计算和输出回归结果，具体使用方法请参考该仓库中的文档。