使用 R 进行蒙特卡洛模拟如何为 R“翻译”Stata 代码

问题描述

我需要从 Stata 到 R 运行相同的练习（蒙特卡罗模拟）。

我在 Stata 中使用的代码是下面的代码。我如何使用 R 做到这一点？（我已经搜索了很多教程，但我仍然没有设法在 R 中做到这一点）。

* Simulations (10,100 and 1000 sample replications/iterations)

clear
drop _all
set obs 100
set seed 54231
gen x = ((rnormal()))*10 + 40

* Generating true_y,considering Beta = 0,035
 
gen true_y = 5+0.03500*x
save truth,replace
twoway scatter true_y x

program hetero1
version 13
args c
use truth,clear
gen y = true_y + rnormal()
regress y x
end


foreach i in 10 100 1000 {
simulate _b,reps (`i'): hetero1
sum _b_x
twoway histogram _b_x,fraction xline(+0.03500,lcolor(red)) xline(`r(mean)',lcolor(green)) fcolor(none) lcolor(gs0) legend(off) title(`i' Repetições)
graph save graf`i'.gph,replace
}
 gr combine graf10.gph graf100.gph graf1000.gph
 graph export "graf.png",as(png) replace

最后，考虑到 10、100 和 1000 个样本重复，我需要获得这 3 个直方图（估计的 beta/系数）。红线指的是“真实”系数，绿线是估计系数的平均值 - [见链接中的图片]

解决方法

应该这样做：

# Simulations (10,100 and 1000 sample replications/iterations)
library(ggplot2)
library(dplyr)
library(gridExtra)
n <- 100
set.seed(54231)
x <- rnorm(n)*10 + 40

# Generating true_y,considering Beta = 0,035

true_y <- 5+0.03500*x

plot(x,true_y)

b <- t(replicate(1110,coef(lm(true_y + rnorm(n) ~ x))))

b <- as.data.frame(b) %>% 
  rename("a" = "(Intercept)","b" = "x") %>% 
  mutate(
    obs = 1:n(),n = case_when(
      obs %in% 1:10 ~ "N = 10",obs %in% 11:110 ~ "N = 100",TRUE ~ "N = 1000"),n = factor(n,levels=c("N = 10","N = 100","N = 1000")))
  

b10 <- b %>% filter(n == "N = 10") 
g1 <-   ggplot() + 
  geom_histogram(data = b10,aes(x=b),bins=3,col="white") + 
  geom_vline(xintercept = 0.03500,col="red") + 
  geom_vline(data = b10 %>% summarise(b=mean(b)),aes(xintercept = b),col="green") + 
  facet_wrap(~n) + 
  theme_bw()

b100 <- b %>% filter(n == "N = 100") 
g2 <-   ggplot() + 
  geom_histogram(data = b100,bins=10,col="red") + 
  geom_vline(data = b100 %>% summarise(b=mean(b)),col="green") + 
  facet_wrap(~n) + 
  theme_bw()

b1000 <- b %>% filter(n == "N = 1000") 
g3 <-   ggplot() + 
  geom_histogram(data = b1000,bins=25,col="red") + 
  geom_vline(data = b1000 %>% summarise(b=mean(b)),col="green") + 
  facet_wrap(~n) + 
  theme_bw()

library(gridExtra)
grid.arrange(g1,g2,g3,nrow=2)

montecarlo r r stata