R中Arima()与Python SARIMAX多元预测结果不一致问题求助

我尝试用Python复现R中Arima()的多元预测结果，已为Python的SARIMAX设置与R中Arima()相同的p、d、q阶数，且将SARIMAX的优化方法改为R中Arima默认的bfgs（原SARIMAX默认lbfgs），但仍无法得到一致结果。恳请帮助排查原因，或告知可复现R中Arima()结果的Python函数。

R代码实现 #

h1 <- dplyr::select(sample_data_dev, x1,x2,x3) %>% as.matrix() 
ts_h <-  ts(sample_data_dev$y,
            start=c(year_start_dev,quarter_start_dev),
            end=c(year_end_dev,quarter_end_dev),
            frequency = 4)

h_model1 <- Arima(y = ts_h, order = c(1,1,1),season = c(0,0,0),xreg = h1) 
h_pred1 <- as_tibble((forecast::forecast(h_model1,
                   xreg=dplyr::select(out_data_forecast[1:h,], x1,x2,x3) %>% as.matrix()))$mean)

Python代码实现 #

import statsmodels.api as sm
import pandas as pd

dev_start = '1992-06-30'  
dev_end = '2017-12-31'     

date_range = pd.date_range(start= dev_start,end= dev_end,freq= 'Q')

h_data = pd.Series(sample_data_dev['y'],index=date_range)

from statsmodels.tsa.statespace.sarimax import SARIMAX

# 注意：原代码参数名错误，应为seasonal_order而非season
h_model = SARIMAX(endog = h_data, exog = hel , method = 'bfgs',order = (1,1,1), seasonal_order = (0,0,0,4))

h_results = h_model.fit()

forecast_steps = h
forecast_x =  out_data_forecast.loc[0:h-1, ['x1','x2','x3']]
 
forcst = h_results.get_forecast(steps = h, exog=forecast_x)
forecaset_mean = forcst.predicted_mean

结果对比 #

排查与解决要点 #

• 参数名修正：Python中SARIMAX的季节性参数是seasonal_order=(P,D,Q,s)，其中s是时间序列频率（此处为4，对应季度数据），原代码错误使用了season参数，需改为seasonal_order=(0,0,0,4)以匹配R的season=c(0,0,0)设置。
• 自变量中心化：R的Arima()默认会对xreg中的自变量进行中心化（减去均值），而Statsmodels的SARIMAX默认不执行此操作。需手动对训练集exog和预测集forecast_x做中心化处理：

  # 对训练集自变量中心化
  exog_mean = hel.mean()
  hel_centered = hel - exog_mean
  # 对预测集自变量用相同均值中心化
  forecast_x_centered = forecast_x - exog_mean
  

• 初始化方式匹配：R的Arima()使用CSS-ML混合方法初始化参数，Statsmodels可通过fit(start_params='approximate')模拟该逻辑：

  h_results = h_model.fit(start_params='approximate')
  

• 数据一致性检查：确保R与Python中的输入数据完全一致，包括时间戳对齐、缺失值处理（R的ts与Python的pd.Series对缺失值的处理逻辑需统一）。
• 替代库推荐：若上述调整后仍有差异，可尝试使用pmdarima库，其auto_arima函数更贴近R的forecast::Arima实现，支持设置method='bfgs'、seasonal=True等参数来匹配R的行为。

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