I am attempting to run STATA code from Abbott (2009) (cited below). I am attempting to run this to analyse married couple's labour particpiation decisions, so these are the two main base equations I am analysing,

\[y_m = x_b + \alpha y_w + \epsilon^m\]
\[y_w = x_b + \alpha y_m + \epsilon^w\]

where y_m is a dummy variable for husband decisions to work and vice versa for y_w.

Here is a sample of the code from the paper.

Code:
program define nash_lf
              version 10.0
              args lnf xb1 xb2 alpha1 alpha2 rho
quietly replace `lnf' = ln( binormal( -`xb1', -`xb2', `rho') –
0.5 * ( binormal( -`xb1', -`xb2', `rho' ) –
binormal( -`xb1', -`xb2' -`alpha2', `rho' ) –
binormal( -`xb1' - `alpha1', -`xb2', `rho' ) +
                            binormal( -`xb1' - `alpha1', -`xb2' - `alpha2', `rho' ) ) )
if `alpha1' >= 0 & `alpha2' >= 0 & $ML_y1==0 & $ML_y2==0 ;
The problem I am facing is the actual running of the program, i.e how I write the equations. The other control variables I am including would all be estimated with the xb part of the code but the problem comes in telling STATA that the alpha coeeficient is for a specific factor (whihc later becomes the dependent variable in the second equation). I am unsure on the format/syntax of the ml model lf command for this purpose.


Source: Abbott, m.r. (2009), modelling decisions to volunteer at a household level, b.s. honors monograph, university of new south wales- school of economics