Chương 7 Kiểm định thống kê và Khoảng tin cậy với Hồi quy đơn giản

Trong chương này ta tiếp tục tìm hiểu về mô hình hồi quy đơn giản. Cụ thể, ta sẽ biết cách vận dụng kiến thức về phân phối lấy mẫu của OLS estimator để đưa ra nhận định về sự không chắc chắn.

Trong R ta dùng các gói sau: AER (Kleiber & Zeileis, 2018)3scales (Wickham, 2018)4.

7.1 Kiểm định hai đuôi liên quan đến Slope

Ở chương trước ta biết rằng hệ số \(\hat{\beta}_1\) tuân theo quy luật xấp xỉ phân phối Normal đối với quy mô mẫu đủ lớn, kiểm định thống kê về giá trị đúng của \(\beta_1\) có thể thực hiện như chương . Theo đó thống kê \(t\):

\[t = \frac{\text{Giá trị ước lượng}- \text{Giá trị giả thuyết}}{\text{Sai số chuẩn của estimator}}\]

trong đó, sai số chuẩn (\(SE\)) của \(\hat{\beta}_1\) là:

\[SE(\hat{\beta}_1) = \sqrt{\hat{\sigma}_{\hat{\beta}_1}} = \sqrt{\frac{1}{n} \frac{\frac{1}{n-2}\sum_{i=1}^n (X_i - \bar{X})^2\hat{u}_i^2}{\left[\frac{1}{n}\sum_{i=1}^n (X_i-\bar{X})^2\right]^2}}\]

Giả sử mô hình hồi quy của chúng ta như sau.

\[\hat{\text{Score}} = 698.9 -2.28\times \hat{STR}\]

Kết quả hồi quy được thể hiện trong bảng sau.

Fitting linear model: score ~ STR
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 698.9 9.467 73.82 6.57e-242
STR -2.28 0.4798 -4.751 2.783e-06

Theo đó, kết quả thống kê \(t\) cho hệ số hồi quy của biến \(STR\) bằng \(-4.751\). Ta hoàn toàn có thể tính được giá tri này nếu giá trị kiểm định của hệ số này là \(0\) như sau:

\[t = \frac{-2.28 - 0}{0.4798} = -4.751\]

Chú ý rằng R cung cấp cho ta giá trị p-value thông qua cột Pr(>|t|) nhưng giá trị này không được lấy từ phân phối cận Normal mà lấy từ phân phối Student, với bậc tự do \(DF = n-k-1\) với \(n\) là số quan sát, \(k\) là số biến \(X\). Giá trị p-value đối với phân phối cận Normal nên là -2.020860210^{-6} thay vì 2.783e-06 như kết quả hồi quy theo R. Ta thấy rằng sự khác biệt không quá lớn, bởi lẽ \(n\) đủ lớn.

Kết quả kiểm định cho thấy, nếu \(H_0: \beta_1= 0\) là đúng và ta khi lặp đi lặp lại quá trình thu thập quan sát và ước lượng mô hình, thì khả năng thấy \(\hat{\beta}_1 \ge |-2.28|\) là rất thấp. Để hình dung, ta có thể nhìn vào đồ thị phân phối dưới đây.

Giá trị p-value là diện tích nằm dưới đường cong bên phải \(4.75\) và bên trái \(-4.75\), và như ta đã tính toán, nó rất nhỏ.

7.2 Khoảng tin cậy

Ta không bao giờ ước lượng chính xác giá trị đúng của các tham số hồi quy từ dữ liệu mẫu, tuy nhiên ta hoàn toàn có thể xây dựng khoảng tin cậy cho chúng. Khoảng tin cậy \(95\%\) cho hệ số \(\beta_1\) được tính như sau:

\[\text{CI}_{95\%}^{\beta_1} = \left[\hat{\beta}_1-1.96\times SE(\hat{\beta}_1),\hat{\beta}_1+1.96\times SE(\hat{\beta}_1)\right]\]

