Python中含随机NaN的滞后时间序列自相关值大于1问题排查

我正在对一组含随机NaN的N×N数组做滞后自相关分析，其中eke_array[i,i]对应滞后0，eke_array[i,i+n]对应滞后n。但分析结果在滞后40阶后出现自相关值大于1的异常，仅前几阶结果合理。以下是我的实现函数：

def auto_corr(eke_array):

    nlag = np.shape(eke_array)[0]
    eke_lag_0 = []
    auto_store = []
    std_store = []
    for mylag in range(0,80):
        autocorr = 0
        numvalid = 0
        norm = 0
        std = 0
        for i in range(0, 105):
            if i+mylag > 104:
                break
            lag0 = eke_array[i,i]
            lag_itt = eke_array[i,i+mylag]
            if not np.isnan(lag_itt):
                if not np.isnan(lag0):
                    numvalid += 1

                    diff = lag0 - lag_itt
                    std = std + diff**2
                    autocorr = autocorr + lag0*lag_itt
        print(numvalid, mylag,i )
        auto_store.append((autocorr/numvalid))
        std_store.append(std / numvalid)

    for k in range(105):
        eke_lag_0.append(eke_array[k,k])
    eke_lag_0 = np.asarray(eke_lag_0)
    eke_lag_0 = eke_lag_0[~np.isnan(eke_lag_0)]
    norm = np.sum(eke_lag_0 ** 2) / np.size(eke_lag_0)

    std_store = np.sqrt(np.asarray(std_store))
    auto_store = auto_store/norm
    return std_store, auto_store, norm

问题根源 #

• 归一化逻辑偏差：代码用全局滞后0阶的均值平方（norm）归一化所有滞后阶的协均值，但正确的自相关计算需要基于对应滞后阶的有效配对样本计算协方差与方差的比值。全局归一化忽略了高滞后阶样本量骤减带来的波动，导致结果超出[-1,1]范围。
• 样本量不足引发波动：滞后40阶时，有效配对样本数（numvalid）大幅减少，再加上NaN过滤，样本量过小会让协均值计算结果波动剧烈，甚至出现绝对值大于全局norm的情况，最终归一化后超过1。
• 标准差计算逻辑错误：当前std_store计算的是均方误差的平方根，并非序列标准差，用于置信区间等场景会产生偏差。

修正方案 #

正确自相关计算逻辑 #

自相关系数的标准公式为：
$$
r_k = \frac{E[(X_t - \mu)(X_{t+k} - \mu)]}{\sigma^2}
$$
若序列为零均值，可简化为：
$$
r_k = \frac{E[X_t X_{t+k}]}{E[X_t^2]}
$$
核心是保证协方差与方差的计算基于同一组有效配对样本，而非全局样本。

修正后的代码 #

import numpy as np

def auto_corr(eke_array):
    n = eke_array.shape[0]
    max_lag = 80
    auto_store = []
    std_store = []

    # 计算全局滞后0阶的均值与方差
    lag0_vals = eke_array.diagonal()
    valid_lag0 = lag0_vals[~np.isnan(lag0_vals)]
    mu0 = np.mean(valid_lag0)
    var0 = np.var(valid_lag0, ddof=0)

    for mylag in range(max_lag):
        x_t = []
        x_tk = []
        # 收集当前滞后阶的有效配对样本
        for i in range(n - mylag):
            val0 = eke_array[i, i]
            valk = eke_array[i, i + mylag]
            if not np.isnan(val0) and not np.isnan(valk):
                x_t.append(val0)
                x_tk.append(valk)
        
        x_t = np.array(x_t)
        x_tk = np.array(x_tk)
        numvalid = len(x_t)
        
        if numvalid == 0:
            auto_store.append(np.nan)
            std_store.append(np.nan)
            continue
        
        # 计算去均值后的协方差
        cov = np.mean((x_t - mu0) * (x_tk - mu0))
        # 自相关系数 = 协方差 / 滞后0阶方差
        autocorr = cov / var0
        auto_store.append(autocorr)
        
        # 计算当前配对样本的标准差
        std = np.sqrt(np.var(x_tk - x_t, ddof=0))
        std_store.append(std)
        
        print(numvalid, mylag)

    auto_store = np.array(auto_store)
    std_store = np.array(std_store)
    return std_store, auto_store, var0

关键改进点 #

• 对每个滞后阶，仅使用同时有效的$(X_t, X_{t+k})$配对样本计算协方差，避免样本不匹配导致的偏差。
• 采用去均值后的协方差计算自相关，符合标准定义，且根据柯西-施瓦茨不等式，结果会被限制在[-1,1]范围内。
• 有效样本数为0时填充NaN，避免除以0错误。
• 样本量不足时结果会有波动，但不会超出合理范围。

额外建议 #

• 可添加置信区间计算，比如用$\pm 2/\sqrt{N}$（N为有效样本数）判断自相关是否显著。
• 高滞后阶有效样本数过少时，结果可靠性低，可考虑截断或标注结果不可靠。

内容的提问来源于stack exchange，提问作者drwnthembirds