#R-function to do ordinary regression and get the standard errors out
# x is a matrix of regressors 
Myols=function(y,x)
{
# create a column of ones for the intercept

n=length(y)
ane=rep(1,n)
bigx=cbind(ane,x)
xtx=t(bigx)%*%bigx
xty=t(bigx)%*%y
bb=solve(xtx)%*%xty  #regr coefficients from (x’x)inverse x’y
r1=y-bigx%*%bb  # GET residuals
p=dim(bigx)[2]  #needed for degrees of freedom

s2=(t(r1)%*%r1)/(n-p)#variance or residuals
xtxinv=solve(xtx)
se=sqrt(s2*diag(xtxinv))
list(bb=bb, se=se, r1=r1)
}

#example
y=c(2,6,8,10)
x=c(11,3,31,45)
z=c(7,8,5,55)
xx=cbind(x,z)
reg2=Myols(y,xx)
summary(reg2)#this is useless
reg2$bb # regr coefficients
reg2$se #regr standard errors
reg2$r1  #residuals
reg1=lm(y~x+z)
summary(reg1)



#the following is for computing ridge regression
#You must standardize the data before using ridge
# if n=length(x)
# standardized x is xs=(x-mean(x))/(sqrt(n-1)*sd(x))
# k is usually a small number to be given by the user
ridg=function(y,bigx,k)
{
xtx=t(bigx)%*%bigx
xty=t(bigx)%*%y
n=dim(bigx)[2]
eye=diag(rep(1,n))
ki=k*eye
bb=solve(xtx+ki)%*%xty
r1=y-bigx%*%bb
s2=(t(r1)%*%r1)/2
xtxinv=solve(xtx+ki)
se=sqrt(s2*diag(xtxinv))
list(bb=bb, se=se)
}
#example
y=c(2,6,8,10)
x=c(11,3,31,45)
z=c(7,8,5,55)
n=length(y)
cor(cbind(x,z)) #correlation matrix for regressors
xs=(x-mean(x))/(sqrt(n-1)*sd(x))
zs=(z-mean(z))/(sqrt(n-1)*sd(z))
bigx=cbind(xs,zs)
t(bigx)%*%bigx  #this should be correlation matrix
#if it is not the correlation matrix, standardization is wrong
reg2=ridge(y,bigx,0.2)
reg2
reg2$bb
reg2$se