图片由作者和GPT-4提供。

这很明显，我们未来大部分的工作将会被AI自动化。这将实现是因为许多研究人员和专业人士正在努力将他们的工作数字化或在线化。这些贡献不仅帮助我们理解基本概念，还促进了人工智能模型的优化，最终为我们腾出时间去关注其他任务。

然而，有一个概念即使在专家中也经常被误解，那就是时间序列分析中的伪回归。这种问题出现在回归模型显示变量之间存在强烈关系时，即使实际上并不存在这样的关系。这种情况通常出现在时间序列回归方程看起来拟合度很高——表现为R²（多重相关系数）很高——但Durbin-Watson统计量(d)却很低，这表明误差项存在强烈的自相关。

更令人惊讶的是，几乎所有计量经济学教材都会提醒自相关误差的风险，但这个问题在很多已发表的研究中依然存在。Granger和Newbold（1974）发现了几个例子。例如，他们发现了一些R² 为0.997的方程，杜宾-沃森统计量（d）为0.53。最极端的例子是一个R² 为0.999，d为0.093的方程。

在经济学和金融领域，这种情况尤其棘手，这些领域中，许多关键变量表现出自相关性或者相邻值之间的序列相关性，特别是在采样间隔较短，如一周或一个月时，如果处理不当，可能会导致错误的结论。例如，今天的GDP与上一季度的GDP密切相关。我们在文章中详细解释了Granger和Newbold（1974）的研究发现，并通过Python模拟（详见第7节）再现了他们文章中的关键结果。

无论你是哪一类专业人员，如经济学家、数据科学家，还是分析师，在处理时间序列数据时，理解这一问题对于确保你所使用的模型产生有意义的成果至关重要。

为了带大家了解这篇论文，下一节将介绍随机游走和ARIMA(0,1,1)过程。在第三小节，我们将解释Granger和Newbold（1974）是如何描述无意义回归的产生的，并在第四部分通过示例来解释。最后，我们来谈谈如何避免在处理时间序列数据时出现无意义回归。

设 𝐗ₜ 是一个时间序列。我们说 𝐗ₜ 遵循随机游走，如果它的形式为：

𝐗ₜ = 𝐗ₜ₋₁ + 𝜖ₜ。此公式表示当前时间点的状态是由前一时间点的状态加上一个误差项构成的。（1）

其中 𝜖ₜ 是白噪音。它可以表示为多个白噪音的总和，这对于模拟很有帮助。它是一个非平稳时间序列，因为它的时间方差会随着 t 的变化而变化。

2.2 ARIMA(0,1,1) 模型

这个 ARIMA(0,1,1) 过程是：

[𝐗ₜ = 𝐗ₜ₋₁ + 𝜖ₜ − 𝜃 𝜖ₜ₋₁.] （2）

其中，(\epsilon_t) 表示当前时刻的误差项，(\theta) 是一个参数，而 (𝐗ₜ) 和 (𝐗ₜ₋₁) 分别表示当前时刻和前一时刻的状态。

其中 ( \epsilon_t ) 是白噪声，。ARIMA(0,1,1) 模型是非平稳的。它可以写成一个随机游走和白噪声的和：

$ yt = y{t-1} + \theta \epsilon_{t-1} + \epsilon_t $。

这里随机游走项 ( y{t-1} ) 和白噪声项 ( \theta \epsilon{t-1} ) 是独立的。

[ X_t = X_0 + \text{随机游走过程} + \text{加白噪声。} \ (3)]这种形式在模拟中很有用。

那些非平稳的时间序列常被用来作为其他模型预测性能的参考标准。

首先，我们先来回顾一下线性回归模型。线性回归模型可以表示如下：

Y = Xβ + ϵ。（4）

其中 (𝐘) 是因变量的 (T \times 1) 向量，(𝛽) 是系数的 (K \times 1) 向量，(𝐗) 是一个 (T \times K) 矩阵，包含一列全1和 (K-1) 列独立变量，每列有 (T) 个观测值。这些独立变量（即 (K-1) 列）是随机的，但与误差项的 (T \times 1) 向量 (𝜖) 互相独立。通常假定：

𝐄(𝜖) = 0 (5)

需要更多上下文来准确翻译。

E(𝜖𝜖′) = σ²I。公式(6)

其中 (I) 就是单位矩阵（Identity Matrix）。

F检验是用来测试独立变量对因变量解释的贡献的。它的零假设是：

H₀: 假设这些β值等于零（即β₁ = β₂ = ⋯ = βₖ减1 = 0），（方程7）

