x2<-c(1,2,3,2)
x3<-c(2,1,1,2)
y0<-c(1,1,2,3)
xls<-matrix(c(1,2,3,2,2,1,1,2), nrow=4, ncol=2)
x<-matrix(c(1,1,1,1,1,2,3,2,2,1,1,2), nrow=4, ncol=3)
y<-matrix(c(1,1,2,3), nrow=4)
n<-length(y)

x
xpx=t(x)%*%x
xpx
xpxinv=solve(xpx)
xpxinv
xpy=t(x)%*%y
xpy
beta<-xpxinv%*%xpy
beta

ypy=t(y)%*%y
ypy
bpxpxb=t(beta)%*%t(x)%*%x%*%beta
bpxpxb
epe=ypy-bpxpxb
epe
sigusq=epe/(n-3)
sigusq

rsq<-(bpxpxb-n*mean(y)^2)/(ypy-n*mean(y)^2)
rsq

varcov<-0.25*xpxinv
varcov
se<-sqrt(diag(varcov))
se

x0<-matrix(c(1,3,2))
yhat<-t(x0)%*%beta
x0
yhat

sigisq<-0.25*(1+t(x0)%*%xpxinv%*%x0)
sigisq
sigi<-sqrt(sigisq)
sigi
yhat+c(-1,1)*qt(.975, 1)*sigi

sigesq<-0.25*(t(x0)%*%xpxinv%*%x0)
sigesq
sige<-sqrt(sigesq)
sige
yhat+c(-1,1)*qt(.975, 1)*sige