部分线性模型模拟研究（现代统计模型）

估计方法

library(magrittr)

#样本量 

n=100

#循环次数
MN = 500

## 设置核函数，EV核
K <- function(x){
  0.75*(1-x^2)*(abs(x) <= 1) 
}
g <- function(U) {
  50*(U-0.5)^2
}
U <- runif(n, 0, 1)
gu=g(U)

## 设置结果矩阵
beta_NW <- matrix(NA, nrow = MN, ncol = 2)
gu_NW <- matrix(NA, nrow = MN, ncol = n)
MSE_NW <- c()
beta_LL <- matrix(NA, nrow = MN, ncol = 2)
gu_LL <- matrix(NA, nrow = MN, ncol = n)
MSE_LL <- c()

for (a in 1:MN) {
  X <- matrix(rnorm(2 * n), nrow = n, ncol = 2)
  beta <- c(0.8,0.6)
  err <- rnorm(n, 0, 0.6)
  Y = X%*%beta+gu+err
  
  ## 设置带宽
  h = sd(U)*n^(-1/5)
  
  ## 方法一核权，估计g1=E(Y|U=u)和g2=E(X|U=u)
  g1_NW <- c()
  for (i in 1:n) {
    u1 <- U[abs(U-U[i])<=h]
    y1 <- Y[abs(U-U[i])<=h]
    Kvec <- K((u1-U[i])/h)
    length(y1)=length(Kvec)
    S <- sum(Kvec)
    if(S == 0){S = S+0.0001}
    g1_NW[i] <- sum(Kvec*y1)/S
  }
  
  g2_NW<-matrix(NA,nrow = n, ncol = 2)
  for (i in 1:n) {
    u1 <- U[abs(U-U[i])<=h]
    x2 <- X[abs(U-U[i])<=h,]
    Kvec <- K((u1-U[i])/h)
    S <-sum(Kvec)
    if(S == 0){S = S+0.0001}
    g2_NW[i,] <- colSums(Kvec*x2)/S
  }
  
  ## 估计beta
  beta_NW[a,] <- solve(t(X-g2_NW)%*%(X-g2_NW))%*%t(X-g2_NW)%*%(Y-g1_NW)
  
  ## g(u)的最终估计量
  gu_NW[a,] <- g1_NW-g2_NW%*%beta_NW[a,]
  
  ## 非参数估计的拟合结果度量
  MSE_NW[a] <- mean((g(sort(U))-gu_NW[a,][order(U)])^2)
  
  ## 方法二局部线性权，估计g1=E(Y|U=u)和g2=E(X|U=u)
  g1_LL <- c()
  for (i in 1:n) {
    u1 <- U[abs(U-U[i])<=h]
    y1 <- Y[abs(U-U[i])<=h]
    Kvec <- K((u1-U[i])/h)
    S_n_0 <- mean(K((u1-U[i])/h))
    S_n_1 <- mean((u1-U[i])*K((u1-U[i])/h))
    S_n_2 <- mean((u1-U[i])^2*K((u1-U[i])/h))
    if(S_n_0*S_n_2-S_n_1^2 == 0){S_n_0*S_n_2-S_n_1^2 = S_n_0*S_n_2-S_n_1^2+0.0001}
    g1_LL[i] <- mean(Kvec*(S_n_2-(u1-U[i])*S_n_1)*y1)/(S_n_0*S_n_2-S_n_1^2)
  }
  
  g2_LL<-matrix(NA,nrow = n, ncol = 2)
  for (i in 1:n) {
    u1 <- U[abs(U-U[i])<=h]
    x1 <- X[abs(U-U[i])<=h,]
    Kvec <- K((u1-U[i])/h)
    S_n_0 <- mean(K((u1-U[i])/h))
    S_n_1 <- mean((u1-U[i])*K((u1-U[i])/h))
    S_n_2 <- mean((u1-U[i])^2*K((u1-U[i])/h))
    if(S_n_0*S_n_2-S_n_1^2 == 0){S_n_0*S_n_2-S_n_1^2 = S_n_0*S_n_2-S_n_1^2+0.0001}
    g2_LL[i,] <- colMeans(Kvec*(S_n_2-(u1-U[i])*S_n_1)*x1)/(S_n_0*S_n_2-S_n_1^2)
  }
  
  ## 估计beta
  beta_LL[a,] <- solve(t(X-g2_LL)%*%(X-g2_LL))%*%t(X-g2_LL)%*%(Y-g1_LL)
  
  ## g(u)的最终估计量
  gu_LL[a,] <- g1_LL-g2_LL%*%beta_LL[a,]
  
  ## 非参数估计的拟合结果度量
  MSE_LL[a] <- mean((g(sort(U))-gu_LL[a,][order(U)])^2)
}

## NW参数部分结果
result_NW <- as.matrix(cbind(colMeans(beta_NW)-beta, 
                          c(sd(beta_NW[,1]),sd(beta_NW[,2])),
                          colMeans((beta_NW-beta)^2)))
colnames(result_NW) <- c("Bias","SD","MSE")
rownames(result_NW) <- c("beta1_NW","beta2_NW")

## LL参数部分结果
result_LL <- as.matrix(cbind(colMeans(beta_LL)-beta, 
                             c(sd(beta_LL[,1]),sd(beta_LL[,2])),
                             colMeans((beta_LL-beta)^2)))
colnames(result_LL) <- c("Bias","SD","MSE")
rownames(result_LL) <- c("beta1_LL","beta2_LL")
result_combined <- rbind(result_NW, result_LL)

## 非参数部分拟合结果图形展示
op=par(mfrow = c(1,2))
plot(sort(U), g(sort(U)), # 图1：曲线拟合图 
     type = "l", 
     xlab = "u", 
     ylab = "g(u)", 
     xlim = c(0,1), 
     ylim = c(min(g(U))-1, max(g(U))+3))
lines(sort(U), colMeans(gu_NW)[order(U)],type = "l", col = "blue", lty = 2, lwd = 2)#NW
lines(sort(U), colMeans(gu_LL)[order(U)],type = "l", col = "red", lty = 3, lwd = 2)#LL
grid()
legend("topright", 
       legend = c(expression(g(u)),expression(hat(g)(u)("NW")), expression(hat(g)(u)("LL"))), 
       col = c("black","blue","red"),
       lty = c(1,2,3),
       lwd = 2,
       cex = 0.6)
title(sub="(a)")

boxplot(data.frame(MSE_NW,MSE_LL),xlab = "",ylab = "MSE",
        range=3.5, names=c(expression(MSE[1]-NW),expression(MSE[2]-LL)))
grid()
title(sub="(b)")

参数和非参数估计结果