Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: re: maximum likelihood problem


From   Kit Baum <baum@bc.edu>
To   statalist@hsphsun2.harvard.edu
Subject   st: re: maximum likelihood problem
Date   Sat, 8 Dec 2007 11:19:55 -0500

Ben wrote

The program I've written in Stata is the following
(which, for obvious reasons, isn't working):

program mymlprog
version 10.0
args y x

quietly replace `y' = ln(1/sqrt(2*_pi*Sigma_hat^2))-
((1/(2*Sigma_hat^2))*(`y'+Alpha_1-Beta_hat*`x')^2 if $ML_y1 <0

quietly replace `y' =
ln(normal((Alpha_2-Beta_hat*`x')/Sigma_hat)-normal((Alpha_1- Beta_hat*`x')/Sigma_hat))
if $ML_y1 == 0

quietly replace `y' =
ln[1/sqrt(2*_pi*Sigma_hat^2)]-(1/(2*Sigma_hat^2))*(`y'+Alpha_2- Beta_hat*`x')^2
if $ML_y1 > 0

end


I think this will do it. Note that there are some fundamental conceptual flaws in what Ben has written above. This can be viewed as a linear-form (lf) model of a rather strange sort (given that there are two parameter vectors that share an element). I think that could be handled by applying constraints that the two slope coefficients should be the same.

--------mymlprog.ado--------
program mymlprog
version 10.0
args lnf mu1 mu2

quietly replace `lnf' = log(1/sqrt(2*_pi*$Sigma_hat^2))- ///
(1/(2*$Sigma_hat^2))*($ML_y1+`mu1')^2 if $ML_y1 <0

quietly replace `lnf' =log(normal((`mu2')/$Sigma_hat)- ///
normal((`mu1')/$Sigma_hat)) if $ML_y1 == 0

quietly replace `lnf' =log(1/sqrt(2*_pi*$Sigma_hat^2))- ///
(1/(2*$Sigma_hat^2))*($ML_y1+`mu2')^2 ///
if $ML_y1 > 0

end
---------end of program------

to run, you must first set the global value of Sigma_hat to some value (as that does not seem to be a parameter to be estimated)

global Sigma_hat 10

I created a dataset with the desired characteristics with

webuse auto
g mpgzap = mpg - 25
ml model lf mymlprog (mpgzap = price) (mpgzap=price)
ml check
ml maximize


Number of obs = 74
Wald chi2(1) = 4.76
Log likelihood = -243.23578 Prob > chi2 = 0.0292

------------------------------------------------------------------------ ------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------- +----------------------------------------------------------------
eq1 |
price | .0008808 .0004039 2.18 0.029 . 0000891 .0016725
_cons | -2.170081 2.78924 -0.78 0.437 -7.63689 3.296728
------------- +----------------------------------------------------------------
eq2 |
price | .0050802 .0023531 2.16 0.031 . 0004683 .0096922
_cons | -19.44901 10.88796 -1.79 0.074 -40.78903 1.891009
------------------------------------------------------------------------ ------

If we now do

cons 1 [eq1]price = [eq2]price
ml model lf mymlprog (mpgzap = price) (mpgzap=price) , constraints(1)
ml maximize

the two slope parameters will be equated;

Number of obs = 74
Wald chi2(0) = .
Log likelihood = -246.37235 Prob > chi2 = .

( 1) [eq1]price - [eq2]price = 0
------------------------------------------------------------------------ ------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
------------- +----------------------------------------------------------------
eq1 |
price | .0010458 .0004001 2.61 0.009 . 0002615 .0018301
_cons | -3.476855 2.776426 -1.25 0.210 -8.918549 1.964839
------------- +----------------------------------------------------------------
eq2 |
price | .0010458 .0004001 2.61 0.009 . 0002615 .0018301
_cons | -.181225 2.8025 -0.06 0.948 -5.674024 5.311574
------------------------------------------------------------------------ ------

I think this is estimating the model that Ben has specified. As he has an independent
estimate of the model via Excel Solver, it should be easy to check that out.

Kit

Kit Baum, Boston College Economics and DIW Berlin
http://ideas.repec.org/e/pba1.html
An Introduction to Modern Econometrics Using Stata:
http://www.stata-press.com/books/imeus.html


*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/




© Copyright 1996–2022 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index