Điều này tương đương với cách nói, đây là khoảng giá trị mà bộ giả thuyết \(H_0\) cho kiểm định hai đuôi \(5\%\) sẽ không bị bác bỏ.

Để hiểu rõ hơn về khoảng tin cậy, ta thực hiện giả lập thí nghiệm như sau.

7.2.1 Giả lập thí nghiệm

Giả sử có mẫu \(n=100\) quan sát gồm biến \(Y\) nào đó với \(Y \stackrel{iid}{\sim} N(5,25)\). Khảo sát đồ thị điểm rải như sau.

# set seed for reproducibility
set.seed(4)

# generate and plot the sample data
Y <- rnorm(n = 100, 
           mean = 5, 
           sd = 5)

plot(Y, 
     pch = 19, 
     col = "steelblue")

Giả sử dữ liệu được tạo bởi mô hình:

\[Y_i = \mu + \varepsilon_i\]

trong đó \(\mu\) là một hằng số chưa biết và ta biết rằng \(\varepsilon \stackrel{iid}{\sim} N(0,25)\). Trong mô hình này, OLS estimator cho \(\mu\) như sau:

\[\mu = \bar{Y} = \frac{1}{n}\sum_{i=1}^n Y_i\]

Khi đó, theo chương :

\[SE(\hat{\mu}) = \frac{\sigma_{\varepsilon}}{\sqrt{n}} = \frac{5}{\sqrt{100}} = 0.5\]

Như vậy, khoảng tin cậy cho \(\mu\) là:

\[\text{CI}_{95\%}^{\mu} = \left[\hat{\mu}-1.96\times 0.5, \hat{\mu}+1.96\times 0.5\right]\]

Trong R ta có thể tính giá trị này như sau.

cbind(CIlower = mean(Y) - 1.96 * 5 / 10, CIupper = mean(Y) + 1.96 * 5 / 10)
##       CIlower  CIupper
## [1,] 4.502625 6.462625

Ta biết rằng \(\mu=5\) nên rõ ràng khoảng tin cậy ở trên có bao gồm giá trị đúng của \(\mu\).

Sâu hơn, ta có thể dùng phương pháp lấy mẫu thay thế để tạo lập các khoảng tin cậy nhằm xem xét liệu CI estimator này có tốt hay không. Tiến trình thực hiện như sau:

  • lấy mẫu \(100\) quan sát, lặp lại \(10000\) lần.
  • giả lập \(10000\) khoảng tin cậy.
# set seed
set.seed(1)

# initialize vectors of lower and upper interval boundaries
lower <- numeric(10000)
upper <- numeric(10000)

# loop sampling / estimation / CI
for(i in 1:10000) {
  
  Y <- rnorm(100, mean = 5, sd = 5)
  lower[i] <- mean(Y) - 1.96 * 5 / 10
  upper[i] <- mean(Y) + 1.96 * 5 / 10
  
}

# join vectors of interval bounds in a matrix
CIs <- cbind(lower, upper)

Ta kỳ vọng xác suất của \(10000\) khoảng tin cậy giả lập có chứa giá trị đúng \(\mu=5\) khoảng \(95\%\).

mean(CIs[, 1] <= 5 & 5 <= CIs[, 2])
## [1] 0.9487

Ta thấy rằng 0.9487 rất gần với giá trị ta kỳ vọng. Hình hoạ hoá \(100\) khoảng tin cậy đầu tiên sẽ cho ta cái nhìn tốt hơn.

# identify intervals not covering mu
# (4 intervals out of 100)
ID <- which(!(CIs[1:100, 1] <= 5 & 5 <= CIs[1:100, 2]))

# initialize the plot
plot(0, 
     xlim = c(3, 7), 
     ylim = c(1, 100), 
     ylab = "Sample", 
     xlab = expression(mu), 
     main = "Confidence Intervals")

# set up color vector
colors <- rep(gray(0.6), 100)
colors[ID] <- "red"

# draw reference line at mu=5
abline(v = 5, lty = 2)

