1 Их тооны хууль ба төвийн хязгаарын теорем
- 1.1 Их тооны хууль
- 1.2 Төвийн хязгаарын теорем
2 OLS t-тестийн асимптот тархалт
3 OLS \(t\)-тестийн чадал

1 Их тооны хууль ба төвийн хязгаарын теорем

1.1 Их тооны хууль

Их тооны хууль (Laws of Large Numbers-LLN) \(E(x)=\mu\) байх санамсаргүй хэмжигдэхүүнээс \(n\) хэмжээтэй түүвэр авахад \[\frac{1}{n}\sum_{i=1}^nx_i\to_p \mu\] байна. Өөрөөр хэлбэл түүврийн дундаж нь эх олонлогын дундаж руугаа магадлалт нийлнэ.

1.1.1 Шоо хаях

Шоо хаяхад хүлээгдэх утга нь 3.5 байна

set.seed(1)
n <- 1000
out <- sample(1:6, n, replace = T) 
dundaj <- cumsum(out)/1:n

## Warning: Removed 2 rows containing missing values (geom_path).

1.2 Төвийн хязгаарын теорем

Төвийн хязгаарын теорем (Central Limit Theorem-CLT) \(E(X)=\mu\), \(Var(Х)=\sigma^2\) байх санамсаргүй хэмжигдэхүүнээс \(n\) хэмжээтэй түүвэр авахад \[\sqrt{n}(\bar x -\mu)\to_d N(0, \sigma^2)\] биелнэ.

Дээрхи шоо хаях жишээн дээр CLT биелэхийг шалгая.

nsim <- 5000
n <- 100
out <- rep(NA, nsim)
for(i in 1:nsim){
  out[i] <- mean(sample(1:6, n, replace = TRUE))
}
res <- sqrt(n)*(out-3.5)/sqrt(2.916667)

ggplot2ашиглан хистограммыг зурвал

clr <- RColorBrewer::brewer.pal(8, "Set3")
df <- data.frame(x=res)
ggplot(df, aes(x=x)) + 
  geom_histogram(aes(y = ..density..), binwidth=0.4,  colour="black", fill=clr[5]) +
  stat_function(fun = dnorm, colour="red", size=1.2, alpha=0.7)

2 OLS t-тестийн асимптот тархалт

\[\boldsymbol y=\beta_0+\beta_1 \boldsymbol x +\boldsymbol u\] загварын хувьд \(H:\beta_1=0\) таамаглалын \(t\)-статистикийн асимптот тархалтыг олъё.

