У меня возникли проблемы с моей пользовательской функцией. Основная идея заключается в том, что я хочу показать, что бета-оценка приближается к 1 по мере увеличения количества наблюдений. Поэтому я написал эту пользовательскую функцию и попробовал ее для n = 200 и n = 1000. Проблема в том, что изменение количества наблюдений ничего не дает, насколько я могу судить. Любой вклад высоко ценится!
Я предполагаю, что я плохо закодировал свою функцию (в отличие от некоторой статистической ошибки).
Вот функция:
func <- function(iterations, observations){
set.seed(1001)
betas <-c()
for (i in 1:iterations) {
c <- 1:observations
x <- c/5
e <- rnorm(1)*x # Remove x if you want homoskedastic errors
y <- x - 10 + e
#y[y < 0] <- 0 # censoring
model <- lm(y ~ 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 = 200)
func(iterations = 2000, observations = 1000)
mean(obs_200);mean(obs_1000) # Compare means
Данные графика
library(ggplot2)
# Make individual data frames
a <- data.frame(group = "a", value = obs_200)
b <- data.frame(group = "b", value = obs_1000)
# Combine into one long data frame
plot_data <- rbind(a, b)
# plot
ggplot(plot_data) +
geom_boxplot(aes(x = group,
y = value,
fill = group))
РЕДАКТИРОВАТЬ
library(tidyverse)
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)
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 = 200)
func(iterations = 2000, observations = 400)
func(iterations = 2000, observations = 600)
func(iterations = 2000, observations = 800)
func(iterations = 2000, observations = 1000)
mean(obs_200);mean(obs_1000) # compare means
# Make individual data frames
plot_data <- tibble(a = obs_200,
b = obs_400,
c = obs_600,
d = obs_800,
e = obs_1000) |>
pivot_longer(cols = everything(),
names_to = "n")
# plot
ggplot(plot_data) +
geom_boxplot(aes(x = n,
y = value,
fill = n)) +
scale_fill_grey(start = 0.5, labels = c("200", "400", "600", "800", "1000")) +
theme_bw() +
labs(x = NULL,
y = NULL,
fill = "Observations",
title = "Monte Carlo simulation using OLS",
subtitle = "(Heterosedastic errors)")
@Роланд. Я думаю, это все. Вернусь, когда буду дома...
Примечание: избегайте использования «c» в качестве имени переменной. Хотя это обычно не вызывает «путаницы» с функцией c
, могут произойти плохие вещи. Я бы тоже не стал использовать "e", так как когда-нибудь это может испортить что-то вроде 1.2e5
Несколько вопросов:
e <- rnorm(1)*x
означает, что все ошибки коррелированы. Вместо этого используйте e <- rnorm(observations)*x
.Код:
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
Спасибо. Я обновил свой оригинальный пост. Я не добавлял веса, так как хотел проверить непротиворечивость оценки без дополнений (но все равно спасибо!). Кажется, что переход с 200 на 1000 делает совсем немного! Я отмечу ваш ответ как решение, если вы не хотите что-то добавить к моему редактированию?
Выглядит неплохо! Мне нравится созданный вами сюжет.
Я думаю, вы хотите запустить симуляцию Монте-Карло с разными значениями. Я присылаю сюда код, который я использовал для своей курсовой работы по эконометрике, вы можете изменить его с помощью параметра, который вам нужен, очень легко. Извините, если я не буду вдаваться в подробности вашего кода, но я думаю, вы все еще пытаетесь понять, что вы хотите с ним делать, поэтому вам может быть полезно начать с базового кода, подобного моему.
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()
Я предполагаю, что вы хотите
e <- rnorm(observations) * x
. (Кроме того, избегайте побочных эффектов, т. е.assign
. Пусть ваша функция возвращает результат и выполняетobs_200 <- func(iterations = 2000, observations = 200)
.)