# add horizontal bars representing the CIs
for(j in 1:100) {
  
  lines(c(CIs[j, 1], CIs[j, 2]), 
        c(j, j), 
        col = colors[j], 
        lwd = 2)
  
}

Ta thấy rằng có \(4\) khoảng tin cậy không chứa giá trị \(\mu=5\), được tô màu đỏ, trong các trường hợp này, \(H_0\) sẽ bị bác bỏ.

Quay trở lại với mô hình hồi quy đơn giản, ta có thể tính toán khoảng tin cậy cho các hệ số hồi quy như sau.

# compute 95% confidence interval for coefficients in 'linear_model'
confint(linear_model)
##                 2.5 %     97.5 %
## (Intercept) 680.32312 717.542775
## STR          -3.22298  -1.336636

Ta có thể kiểm tra lại kết quả từ R và công thức chúng ta tự tính toán:

# compute 95% confidence interval for coefficients in 'linear_model' by hand
lm_summ <- summary(linear_model)

c("lower" = lm_summ$coef[2,1] - qt(0.975, df = lm_summ$df[2]) * lm_summ$coef[2, 2],
  "upper" = lm_summ$coef[2,1] + qt(0.975, df = lm_summ$df[2]) * lm_summ$coef[2, 2])
##     lower     upper 
## -3.222980 -1.336636

7.3 Hồi quy khi \(X\) là biến nhị phân

Thay vì biến \(X\) liên tục, trong một vài trường hợp ta quan tâm đến mô hình hồi quy sau:

\[Y_i = \beta_0 + \beta_1 D_i + u_i\]

trong đó \(D_i\) là biến nhị phân, hay còn gọi là dummy variable. Chẳng hạn, \(D_i\) được định nghĩa:

\[D_i =\begin{cases} 1 \text{ nếu } STR <20 \\ 0 \text{ nếu } STR \ge 20\end{cases}\]

Khi đó đồ thị điểm rải như sau.

# Create the dummy variable as defined above
CASchools$D <- CASchools$STR < 20

# Plot the data
plot(CASchools$D, CASchools$score,            # provide the data to be plotted
     pch = 20,                                # use filled circles as plot symbols
     cex = 0.5,                               # set size of plot symbols to 0.5
     col = "Steelblue",                       # set the symbols' color to "Steelblue"
     xlab = expression(D[i]),                 # Set title and axis names
     ylab = "Test Score",
     main = "Dummy Regression")

Rõ ràng ta không thể nào xây được đường tuyến tính như ở các ví dụ trước, thay vào đó, kết quả hồi quy sẽ được giải thích như sau:

  • \(E(Y_i|D_i=0) = \beta_0\), do đó \(\beta_0\) là điểm kỳ vọng khi \(D_i=0\) hay \(STR\) nhỏ hơn 20.
  • \(E(Y_i|D_i=0) = \beta_0 + \beta_1\), do đó \(\beta_1\) là sự chênh lệch của điểm thi kỳ vọng giữa những lớp học \(STR <20\)\(STR \ge 20\).

Để ước lượng, ta dùng R như sau.

# estimate the dummy regression model
dummy_model <- lm(score ~ D, data = CASchools)
summary(dummy_model)
## 
## Call:
## lm(formula = score ~ D, data = CASchools)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -50.496 -14.029  -0.346  12.884  49.504 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  650.077      1.393 466.666  < 2e-16 ***
## DTRUE          7.169      1.847   3.882  0.00012 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.74 on 418 degrees of freedom
## Multiple R-squared:  0.0348, Adjusted R-squared:  0.0325 
## F-statistic: 15.07 on 1 and 418 DF,  p-value: 0.0001202

Ta có thể thấy các giá trị kỳ vọng của điểm thi ở cả hai trường hợp trên đồ thị điểm rải như sau.

plot(CASchools$D, CASchools$score,            
     pch = 20,                                
     cex = 0.5,                               
     col = "Steelblue",                       
     xlab = expression(D[i]),                 
     ylab = "Test Score",
     main = "Dummy Regression")
