# Monte Carlo - Confidence Intervalls # 13.04.2014 rm(list=ls(all=TRUE)) n <- 50 # sample size Reps <- 1000 # no of replications alpha <- 0.05 # Significance level sigma <- 100 # standard deviation of eps x <- seq(n) # Trend CI <- matrix(rep(NA, 3 * Reps), ncol = 3) inCI <- 0 for (r in 1:Reps) { # true model eps <- rnorm(n, mean = 0, sd = sigma) y <- 5 + 5*x + eps # ß1 = 5, ß2 = 5 # estimation eq1 <- lm(y ~ x) b2 <- summary(eq1)$coefficients[2,1] # coef. se.b2 <- summary(eq1)$coefficients[2,2] # std. error coef. t.crit <- qt((1-alpha/2),summary(eq1)$df[2]) # critical value of t-distr. CI[r,1] <- b2 CI[r,2] <- b2 - t.crit * se.b2 # Conf. Int. lower CI[r,3] <- b2 + t.crit * se.b2 # Conf. Int. upper if (CI[r,2] <= 5 & 5 <= CI[r,3]) { inCI <- inCI + 1 } } # end loop Reps # Show result head(CI) tail(CI) cat("True value is in", inCI/Reps*100, "percent of replications in CI")