y=c(-1,0,1,1,0,-1)
 one=rep(1,6)
 x=c(1,1,1,-1,-1,-1)
 X=cbind(one,x)
 XpX=t(X)%*%X
 XpXi=solve(XpX)
 beta=XpXi%*%t(X)%*%y
 yhat=X%*%beta
 SSE=sum((y-yhat)^2)
 df=(6-2)
 MSE=SSE/(6-2)
 t=qt(1-(.05/2)/2,df)
CI=0 + c(-1,1)*t*sqrt(MSE)*sqrt(1/6)
plot(x,y,type="n",xlim=c(-2,2),ylim=c(-2,2),main="square - simultaneous intervals, circle - simultaneous region",
 xlab="beta 0",ylab="beta 1")
lines(c(CI[1],CI[1]),c(CI[1],CI[2]))
lines(c(CI[2],CI[1]),c(CI[2],CI[2]))
lines(c(CI[2],CI[1]),c(CI[1],CI[1]))
lines(c(CI[2],CI[2]),c(CI[1],CI[2]))
Area1=(CI[2]-CI[1])^2
#simultanous confidence region is a circle centered at 0, with radius squared
#equal to (1/6)(2)(MSE)(qf(.95,2,df))
r2=(1/6)*2*MSE*qf(.95,2,df)
Area2=pi*r2
xl=-sqrt(r2); xh=abs(xl)
xplot=seq(xl,xh,length.out=1000)
yh=sqrt(r2-xplot^2)
yl=-yh
lines(xplot,yh); lines(xplot,yl)


lmout=lm(y~x)
summary(lmout)