# add group specific predictions to the plot
points(x = CASchools$D, 
       y = predict(dummy_model), 
       col = "red", 
       pch = 20)

7.4 Phương sai đồng nhất và Phương sai thay đổi

Mọi suy luận thống kê ở các phần trước đều dựa trên một giả định quan trọng: phương sai sai số không thay đổi khi giá trị của \(X\) thay đổi. Tuy nhiên các vấn đề trong thực tế đều không phải lúc nào cũng thoả mãn giả định này.

Phương sai của \(u_i\) được cho là đồng nhất( homoskedastic ) nếu phương sai của phân phối điều kiện của \(u_i\) cho trước \(X_i\) là hằng số:

\[\text{Var}(u_i|X_i = x) = \sigma^2 \text{ } \forall i=1,\dots,n\]

Nếu như có sự phụ thuộc của phương sai phân phối này, ta nói phương sai sai số thay đổi ( heteroskedasticity).

Để hình dung hiện tượng này, ta khảo sát hai biến như sau.

# generate some heteroskedastic data:

# set seed for reproducibility
set.seed(123) 

# set up vector of x coordinates
x <- rep(c(10, 15, 20, 25), each = 25)

# initialize vector of errors
e <- c()

# sample 100 errors such that the variance increases with x
e[1:25] <- rnorm(25, sd = 10)
e[26:50] <- rnorm(25, sd = 15)
e[51:75] <- rnorm(25, sd = 20)
e[76:100] <- rnorm(25, sd = 25)

# set up y
y <- 720 - 3.3 * x + e

# Estimate the model 
mod <- lm(y ~ x)

# Plot the data
plot(x = x, 
     y = y, 
     main = "Heteroskedasticity",
     xlab = "Student-Teacher Ratio",
     ylab = "Điểm thi",
     cex = 0.5, 
     pch = 19, 
     xlim = c(8, 27), 
     ylim = c(600, 710))

# Add the regression line to the plot
abline(mod, col = "darkred")

# Add boxplots to the plot
boxplot(formula = y ~ x, 
        add = TRUE, 
        at = c(10, 15, 20, 25), 
        col = alpha("gray", 0.4), 
        border = "black")

Ta đã dùng y~x trong hàm boxplot() để chỉ rõ rằng vector y được chia thành các nhóm dựa vào x. Nhìn vào đồ thị hình hộp, ta có thể thấy phương sai của Điểm thi (và phương sai của error term) tăng cùng với \(STR\).

7.4.1 Ví dụ thực tế

Ta xem xét về giá trị kinh tế của giáo dục. Nếu giáo dục đại học không tạo được giá trị kinh tế tăng thêm, rõ ràng sẽ chẳng có sinh viên nào thèm đi học. Do đó ta liên tưởng đến mối quan hệ giữa tiền lương và trình độ giáo dục. Ta kỳ vọng rằng, trình độ càng cao, tiền lương càng lớn, có nghĩa một cách nào đó, ta khảo sát quan hệ tuyến tính giữa hai biến theo mô hình sau:

\[\text{Tiền lương}_i = \beta_0 + \beta_1 \times \text{Trình độ}_i + u_i\]

Rõ ràng, ta kỳ vọng \(\beta_1 >0\), đồng thời, thu nhập của những người có giáo dục cao hơn sẽ có sự dàn trải rộng hơn so với những ngiowif có trình độ thấp. Cụ thể, một công nhân dệt may tốt nghiệp cấp 3 mức lương vào khoảng 3 triệu đến 5 triệu một tháng, tuy nhiên một cử nhân vừa tốt nghiệp đại học ngành dệt may có thể có mức lương dao động trong khoảng 8 triệu đến 16 triệu một tháng. Đơn giản là vì giáo dục cao chưa chắc làm việc tốt hoặc họ nhận thức được “sứ mệnh” của họ ở một số lĩnh vực khác ngành đào tạo, do đó họ chấp nhận các mức lương thấp. Có vô số cách giải thích cho hiện tượng: giáo dục cao nhưng mức lương thấp. Tuy nhiên, với một công việc mức lương cao, điều kiện cần là giáo dục phải cao tương đối. Nói cách khác, những người có bằng cấp cao sẽ có nhiều cơ hội tiếp cận công việc mức lương cao hơn so với những người có bằng cấp thấp hơn.

