y = kids
x = cbind(1, educ)
Estimate_Beta = solve(t(x)%*%x)%*%(t(x)%*%y) # the estimate of B0 and B1
y_hat = x%*%Estimate_Beta
Estimate_sigmasquare = sum ((y - y_hat)^2)/(length(y) - 2) # the estimate of variance
Var_Beta = Estima_sigmasquare*solve(t(x)%*%x) #the diagonal elements are the estimates of variances for Betas