*! version 1.0.0  12apr2004 M. Lokshin, Z. Sajaia              (SJ4-3: st0000)
#delim ;

cap program drop mspredict;
program define mspredict;
        version 8;
	syntax anything   [if] [in] [, PSel xb1 xb2 yc1_1 yc1_2 yc2_1 yc2_2 mills1 mills2];

	marksample touse;

        syntax newvarname [if] [in] [, PSel xb1 xb2 yc1_1 yc1_2 yc2_1 yc2_2 mills1 mills2];

        if ("`e(cmd)'" ~= "movestay") error 301;
        local check="`psel' `xb1' `xb2' `yc1_1' `yc1_2' `yc2_1' `yc2_2' `mills1' `mills2'";
        if (wordcount("`check'")>1)  {; disp as error "Only one statistic is allowed" ;  error;};
	if (wordcount("`check'")==0) local nooptions = 1; else local nooptions = 0;

        quietly {;

        local r1x1  : word 1 of `e(ivar_reg1)';
        local r2x1  : word 1 of `e(ivar_reg2)';
        local px1   : word 1 of `e(ivar_prob)';
        local y_reg1: word 1 of `e(depvar)';
        local y_reg2: word 2 of `e(depvar)';
        local y_prob: word 3 of `e(depvar)';

        tempname b TMP sigma1 sigma2 rho1 rho2 reg1 reg2 prob;
             matrix `b'   = e(b);
             matrix `TMP' = `b'[1,"lns1:_cons"];                        scalar `sigma1' =  exp(`TMP'[1,1]);
             matrix `TMP' = `b'[1,"lns2:_cons"];                        scalar `sigma2' =  exp(`TMP'[1,1]);
             matrix `TMP' = `b'[1,"r1:_cons"]  ;                        scalar `rho1'   = tanh(`TMP'[1,1]);
             matrix `TMP' = `b'[1,"r2:_cons"]  ;                        scalar `rho2'   = tanh(`TMP'[1,1]);
             matrix `TMP' = `b'[1,1.."`y_reg1':_cons"];		      	matrix `reg1'   = `TMP';
             matrix `TMP' = `b'[1,"`y_reg2':`r2x1'".."`y_reg2':_cons"]; matrix `reg2'   = `TMP';
             matrix `TMP' = `b'[1,"`y_prob':`px1'".."`y_prob':_cons"];  matrix `prob'   = `TMP';

        tempvar xb_prob x_b1 x_b2;
        generate `xb_prob' = 0 if `touse';
        local i = 1; local param = "`e(ivar_prob)' 1";
        while (`prob'[1,`i']~=.) {;
              local x : word `i' of `param';
              replace `xb_prob' = `xb_prob' + `prob'[1,`i']*`x' if `touse';
              local i = `i' + 1;
              }; /* end while */

        if (("`psel'"~="") | (`nooptions')) {;  	    					// psel or default
           noisily display as text "(Option psel assumed; Pr(`y_prob'))";
           generate `varlist' = norm(`xb_prob') if `touse';
           exit;
           }; // endif

        if ("`mills1'"~="") {;
		generate `varlist' = normden(`xb_prob')/norm(`xb_prob')     if `touse';
		exit; };
        if ("`mills2'"~="") {;
		generate `varlist' = normden(`xb_prob')/(1-norm(`xb_prob')) if `touse';
		exit; };

        generate `x_b1' = 0 if `touse'; generate `x_b2' = 0 if `touse';
        local i = 1; local param = "`e(ivar_reg1)' 1";
        while (`reg1'[1,`i']~=.) {;
              local x : word `i' of `param';
              replace `x_b1' = `x_b1' + `reg1'[1,`i']*`x' if `touse';
              replace `x_b2' = `x_b2' + `reg2'[1,`i']*`x' if `touse';
              local i = `i' + 1;
              }; /* end while */

        if ("`yc1_1'"~="") generate `varlist' = `x_b1'+`sigma1'*`rho1'*normden(`xb_prob')/norm(`xb_prob')     if `y_prob'==1;
        if ("`yc1_2'"~="") generate `varlist' = `x_b2'-`sigma1'*`rho1'*normden(`xb_prob')/(1-norm(`xb_prob')) if `y_prob'==1;
        if ("`xb1'"  ~="") generate `varlist' = `x_b1';

        if ("`yc2_1'"~="") generate `varlist' = `x_b1'+`sigma2'*`rho2'*normden(`xb_prob')/norm(`xb_prob')     if `y_prob'==0;
        if ("`yc2_2'"~="") generate `varlist' = `x_b2'-`sigma2'*`rho2'*normden(`xb_prob')/(1-norm(`xb_prob')) if `y_prob'==0;
        if ("`xb2'"  ~="") generate `varlist' = `x_b2';

	replace `varlist' = . if `touse'==0;

        }; // end quietly

end; // end mspredict