Để chứng thực nhận định thực tế này, ta có thể dùng dữ liệu của Khảo sát dân số của Cục Thống kê Lao động của Mỹ ( Bureau of Labor Statistics) và khảo sát sự phụ thuộc của tiền lương theo giờ đối với số năm được đào tạo.

library(AER)
data("CPSSWEducation")
attach(CPSSWEducation)
## The following objects are masked from CPSSWEducation (pos = 5):
## 
##     age, earnings, education, gender
## The following objects are masked from CPSSWEducation (pos = 8):
## 
##     age, earnings, education, gender
## The following objects are masked from CPSSWEducation (pos = 11):
## 
##     age, earnings, education, gender
## The following objects are masked from CPSSWEducation (pos = 19):
## 
##     age, earnings, education, gender
summary(CPSSWEducation)
##       age          gender        earnings        education    
##  Min.   :29.0   female:1202   Min.   : 2.137   Min.   : 6.00  
##  1st Qu.:29.0   male  :1748   1st Qu.:10.577   1st Qu.:12.00  
##  Median :29.0                 Median :14.615   Median :13.00  
##  Mean   :29.5                 Mean   :16.743   Mean   :13.55  
##  3rd Qu.:30.0                 3rd Qu.:20.192   3rd Qu.:16.00  
##  Max.   :30.0                 Max.   :97.500   Max.   :18.00
# estimate a simple regression model
labor_model <- lm(earnings ~ education)

# plot observations and add the regression line
plot(education, 
     earnings, 
     ylim = c(0, 150))

abline(labor_model, 
       col = "steelblue", 
       lwd = 2)

Đồ thị cho thấy trung bình của phân phối thu nhập tăng cùng với trình độ đào tạo, xác nhận kỳ vọng hệ số \(\beta_1>0\) là có cơ sở. Khoảng tin cậy của hệ số này như sau.

# compute a 95% confidence interval for the coefficients in the model
confint(labor_model)
##                 2.5 %    97.5 %
## (Intercept) -5.015248 -1.253495
## education    1.330098  1.603753

Với kết quả trên, ta có thể bác bỏ giả thuyết hệ số \(\beta_1=0\) ở mức ý nghĩa \(5\%\).

Tuy nhiên nhìn vào đồ thị một lần nữa, ta thấy phương sai của tiền lương (mức độ phân tán của biến \(Y\)) tăng cùng với giá trị của \(X\), do đó rõ ràng có hiện tượng phương sai sai số thay đổi.

7.4.2 Hậu quả của Heteroskedasticity

Trong trường hợp phương sai sai số đồng nhất, ta biết ở chương trước rằng:

\[\sigma^2_{\hat{\beta_1}} = \frac{\sigma^2_u}{n\sigma^2_X}\]

Đây là một phiên bản đơn giản hoá vì ta bỏ qua bậc tự do của các biến. Công thức chính xác và đầy đủ được thể hiện ở . Phần mềm R ước lượng phương sai này bằng estimator như sau:

\[\tilde{\sigma}^2_{\hat{\beta_1}} = \frac{SER^2}{\sum_i (X_i-\bar{X})^2}\]

trong đó \(SER = \frac{1}{2}\sum_i \hat{u}_i^2\).

Tuy nhiên trong trường hợp phương sai sai số thay đổi, ước lượng này không còn consistent. Hậu quả là thống kê \(t\) sẽ không còn tuân theo phân phối cận Normal khi \(n\) đủ lớn. Do đó bất kỳ kiểm định, hay khoảng tin cậy nào được xây dựng đều không còn tin cậy.