\((\beta_0, \beta_1)\) коэффициентийн OLSE \[\left(\begin{array}{c} \hat \beta_0\\ \hat\beta_1 \\ \end{array}\right)=(X'X)^{-1}X'y\] энд \(X=\left(\begin{array}{c c} \boldsymbol 1& \boldsymbol x\\ \end{array}\right)\)

Иймд дээрх таамаглалын \(t\)-статистик нь \[t_{\hat{\beta_1}}=\frac{\hat\beta_1}{s_{\hat\beta_1}}\to_d N(0, 1)\] \(s_{\hat\beta_1}^2=\hat\Omega_{(2,2)}\)

\(s^2=\frac{\sum_{i=1}^n(y_i-\hat\beta_0-\hat\beta_1 x_i)^2}{n-1}\)
\(\hat\Omega=s^2(X'X)^{-1}\)

Бид lm функц ашиглан олж болох боловч симуляци хийхэд хурдан биш учир учир дээрх томъёог ашиглан \(t\)-утгуудыг бодно. Түүврийн хэмжээ n=1000-ийн хувьд nsim=5000 утгыг симуляци хийн хистограммыг зурж үзүүлбэл том түүвэрт \(t\)-утгууд стандарт нормаль тархалттай байдгийг ерөнхийдөө харуулж байна.

nsim <- 5000
n <- 1000
tstat <- rep(NA, nsim)
for(i in 1:nsim){
  y <- 1+rnorm(1000)
  x <- 5*rnorm(1000)
  X <- cbind(1, x)
  b <- solve(crossprod(X), crossprod(X, y)) # better, faster expression for solve(t(X)%*%X)%*%t(X)%*%y
  s2 <- sum((y - X%*%b)^2)/(n-2)
  Omega <- s2*solve(crossprod(X))
  tstat[i] <- b[2]/sqrt(Omega[2, 2])
}

Base R-ийн график функц ашиглан хистограммыг үзүүлбэл \(N(0, 1)\) тархалтад ойр буйг шалгаж болно

hist(tstat, density=20, breaks=20, prob=TRUE, 
         xlab="x", xlim = c(-4, 4), ylim = c(0, 0.5),
         main="Histogram of Simulated t-values")
curve(dnorm(x), 
          col="darkblue", lwd=2, add=TRUE)

\(\alpha=(0.01, 0.025, 0.05, 0.9, 0.95, 0.975, 0.99)\) цэгүүдийн хувьд квантилуудыг харьцуулбал

alpha <- c(0.01, 0.025, 0.05, 0.9, 0.95, 0.975, 0.99)
quantile(tstat, alpha)

##        1%      2.5%        5%       90%       95%     97.5%       99% 
## -2.409773 -1.981000 -1.643014  1.259689  1.639269  1.973854  2.426397

qnorm(alpha)

## [1] -2.326348 -1.959964 -1.644854  1.281552  1.644854  1.959964  2.326348

3 OLS \(t\)-тестийн чадал

Тэг таамаглал \(H_0\) худлаа үед \(H_0\)-г няцаах магадлалыг тестийн чадал гэдэг. 2-талт \(t\)-тестээр хэрэв \(|t_{\hat{\beta}}|>C_{1-\alpha/2}\) бол \(\alpha\) түвшинд \(H_0\) няцаагдана.

\[\boldsymbol y=\beta_0+\beta_1 \boldsymbol x +\boldsymbol u\] загварын хувьд \(\beta=[-5, 5]\) завсарт \(t\)-тестийн чадлыг \(\alpha=0.05\) түвшинд тооцъё (түүврийн хэмжээ \(n=20\)) . \(t\) тестийн хувьд онолын чадалыг мэдэх учир шууд харьцуулан, ойролцоолол хир сайн байгааг шалгаж болох юм.

Доорхи програмд tvals-д \(\beta_1\)-ийн утгад харгалзах \(t\)-утга, pow нь тухайн цэг дээрхи онолын чадалыг оноож байна. tpower функц нь nsim- оролттой, beta, n зэрэг функцаас гадуур тодорхойлогдсон утгыг global environment-с оруулан, beta-д харгалзах онолын болон симуляцийн утгыг data.frame байдалтай гаргана.

beta <- seq(-5, 5, 0.1) # -200:200/100
n <- 20
alpha <- 0.05
crit <- qt(1-alpha/2, df=n-2)  
x <- runif(n)
X <- cbind(1, x)

tpower <- function(nsim){
  #This function computes power of t test
  #in regression y=b0+b1*x+u 
  l <- length(beta)
  tval <- matrix(0, nsim, l)
  pow <- matrix(0, nsim, l)
  for (i in 1:nsim){
    u <- rnorm(n)
    for (j in 1:l){
      b = beta[j]
      y = 0.25 + x*b + u # true model
      betahat =  solve(crossprod(X), crossprod(X, y))
      residuals <- y - X %*% betahat
      s2 <- sum(residuals^2)/(n-2)
      S <- solve(crossprod(X))
      se <- sqrt(s2*S[2, 2])
      tval[i, j] <- betahat[2] / se
      pow[i, j] <- pt(-crit, df=n-2, ncp = b/se) + (1-pt(crit, df=n-2, ncp = b/se))
    }
  }        
  data.frame(x=beta,
             simulated =colMeans(abs(tval)>crit),
             theoretical = colMeans(pow))
}

power <- tpower(1000)

Зөвхөн 1000 симуляци хийсэн боловч онолын утгатайгаа ойрхон байгааг харж болно.