测试的统计信息如下所示：

(F = \frac{R^2 / (K-1)}{(1-R^2) / (T-K)}). (8)

其中 R² 表示决定系数（R²）。

如果我们想构建检验统计量，假设零假设成立，尝试用形式如（方程4）的方式对经济时间序列进行回归。假设这些序列不是平稳的或高度自相关的。在这种情况下，这种检验过程是无效的，因为在零假设下（方程7），（方程8）中的𝐹不遵循F分布。事实上，在零假设下，（方程4）中的误差或残差如下：

𝜖ₜ = 𝐘ₜ − 𝐗𝛽₀ ; t = 1, 2, …, T. (9)

并将和原始序列 𝐲 有相同的自相关结构。

在某种情况下，比如，就能大致了解分配问题：

Yₜ = β₀ + Xₜβ₁ + εₜ. (10)
其中，𝐘ₜ 表示时间t的因变量, β₀ 是截距项, 𝐗ₜ 是时间t的自变量, β₁ 是𝐗ₜ的系数, 而𝜖ₜ 是误差项。

其中，𝑌𝑡 和 𝑋𝑡 遵循一阶自回归模型：

( yt = \rho y{t-1} + \eta_t ) 和 ( xt = \rho^* x{t-1} + \nu_t )。（11）

其中 𝜂ₜ 和 𝜈ₜ 是白噪声（white noise）。

我们知道在这种情况下，𝑅²是𝑌𝑡和𝑋𝑡之间相关系数的平方。他们使用了Kendall在Knowles（1954）文章中的结果，该结果表达了𝑅的方差的表达式。

𝐕𝐚𝐫(𝐑) = (1/T) (1 + ρρ) /(1 − ρρ*). (12)

由于 𝐑 的范围被限定在 -1 和 1 之间，其方差若超过 1/3，那么 𝐅的分布就不会在 0 有众数。这意味着，𝜌𝜌* > (T−1) / (T+1)。

因此，例如，如果 T = 20 并且 𝜌 = 𝜌（这里𝜌表示某个特定值），当 𝜌 > 0.86 时，将不会得到在 0 处的单峰分布，而当 𝜁 = 0.9 时，𝑉𝑎𝑟(𝑅) = 0.47。因此 𝔼(𝑅²) 将接近 0.47。

已经表明，当 𝜌 接近 1 时，𝑅² 可以非常高得惊人，这表明 𝐘ₜ 和 𝐗ₜ 之间存在很强的联系。然而，在实际情况中，这两个序列实际上是完全不相关的。当 𝜌 接近 1 时，两个序列的行为就像随机游走。此外，两个序列的高度自相关使得回归的残差也表现出很强的自相关性。因此，Durbin-Watson 统计量 𝐝 将非常低。

这就是为什么在这种情况下，高的 (R^2) 值不应被看作两个序列真正有关系的证据。

为了探讨在两个独立的随机游走过程中进行回归时出现伪回归的可能性，在下一节中，我们将进行一系列由 Granger 和 Newbold 提出的模拟实验。

用 Python 模拟的结果。

在本节中，我们将通过模拟展示使用带有独立随机游走偏差的回归模型会使得系数估计产生偏差，并从而使系数的假设检验变得无效。将在第6节展示用于生成模拟结果的Python代码。