Khi đó, ta cần dùng một estimator khác cho phương sai hệ số hồi quy:

\[\sigma^2_{\hat{\beta_1}} = \frac{1}{n} \frac{\frac{1}{n}\sum_i (X_i-\bar{X})^2\hat{u}^2_i}{\left[\frac{1}{n} \sum_i (X_i-\bar{X})^2\right]^2}\]

Đây là estimator theo Eicker-Huber-White standard errors5. Trong R ta có thể tính giá trị này bằng hàm vcovHC() trong gói sandwich với type = "HC0".

Trong thực tế, người ta hay sử dụng phiên bản estimator của STATA:

\[\sigma^2_{\hat{\beta_1},HC1} = \frac{1}{n} \frac{\frac{1}{n-2}\sum_i (X_i-\bar{X})^2\hat{u}^2_i}{\left[\frac{1}{n} \sum_i (X_i-\bar{X})^2\right]^2}\]

Điểm khác biệt duy nhất nằm ở tử thức. Trong R ta đặt type = "HC1" để tính dạng estimator này. Hàm vcovHC() cho ra một ma trận hiệp phương sai, do vậy, để tính phương sai, ta chỉ cần quan tâm tới diagonal line của ma trận này.

library(sandwich)
# compute heteroskedasticity-robust standard errors
vcov <- vcovHC(labor_model, type = "HC1")
vcov
##             (Intercept)   education
## (Intercept)  0.85726285 -0.06573674
## education   -0.06573674  0.00517569
# compute the square root of the diagonal elements in vcov
robust_se <- sqrt(diag(vcov))
robust_se
## (Intercept)   education 
##  0.92588490  0.07194227

Để tổng kết lại kết quả mô hình hồi quy giữa hai trường hợp phương sai đồng nhất và phương sai thay đổi, ta xem xét đoạn code dưới đây:

# Homoskedasticity
summary(labor_model)
## 
## Call:
## lm(formula = earnings ~ education)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.270  -5.355  -1.513   3.194  77.164 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.13437    0.95925  -3.268   0.0011 ** 
## education    1.46693    0.06978  21.021   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.769 on 2948 degrees of freedom
## Multiple R-squared:  0.1304, Adjusted R-squared:  0.1301 
## F-statistic: 441.9 on 1 and 2948 DF,  p-value: < 2.2e-16
# Hetoroskedasticity
library(lmtest)
# we invoke the function `coeftest()` on our model
coeftest(labor_model, vcov. = vcov)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value  Pr(>|t|)    
## (Intercept) -3.134371   0.925885 -3.3853 0.0007204 ***
## education    1.466925   0.071942 20.3903 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Để phát hiện phương sai sai số thay đổi, bên cạnh các kiểm định phức tạp, người ta thường nhìn vào đồ thị giữa error term và biến \(X\) hoặc \(\hat{Y}\). Nếu có xu hướng gì đó giữa các điểm rải này, thì rõ ràng có hiện tượng phương sai sai số thay đổi. Trong R, ta thực hiện như đoạn code sau và để ý đồ thị nằm góc trên bên trái và ở dưới bên trái.

CPSSWEducation$sigres <- labor_model$residuals
par(mfrow = c(2,2))
plot(labor_model)

7.4.3 Định lý Gauss-Markov

Khi ước lượng mô hình hồi quy, ta biết rằng kết quả sẽ ngẫu nhiên. Tuy nhiên, khi sử dụng các estimators không thiên lệch, ít ra về măt trung bình, ta đang ước lượng tham số đúng. Khi so sánh các estimators khác nhau đều cùng thiên lệch, ta cần một estimator có độ chính xác cao nhất, nghĩa là xác suất mà estimator càng gần giá trị đúng càng cao. Nghĩa là phương sai sai số nhỏ nhất trong các estimator đó.

Nói tóm lại, ta cần lựa chọn estimator nào phải consistent, cho trước sự không thiên lệch. Theo định lý Gauss-Markov, OLS estimator thoả mãn những điều kiện đã nêu, miễn là các giả định của nó được đáp ứng.

