clear
set mem 200m
set more off
*cd "c:\upcd1\work\stata_prg" /* please change to your own working directory */

log using mle_probit_ex.log, replace

*webuse auto.dta, clear 
 /* to use data directly from STATA official website */
*save auto.dta, replace 
 /* to save the data set you just netted from the website to your working directory */

*webuse weib.dta,clear
 /* to use data directly from STATA official website */
*save weib.dta, replace
 /* to save the data set you just netted from the website to your working directory */

***************************
* Overview of ml *
***************************

**********************************************
* one equation one dependent variable *
**********************************************

use auto.dta, clear
capture program drop myprobit

program myprobit 
 /* to use maximum likelihood method to estimate the probit model */
 version 8.0
 args lnf theta1
 quietly replace `lnf' = ln(norm(`theta1')) if $ML_y1==1
 quietly replace `lnf' = ln(norm(-`theta1')) if $ML_y1==0
end

ml model lf myprobit (foreign = mpg weight) 
  /* to estimate the basic probit model using method lf */
ml check 
  /* to verify that the log-likelihood evaluator you have written seems to work */
ml search 
  /* to find a better initial value so that the convergence of log-likelihood takes less time*/
ml maximize
ml graph 
  /* to graph the log-likelihood values against the iteration number to verify the process of convergence */

ml model lf myprobit (foreign=mpg weight, nocons) 
  /* to estimate the restricted model without the intercept */
ml maximize

ml model lf myprobit (foreign=) 
  /* to estimate the restricted model whose only explantory variable is the constant */
ml maximize

constraint define 1 _b[mpg]==0 /* to define a constraint which drops a coefficient */
ml model lf myprobit (foreign = mpg weight), constraint (1) 
  /* to estimate the restricted model using ML*/
ml maximize

ml model lf myprobit (foreign = mpg weight), robust 
  /* to specify the robust variance estimator and report it */
ml maximize

ml model lf myprobit (foreign = mpg weight) if mpg<= 25 
  /* to specify a subsample to estimate the model*/
ml maximize

ml model lf myprobit (foreign = mpg weight)
ml maximize, gradient hessian trace 
  /* to add to the iteration log a report on the current gradient vector,
   current negative Hessian matrix and current parameter vector */

ml model lf myprobit (foreign = mpg weight) 
ml maximize,level(99) 
  /* to change the significance level from 95 to 99*/


***************************************************
* two equations, one dependent variable *
***************************************************

use auto.dta,clear
capture program drop myreg 
  /* to use maximum likelihood method to estimate the linear regression model */
program myreg
version 8.0
args lnf theta1 theta2 
  /* to use theta2 to stand for the standard deviation, so this is a two equation model */
quietly replace `lnf' = ln(normd(($ML_y1-`theta1')/`theta2'))-ln(`theta2') 
end

ml model lf myreg (mpg =weight displ)/sigma 
  /* to use sigma to specify the equation for standard deviation */
ml check
ml search
ml maximize
ml graph

ml model lf myreg (mpg =weight displ) () 
  /* an equivalent specification as the former one which uses sigma */
ml maximize

ml model lf myreg (mpg =weight displ) (sigma:) 
  /* another equivalent way to specify the second equation, and the name  
   of the second quation has been changed to sigma */
ml maximize

ml model lf myreg (mpg =weight displ) (sigma: weight) 
  /* to estimate the conditional heteroskedasticity model, and to
   specify that standard deviation depends on weight and constant */
ml maximize

ml model lf myreg (mpg =weight displ) (sigma: weight, nocons) 
  /* to specify the second equation without intercept */
ml maximize

***********************************************
* two equations, two dependent variables *
*********************************************** 

use weib.dta,clear
capture program drop myweib 
  /* to use maximum likelihood method to estimate the Weibull model */ 

program define myweib
 version 8.0
 args lnf theta1 theta2
 tempvar p M R 
   /* to make programming easier by introducing temporary variabls*/
 quietly gen double `p' = exp(`theta2') 
   /* temporary variables should be generated as doubles */
 quietly gen double `M' = ($ML_y1*exp(-`theta1'))^`p'
 quietly gen double `R' = ln($ML_y1)-`theta1'
 quietly replace `lnf' = -`M' + $ML_y2*(`theta2'-`theta1'+(`p'-1)*`R')
end

ml model lf myweib (studytime died =drug2 drug3 age) () 
  /* the () corresponds to the Weibull shape parameter s*/
ml check
ml search 
ml maximize
ml graph

ml model lf myweib (studytime died =drug2 drug3 age) /sigma 
  /* An alternative way to specify the second equation*/
ml maximize

ml model lf myweib (studytime died =drug2 drug3 age) /s 
  /* Another equivalent way to specify the second equation*/
ml maximize

**********************************************************************
* another example of two equations, two dependent variables model *
**********************************************************************

use auto.dta,clear
capture program drop myreg1 
program myreg1
version 8.0
args lnf theta1 theta2 
quietly replace `lnf' = ln(normd(($ML_y1-`theta1')/`theta2'))-ln(`theta2') 
end

ml model lf myreg1 (mpg =weight displ) (price=weight) 
 /* the second dependent variable is price which depends on weight and constant*/ 
ml maximize

ml model lf myreg1 (mpg price=weight displ) (weight) 
 /* an equivalent specification*/
ml maximize

ml model lf myreg1 (mpg =weight displ) (price=weight, nocons) 
 /* the second dependent variable depends only on weight*/ 
ml maximize

ml model lf myreg1 (mpg price=weight displ) (weight, nocons) 
 /* an equivalent specification of the restricted model*/
ml maximize

ml model lf myreg1 (mpg =weight displ) (price=) 
 /* the second dependent varible is a constant, in this case
  the model degenerates into one dependent variable model*/
ml maximize

ml model lf myreg1 (mpg =weight displ)/s 
  /* an equivalent specification of the one dependent variable model*/
ml maximize