OLS optimization in R (optim)
In this post, we will discuss how to perform a simple optimization in R. I will show how to perform an OLS estimation using the optim
function in R
. Let’s get started.
You have the following regression equation, where \(y\) is linearly dependent on four independent variables (\(x\)).
\begin{equation} y_{i} = 2 x_{1i} + 3 x_{2i} + x_{3i} + 2 x_{4i} \end{equation}
Let’s say that you only observe the variables \(x_{1i}, x_{2i}\). You want to estimate \(\beta_{1}, \beta_{2}\) from your data.
\begin{equation} y_{i} = \beta_{1} x_{1i} + \beta_{2} x_{2i} + \epsilon_{i} \end{equation}
The OLS regression is an optimization problem that estimates the values of the betas (\(\beta\)) to minimize the sum of squared errors, which are the differences between the predicted values from your model and the observed values of \(y\).
This is called the “objective function”
\begin{equation} \sum_{i=1}^{n}(y - \underbrace{ \beta_{1} x_{1i} + \beta_{2} x_{2i}}_{\text{predicted values}} )^{2} \end{equation}
With R
, it is very simple the estimates the parameters using an optimization strategy. The function is called optim
.
First, let’s simulate the data with
# our beta values
betas = c(2, 3, 1, 2)
# sample size
n=1000
X = matrix(rnorm(n), ncol = 4)
# predict y as linear function of Xb
y = X%*%betas
X_obs = X[,1:2]
Now let’s write the objective function. We write a function in R
f_min_SS = function(y, X_obs, par){
SSQ = sum( (y - (X_obs %*% par))^2)
return(SSQ)
}
We can test different values for the parameters and see what the sum of squares is.
f_min_SS(y, X_obs, par = c(2,2))
f_min_SS(y, X_obs, par = c(2,3))
Now we can use optim
. We feed the function starting values, here I put 1, 1
. Then add the inputs of the function, the \(y\) and the \(X\). In the fn
, we input the objective function, fn = f_min_SS
. Finally, we can add some lower and upper bounds for the betas values.
This is the command
optim(par = c(1,1), y = y, X_obs = X_obs,
fn = f_min_SS,
lower = c(-5, -5), upper = c(5,5),
method = "L-BFGS-B")
We can compare the results from a simple OLS with
summary(lm(y ~ X_obs))