Симуляции в r (ols) - PullRequest
Новая
       76

Симуляции в r (ols)

0 голосов
/ 15 января 2020

Я пытаюсь запустить ols-симуляции в R, основанные на сценариях Carsey & Harden, с разными размерами выборок и уровнями корреляций (размер выборок: 26, 50, 200, 1000 и уровни корреляции .1, 0,2, 0,3, 0,4, 0,5, 0,6, 0,7, 0,8, 0,9). Затем я использую функцию покрытия (ту же, что указана в книге), чтобы вычислить вероятности покрытия, и она работает на выборках размером 26 и 50, но затем у 200, 1000 и более крупных она начинает возвращать только 0 или 0,001. Я проверил коэффициенты и в матрице par.est) Я вижу, что симуляция фактически не захватывает какое-либо значение внутри коэффициента интервала, но это странно, потому что 1000 - довольно хороший размер выборки и, по крайней мере, на некотором уровне корреляции, он должен вернуть некоторую вероятность больше, чем 0,001. Также это так странно, что я работаю в 26 и 50, но не больше этого. Что-нибудь может быть не так с моим кодом? Или это возможно, что R просто сходит с ума после выполнения некоторого количества симуляций?

Любой ответ будет очень признателен

это код для запуска моделей: 

par.est<-matrix(NA,nrow=reps,ncol=4) 
b1<-.4
b2<-.2
n<-26 #sample size
t1<-.1
C <- rnorm(n,0,1) 
V1<-C*t1^0.5 + rnorm(n, 0, 1)*(1-t1)^0.5 #Correlated variables
V2<-C*t1^0.5 + rnorm(n, 0, 1)*(1-t1)^0.5

for(i in 1:reps){ 
  Y<-V1*b1+V2*b2+rnorm(n,0,1) #The true DGP, with N(0,1) error
  model1<-lm(Y~V1+V2) 
  vcv1<-vcov(model1)
  par.est[i,1]<-model1$coef[1] 
  par.est[i,2]<-model1$coef[2]  
  par.est[i,3]<-sqrt(diag(vcv1)[1]) 
  par.est[i,4]<-sqrt(diag(vcv1)[2])

}


the coverage function:

coverage <- function(b, se, true, level = .95, df = Inf){ # Estimate, 
                                                          # standard error,
                                                          # true parameter, 
                                                          # confidence level, 
                                                          # and df  
qtile <- level + (1 - level)/2 # Compute the proper quantile
lower.bound <- b - qt(qtile, df = df)*se # Lower bound
upper.bound <- b + qt(qtile, df = df)*se # Upper bound 
# Is the true parameter in the confidence interval? (yes = 1)
true.in.ci <- ifelse(true >= lower.bound & true <= upper.bound, 1, 0)
cp <- mean(true.in.ci) # The coverage probability
mc.lower.bound <- cp - 1.96*sqrt((cp*(1 - cp))/length(b)) # Monte Carlo error  
mc.upper.bound <- cp + 1.96*sqrt((cp*(1 - cp))/length(b))  
return(list(coverage.probability = cp, # Return results
            true.in.ci = true.in.ci,
            ci = cbind(lower.bound, upper.bound),
            mc.eb = c(mc.lower.bound, mc.upper.bound)))
}

then I see the coverage:
cp.beta1_c1<-coverage(par.est[,1],par.est[,3],b1,df=n-model1$rank) 
cp.beta2_c1<-coverage(par.est[,2],par.est[,4],b2,df=n-model1$rank)
cp.beta1_c1$coverage.probability
cp.beta2_c1$coverage.probability
...