# Now we embark on one numerical experiment
set.seed(101)
X <- cbind(1, runif(1000))
theta.true <- c(2,4,6) # error variance = 2, intercept = 4, slope = 6.
cat("True theta = ", theta.true, "\n")
y <- X %*% theta.true[-1] + sqrt(theta.true[1]) * rnorm(1000)

# Estimation by OLS --
d <- summary(lm(y ~ X[,2]))
theta.ols <- c(sigma2 = d$sigma^2, d$coefficients[,1])
cat("OLS theta = ", theta.ols, "\n\n")

# A function that computes the SSE of a given parameter vector w.r.t theta.ols
SSE <- function(x) {e <- x - theta.ols; return(sum(e*e))}