Tobit models

[Jack Serck]

I wrote my own function to fit a Tobit model, but it only works when censoring from below, not above (or both). I can't figure out why.

Suppose you have data censored from above at 1 and below at 0:

sysuse auto, clear
version 15
set more off

// generate fake data: lower limit is zero, upper limit is one
gen y =invnormal(runiform())
replace y = 0 if y < 0
replace y = 1 if y > 1

Here is the canned estimation:

qui tobit y length turn headroom, ul(1)
eststo m1_1

and here is my take:

capture program drop mytobit_lc
program mytobit_lc
args lnf beta sigma
tempvar lnlj
quietly {
gen double `lnlj' = log(1-normal(`beta'/`sigma')) if $ML_y1==1 
replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1>0
replace `lnf' = `lnlj'
}
end
ml model lf mytobit_lc (y = length turn headroom) /sigma
qui ml maximize
eststo m1_2

The bold line in the above chunk is the likelihood for censored observations. If I change that to capture left-censored observations, I get identical results as

 qui tobit y length turn headroom, ll(0)

. But it doesn't work for right-censored observations. What is going on?

Here are the full results:

. esttab m1_1 m1_2 m2_1 m2_2, mtitles("canned LC" "user LC" "canned RC" "user RC")

----------------------------------------------------------------------------
                      (1)             (2)             (3)             (4)  
                canned LC         user LC       canned RC         user RC  
----------------------------------------------------------------------------
main                                                                        
length           -0.00107        -0.00107         0.00203        -0.00265  
                  (-0.13)         (-0.13)          (0.42)         (-0.81)  

turn             -0.00137        -0.00137         -0.0126          0.0116  
                  (-0.04)         (-0.04)         (-0.54)          (0.73)  

headroom         -0.00454        -0.00454         -0.0203          0.0287  
                  (-0.04)         (-0.04)         (-0.29)          (0.61)  

_cons               0.286           0.286           0.494          0.0838  
                   (0.35)          (0.35)          (1.07)          (0.27)  
----------------------------------------------------------------------------
/                                                                          
var(e.y)            0.453***                        0.183***                
                   (3.74)                          (5.41)                  

sigma                               0.673***                        0.284***
                                   (7.47)                         (11.12)  
----------------------------------------------------------------------------
N                      74              74              74              74  
----------------------------------------------------------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

[Phil Bromiley]

This question was posed and responded to previously. The previous posting had a fuller set of code.

[Jack Serck]

The question was re-posted with a shorter, minimal example. The previous question was not answered (link). Phil Bromiley simply wrote a comment, saying the code was too long, and asking why anyone would want to code their own Tobit model. Well, what if you want to code a Tobit with inequality constraints? That requires hard-coding the log-likelihood function, which is the point of this exercise. Phil Bromiley did not respond to that comment.

[Nick Cox]

If repeated questioning doesn't explain my puzzlement, I sometimes quietly slip away from a thread. Whether that is happening here I can't say.

I didn't comment earlier because I too was puzzled at what you were doing. But here goes.

You set up an example with limits 0 and 1. However,

1. I don't understand why you didn't call tobit with both ul() and ll() options. That's the reference case.

2. Your own code looks buggy.

gen double `lnlj' = log(1-normal(`beta'/`sigma')) if $ML_y1==1 
replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1>0

The second statement will undo the first as 1 > 0.

There is no code for data that are exactly 0, which you need for your example.

3. The rationale for all this is still mysterious as tobit will cope with responses [0, 1]. You don't seem to have a problem that is different.

[Jack Serck]

Thanks. All I'm trying to do is learn how to code a Tobit model by hand in Stata. Why? Two reasons. First, it is useful practice. Second, it has a practical implication. What if you want to run a Tobit model with an inequality constraint? According to Stata, non-linear models with interval constraints must be coded by hand. If for example one wanted to constrain coefficients to be greater than or equal to zero, with this code once could simply transform the coefficients with

exp()

.

Sorry for the bug. The line you referred to should read:

gen double `lnlj' = log(1-normal(`beta'/`sigma')) if $ML_y1==0
replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1>0

and it reproduces the results from

 tobit ..., ll(0)

.

For the right-censored Tobit, I change the lines to

gen double `lnlj' = log(1-normal(`beta'/`sigma')) if $ML_y1==1
replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1<1

so that once again censored observations are modeled by the CDF and uncensored observations by the PDF. But this does not recreate

 tobit ..., ul(1)

.

[Christophe Kolodziejczyk]

The first and second lines of your likelihood should be

 gen double `lnlj' = log(normal((`beta'-1)/`sigma')) if $ML_y1==1 
replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1<1

and generally if `ll' and `ul' are respectively the left and right-censoring limits

gen double `lnlj' = log(normal((`beta'-`ul')/`sigma')) if $ML_y1==`ul'
gen double `lnlj' = log(normal((`ll'-`beta')/`sigma')) if $ML_y1==`ll'
replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1<`ul'&$ML_y1>`ll'

[Jack Serck]

Thanks, but your code throws an error.

Code (working example):

sysuse auto, clear
gen y =invnormal(runiform())
replace y = 0 if y < 0
replace y = 1 if y > 1
capture program drop mytobit_twolimit
program mytobit_twolimit
    args lnf beta sigma
    tempvar lnlj
    quietly {
        gen double `lnlj' = log(normal((`beta' - 1 ) / `sigma')) if $ML_y1 == 1
        gen double `lnlj' = log(normal((0 - `beta')/`sigma')) if $ML_y1 == 0
        replace `lnlj' = log((1/`sigma')*normalden(($ML_y1-`beta')/`sigma')) if $ML_y1 < 1 & $ML_y1 > 0
        replace `lnf' = `lnlj' 
    }
end
ml model lf mytobit_twolimit (y = length turn headroom) /sigma
qui ml maximize, nolog

Error:

variable __000009 already defined
r(110);

Is the error saying the program is trying to redefine one of the macros when it shouldn't?

[Jack Serck]

Ignore my post. The error comes from the local macro `lnlj' being generated twice. Replacing the second call to `gen double` with `replace` resolves the error. My mistake.

[Vasudha Chopra]

Hey I modified your code and this works for right censoring see below,
capture program drop tobit_ul // works perfectly
program tobit_ul
args lnf xb sigma
tempvar lnlj
quietly {
gen double `lnlj'= ln(1-normal(($ML_y1 -`xb')/`sigma')) if $ML_y1==3.5
replace `lnlj'= ln((1/`sigma')*normalden(($ML_y1 -`xb')/`sigma')) if $ML_y1<3.5
replace `lnf' = `lnlj'
}
end
ml model lf tobit_ul (b: colgpa1 = tothrs athlete female black hsize hsperc)(sigma
ml maximize, difficult tolerance(1e-6) ltolerance(1e-7)


Btw, thanks a lot for your code, super helpful.

[Vasudha Chopra]

There isn't supposed to be a smiley btw it is : and ) written together after sigma.