---
title: "Lab 3: Practice Key"
subtitle: "Ordinary Least Squares and Maximum Likelihood Estimation"
author: "Minji Jeong"
date: \today
output: pdf_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

# Good practice to remove all objects from the workspace
rm(list = ls())

# Use library() for packages you need, or source() for other R files.
library(tidyverse)

```

# Optimization algorithms


The purpose of this short coding exercise is to start practicing with the `optim()` function. This function allows us to estimate likelihood parameters from any Probability Density Function (PDF) using Maximum Likelihood Estimation (MLE). Recall that `optim()` has several key arguments:

  - **par**:  
  The initial guess or starting values for the parameters to be optimized. You must provide as many starting values as there are parameters to estimate.
  
  - **fn**:  
  The objective function to be optimized. For OLS, we can provide the Sum of Squared Residuals (SSR). In contrast, for MLE, we should provide a function that characterizes the log-likelihood of the data-generating process (DGP) of interest.
    - **IMPORTANT**: `optim()` is a minimizer by default. Therefore, when estimating likelihoods, you should multiply the log-likelihood function by $-1$, so you are maximizing the likelihood instead of minimizing it.

  - **method**:  
  Specifies the optimization algorithm to be used (e.g., BFGS, Nelder-Mead, etc.).
  
  - **hessian**: whether the Hessian matrix should be returned. This argument is required to compute the standard errors for ML estimates, but not for OLS.
  
  - **Other arguments**:  
  Additional arguments for the objective function can be passed to `optim()`. For example, you might pass the predictor variables ($X_1$, $X_2$) and response ($Y$) as additional arguments.


# Data

## Simulation example

Let's estimate a linear model with three parameters: $\beta_0$, $\beta_1$, and $\beta_2$, by minimizing the sum of squared residuals (SSR) using the `optim()` function. The model you want to estimate is:

$$
Y \sim \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon
$$

where:

  - $Y$ is normally distributed, and $\varepsilon \sim \mathcal{N}(0,1)$.
  - $X_1 \sim \mathcal{N}(10,2)$ and $X_2 \sim \mathcal{N}(5,1)$.
  - The true estimands are $\beta_0 = 5$, $\beta_1 = 2$, and  $\beta_2 = -1.5$.


```{r, fig.align='center', fig.width=5, fig.height=3.5}
# Set simulation parameters
set.seed(123)  # for reproducibility

n <- 100       # Sample size

b0 <- 5    # True intercept
b1 <- 2    # True coefficient for X1
b2 <- -1.5 # True coefficient for X2

# Simulate random variables

## Note: the square root in the sd argument of rnorm()

## Predictors
x1 <- rnorm(n, mean = 10, sd = sqrt(2))
x2 <- rnorm(n, mean = 5, sd = sqrt(1))

## Residuals
epsilon <- rnorm(n, mean = 0, sd = sqrt(1)) 

# DGP of the outcome variable
Y <- b0 + b1 * x1 + b2 * x2 + epsilon

# Visualize the distribution of Y
hist(Y, main = "Histogram of Y", xlab = "Y")


```



## Salaries dataset

**Description:** This dataset contains information on the 2008-09 nine-month academic salaries for Assistant Professors, Associate Professors, and Professors at a U.S. college. The data were collected as part of the college’s ongoing efforts to monitor salary differences between male and female faculty members.


Below, I am cleaning the data for you.

```{r}

# install.packages("carData")

# Load the data
data("Salaries", package = "carData")

# Prepare data

df <- 
  Salaries |>
  # optional: turn it into a tibble
  as_tibble() |>
  # re-scale salary
  mutate(salary=salary/1000) |> 
  # select predictors for analysis
  select(salary,sex,yrs.since.phd)

```




# Ordinary Least Squares


In this section, we will use the built-in dataset *Salaries* from the `car` / `carData` package, that we used in the previous lab code practice. However, before we will fit a toy example with simulated data. The goal is to fit least squares estimates using `optim()` instead of `lm()`.


## 1. Ordinary Least Squares estimation using OLS

Create a function called `SSR` that 

  1- Maps a parameter argument `params=` with each parameter to be estimated ($\beta_0, \beta_1$, and $\beta_2$).
  
  2- Provides an objective function to estimate the linear projection of $Y$. That is $\hat{Y} = \hat{\beta_0} + \hat{\beta_1} x1 + \hat{\beta_2} X_2$. Do not forget to include an intercept.
  
  3- An objective function to estimate the Sum of Squared Residuals (SSR) which will return an output. This output is the one `optim()` will optimize until finding the global minimum.

Once you have created `SSR()`, use `optim()` to estimate $\hat{\beta_0}$, $\hat{\beta_1}$, and $\hat{\beta_2}$. Double-check your results with `lm()`.


```{r}

# Define the residual sum of squares (SSR) function

