如何手动计算自回归(AR)模型的残差？

要手动复现stats::ar输出的残差，核心是对齐该函数的模型设定逻辑，以下是具体实现步骤和代码：

模型逻辑说明 #

stats::ar默认开启demean=TRUE，即先对时间序列去均值，再拟合AR模型。对于多变量AR(1)模型，残差的计算逻辑为：

1. 对原始序列中心化：$\tilde{y}_t = y_t - \mu$，其中$\mu$是序列的均值（对应模型的x.mean属性）
2. 残差 = 当前中心化值 - AR系数矩阵与前一期中心化值的乘积：$\epsilon_t = \tilde{y}t - \Phi \tilde{y}{t-1}$，其中$\Phi$是AR系数矩阵（对应模型的ar属性）

手动计算残差的函数 #

compute_ar_residuals <- function(ts_data, ar_coef, x_mean) {
  n <- nrow(ts_data)
  n_var <- ncol(ts_data)
  # 初始化残差矩阵，第一行设为NA（无前期观测）
  resid_mat <- matrix(NA, nrow = n, ncol = n_var)
  colnames(resid_mat) <- colnames(ts_data)
  rownames(resid_mat) <- rownames(ts_data)
  
  # 生成中心化序列
  ts_centered <- ts_data - matrix(x_mean, nrow = n, ncol = n_var, byrow = TRUE)
  
  # 逐期计算残差
  for (t in 2:n) {
    # 计算AR模型的中心化拟合值
    fit_centered <- ar_coef %*% as.matrix(ts_centered[t-1, ])
    # 计算并赋值残差
    resid_mat[t, ] <- ts_centered[t, ] - t(fit_centered)
  }
  
  # 转换为时间序列对象，保留原始时序属性
  ts(resid_mat, start = start(ts_data), frequency = frequency(ts_data))
}

验证与使用示例 #

复现usconsumption的残差 #

library(fpp)
# 拟合AR(1)模型
ar_model <- ar(usconsumption, aic=FALSE, order.max=1)
# 手动计算残差
manual_resid <- compute_ar_residuals(usconsumption, ar_model$ar, ar_model$x.mean)

# 对比结果（应返回TRUE）
all.equal(manual_resid, ar_model$resid, check.attributes = FALSE)
# 查看前几行结果
head(manual_resid)

跨数据集应用AR系数 #

如果要将已估计的AR系数应用到另一个同维度的时间序列new_ts，直接传入参数即可：

# 假设new_ts是与usconsumption同维度的时间序列
new_resid <- compute_ar_residuals(new_ts, ar_model$ar, ar_model$x.mean)

使用自定义AR系数 #

若要替换为其他AR系数（比如调整后的系数），只需修改ar_coef参数：

# 自定义1×2×2的AR系数矩阵
custom_coef <- array(c(0.3, 0.5, 0.1, -0.2), dim = c(1,2,2))
custom_resid <- compute_ar_residuals(usconsumption, custom_coef, ar_model$x.mean)

关于timsac::mfilter的问题 #

timsac::mfilter未得到一致结果，大概率是因为该函数的模型设定与stats::ar不同——比如默认是否去均值、AR模型的拟合方向或残差计算逻辑存在差异，需要对齐模型设定才能匹配结果。

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