7.4.4 Chứng minh công thức tính phương sai của \(\beta_1\)

Estimator của \(\beta_1\):

\[\hat{\beta_1} = \frac{\sum_i (X_i - \bar{X})(Y_i - \bar{Y})}{\sum_i (X_i - \bar{X})^2}\]

Theo đó, tử thức được biến đổi như sau:

\[\begin{aligned} \sum_i (X_i - \bar{X})(Y_i - \bar{Y}) &= \sum_i (X_i - \bar{X})Y_i - \sum_i(X_i-\bar{X})\bar{Y} \\ &= \sum_i (X_i - \bar{X})Y_i - \bar{Y}\sum_i (X_i - \bar{X}) \\ &= \sum_i(X_i - \bar{X})Y_i \\ &= \sum_i(X_i - \bar{X})(\beta_0+\beta_1X_i+u_i) \\ &=\beta_0\sum_i(X_i-\bar{X}) +\beta_1\sum_i (X_i - \bar{X})X_i +\sum_i (X_i - \bar{X})u_i \\ &= \beta_0 \times 0 + \beta_1 \sum_i(X_i - \bar{X})^2 + \sum_i(X_i - \bar{X})u_i \\ &= \beta_1\sum_i(X_i - \bar{X})^2 + \sum_i(X_i - \bar{X})u_i \end{aligned}\]

Chú ý rằng: \(\sum_i(X_i - \bar{X})(Y_i-\bar{Y}) = \sum_i(X_i - \bar{X})Y_i\), nếu ta thay \(Y=X\) thì \(\sum_i(X_i - \bar{X})X_i = \sum_i(X_i - \bar{X})\times (X_i - \bar{X}) = \sum_i(X_i - \bar{X})^2\).

Khi đó:

\[\begin{aligned}\hat{\beta_1} &= \frac{\beta_1\sum_i(X_i - \bar{X})^2 + \sum_i(X_i - \bar{X})u_i}{\sum_i (X_i - \bar{X})^2} \\ &= \beta_1 + \frac{\sum_i (X_i - \bar{X})u_i}{\sum_i (X_i - \bar{X})^2} \end{aligned}\]

Như vậy:

\[\begin{aligned}\text{Var}(\hat{\beta}_1) &= \text{Var}\left(\frac{\sum_i (X_i - \bar{X})u_i}{\sum_i (X_i - \bar{X})^2}\right) \\ &= \frac{1}{\left(\sum_i(X_i-\bar{X})^2\right)^2} \text{Var}(\sum_i(X_i-\bar{X})u_i) \\ &= \frac{1}{\left(\sum_i(X_i-\bar{X})^2\right)^2} E\left[ \left(\sum_i(X_i-\bar{X})u_i\right)^2\right] \\ &= \frac{1}{\left(\sum_i(X_i-\bar{X})^2\right)^2} E\left[ \left(\sum_i(X_i-\bar{X})^2u_i^2 +2\sum_{i\ne j}(X_i-\bar{X})(X_j-\bar{X})u_iu_j\right)\right] \\ &= \frac{1}{\left(\sum_i(X_i-\bar{X})^2\right)^2}E\left[\sum_i(X_i -\bar{X})^2u_i^2\right] \\ &= \frac{1}{\left(\sum_i(X_i-\bar{X})^2\right)^2} \sum_i(X_i -\bar{X})^2 E(u_i^2) \\ &= \frac{\sigma^2}{\sum_i(X_i -\bar{X})^2} \end{aligned}\]

Trong đó: \(\sigma^2 = \text{Var}(u_i) = E(u_i^2) - (E(u_i))^2 = E(u_i^2)\) do \(E(u_i)=0\).


  1. Kleiber, C., & Zeileis, A. (2008). Applied Econometrics with R. Springer.

  2. Wickham, H. (2018). scales: Scale Functions for Visualization.

  3. https://en.wikipedia.org/wiki/Heteroscedasticity-consistent_standard_errors