SSR <- function(params, y, x1, x2) {
  
  # 1- Map a parameter arguments
  b0 <- params[1]  # Intercept
  b1 <- params[2]  # Coefficient for X1
  b2 <- params[3]  # Coefficient for X2
  
  # 2- Linear projection model of y (aka model to predict y)
  # rep(1,n) is constant
  y_hat <- b0 * rep(1, n) + b1 * x1 + b2 * x2
  
  # Objective funtion to estimate SSR
  SSR <- sum((y - y_hat) ^ 2)
  
  return(SSR)
}

```



### 2. Use `SSR()` and `optim()` to estimate the following model:  $Y_i = \alpha + \beta_1  X_{1i}  + \beta_2  X_{2i} + \epsilon$ using the *Salaries* dataset. Compare `optim()` results with `lm()`.

where:

  - **salary** ($Y$): Outcome variable representing the nine-month salary in dollars.
  - **sex** ($x1$): Binary indicator for the respondent’s reported sex (Female or Male).
  - **yrs.since.phd** ($X_2$): Count variable indicating the number of years since the respondent received their PhD.
    - Note that input the variable *sex* as **numeric** class.


```{r}

# Set vectors/matrix of covariates
y <- df$salary
x1 <- as.numeric(df$sex)
x2 <- df$yrs.since.phd

# adjust n for the intercept

n <- length(Y)

# Initial parameter guesses for b0, b1, and b2
# optim requrires initial parameters
stval <- c(b0 = 0, b1 = 0, b2 = 0)

# Use optim to minimize SSR
optim_ols <- optim(par = stval,
                   fn = SSR,
                   method = "BFGS",
                   y = y,
                   x1 = x1,
                   x2 = x2)


# Compare with lm()
lm_ols <- lm(y ~ x1 + x2)


# Display results:
cat("Minimum SSR:", optim_ols$value, "\n")
cat("Estimated parameters with optim (b0, b1, b2):", coef(lm_ols), "\n")
cat("Estimated parameters with lm (b0, b1, b2):", optim_ols$par, "\n")

```




# Maximum Likelihood Estimation

From `lab2_practice.Rmd`, we know that model 3 (the one estimated in the previous section) violates the assumption of homoskedasticity. Therefore, in this section, we will re-estimate the model using a heteroskedastic normal model via Maximum Likelihood Estimation (MLE).

## 3. Estimating a Heteroskedastic Normal Model Using MLE

To estimate the heteroskedastic normal model, you will first need to create a function that connects the parameters to the objective function for estimation.

  - You can modify the function we covered during the lab session or create your own. This time, ensure that the covariates/predictors are set into matrices for each parameter of the normal distribution that we want to model ($\beta$ and $\gamma$).



```{r}


# Set vectors/matrix of covariates

## vector of outcomes
y <- df$salary

## matrix of covariates for mean and variance parameters
xcov <- cbind(as.numeric(df$sex),df$yrs.since.phd)
zcov <- cbind(as.numeric(df$sex),df$yrs.since.phd)


llk.hetnorm <- function(param, y, x, z) {
  
  # 1- Map a parameter arguments, do not forget to include an intercept!
  x <- as.matrix(x)		             # x (some covariates) as a matrix for mean
  z <- as.matrix(z)		             # z (some covariates) as a matrix for variance
  
  os <- rep(1, nrow(x)) 	         # Set the intercept as 1 (constant)
  
  x <- cbind(os, x)		             # Add intercept to covariates x
  z <- cbind(os, z)                # Add intercept to covariates z
  
  # 2- Define vector of parameters (b and g)
  b <- param[1 : ncol(x)]                         # Parameters for x
  g <- param[(ncol(x) + 1) : (ncol(x) + ncol(z))] # Parameters for z
  
  # 3- Define systematic components
  # Systematic components for mean
  xb <- x %*% b
  
  # Systematic components for variance
  s2 <- exp(z %*% g)
  
  # 4- Find Objective function
  sum(0.5 * (log(s2) + (y - xb)^2 / s2)) # Likelihood we want to maximize 
}


# Initial guess, how many parameters we are going to estimate?
stval <- c(0, 0, 0, 0, 0, 0)  

# inverse the "value" (multiply -1)
# Run ML, and get the output we need
opmtim_mle <- optim(par = stval, # Initial guesses
                    fn= llk.hetnorm, # Likelihood function
                    method = "BFGS", # Gradient method
                    hessian = TRUE,  # Return Hessian matrix
                    y = y,           # Outcome variable
                    x = xcov,        # Covariates x (w/o constant)
                    z = zcov         # Covariates z (w/o constant)
) 

# Extract point estimates and standard errors

pe <- opmtim_mle$par            # Point estimates
vc <- solve(opmtim_mle$hessian) # Var-cov matrix (for computing s.e.)
se <- sqrt(diag(vc))            # Standard errors (s.e.)


# Compare with mle and lm() results

#options(tibble.width = Inf)

tibble(
  "Variables" = c("(Intercept)", "x1", "x2"),
  "betas pe" = pe[1:3],
  "betas se" = se[1:3],
  "gamma pe" = pe[4:6],
  "gamma se" = se[4:6]
)

round(coef(summary(lm_ols))[, 1:2],3)

# MLE-heteroskedastic normal is better estimator when we have heteroscedasticity???

```














