I would like to implement multinomial logit model using maximum likelihood command to generate the same result as mlogit.
Pr(y=1) = exp(beta1*x)/(1+beta1*x+beta2*x)
Pr(y=1) = exp(beta2*x)/(1+beta1*x+beta2*x)
Pr(y=1) = 1/(1+beta1*x+beta2*x)

The following code did not work (not generating anything), but I have no idea how to fix it.
Code:
program define mymlogit
    args lnf Xb1 Xb2
    quietly replace `lnf' = -`Xb1' - ln(1+exp(-`Xb1')+exp(-`Xb2')) if $ML_y1==1
    quietly replace `lnf' = -`Xb2' - ln(1+exp(-`Xb1')+exp(-`Xb2')) if $ML_y1==2
    quietly replace `lnf' = -ln(1+exp(-`Xb1')+exp(-`Xb2')) if $ML_y1==3
end

ml model lf mymlogit (y= x1 x2)
It would be awesome if anyone could also explain how does ml model return values to the program.