Как создать функцию Монте-Карло, которая имитирует непротиворечивость бета-оценщика?

Несколько вопросов:

1. e <- rnorm(1)*x означает, что все ошибки коррелированы. Вместо этого используйте e <- rnorm(observations)*x.
2. Переменные должны быть сохранены вне цикла for.
3. Так как член ошибки не является постоянным. Взвешенный метод наименьших квадратов больше подходит.
4. n=200 и 1000 не будут иметь большого значения. Вместо этого попробуйте большую разницу.

Код:

func <- function(iterations, observations){
  set.seed(1001)
  
  betas <-c()
  for (i in 1:iterations) {
    c <- 1:observations
    x <- c/5
    e <- rnorm(observations)*x # Remove x if you want homoskedastic errors
    y <- x - 10 + e
    #y[y < 0] <- 0 # censoring
    model <- lm(y ~ x,weights=x)
    betas[i] <- model$coefficients[2]
    
    
    
  }
  #name:
  arg_name <- deparse(substitute(observations))
  var_name <- paste("obs", arg_name, sep = "_") # Construct the name
  assign(var_name, betas, env=.GlobalEnv) # Assign values to variable
  # variable will be created in .GlobalEnv 
}

Тестирование:

func(iterations = 2000, observations = 1000)
func(iterations = 2000, observations = 10000)
mean(obs_1000)
[1] 1.00306
mean(obs_10000)
[1] 0.9989306
var(obs_1000)
[1] 0.01282296
var(obs_10000)
[1] 0.00123895

Я думаю, вы хотите запустить симуляцию Монте-Карло с разными значениями. Я присылаю сюда код, который я использовал для своей курсовой работы по эконометрике, вы можете изменить его с помощью параметра, который вам нужен, очень легко. Извините, если я не буду вдаваться в подробности вашего кода, но я думаю, вы все еще пытаетесь понять, что вы хотите с ним делать, поэтому вам может быть полезно начать с базового кода, подобного моему.

library(MASS)
set.seed(123)

n_vec <- c(100, 500, 1000)
rho_vec <- c(0.1, 0.5, 0.9)
beta2_vec <- c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3)

results <- list()

for (n in n_vec) {
  for (rho in rho_vec) {
    for (beta2 in beta2_vec) {
      
      # set other parameters
      beta1 <- 0
      cov_mat <- matrix(c(1, rho, rho, 1), nrow = 2)
      k <- 1000
      
      # placeholder for current iteration results
      iter_result <- list()
      
      # run Monte Carlo simulations
      b1 <- numeric(k)
      b2 <- numeric(k)
      short_b1 <- numeric(k)
      test1 <- numeric(k)
      test2 <- numeric(k)
      test3 <- numeric(k)
      for (i in 1:k) {
        data <- as.data.frame(mvrnorm(n, mu = c(0, 0), Sigma = cov_mat))
        names(data) <- c("X1", "X2")
        data$Y <- beta1 * data$X1 + beta2 * data$X2 + rnorm(n, mean = 0, sd = 1)
        
        # long model
        reg <- lm(Y ~ X1 + X2, data = data)
        b1[i] <- coef(reg)["X1"]
        b2[i] <- coef(reg)["X2"]
        
        # short model
        short_reg <- lm(Y ~ X1, data = data)
        short_b1[i] <- coef(short_reg)["X1"]
        
        # compute t-tests
        se_b1 <- summary(reg)$coefficients["X1", "Std. Error"]
        se_b2 <- summary(reg)$coefficients["X2", "Std. Error"]
        t1 <- sqrt(n) * (b1[i] / se_b1)
        t2 <- sqrt(n) * (b2[i] / se_b2)
        t3 <- sqrt(n) * (short_b1[i] / summary(short_reg)$coefficients["X1", "Std. Error"])
        
        test1[i] <- as.integer(t1 < sqrt(log(n)))
        test2[i] <- as.integer(t2 < sqrt(log(n)))
        test3[i] <- as.integer(t3 < sqrt(log(n)))
      }
      
      # store results for current iteration
      iter_result$n <- n
      iter_result$rho <- rho
      iter_result$beta1 <- beta1
      iter_result$beta2 <- beta2
      iter_result$b1_summary <- summary(b1)
      iter_result$b2_summary <- summary(b2)
      iter_result$short_b1_summary <- summary(short_b1)
      iter_result$test1_summary <- summary(test1)
      iter_result$test2_summary <- summary(test2)
      iter_result$test3_summary <- summary(test3)
      
      # append current iteration results to overall results list
      results <- c(results, iter_result)
      
    }
  }
}

sink("output.txt")

cat("My list:\n")
print(results)

sink()