4
$\begingroup$

I am considering the following censored regression model. I have a question about my R code.

enter image description here

Then I obtain likelihood function as follow. [
logL(\beta,\sigma_u) = \sum_{i\in I} \qty[\mathbbm{1}{y_i > 0}\qty[log\phi\qty(\frac{y_i-x_i^T\beta}{\sigma_u})-log(\sigma_u)] + \mathbbm{1}{y_i = 0}log\qty[1-\Phi\qty(\frac{x_i^T\beta}{\sigma_u})]]
]

Then i obtain log-likelihood function. enter image description here

##Data Generation#################################################
set.seed(123)
n <- 1000; u <- rnorm(n);x<- rnorm(n,0,sqrt(2))
y_star <- 1 + x + u
y <- ifelse(y_star > 0, y_star, 0)
data <- data.frame(x,y)

##Censored Regression##############################################
Censored_reg <- function(initial_value,y_in, x_in) {
  f_censored <- function(input, y_in, x_in) {
    sigma <- input[1]
    beta <- input[2:length(input)]
    xb <- cbind(1, x_in) %*% beta
    z <- (y_in - xb) / sigma
    y.indic <- ifelse(y_in > 0, 1, 0)
    
    log_lik <- sum(y.indic*(log(dnorm(z))-log(sigma)) + (1 - y.indic)*log(1 - pnorm(z)))
    return(-log_lik)
  }
  
  result <- optim(par = initial_value, fn = f_censored, y_in = y_in, x_in = x_in)
  return(result)
}
#log_lik <- -sum(y.indic*(log(dnorm(z))-log(sigma)) + (1 - y.indic)*log(1 - pnorm(z)))

# Example usage:
result <- Censored_reg(initial_value = c(2,2,2),y_in = data$y, x_in = data$x)
print(result)

I thought it works well, however, the result does not coinside with the results from censReg package.

library(quantreg)

rq <- rq(y ~ 1 + x, tau = 0.5, data = data)
screenreg(rq) #library(texreg) # similar to "stargazer"

I am very appriciate if you find my error. Thank you in advance.

Chat GPT does not understand Censored Reg, so pls help me.

$\endgroup$
2
  • $\begingroup$ LaTeX is not built into SO. Perhaps it'd be better if you could render those equations and pasting them as images into SO. $\endgroup$
    – duffymo
    Commented Sep 6, 2023 at 11:17
  • $\begingroup$ pnorm(z) does not correspond to the formula of the loglik. Should be pnorm(xb / sigma). $\endgroup$ Commented Sep 6, 2023 at 11:32

1 Answer 1

5
$\begingroup$

After replacing pnorm(z) with pnorm(xb/sigma), I find:

> print(result)
$par
[1] 0.9889646 1.0033864 1.0710261

censReg gives:

> censReg(y ~ 1 + x, data = dat)

Call:
censReg(formula = y ~ 1 + x, data = dat)

Coefficients:
(Intercept)           x    logSigma 
    1.00337     1.07109    -0.01101 

The results are almost the same.

$\endgroup$
1
  • 1
    $\begingroup$ Dear Stéphane Laurent Thank you for finding my mistake!! I considered more than an hour. sincerely appreciate your support : ) $\endgroup$
    – Hiroki
    Commented Sep 6, 2023 at 12:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.