for(i in 1:1000){
b0 <- 5
b1 <- 0.2
b2 <- 0.2
x1 <- runif(50, 0, 10)
x2 <- runif(50, 0, 10)
e_var <- rnorm(50, 0, 3)

y <- b0 + x1 * b1 + x2 * b2 + e_var

mat[i, 1] <- summary(lm(y ~ x1 + x2))$coefficients[11] < 0.05
mat[i, 2] <- summary(lm(y ~ x1 + x2))$coefficients[12] < 0.05
}

mat$pw <-  mat[,1] * mat[,2]
sum(mat[,1]) / 1000
sum(mat[,2])  / 1000
sum(mat[,3]) / 1000