{
#Suppose you have your data in X and Y
#Y being the vector dependent variables and X being the vector (or matrix)
#of independent variables.

Y=scan("UBN_precip.txt")
X=matrix(scan("UBN_predictors.txt"),ncol=7,byrow=T)

n=length(Y)	# of data points..
X1=cbind(rep(1,n),X)	#add the first column of 1s

betas=solve(t(X1) %*% X1)%*%t(X1)%*%Y	#betas = (X^TX)^-1 X^T Y

#command solve gives the inverse of a matrix
#t(X) gives the transpose of a matrix
#%*% is the command for matrix multiplication.

#you can check the betas obtained as above from lsfit..

zz=lsfit(X,Y)
#zz$coef  should be exactly equal to betas

#to get the variances of the betas

errvar=var(zz$resid)	#error variances..

varbetas = solve(t(X) %*% X)errvar #this gives a (k+1 X k+1) matrix of
				#variance-covariance..
#diagonal elements are the variances of the betas.

#to get the confidence intervals for beta
dofs = n-(k+1)-1 #degrees of freedom..
#k is the number of variables in the X matrix
#95% confidence interval
betas[i] + qt(0.025, dofs)*sqrt(varbetas[i,i])
betas[i] - qt(0.025, dofs)*sqrt(varbetas[i,i])

#You can also use this to get the t-statistic for the betas
#and test their significances

#confidence interval of the regression line

#for a given point i of the original data.. the 95% confidence intervals are

sum(X1[i,] * betas) + qt(0.025,dofs)*sqrt(t(X1[i,])%*%solve(t(X1)%*%X1)%*%X1[i,]*errvar)
sum(X1[i,] * betas) - qt(0.025,dofs)*sqrt(t(X1[i,])%*%solve(t(X1)%*%X1)%*%X1[i,]*errvar)


#For a new point xo it is a prediction interval..
sum(1[i,] * betas) + qt(0.025,dofs)*sqrt((1+t(X1[i,])%*%solve(t(X1)%*%X1)%*%X1[i,])*errvar)
sum(1[i,] * betas) - qt(0.025,dofs)*sqrt((1+t(X1[i,])%*%solve(t(X1)%*%X1)%*%X1[i,])*errvar)

# to get confidence and prediction intervals.. you can use lm
# lm and lsfit provide the same answers..
# but lm is easy to work with in terms of prediction etc..

zz=lm(Y ~ X)

predict.lm(zz,interval="confidence")

#this will give you the confidence interval at the original data points that went into
#fitting of the model.

# SAy your point of interest is in the vector xo
Xtrue=X #save the original data..

#replace X with xo
X=matrix(xo,ncol=k,byrow=T)

predict(zz,data.frame(X),interfval="confidence")

#predict.lm requires the data to be in the same name as the variable that went into
# the lm(..) command - here it is X

}