在R语言中手动实现OLS高阶矩估计的问题求助

我尝试在R语言中手动实现一个自定义OLS估计量，目的是研究epsilon的高阶矩计算带来的影响，因此采用递归方式完成实验，使用moment函数计算矩。我的代码如下：

library(moments)
x=list()
epsilon=list()
y=list()
eta=list()

for (i in 1:1000) {
  x[[i]]=rnorm(100, mean=1, sd=1)
  epsilon[[i]]=rnorm(100, mean=0, sd=1)
}
i=c(1,2,4,8)
for (k in i) {
  eta[[k]]=lapply(epsilon[[k]], moment, order=k, central=T)
}
eta=eta %>% keep( ~ !is.null(.) )

for (k in i) {
  y[[k]]=mapply(function(x, eta) x[[k]]+eta[[k]], x, eta)
  
}
y=y %>% keep( ~ !is.null(.) )

}
estimator=list()
for (k in i) {
  estimator[[k]]=mapply(function(x,y) inv(t(x[[k]])%*%x[[k]])%*%(t(x[[k]])%*%y[[k]]), x, y )
}

当前代码存在两个核心问题：

• 执行t(x[[1]]%*%x[[1]])时得到的是标量而非100×100的矩阵
• y应该是包含4个元素的列表，每个元素对应一组样本，每组样本由1000个长度为100的向量组成，以保证OLS估计量计算时维度匹配

修正方案与代码 #

核心问题解析

1. 矩阵维度错误：x[[i]]是单个100维观测向量，直接做x[[i]] %*% x[[i]]会计算向量内积得到标量。OLS计算需要的是所有观测组成的设计矩阵的交叉乘积t(X) %*% X，其中X是1000×100的矩阵（1000个观测，每个观测100维）。
2. y的构造逻辑错误：原代码中对单个epsilon元素计算矩，而非对每组epsilon整体计算k阶中心矩，导致eta和y的维度完全不匹配。

修正后的代码

library(moments)
library(purrr)

# 初始化列表，明确长度与命名
x <- vector("list", 1000)
epsilon <- vector("list", 1000)
k_values <- c(1, 2, 4, 8)
eta <- setNames(vector("list", 4), k_values)
y <- setNames(vector("list", 4), k_values)
estimator <- setNames(vector("list", 4), k_values)

# 生成1000组100维的观测向量
for (i in 1:1000) {
  x[[i]] <- rnorm(100, mean = 1, sd = 1)
  epsilon[[i]] <- rnorm(100, mean = 0, sd = 1)
}

# 计算每组epsilon的k阶中心矩
for (k in k_values) {
  eta[[as.character(k)]] <- lapply(epsilon, function(eps) moment(eps, order = k, central = TRUE))
}

# 构造符合维度要求的y：每个k对应1000个100维向量
for (k in k_values) {
  y[[as.character(k)]] <- mapply(function(x_i, eta_i) x_i + eta_i, 
                                 x, eta[[as.character(k)]], 
                                 SIMPLIFY = FALSE)
}

# 计算OLS估计量
for (k in k_values) {
  # 将观测列表转为设计矩阵（1000×100）
  X <- do.call(rbind, x)
  # 将y列表转为响应向量（1000×1）
  y_vec <- do.call(c, y[[as.character(k)]])
  
  # OLS核心计算：beta = (X'X)^(-1)X'y
  XtX <- t(X) %*% X
  Xty <- t(X) %*% y_vec
  estimator[[as.character(k)]] <- solve(XtX) %*% Xty
}

关键修改点

• 设计矩阵构造：用do.call(rbind, x)将1000个100维向量组合成1000×100的设计矩阵X，此时t(X) %*% X自然得到100×100的交叉乘积矩阵。
• eta与y的维度修正：对每组epsilon整体计算k阶中心矩，得到1000个标量，再与对应组的x向量相加，生成的y[[k]]是包含1000个100维向量的列表，完全匹配OLS计算的维度要求。
• 简化代码逻辑：提前初始化列表并命名，避免冗余的keep操作；使用R原生的solve()进行矩阵求逆，无需额外依赖。

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