rm(list=ls()) # Clear memory set.seed(123456) # For reproducible random numbers library(MASS) # Load packages library(simcf) library(tidyverse) n <- 1500 # Generate 1500 observations w0 <- rep(1, n) # Create the constant w1 <- runif(n, min = 0, max = 1) # Create two covariates w2 <- runif(n, min = 0, max = 1) x <- cbind(w0, w1, w2) # Create a matrix of the covariates z <- x # i.e., same covariates affect mu and sigma beta <- c(0, 5, 15) # Set a parameter vector for the mean gamma <- c(1, 0, 3) # Set a parameter vector for the variance mu <- x %*% beta # Create systemtic component for the mean sigma2 <- exp(z %*% gamma) # Create systematic component for the variance y <- rnorm(n = n, # Create the outcome variable mean = mu, # Think about the stochastic component! sd = sqrt(sigma2) ) data <- cbind(y, w1, w2) # Save the data to a data frame data <- as.data.frame(data) names(data) <- c("y", "w1", "w2") ls.result <- lm(y ~ w1 + w2, data = data) # Fit a linear model using lm() # A likelihood function for ML heteroskedastic Normal llk.hetnormlin <- function(param, y, x, z) { x <- as.matrix(x) # x (some covariates) as a matrix z <- as.matrix(z) # z (some covariates) as a matrix 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 b <- param[1 : ncol(x)] # Parameters for x g <- param[(ncol(x) + 1) : (ncol(x) + ncol(z))] # Parameters for z xb <- x %*% b # Systematic components for mean s2 <- exp(z %*% g) # Systematic components for variance sum(0.5 * (log(s2) + (y - xb)^2 / s2)) # Likelihood we want to maximize } # Create input matrices xcovariates <- cbind(w1, w2) zcovariates <- cbind(w1, w2) # Initial guesses of beta0, beta1, ..., gamma0, gamma1, ... stval <- c(0, 0, 0, 0, 0, 0) # Run ML, and get the output we need hetnorm.result <- optim(stval, # Initial guesses llk.hetnormlin, # Likelihood function method = "BFGS", # Gradient method hessian = TRUE, # Return Hessian matrix y = y, # Outcome variable x = xcovariates, # Covariates x (w/o constant) z = zcovariates # Covariates z (w/o constant) ) pe <- hetnorm.result\$par # Point estimates vc <- solve(hetnorm.result\$hessian) # Var-cov matrix (for computing s.e.) se <- sqrt(diag(vc)) # Standard errors (s.e.)