一个由格兰杰和纽伯尔德提出的回归方程（参见（1974）[https://jumbong.github.io/personal-website/Others/spurious_reg.html#ref-granger1974spurious]）为：

𝐘ₜ = 𝛽₀ + 𝐗ₜ𝛽₁ + 𝜖ₜ

当 𝐘ₜ 和 𝐗ₜ 作为独立的随机游走被生成时，每个序列长度为 50。对于 100 次模拟，表示检验 𝛽₁ 是否显著的统计量 𝐒 = |𝛽̂₁| / √(𝐒𝐄̂(𝛽̂₁)) 的值将在下面的表格中汇报。

表1：对两个独立随机漫步的回归分析

如果 𝐒 > 2，则在5%的显著性水平上拒绝原假设（即 𝐘ₜ 和 𝐗ₜ 之间无关系）。此表显示，原假设（𝛽 = 0）大约有四分之一（71次）的情况下被错误拒绝了。这种情况相当尴尬，因为这两个变量是独立的随机游走，实际上没有实际的关系。让我们来看看为什么会这样。

如果 𝛽̂₁ / 𝐒𝐄̂ 服从 𝐍(0,1)，则 𝐒 的期望绝对值应为 √2 / π ≈ 0.8（√2/π 是标准正态分布绝对值的均值，即约为 0.8）。然而，模拟结果的平均值为 4.59，这意味着估计的 𝐒 被低估了大约 5.74 倍（4.59 / 0.8 ≈ 5.74）。

4.59 / 0.8 = 5.7 等于 5.7，结果是 5.7

在经典统计学中，我们通常使用大约2的t检验阈值来检查系数是否显著。然而，这些结果显示，在这种特定情况下，您需要使用11.4这个阈值才能正确检验显著性。

2 × (4.59 / 0.8)，这等于 11.4

我们刚刚展示，如果模型中包括不属于它的变量——特别是随机游走过程——可能会导致对这些系数的显著性检验完全无效。

为了使他们的模拟更加清楚明了，Granger 和 Newbold (1974) 进行了一系列回归分析，使用的是遵循随机游走或 ARIMA(0,1,1) 模型过程的变量。

他们就是这样设置模拟实验的：

他们将依赖序列 𝐘ₜ 对 m 个独立序列 𝐗ⱼ,ₜ（其中 j = 1, 2, …, m）进行了回归分析，将 m 的取值从 1 到 5 进行变化。依赖序列 𝐘ₜ 和独立序列 𝐗ⱼ,ₜ 遵循相同的进程类型，他们测试了四种不同的情况：

• 情况1（水平值）： 𝐘ₜ 和 𝐗ⱼ,ₜ 遵循随机游走
• 情况2（差分）：它们使用随机游走的第一差分，这些差分是平稳的时间序列
• 情况3（水平值）： 𝐘ₜ 和 𝐗ⱼ,ₜ 遵循 ARIMA(0,1,1) 模型
• 情况4（差分）：它们使用先前 ARIMA(0,1,1) 模型的第一差分

每个系列有50个观察值，每种情况下他们都进行了100次模拟。

所有误差项都符合𝑁(0,1)分布，而ARIMA(0,1,1)序列是由随机漫步与独立白噪声相加得到的。基于100次重复，每次序列长度为50的模拟结果将在下表中总结。

表2：针对m个独立的“解释”变量的回归分析。

结果的解读：

• 可以看出，当使用随机游走序列（rw-levels）进行回归时，当 ( m \geq 3 ) 时，不拒绝 ( Yt ) 和 ( X{j,t} ) 之间无关系的原假设的概率变得非常小。此时 ( R^2 ) 和平均Durbin-Watson值会增加。类似地，当使用ARIMA(0,1,1)序列（arima-levels）进行回归时，会得到类似的结果。
• 当使用白噪声序列（rw-diffs）时，经典的回归分析是有效的，因为误差序列将是白噪声，最小二乘法将是最有效的。
• 然而，当使用ARIMA(0,1,1)序列的差分（arima-diffs）或一阶移动平均序列MA(1)进行回归时，平均来说，会拒绝原假设：

这个计算的结果是 (10 + 16 + 5 + 6 + 6) / 5 等于 8.6

哪一项超过了5%的时间？

如果你的变量接近或符合随机游走，并且在回归分析中加入了不必要的变量，你容易得到错误的结果。高𝑅²和低Durbin-Watson值并不意味着存在真实的关系，而是可能表明一种虚假关系。

5. 如何避免时间序列中的虚假回归

确实很不容易列出所有避免虚假回归的途径或方法。不过，有一些好的做法可以遵循，以尽可能地减少这种风险。

如果当用时间序列数据做回归分析时，发现残差自相关性很强，则解释方程系数时会遇到严重问题。为了检查残差自相关性，可以使用杜宾-瓦森检验或波特曼检验。

基于上述研究，我们可以得出结论，使用经济变量进行的回归分析会产生高度自相关的残差，这表示杜宾-沃森统计量较低。无论观察到的 R² 值是多少，分析的结果很可能是错误的或者是无意义的。

在这种情况下，了解这种偏差的来源很重要。根据文献，模型设定误差通常分为三类：（i）遗漏重要变量，（ii）引入无关变量，或（iii）误差项的自相关性。大多数时候，这种错误是这三种因素混合的结果。

为了避免时间序列中的伪回归，这里有一些有用的建议：

• 第一条建议是选择合适的宏观经济变量来解释因变量。这可以通过查阅相关文献或咨询该领域的专家来完成。
• 第二条建议是使序列平稳化，通常通过取一阶差分来实现。在处理宏观经济数据时，强烈建议对序列进行一次差分以减少残差的自相关性，特别是在当样本量较小时。这些变量有时确实会观察到强烈的自相关性。简单的计算表明，一阶差分几乎总是比原始序列自相关性更小。
• 第三条建议是使用博克斯-詹金斯法分别对每个宏观经济变量建模，然后通过比较每个单独模型的残差来寻找序列之间的关系。这里的想法是，博克斯-詹金斯过程提取了序列可解释的部分，留下仅包含无法由序列自身解释的未解释部分。这使得更容易检查这些未解释的部分在不同变量之间是否存在相关性。

许多计量经济学教科书都警告回归模型中的设定错误问题，但这个问题仍然出现在许多已发表的论文中。Granger和Newbold (1974) 强调了虚假回归的风险，即你可能会得到一个很高的R方值，但同时却非常低的Durbin-Watson统计量。

通过Python模拟，我们展示了这些虚假回归现象的一些主要原因，特别是包括不应在模型中的高度自相关的变量。我们还展示了这些问题是如何完全扭曲了对系数的假设检验结果，使句子意思更完整。

希望能帮助大家在未来做计量经济学分析时避免遇到伪回归的问题。

表1的模拟代码段.

======================

    import numpy as np  
    import pandas as pd  
    import statsmodels.api as sm  
    import matplotlib.pyplot as plt  

    np.random.seed(123)  
    M = 100   
    n = 50  
    S = np.zeros(M)  
    for i in range(M):  
    #---------------------------------------------------------------  
    # 生成数据部分  
    #---------------------------------------------------------------  
        epsilon_y = np.random.normal(0, 1, n)  
        epsilon_x = np.random.normal(0, 1, n)  
        Y = np.cumsum(epsilon_y)  
        X = np.cumsum(epsilon_x)  
    #---------------------------------------------------------------  
    # 拟合线性回归模型  
    #---------------------------------------------------------------  
        X = sm.add_constant(X)  
        model = sm.OLS(Y, X).fit()  
    #---------------------------------------------------------------  
    # 计算统计值  
    #------------------------------------------------------  
        S[i] = np.abs(model.params[1])/model.bse[1]  

    #------------------------------------------------------   
    # S的最大值是多少  
    #------------------------------------------------------  
    S_max = int(np.ceil(max(S)))  
    #------------------------------------------------------   
    # 创建区间（bins）  
    #------------------------------------------------------  
    bins = np.arange(0, S_max + 2, 1)    
    #------------------------------------------------------  
    # 计算S的直方图  
    #------------------------------------------------------  
    frequency, bin_edges = np.histogram(S, bins=bins)  
    #------------------------------------------------------  
    # 创建数据框（dataframe）  
    #------------------------------------------------------  
    df = pd.DataFrame({  
        "S区间": [f"{int(bin_edges[i])}-{int(bin_edges[i+1])}" for i in range(len(bin_edges)-1)],  
        "频次": frequency  
    })  
    print(df)  
    print(np.mean(S))

分割线

    import numpy as np  
    import pandas as pd  
    import statsmodels.api as sm  
    from statsmodels.stats.stattools import durbin_watson  
    from tabulate import tabulate  

    np.random.seed(1)  # 使得结果可重复  
    #------------------------------------------------------  
    # 函数定义  
    #------------------------------------------------------  
    def generate_random_walk(T):  
        """  
        生成长度为T的随机游走序列：  
            Y_t = Y_{t-1} + e_t，  
        其中 e_t ~ N(0,1)。  
        """  
        e = np.random.normal(0, 1, size=T)  
        return np.cumsum(e)  
    def generate_arima_0_1_1(T):  
        """  
        生成一个ARIMA(0,1,1)模型的序列，按照Granger和Newbold的方法：  
        序列由一个随机游走序列和一个独立的白噪声序列相加得到。  
        """  
        rw = generate_random_walk(T)  
        wn = np.random.normal(0, 1, size=T)  
        return rw + wn  
    def difference(series):  
        """  
        计算一维序列的一阶差。  
        返回长度为T-1的序列。  
        """  
        return np.diff(series)  
    #------------------------------------------------------  
    # 参数  
    #------------------------------------------------------  
    T = 50           # 每个序列的长度  
    n_sims = 100     # 蒙特卡洛模拟次数  
    alpha = 0.05     # 显著性水平  
    #------------------------------------------------------  
    # 模拟函数定义  
    #------------------------------------------------------  
    def run_simulation_case(case_name, m_values=[1,2,3,4,5]):  
        """  
        case_name : 生成类型标识符：  
            - 'rw-levels' : 随机游走（原序列）  
            - 'rw-diffs' : 随机游走的一阶差（白噪声）  
            - 'arima-levels' : ARIMA(0,1,1)模型的序列（原序列）  
            - 'arima-diffs' : ARIMA(0,1,1)模型的序列的一阶差（MA(1)）  

        m_values : 自变量数量列表。  

        返回一个DataFrame，包含每个m值下的：  
            - 拒绝零假设的百分比  
            - 平均Durbin-Watson统计量  
            - 平均调整的R方  
            - R方大于0.7的百分比  
        """  
        results = []  

        for m in m_values:  
            count_reject = 0  
            dw_list = []  
            r2_adjusted_list = []  

            for _ in range(n_sims):  
    #--------------------------------------  
    # 1) 生成Y_t和X_{j,t}的独立随机序列。  
    #----------------------------------------  
                if case_name == 'rw-levels':  
                    Y = generate_random_walk(T)  
                    Xs = [generate_random_walk(T) for __ in range(m)]  

                elif case_name == 'rw-diffs':  
                    # Y和X是一阶差的随机游走序列，即～白噪声  
                    Y_rw = generate_random_walk(T)  
                    Y = difference(Y_rw)  
                    Xs = []  
                    for __ in range(m):  
                        X_rw = generate_random_walk(T)  
                        Xs.append(difference(X_rw))  
                    # 注意：现在Y和Xs的长度为T-1  
                    # => 实际上我们将使用T-1个点进行回归  

                elif case_name == 'arima-levels':  
                    Y = generate_arima_0_1_1(T)  
                    Xs = [generate_arima_0_1_1(T) for __ in range(m)]  

                elif case_name == 'arima-diffs':  
                    # ARIMA(0,1,1)模型的序列的一阶差=> MA(1)  
                    Y_arima = generate_arima_0_1_1(T)  
                    Y = difference(Y_arima)  
                    Xs = []  
                    for __ in range(m):  
                        X_arima = generate_arima_0_1_1(T)  
                        Xs.append(difference(X_arima))  

                # 2) 准备回归数据  
                #    根据情况，长度为T或T-1  
                if case_name in ['rw-levels','arima-levels']:  
                    Y_reg = Y  
                    X_reg = np.column_stack(Xs) if m>0 else np.array([])  
                else:  
                    # 在差分的情况下，长度为T-1  
                    Y_reg = Y  
                    X_reg = np.column_stack(Xs) if m>0 else np.array([])  

                # 3) OLS回归  
                X_with_const = sm.add_constant(X_reg)  # 添加截距  
                model = sm.OLS(Y_reg, X_with_const).fit()  

                # 4) F检验：H0：所有beta_j=0  
                #    查看p-value<alpha  
                if model.f_pvalue is not None and model.f_pvalue < alpha:  
                    count_reject += 1  

                # 5) R^2, Durbin-Watson统计量  
                r2_adjusted_list.append(model.rsquared_adj)  

                dw_list.append(durbin_watson(model.resid))  

            # n_sims次模拟的统计量  
            reject_percent = 100 * count_reject / n_sims  
            dw_mean = np.mean(dw_list)  
            r2_mean = np.mean(r2_adjusted_list)  
            r2_above_0_7_percent = 100 * np.mean(np.array(r2_adjusted_list) > 0.7)  

            results.append({  
                'm': m,  
                'Reject %': reject_percent,  
                'Mean DW': dw_mean,  
                'Mean R^2': r2_mean,  
                '% R^2_adj>0.7': r2_above_0_7_percent  
            })  

        return pd.DataFrame(results)  

    #------------------------------------------------------  
    # 应用模拟  
    #------------------------------------------------------         
    cases = ['rw-levels', 'rw-diffs', 'arima-levels', 'arima-diffs']  
    all_results = {}  
    for c in cases:  
        df_res = run_simulation_case(c, m_values=[1,2,3,4,5])  
        all_results[c] = df_res  
    #------------------------------------------------------  
    # 存储结果到表格中  
    #------------------------------------------------------  
    for case, df_res in all_results.items():  
        print(f"\n\n{case}：“")  
        print(tabulate(df_res, headers='keys', tablefmt='fancy_grid'))

Granger, C. W. J., 和 Paul Newbold. 1974. “计量经济学中的虚假回归问题.” Journal of Econometrics 卷2 (期2): 111–20.

Knowles, EAG. 1954. “理论统计练习集。”牛津大学出版社。