library(MASS) 
library(car) 
set.seed(1000)

#High correlation
mu <- c(0,0)
Sigma <- matrix(c(1,0.9999,0.9999,1),2,2)
x <- mvrnorm(n=10,mu,Sigma)
x1 <- x[,1]
x2 <- x[,2]
x3 <- rnorm(10,0,1)
X <- cbind(x1,x2,x3)
cor(X)

XTX<-t(X) %*% X
evec1<-eigen(XTX)$values
max(evec1)/min(evec1)

y <- x1+x2+x3+rnorm(10,0,0.1)
 t1 <- lm(y~x1+x2+x3)
summary(t1)

#Variance Inflation factor
vif(t1)

#Independent case
Sigma <- matrix(c(1,0,0,1),2,2)
x <- mvrnorm(n=10,mu,Sigma)
x1 <- x[,1]
x2 <- x[,2]
x3 <- rnorm(10,0,1)
X <- cbind(x1,x2,x3)
cor(X)

XTX<-t(X) %*% X
evec2<-eigen(XTX)$values
max(evec2)/min(evec2)


y <- x1+x2+x3+rnorm(10,0,0.1)
 t2 <- lm(y~x1+x2+x3)
summary(t2)
vif(t2)

#Medium correlation
#High correlation
mu <- c(0,0)
Sigma <- matrix(c(1,0.8,0.8,1),2,2)
x <- mvrnorm(n=10,mu,Sigma)
x1 <- x[,1]
x2 <- x[,2]
x3 <- rnorm(10,0,1)
X <- cbind(x1,x2,x3)
cor(X)

XTX<-t(X) %*% X
evec3<-eigen(XTX)$values
max(evec3)/min(evec3)

y <- x1+x2+x3+rnorm(10,0,0.1)
 t3 <- lm(y~x1+x2+x3)
summary(t3)

#####################################################
#Ridge

#High correlation
set.seed(3274)
mu <- c(0,0)
Sigma <- matrix(c(1,0.999,0.999,1),2,2)
x <- mvrnorm(n=10,mu,Sigma)
x1 <- x[,1]
x2 <- x[,2]
x3 <- rnorm(10,0,1)
X <- cbind(x1,x2,x3)
y <- x1+x2+x3+rnorm(10,0,0.1)
(t.ridge1 <- lm.ridge(y~x1+x2+x3,lambda=0))

select(lm.ridge(y~x1+x2+x3,lambda=0:100/1000))

(t.ridge2 <- lm.ridge(y~x1+x2+x3,lambda=0.006))

####################################################
nreps<-1000
set.seed(3274)
mu <- c(0,0)
Sigma <- matrix(c(1,0.999,0.999,1),2,2)
ests.lm<-ests.ridge<-matrix(0,nrow=nreps,3)
lam.best.vec<-rep(0,nreps)
for(irep in 1:nreps){
	x <- mvrnorm(n=10,mu,Sigma)
	x1 <- x[,1]
	x2 <- x[,2]
	x3 <- rnorm(10,0,1)
	X <- cbind(x1,x2,x3)
	y <- x1+x2+x3+rnorm(10,0,0.1)
	t.ridge1 <- lm.ridge(y~x1+x2+x3,lambda=0)
	lam.seq<-seq(0,0.2,by=0.001)
	lam.best<-which.min(lm.ridge(y~x1+x2+x3,lambda=lam.seq)$GCV)
	lam.best<-as.numeric(names(lam.best))
	lam.best.vec[irep]<-lam.best
	t.ridge2 <- lm.ridge(y~x1+x2+x3,lambda=lam.best)
	ests.lm[irep,1:3]<-coef(t.ridge1)[2:4]
	ests.ridge[irep,1:3]<-coef(t.ridge2)[2:4]
}
hist(lam.best.vec)

par(mfrow=c(1,2),mar=c(4,2,4,1))
boxplot(ests.lm,main="OLS")
boxplot(ests.ridge,main="Ridge")
par(mfrow=c(1,1))

br12<-seq(-7.5,7.5,by=0.5)
br3<-seq(0.5,1.5,by=0.05)
par(mfrow=c(2,3),oma=c(0,0,2,0))
hist(ests.lm[,1],xlab=expression(beta[1]),main=expression(beta[1]),breaks=br12,ylim=range(0,500))
abline(v=1,col="red",lwd=2)
hist(ests.lm[,2],xlab=expression(beta[2]),main=expression(beta[2]),breaks=br12,ylim=range(0,500))
abline(v=1,col="red",lwd=2)
hist(ests.lm[,3],xlab=expression(beta[3]),main=expression(beta[3]),breaks=br3,ylim=range(0,500))
abline(v=1,col="red",lwd=2)
hist(ests.ridge[,1],xlab=expression(beta[1]),main=expression(beta[1]),breaks=br12,ylim=range(0,500))
abline(v=1,col="red",lwd=2)
hist(ests.ridge[,2],xlab=expression(beta[2]),main=expression(beta[2]),breaks=br12,ylim=range(0,500))
abline(v=1,col="red",lwd=2)
hist(ests.ridge[,3],xlab=expression(beta[3]),main=expression(beta[3]),breaks=br3,ylim=range(0,500))
abline(v=1,col="red",lwd=2)
title("Comparison of OLS and Ridge estimators", outer = TRUE )

apply(ests.lm,2,mean)
apply(ests.lm,2,var)

apply(ests.ridge,2,mean)
apply(ests.ridge,2,var)