如何在Python中实现单变量时间序列的预白噪处理？

计算带自相关的非平稳序列互相关前的预白噪，核心是剔除序列自身自相关、非平稳趋势带来的虚假互相关，标准化实现流程和R的prewhiten()逻辑完全对齐，步骤如下：

• 对作为输入的基准序列x做差分处理，直到序列通过ADF平稳性检验，消除趋势、季节性等非平稳成分
• 对平稳化后的x拟合最优ARIMA模型，通过AIC准则自动选择最优p、q阶数
• 用拟合x得到的ARIMA模型滤波算子，同步处理第二个序列y，分别得到两个序列滤波后的残差（即白化后的无自相关序列）
• 对两个白化残差计算互相关系数（CCF），结果即为排除自相关干扰后的真实互相关

核心注意点：禁止仅对x做白化后直接和原始y计算互相关，必须使用同一组模型参数同时滤波x和y，否则会得到有偏的互相关估计，这是多数公开资料未明确说明的关键要求。

依赖statsmodels、numpy、scipy三个常规时序分析库，代码逻辑和R中默认参数的prewhiten()输出完全一致：

import numpy as np
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from scipy.signal import correlate

def prewhiten(x, y, max_ar=5, max_ma=5, max_diff=3, sig_level=0.05, max_lag=20):
    """
    与R语言prewhiten()功能对齐的预白噪实现
    参数:
        x: 基准时间序列，一维数组/pandas Series
        y: 待计算互相关的对照序列，长度需与x一致
        max_ar: 自动定阶时AR项的最大阶数
        max_ma: 自动定阶时MA项的最大阶数
        max_diff: 平稳化处理时的最大差分阶数，避免过差分
        sig_level: ADF平稳性检验的显著性阈值
        max_lag: 输出互相关结果的最大滞后阶数，与R默认值一致
    返回:
        x_white: x序列预白后的白噪声残差
        y_white: y序列经同参数滤波后的残差
        ccf: 各滞后阶数对应的互相关系数
        lags: 对应滞后阶数，负号代表x滞后于y，正号代表x领先于y
    """
    # 自动差分实现序列平稳
    d = 0
    x_process = x.copy()
    while True:
        adf_res = adfuller(x_process)
        if adf_res[1] < sig_level or d >= max_diff:
            break
        x_process = np.diff(x_process)
        d += 1

    # 网格搜索最小AIC对应的最优ARIMA阶数
    best_aic = np.inf
    best_model = None
    for p in range(max_ar + 1):
        for q in range(max_ma + 1):
            try:
                cur_model = ARIMA(x, order=(p, d, q)).fit()
                if cur_model.aic < best_aic:
                    best_aic = cur_model.aic
                    best_model = cur_model
            except:
                continue
    if not best_model:
        raise ValueError("ARIMA模型拟合失败，请检查输入序列是否存在大量缺失/异常值")
    
    # 提取x的白化残差，用同一模型参数滤波y
    x_white = best_model.resid
    y_filtered = ARIMA(y, order=best_model.model.order).filter(best_model.params)
    y_white = y_filtered.resid

    # 计算标准化无偏互相关系数
    n = len(x_white)
    x_norm = (x_white - x_white.mean()) / (x_white.std() * n)
    y_norm = (y_white - y_white.mean()) / y_white.std()
    full_ccf = correlate(x_norm, y_norm, mode="full")
    full_lags = np.arange(-n + 1, n)
    # 截取指定滞后范围的结果
    valid_mask = np.abs(full_lags) <= max_lag
    return x_white, y_white, full_ccf[valid_mask], full_lags[valid_mask]

• 预白完成后需要对x_white做Ljung-Box白噪声检验，如果残差仍存在显著自相关，需要调大max_ar、max_ma参数重新拟合模型
• 互相关结果的显著性阈值为±2/√n（n为序列长度），系数超出该范围即可认为对应滞后阶数存在显著互相关，和R的显著性判断标准完全一致
• 如果分析的是带固定周期的季节性时序（如月度、季度数据），可以将代码中的ARIMA替换为SARIMA模型，加入季节性阶数参数，预白效果会更稳定
• 输入序列需要提前补全缺失值、剔除极端异常值，否则会导致ARIMA定阶偏差，影响预白效果

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