/* 
** THIS IS FREESHARE GAUSS CODE: THE USUAL CAVEATS APPLY                                        **
** ESTIMATES:                                                                                   **
** 1. MNL MODEL                                                                                 **
** 2. (STANDARD) OGEV MODEL: Small, A. K. (1987). "A discrete Choice Model for Ordered          **
** Alternatives", Econometrica, 55(2): 409-424.                                                 **
** 3. DOGIT MODEL: Gaudry, M., and M. Dagenais (1979) "The Dogit Model,"                        **
** Transportation Research-B, 13B:105-112.                                                      **
** 4. DOGEV MODEL: FRY, T.R.L. and HARRIS, M.N. (2005); The Dogit Ordered Generalised Extreme   **
** Value Model, Australian and New Zealand Journal of Statistics, Volume 47 No. 4 December,     **
** pp. 531 - 542.                                                                               **
** 5. ORDERED PROBIT MODEL                                                                      **
** AUTHORS USING THIS CODE SHOULD CITE FRY AND HARRIS (2005)                                    **
** ADDITIONALLY AUTHORS USING THE OGEV CODE SHOULD CITE                                         **
** HARRIS, M.N., RAMFULL, P. and ZHAO, X. (2006); An Ordered Generalised Extreme Value Model    **
** with Application to Alcohol Consumption in Australia, Journal of Health Economics,           **
** Volume 25, pp. 782 - 801.                                                                    **
*/


			/* GEVS.PRC: MULTINOMIAL CHOICE MODELS */

/* 
** Estimates: Multinomial Logit, Ordered Extreme Value, DOGIT and DOGEV models
** Takes form:
**{reg, logit, vcovlog, logllgt, plgt, retlgt,
	dogit, vcovdog, logldgt, pdgt, retdgt,
	ogev, vcovogv, loglogv, pogv, retogv,
	dogev, vcovdgv, logldgv, pdgv, retdgv,
	ordprb, vcovop, loglopb, pordp, retopb} = mnmodels("alc_dat",model,thetas,
										startmnl,startdgt,startogv,startdgv,startopb,res); 

** Input: 	dataset 	- first column is y
**		model 	- string - request models: logit, dogit, ogev, dogev and ordered
**		thetas 	- a j by 1 vector, j'th row = 0 if theta_j is restricted to 0
**		startmnl 	- optional - user supplied starting values, a 0 takes defaults
**		startdgt 	- optional - user supplied starting values, a 0 takes defaults
**		startogv 	- optional - user supplied starting values, a 0 takes defaults
**		startdgv 	- optional - user supplied starting values, a 0 takes defaults
**		startopb 	- optional - user supplied starting values, a 0 takes defaults
**          	res         	- a 1 by k vector = k'th colum = 1 if beta_k is constant across j 
** Output: all model types give model and logit, apart from dogev which gives all variants   	
**			- reg: OLS
**			- logit: parameter estimates and t-stats
**			- logllgt: maximised LogL
**			- plgt: n by J matrix of estimated probablilities for choice j
**			- retlgt: string - return code
**			- dogit: parameter estimates and t-stats (incudes theta's)
**			- logldgt: maximised LogL
**			- pdgt: n by J matrix of estimated probablilities for choice j
**			- ogev: ordered EV parameter estimates and t-stats (incude's rho)
**			- loglogv: maximised LogL
**			- pogv: n by J matrix of estimated probablilities for choice j
**			- dogev: ordered EV DOGIT parameter estimates and t-stats (plus theta's & rho)
**			- logldgv: maximised LogL
**			- pdgv: n by J matrix of estimated probablilities for choice j
**			- ordprb: ordered probit parameter estimates and t-stats (plus mu's)
**			- loglop: maximised LogL
**			- pordp: n by J matrix of estimated probablilities for choice j
**			- NOTE NORMALISATION IS ON FIRST CHOCE
*/


                              /* PROCEDURE TO ESTIMATE ORDER DOGIT MODEL */
cls;
trap 0;
 
proc 36 = mnmodels(dset,model,thetas,startmnl,startdgt,startogv,startdgv,startopb,res,M);
local infile1, fh, varsz, y, z, k, n, zcheck, zchkvec, cdo; 
local zrep, olsparam, vnames, mtype, capj, jrep, yj, btmp; 
local bols, logit, tlogit, maxflgt, ogev, togev, maxfogv, dogev, tdgev, maxfdgv;
local glgt, retlgt, lgtse;
local gogev, retogv, ogevse;
local gdogev, retdgv, dgvse;
local theta, dgtbndsa, dgtbndsb, dgtbnds;
local dgvbndsa, dgvbndsb, dgvbnds; 
local dogit,tdogit,maxfdgt,gdgt,retdgt,dgtse; 
local ordpbt, opbse, topb, maxfopb, gopb, retopb, mu;
local fhthetas, twrite, c;
local d_ij, phat, startmu1, startmu, startrep, mark, startcon;
local betahat, v_ij, plgt;
local rho, jj, v_ijall, p_denom, p_a, p_b, p_ogev, pdgt; 
local pogv, pdgv, pordp, thetahat, zbetahat, muhat, p_ij, jrep2, check;
local noxs, dgtchk, dgvchk;
local CC, DD, thetast;
local vnam,mean,var,std,min,max,valid,mis;
local resbols, betahatr, beta1, betanr, beta, b_rep, row_rep, lgtse1, lgtsenr, se, lgtt1, lgttnr, t;
local beta_1, dgtse_1, tdgt_1, dgtsenr, tdgtnr, sebeta, tbeta;
local ogvse_1,	 togv_1, ogvsenr, togvnr;
local dgvse_1,	 tdgv_1, dgvsenr, tdgvnr;
local covlgt, covogev, covdgev, covopb, covdgt;
local mark0, mark1, mark2;
local zbar, g, mnl_mes, ogv_mes, dgt_mes, dgv_mes, op_mes;
local mestmp, mesall, mnl_mese, ogv_mese, dgt_mese, dgv_mese, op_mese;
local gam_sim;

clearg j, b_me, theta_me, rho_me, mu_me;


/* WHICH MODELS? */

let mtype = logit dogit ogev dogev ordered; 
mtype = indcv(model,mtype); 


/* READ IN FULL DATA SET */

infile1 = dset;
open fh = ^infile1 varindxi;
varsz = getname(infile1);
varsz = varsz[2:rows(varsz)];

z 	= readr(fh,rowsf(fh)); 
y	= z[.,1];
z 	= z[.,2:cols(z)];
k 	= cols(z);
n 	= rows(z);
fh 	= close(fh);

y 	= y - minc(y);						@ run y from 0 to J-1 @
capJ	= maxc(y) - minc(y) + 1;


			/* CHECK FOR SUFFICIENT VARIATION IN Z */

if k gt 1;
	zcheck 	= stdc(z[.,2:k]);
	zchkvec 	= (zcheck .le 0.1);

	zrep=1;
	do until zrep>cols(z)-1;
      		if zchkvec[zrep] eq 1;
        		"Potentially Insufficient Variation in:";$varsz[zrep+1];
        		zrep = zrep + 1;
      		else;
        		zrep =zrep+1;
      		endif;
	endo;
endif;

if sumc(res') eq k;
    if stdc(z[.,1]) eq 0;
        ?";ERROR CONDITIONAL LOGIT: DO NOT INCLUDE OVERALL CONSTANT"; 
    endif;  
endif;


			/* SET DIRECTORIES */

cdo 		= cdir(0);
load path 	= ^cdo;
save path 	= ^cdo;


			/* SAVE THETA'S & RESTRICTION VECTOR */

save thetas, res, capJ; 


			/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

jrep = 0;
d_ij	= zeros(n,1);
do while jrep lt capj;
	d_ij = d_ij~(y .eq jrep);
jrep=jrep+1;
endo;

d_ij = d_ij[.,2:capj+1];

phat = meanc(d_ij);

save phat, d_ij;


			/* ESTIMATE OLS STARTING VALUES - IGNORE BETA_0 */


if rows(startmnl) eq 1 and sumc(res') lt k;

	jrep		= 1;
	bols		= zeros(capJ-1,cols(z));

	do while jrep lt capJ;
		btmp = inv(z'*z)*z'*(y .eq jrep); 
		bols[jrep,.] = btmp';
	jrep=jrep+1;
	endo;

endif;

if sumc(res') eq 0 and rows(startmnl) eq 1;

	startmnl	= vec(bols'); 

elseif sumc(res') ge 1 and rows(startmnl) eq 1;

    if sumc(res') lt k;

	    startmnl	= bols[1,.]';
	    resbols	    = delif(bols',res')';
	    resbols	    = resbols[2:capJ-1,.];
	    startmnl 	= startmnl|vec(resbols');
	
    endif;

endif;

noxs	= rows(startmnl);

save noxs;


if sumc(res') eq k;

    noxs        = k;
    startmnl    = y/z; 
/*
    startmnl    = zeros(capJ-1,1)|startmnl;            @ include choice specific constants @
*/

endif;



			/* ESTIMATE MNL MODEL */

cmlset;
_cml_dirTol 	= 0.0001;
_cml_maxIters 	= 100;			
@_cml_maxIters 	= 1;			@
__output 		= 1;

_cml_covpar 	= 1; 
_cml_algorithm	= 3; 		/* 1 = BFGS, 2=DFP, 3= NEWTON, 4=BHHH */
__title 			= "Multinomial Logit";
if sumc(res') ge 1 and sumc(res') lt k and sumc(res') lt k;
	_cml_parnames 	= varsz|(ones(capJ-2,1) .*. delif(varsz,res'));  
elseif sumc(res') eq 0;
	_cml_parnames 	= ones(capJ-1,1) .*. varsz;  
elseif sumc(res') eq k;

    if stdc(z[.,1]) eq 0;
    
    endif;

/*
    if capJ eq 3;
        varsz = "Choice 2 Con"|"Choice 3 Con"|varsz;
    elseif capJ eq 4;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|varsz;
    elseif capJ eq 5;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|"Choice 5 Con"|varsz;
    elseif capJ eq 6;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|"Choice 5 Con"|"Choice 6 Con"|varsz;
    elseif capJ eq 7;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|"Choice 5 Con"|"Choice 6 Con"|"Choice 7 Con"|varsz;
    elseif capJ eq 8;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|"Choice 5 Con"|"Choice 6 Con"|"Choice 7 Con"|"Choice 8 Con"|varsz;
    elseif capJ eq 9;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|"Choice 5 Con"|"Choice 6 Con"|"Choice 7 Con"|"Choice 8 Con"|"Choice 9 Con"|varsz;
    elseif capJ eq 9;
        varsz = "Choice 2 Con"|"Choice 3 Con"|"Choice 4 Con"|"Choice 5 Con"|"Choice 6 Con"|"Choice 7 Con"|"Choice 8 Con"|"Choice 9 Con"|"Choice 10 Con"|varsz;
    endif;

*/

    _cml_parnames 	= varsz;  

endif;

{ logit,maxflgt,glgt,covlgt,retlgt } = cmlprt(cml(dset,0,&lgt,startmnl));

rho_me      = 1;
theta_me    = zeros(sumc(thetas),1);
b_me        = logit;


if ismiss(logit) eq 1;
	logit 	= ones(rows(startmnl),1).*-999;
	lgtse 	= logit;
	tlogit 	= logit;
	maxflgt 	= -999;
else;
	lgtse 	= sqrt(diag(covlgt));
	lgtse 	= missrv(lgtse,10000);
	if rows(lgtse) eq 1;
		lgtse = lgtse .* ones(rows(startmnl),1);
	endif;
	tlogit 	= abs(logit ./ lgtse);
	maxflgt	= maxflgt*n;
endif;

if sumc(res') ge 1 and sumc(res') lt k;

	beta1 	= logit[1:k];
	lgtse1	= lgtse[1:k];
	lgtt1	= tlogit[1:k];

	betanr	= reshape(logit[k+1:noxs],capJ-2,k-sumc(res'))';
	lgtsenr	= reshape(lgtse[k+1:noxs],capJ-2,k-sumc(res'))';
	lgttnr	= reshape(tlogit[k+1:noxs],capJ-2,k-sumc(res'))';

	beta	= betanr[1,.];
	se		= lgtsenr[1,.];	
	t		= lgttnr[1,.];	

	b_rep 	= 2;
	row_rep = 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta	= beta|betanr[row_rep,.];
			se		= se|lgtsenr[row_rep,.];
			t		= t|lgttnr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta	= beta|(zeros(1,capJ-2) + beta1[b_rep]);
			se		= se|(zeros(1,capJ-2) + lgtse1[b_rep]);
			t		= t|(zeros(1,capJ-2) + lgtt1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	logit 	= vec(beta1~beta);
	lgtse 	= vec(lgtse1~se);	
	tlogit 	= vec(lgtt1~t);    

endif;

if retlgt eq 0;
	let retlgt = normal convergence;
elseif retlgt eq 1;
	let retlgt = forced exit;
elseif retlgt eq 2;
	let retlgt = max iters exceeded;
elseif retlgt eq 3;
	let retlgt = function failed;
elseif retlgt eq 4;
	let retlgt = gradient failed;
elseif retlgt eq 5;
	let retlgt = Hessian failed;
elseif retlgt eq 6;
	let retlgt = line search failed;
elseif retlgt eq 7;
	let retlgt = bad starting values;
elseif retlgt eq 8;
	let retlgt = error with gradient;
elseif retlgt eq 9;
	let retlgt = error with constraints;
elseif retlgt eq 10;
	let retlgt = secant update failed;
elseif retlgt eq 11;
	let retlgt = max time exceeded;
elseif retlgt eq 12;
	let retlgt = error with weights;
elseif retlgt eq 13;
	let retlgt = quadratic program failed;
elseif retlgt eq 20;
	let retlgt = Hessian failed to invert;
elseif retlgt eq 34;
	let retlgt = couldn't open data set;
else;
	let retlgt = unknown;
endif;


			/* CONSTRUCT PROBABLITIES */
 
betahat	= reshape(logit,capJ-1,k)';

v_ij		= exp(z*betahat); 
plgt		= (1/(1 + sumc(v_ij'))) ~ (v_ij ./ (1 + sumc(v_ij')));

check = meanc(sumc(plgt'));
if check gt 1.00001 or check lt (1-0.00001);
	errorlog "MNL PROBABILITIES DON'T SUM TO 1 - AVERAGE PROB IS";check;stop;
endif;


        /* CONSTRUCT MARGINAL EFFECTS AT ZBAR */

zbar = meanc(z);

save capJ, k;

mnl_mes     = zeros(rows(zbar),1);

j = 1;
do while j le capJ;
    g       = gradp(&dogev_p,zbar);
    mnl_mes = mnl_mes~g';
    j=j+1;
endo;

mnl_mes = mnl_mes[.,2:cols(mnl_mes)];

if sumc(sumc(mnl_mes')) gt 0.00001;
    "ERROR: MARGINAL EFFECTS SUM TO";;sumc(sumc(mnl_mes'));
    stop;
endif;


    /* SIMULATED MARGINAL ERRORS FOR ME'S */

if ismiss(covlgt) eq 1;

    mnl_mese = ones(k,capJ) .*-999;

else;

    for i (1,M,1);

        gam_sim = ((rndn(1,cols(covlgt)) * chol(covlgt)) + (logit' .* ones(1,1)))';
        b_me    = gam_sim;


	    if rows(b_me) le 1;
		    "ONLY ONE PARAMETER";stop;
	    endif;

        mestmp     = zeros(rows(zbar),1);

        j = 1;

        do while j le capJ;        
            g       = gradp(&dogev_p,zbar);
            mestmp  = mestmp~g';
            j=j+1;
        endo;

        mestmp = mestmp[.,2:cols(mestmp)];

        if i eq 1;
            mesall = vec(mestmp)';
        else;
            mesall = mesall|vec(mestmp)'; 
        endif;

    endfor;

    mnl_mese = stdc(mesall);
    mnl_mese = reshape(mnl_mese,capJ,rows(mnl_mes))';

endif;

if mtype eq 1;
	clear dogit, dgtse, tdogit, maxfdgt, ogev, ogevse, togev, maxfogv, 
		dogev, dgvse, tdgev, maxfdgv, ordpbt, opbse, topb, maxfopb,
		pdgt, pogv, pdgv, pordp,
		covogev, covdgev, covopb, covdgt, ogv_mes, op_mes, dgt_mes, dgv_mes,
        ogv_mese, op_mese, dgt_mese, dgv_mese;
	let retdgt = not estimated;
	let retogv = not estimated;
	let retdgv = not estimated;
	let retopb = not estimated;
	goto finish;
elseif mtype eq 2 ;
 	clear dogit, dgtse, tdogit, maxfdgt, ogev, ogevse, togev, maxfogv, 
		dogev, dgvse, tdgev, maxfdgv, ordpbt, opbse, topb, maxfopb,
		pogv, pdgv, pordp,
		covogev, covdgev, covopb, ogv_mes, op_mes, dgt_mes, dgv_mes,
        ogv_mese, op_mese, dgt_mese, dgv_mese;
	let retogv = not estimated;
	let retdgv = not estimated;
	let retopb = not estimated;
elseif mtype eq 3 ;
 	clear dogit, dgtse, tdogit, maxfdgt, ogev, ogevse, togev, maxfogv, 
		dogev, dgvse, tdgev, maxfdgv, ordpbt, opbse, topb, maxfopb,
		pdgt, pogv, pdgv, pordp,
		covdgev, covopb, covdgt, ogv_mes, op_mes, dgt_mes, dgv_mes,
        ogv_mese, op_mese, dgt_mese, dgv_mese;
	let retdgt = not estimated;
	let retdgv = not estimated;
	let retopb = not estimated;
	goto ogvest;
elseif mtype eq 5 ;
 	clear dogit, dgtse, tdogit, maxfdgt, ogev, ogevse, togev, maxfogv, 
		dogev, dgvse, tdgev, maxfdgv, ordpbt, opbse, topb, maxfopb,
		pdgt, pogv, pdgv, pordp,
		covogev, covdgev, covdgt, ogv_mes, op_mes, dgt_mes, dgv_mes,
        ogv_mese, op_mese, dgt_mese, dgv_mese;
	let retdgt = not estimated;
	let retogv = not estimated;
	let retdgv = not estimated;
	goto opbest;
else;
	clear dogev, dgvse, tdgev, maxfdgv, pdgv,covdgev; 
endif;


			/* ESTIMATE DOGIT MODEL */

dgtest:

let theta 			= theta;
cmlset;

CC				= zeros(sumc(thetas),noxs)~eye(sumc(thetas));
DD				= (0.000001 .* ones(sumc(thetas),1));
_cml_C			= CC;
_cml_D			= DD;
_cml_dirTol 	= 0.001;
_cml_maxIters	= 100;
__output 		= 1;
__title			= "DOGIT Model";

/*
dgtbndsa 		= ((ones(rows(logit),1) .* -75)~(ones(rows(logit),1) .* 75));
dgtbndsb 		= (0.000001 .* ones(sumc(thetas),1))~(100 .* ones(sumc(thetas),1));
dgtbnds 			= dgtbndsa|dgtbndsb;
_cml_bounds		= dgtbnds; 
*/

if sumc(res') ge 1;
	_cml_parnames 	= varsz|(ones(capJ-2,1) .*. delif(varsz,res'))|(theta .*. ones(sumc(thetas),1));  
else;
	_cml_parnames 	= (ones(capJ-1,1) .*. varsz)|(theta .*. ones(sumc(thetas),1));  
endif;
	
_cml_covpar 	= 1; 
_cml_algorithm	= 3; 		/* 1 = BFGS, 2=DFP, 3= NEWTON, 4=BHHH */
_cml_IneqProc	= &ineqp; 		/* sets upper bounds for thetas */
/* _cml_EqProc		= &ineqp; */

load phat;
thetast = (selif(phat,thetas)) ./ 2;

if rows(startdgt) eq 1 and sumc(res') eq 0;
	startdgt=logit|thetast;
elseif rows(startdgt) eq 1 and sumc(res') ge 1;
	startdgt=logit|thetast;
endif;

{ dogit,maxfdgt,gdgt,covdgt,retdgt } = cmlprt(cml(dset,0,&dgt,startdgt));


if ismiss(dogit) eq 1;
	dogit	= ones(rows(startdgt),1).*-999;
	dgtse	= dogit;
	tdogit	= dogit;
	maxfdgt	= -999;
else;
	dgtse 	= sqrt(abs(diag(covdgt)));
	dgtse 	= missrv(dgtse,10000) ;
	if rows(dgtse) eq 1;
		dgtse = dgtse .* ones(rows(startdgt),1);
	endif;
	tdogit 	= abs(dogit ./ dgtse);
	maxfdgt	= maxfdgt*n;
endif;


			/* PULL OUT APPROPRIATE PARAMETERS TAKING INTO ACCOUNT RESTRICTIONS */

thetahat    = (dogit[rows(dogit)-sumc(thetas)+1:rows(dogit)]);  
betahat	    = dogit[1:rows(dogit)-sumc(thetas)]; 
sebeta	    = dgtse[1:rows(dogit)-sumc(thetas)]; 
tbeta	    = tdogit[1:rows(dogit)-sumc(thetas)]; 

theta_me    = thetahat;
b_me        = betahat;
rho_me      = 1;


if rows(dogit) ne (rows(thetahat)+rows(betahat));
	stop;
endif;

if sumc(res') eq 0; 

	betahat	= reshape(betahat,capJ-1,k)';

elseif sumc(res') ge 1;

	beta_1	= betahat[1:k];
	dgtse_1	= sebeta[1:k];
	tdgt_1	= tbeta[1:k];

	betanr	= reshape(betahat[k+1:noxs],capJ-2,k-sumc(res'))';
	dgtsenr	= reshape(sebeta[k+1:noxs],capJ-2,k-sumc(res'))';
	tdgtnr	= reshape(tbeta[k+1:noxs],capJ-2,k-sumc(res'))';

	beta	= betanr[1,.]; 
	se		= dgtsenr[1,.]; 
	t		= tdgtnr[1,.]; 

	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			se		= se|dgtsenr[row_rep,.];
			t		= t|tdgtnr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
			se		= se|(zeros(1,capJ-2) + dgtse_1[b_rep]);
			t		= t|(zeros(1,capJ-2) + tdgt_1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	betahat = beta_1~beta;
	se		= dgtse_1~se;
	t		= tdgt_1~t;

	dogit	= vec(betahat)|thetahat;
	dgtse	= vec(se)|(dgtse[noxs+1:rows(dgtse)]);  
	tdogit	= vec(t)|(tdogit[noxs+1:rows(tdogit)]);  

endif;


			/* RETURN CODES */

if retdgt eq 0;
	let retdgt = normal convergence;
elseif retdgt eq 1;
	let retdgt = forced exit;
elseif retdgt eq 2;
	let retdgt = max iters exceeded;
elseif retdgt eq 3;
	let retdgt = function failed;
elseif retdgt eq 4;
	let retdgt = gradient failed;
elseif retdgt eq 5;
	let retdgt = Hessian failed;
elseif retdgt eq 6;
	let retdgt = line search failed;
elseif retdgt eq 7;
	let retdgt = bad starting values;
elseif retdgt eq 8;
	let retdgt = error with gradient;
elseif retdgt eq 9;
	let retdgt = error with constraints;
elseif retdgt eq 10;
	let retdgt = secant update failed;
elseif retdgt eq 11;
	let retdgt = max time exceeded;
elseif retdgt eq 12;
	let retdgt = error with weights;
elseif retdgt eq 13;
	let retdgt = quadratic program failed;
elseif retdgt eq 20;
	let retdgt = Hessian failed to invert;
elseif retdgt eq 34;
	let retdgt = couldn't open data set;
else;
	let retdgt = unknown;
endif;


			/* CONSTRUCT DOGIT PROBABLITIES */

betahat	= zeros(k,1)~betahat;
rho		= 1;

theta	=0;
jj		=1;
i		=1;
do while i le capJ;
	if thetas[i] eq 1;
		theta	= theta|thetahat[jj];
		jj		= jj+1;
	else;
		theta	= theta|0;
	endif;
i=i+1;
endo;

theta=theta[2:rows(theta)];

v_ij		= exp((z*betahat) ./ rho); 
v_ijall	= zeros(n,1)~v_ij~zeros(n,1);

p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]) .^rho;
p_a		= (v_ij ./ sumc(p_denom'));  

p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho-1)) + 
			((v_ij + v_ijall[.,3:capJ+2]) .^ (rho-1)); 
			
p_ogev	= p_a .* p_b;

pdgt = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

check = meanc(sumc(pdgt'));
if check gt 1.00001 or check lt (1-0.00001);
	errorlog "DOGIT PROBABILITIES DON'T SUM TO 1: AVERAGE PROB IS";check;stop;
endif;


			/* CHECK FOR UPPER BOUND CONDITION */

load phat;
phat		= selif(phat,thetas);
dgtchk 	= (thetahat ./ (1+sumc(thetahat))) .le (phat + 0.00001);

if sumc(dgtchk) ne rows(dgtchk);
	print;"DOGIT UPPER BOUNDS CONDITION VIOLATED";
	print;(thetahat ./ (1+sumc(thetahat)))~phat;
	stop;
endif;


        /* CONSTRUCT MARGINAL EFFECTS AT ZBAR */

dgt_mes     = zeros(rows(zbar),1);

j = 1;
do while j le capJ;
    g       = gradp(&dogev_p,zbar);
    dgt_mes = dgt_mes~g';
    j=j+1;
endo;

dgt_mes = dgt_mes[.,2:cols(dgt_mes)];

if sumc(sumc(dgt_mes')) gt 0.00001;
    "ERROR: MARGINAL EFFECTS SUM TO";;sumc(sumc(dgt_mes'));
    stop;
endif;


    /* SIMULATED MARGINAL ERRORS FOR ME'S */


if ismiss(covdgt) eq 1;

    dgt_mese = ones(k,capJ) .*-999;

else;

    for i (1,M,1);

        gam_sim     = ((rndn(1,cols(covdgt)) * chol(covdgt)) + (dogit' .* ones(1,1)))';
        theta_me    = (gam_sim[rows(gam_sim)-sumc(thetas)+1:rows(gam_sim)]);  
        b_me        = gam_sim[1:rows(gam_sim)-sumc(thetas)]; 

	    if rows(b_me) le 1;
		    "ONLY ONE PARAMETER";stop;
	    endif;

        mestmp     = zeros(rows(zbar),1);

        j = 1;

        do while j le capJ;
            
            g       = gradp(&dogev_p,zbar);
            mestmp  = mestmp~g';
            j=j+1;
        endo;

        mestmp = mestmp[.,2:cols(mestmp)];

        if i eq 1;
            mesall = vec(mestmp)';
        else;
            mesall = mesall|vec(mestmp)'; 
        endif;

    endfor;

    dgt_mese = stdc(mesall);
    dgt_mese = reshape(dgt_mese,capJ,rows(dgt_mes))';

endif;

if mtype eq 2; 
	goto finish;
endif;


			/* ESTIMATE OGEV MODEL */

ogvest:

cmlset;
_cml_dirTol 	= 0.0001;
_cml_maxIters	= 100;        
_cml_maxIters	= 2;        
__output 		= 1;
__title			= "OGEV Model";

CC				= zeros(2,noxs)~(1|-1);
DD				= 0.01|-1;
_cml_C			= CC;
_cml_D			= DD;

/*
_cml_bounds		= ((ones(rows(logit),1)*-1000)~(ones(rows(logit),1)*1000))|(0.05~1);
*/

if sumc(res') ge 1;
	_cml_parnames 	= varsz|(ones(capJ-2,1) .*. delif(varsz,res'))|"RHO";  
else;
	_cml_parnames 	=(ones(capJ-1,1) .*. varsz)|"RHO"; 
endif;

_cml_covpar 	= 1;   /* 1,2,3, 1 = default others have 3  */
_cml_algorithm	= 1; 		/* 1 = BFGS, 2=DFP, 3= NEWTON, 4=BHHH */

if rows(startogv) eq 1 and sumc(res') eq 0;
	startogv=logit|0.5; 
elseif rows(startogv) eq 1 and sumc(res') ge 1;
	startogv=logit|0.5; 
endif;


{ ogev,maxfogv,gogev,covogev,retogv } = cmlprt(cml(dset,0,&ogv,startogv));

if ismiss(ogev) eq 1;
	ogev 	= ones(rows(startogv),1).*-999;
	ogevse	= ogev;
	togev	= ogev;
	maxfogv	= -999;
else;
	ogevse 	= sqrt(abs(diag(covogev)));
	ogevse 	= missrv(ogevse,10000) .* ones(rows(startogv),1);
	if rows(ogevse) eq 1;
		ogevse 	= ogevse .* ones(rows(startogv),1);
	endif;
	togev 	= abs(ogev ./ ogevse);
	maxfogv	= maxfogv*n;
endif;


			/* PULL OUT APPROPRIATE PARAMETERS TAKING INTO ACCOUNT RESTRICTIONS */

rho 	    = ogev[rows(ogev)]; 
betahat	    = ogev[1:rows(ogev)-1];
sebeta	    = ogevse[1:rows(ogev)-1]; 
tbeta	    = togev[1:rows(ogev)-1]; 

rho_me      = rho;
b_me        = betahat;
theta_me    = zeros(sumc(thetas),1);


if rows(ogev) ne (rows(rho)+rows(betahat));
	stop;
endif;


if sumc(res') eq 0; 

	betahat	= reshape(betahat,capJ-1,k)';

elseif sumc(res') ge 1;

	beta_1	= betahat[1:k];
	ogvse_1	= sebeta[1:k];
	togv_1	= tbeta[1:k];

	betanr	= reshape(betahat[k+1:noxs],capJ-2,k-sumc(res'))';
	ogvsenr	= reshape(sebeta[k+1:noxs],capJ-2,k-sumc(res'))';
	togvnr	= reshape(tbeta[k+1:noxs],capJ-2,k-sumc(res'))';

	beta		= betanr[1,.]; 
	se		= ogvsenr[1,.]; 
	t		= togvnr[1,.]; 

	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			se		= se|ogvsenr[row_rep,.];
			t		= t|togvnr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
			se		= se|(zeros(1,capJ-2) + ogvse_1[b_rep]);
			t		= t|(zeros(1,capJ-2) + togv_1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	betahat 	= beta_1~beta;
	se		= ogvse_1~se;
	t		= togv_1~t;

	ogev	= vec(betahat)|rho;
	ogevse	= vec(se)|(ogevse[rows(ogevse)]);  
	togev	= vec(t)|(togev[rows(togev)]);  

endif;


			/* RETURN CODES */

if retogv eq 0;
	let retogv = normal convergence;
elseif retogv eq 1;
	let retogv = forced exit;
elseif retogv eq 2;
	let retogv = max iters exceeded;
elseif retogv eq 3;
	let retogv = function failed;
elseif retogv eq 4;
	let retogv = gradient failed;
elseif retogv eq 5;
	let retogv = Hessian failed;
elseif retogv eq 6;
	let retogv = line search failed;
elseif retogv eq 7;
	let retogv = bad starting values;
elseif retogv eq 8;
	let retogv = error with gradient;
elseif retogv eq 9;
	let retogv = error with constraints;
elseif retogv eq 10;
	let retogv = secant update failed;
elseif retogv eq 11;
	let retogv = max time exceeded;
elseif retogv eq 12;
	let retogv = error with weights;
elseif retogv eq 13;
	let retogv = quadratic program failed;
elseif retogv eq 20;
	let retogv = Hessian failed to invert;
elseif retogv eq 34;
	let retogv = couldn't open data set;
else;
	let retogv = unknown;
endif;

			/* CONSTRUCT OGEV PROBABLITIES */

betahat	= zeros(k,1)~betahat;
theta	= zeros(capJ,1);

v_ij		= exp((z*betahat) ./ rho); 
v_ijall	= zeros(n,1)~v_ij~zeros(n,1);

p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]).^rho;

p_a		= (v_ij ./ sumc(p_denom'));  

p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho-1)) + 
			((v_ij + v_ijall[.,3:capJ+2]) .^ (rho-1)); 
			
p_ogev	= p_a .* p_b;

pogv = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

check = meanc(sumc(pogv'));
if check gt 1.00001 or check lt (1-0.00001);
	errorlog "OGEV PROBABILITIES DON'T SUM TO 1: AVERAGE PROB IS";check;stop;
endif;


        /* CONSTRUCT MARGINAL EFFECTS AT ZBAR */

ogv_mes     = zeros(rows(zbar),1);

j = 1;
do while j le capJ;
    g       = gradp(&dogev_p,zbar);
    ogv_mes = ogv_mes~g';
    j=j+1;
endo;

ogv_mes = ogv_mes[.,2:cols(ogv_mes)];

if sumc(sumc(ogv_mes')) gt 0.00001;
    "ERROR: MARGINAL EFFECTS SUM TO";;sumc(sumc(ogv_mes'));
    stop;
endif;


    /* SIMULATED MARGINAL ERRORS FOR ME'S */

if ismiss(covogev) eq 1;

    ogv_mese = ones(k,capJ) .*-999;

else;

    for i (1,M,1);

        gam_sim = ((rndn(1,cols(covogev)) * chol(covogev)) + (ogev' .* ones(1,1)))';
        rho_me  = gam_sim[rows(gam_sim)]; 
        b_me	= gam_sim[1:rows(gam_sim)-1];

	    if rows(b_me) le 1;
		    "ONLY ONE PARAMETER";stop;
	    endif;

        mestmp     = zeros(rows(zbar),1);

        j = 1;
        do while j le capJ;           
            g       = gradp(&dogev_p,zbar);
            mestmp  = mestmp~g';
            j=j+1;
        endo;

        mestmp = mestmp[.,2:cols(mestmp)];

        if i eq 1;
            mesall = vec(mestmp)';
        else;
            mesall = mesall|vec(mestmp)'; 
        endif;

    endfor;

    ogv_mese = stdc(mesall);
    ogv_mese = reshape(ogv_mese,capJ,rows(ogv_mes))';

endif;


if mtype eq 3; 
	goto finish;
endif;


			/* ESTIMATE DOGEV MODEL */

dgvest:
let theta 			= theta;

cmlset;


CC				= (zeros(sumc(thetas)+1,noxs)~eye(sumc(thetas)+1))	|
				   (zeros(1,noxs+sumc(thetas))~-1);
DD				= (0.000001 .* ones(sumc(thetas),1)) | 0.01 | -1;
_cml_C			= CC;
_cml_D			= DD;


_cml_dirTol 	= 0.001;
_cml_maxIters 	= 100;
__output 		= 1;
__title			= "DOGEV Model";

if sumc(res') ge 1;
	_cml_parnames 	= varsz|(ones(capJ-2,1) .*. delif(varsz,res'))|(theta .*. ones(sumc(thetas),1))|"RHO"; 
else;
	_cml_parnames 	= (ones(capJ-1,1) .*. varsz)|(theta .*. ones(sumc(thetas),1))|"RHO";
endif;
 
_cml_covpar 	= 1; 
_cml_algorithm	= 1; 			/* 1 = BFGS, 2=DFP, 3= NEWTON, 4=BHHH */
_cml_IneqProc	= &ineqp; 		/* sets upper bounds for thetas and rho */


if rows(startdgv) eq 1;
	if ogev[rows(ogev)] ge 0.95 or ogev[rows(ogev)] lt 0.06;
			if sumc(res') eq 0;
				startdgv = logit|thetast|0.5;
			elseif sumc(res') ge 1;
				startdgv = logit|thetast|0.5;
			endif;
	else;
		if sumc(res') eq 0;
			startdgv = logit|thetast|ogev[rows(ogev)];
		elseif sumc(res') ge 1;
			startdgv = logit|thetast|ogev[rows(ogev)];
		endif;
	endif;
elseif rows(startdgv) ge 2;
	if startdgv[rows(startdgv)] ge 0.95 or startdgv[rows(startdgv)] lt 0.06;
		startdgv[rows(startdgv)] = 0.5;
	endif;
endif;


{ dogev,maxfdgv,gdogev,covdgev,retdgv } = cmlprt(cml(dset,0,&dgv,startdgv));


if ismiss(dogev) eq 1;
	dogev	= ones(rows(startdgv),1).*-999;
	dgvse	= dogev;
	tdgev	= dogev;
	maxfdgv	= -999;
else;
	dgvse 	= sqrt(abs(diag(covdgev)));
	dgvse 	= missrv(dgvse,10000);
	if rows(dgvse) eq 1;
		dgvse 	= dgvse .* ones(rows(startdgv),1);
	endif;
	tdgev 	= abs(dogev ./ dgvse);
	maxfdgv	= maxfdgv*n;
endif;


			/* PULL OUT APPROPRIATE PARAMETERS TAKING INTO ACCOUNT RESTRICTIONS */

betahat	    = dogev[1:rows(dogev)-sumc(thetas)-1]; 
rho 	    = dogev[rows(dogev)]; 
thetahat	= (dogev[rows(dogev)-sumc(thetas):rows(dogev)-1]); 
sebeta	    = dgvse[1:rows(dogev)-sumc(thetas)-1]; 
tbeta	    = tdgev[1:rows(dogev)-sumc(thetas)-1]; 

b_me        = betahat;
theta_me    = thetahat;
rho_me      = rho;


if rows(dogev) ne (rows(thetahat)+rows(rho)+rows(betahat));
	stop;
endif;


if sumc(res') eq 0; 

	betahat	= reshape(betahat,capJ-1,k)';
	
elseif sumc(res') ge 1;

	beta_1	= betahat[1:k];
	dgvse_1	= sebeta[1:k];
	tdgv_1	= tbeta[1:k];

	betanr	= reshape(betahat[k+1:noxs],capJ-2,k-sumc(res'))';
	dgvsenr	= reshape(sebeta[k+1:noxs],capJ-2,k-sumc(res'))';
	tdgvnr	= reshape(tbeta[k+1:noxs],capJ-2,k-sumc(res'))';

	beta		= betanr[1,.]; 
	se		= dgvsenr[1,.]; 
	t		= tdgvnr[1,.]; 

	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			se		= se|dgvsenr[row_rep,.];
			t		= t|tdgvnr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
			se		= se|(zeros(1,capJ-2) + dgvse_1[b_rep]);
			t		= t|(zeros(1,capJ-2) + tdgv_1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	betahat 	= beta_1~beta;
	se		= dgvse_1~se;
	t		= tdgv_1~t;

	dogev	= vec(betahat)|thetahat|rho;
	dgvse	= vec(se)|(dgvse[noxs+1:rows(dgvse)]);  
	tdgev	= vec(t)|(tdgev[noxs+1:rows(tdgev)]);  

endif;


			/* RETURN OUTPUTS */

if retdgv eq 0;
	let retdgv = normal convergence;
elseif retdgv eq 1;
	let retdgv = forced exit;
elseif retdgv eq 2;
	let retdgv = max iters exceeded;
elseif retdgv eq 3;
	let retdgv = function failed;
elseif retdgv eq 4;
	let retdgv = gradient failed;
elseif retdgv eq 5;
	let retdgv = Hessian failed;
elseif retdgv eq 6;
	let retdgv = line search failed;
elseif retdgv eq 7;
	let retdgv = bad starting values;
elseif retdgv eq 8;
	let retdgv = error with gradient;
elseif retdgv eq 9;
	let retdgv = error with constraints;
elseif retdgv eq 10;
	let retdgv = secant update failed;
elseif retdgv eq 11;
	let retdgv = max time exceeded;
elseif retdgv eq 12;
	let retdgv = error with weights;
elseif retdgv eq 13;
	let retdgv = quadratic program failed;
elseif retdgv eq 20;
	let retdgv = Hessian failed to invert;
elseif retdgv eq 34;
	let retdgv = couldn't open data set;
else;
	let retdgv = unknown;
endif;


			/* CONSTRUCT DOGEV PROBABLITIES */

betahat	= zeros(k,1)~betahat;

theta	=0;
jj		=1;
i		=1;
do while i le capJ;
	if thetas[i] eq 1;
		theta	= theta|thetahat[jj];
		jj		= jj+1;
	else;
		theta	= theta|0;
	endif;
i=i+1;
endo;

theta=theta[2:rows(theta)];

v_ij		= exp((z*betahat) ./ rho); 
v_ijall	= zeros(n,1)~v_ij~zeros(n,1);

p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]).^rho;

p_a		= (v_ij ./ sumc(p_denom'));  

p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho-1)) + 
			((v_ij + v_ijall[.,3:capJ+2]) .^ (rho-1)); 
			
p_ogev	= p_a .* p_b;

pdgv = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

check = meanc(sumc(pdgv'));
if check gt 1.00001 or check lt (1-0.00001);
	errorlog "DOGEV PROBABILITIES DON'T SUM TO 1: AVERAGE PROB IS";check;stop;
endif;


			/* CHECK FOR UPPER BOUND CONDITION */

load phat;
phat		= selif(phat,thetas);
dgvchk = (thetahat ./ (1+sumc(thetahat))) .le phat;

if sumc(dgvchk) ne rows(dgvchk);
	print;"DOGEV UPPER BOUNDS CONDITION VIOLATED";
	print;(thetahat ./ (1+sumc(thetahat)))~phat;
	stop;
endif;


        /* CONSTRUCT MARGINAL EFFECTS AT ZBAR */

dgv_mes     = zeros(rows(zbar),1);

j = 1;
do while j le capJ;
    g       = gradp(&dogev_p,zbar);
    dgv_mes = dgv_mes~g';
    j=j+1;
endo;

dgv_mes = dgv_mes[.,2:cols(dgv_mes)];

if sumc(sumc(dgv_mes')) gt 0.00001;
    "ERROR: MARGINAL EFFECTS SUM TO";;sumc(sumc(dgv_mes'));
    stop;
endif;


    /* SIMULATED MARGINAL ERRORS FOR ME'S */

if ismiss(covdgev) eq 1;

    dgv_mese = ones(k,capJ) .*-999;

else;


    for i (1,M,1);

        gam_sim     = ((rndn(1,cols(covdgev)) * chol(covdgev)) + (dogev' .* ones(1,1)))';
        b_me	    = gam_sim[1:rows(gam_sim)-sumc(thetas)-1]; 
        rho_me      = gam_sim[rows(gam_sim)]; 
        theta_me    = (gam_sim[rows(gam_sim)-sumc(thetas):rows(gam_sim)-1]); 

	    if rows(b_me) le 1;
		    "ONLY ONE PARAMETER";stop;
	    endif;

        mestmp     = zeros(rows(zbar),1);

        j = 1;

        do while j le capJ;    
            g       = gradp(&dogev_p,zbar);
            mestmp  = mestmp~g';
            j=j+1;
        endo;

        mestmp = mestmp[.,2:cols(mestmp)];

        if i eq 1;
            mesall = vec(mestmp)';
        else;
            mesall = mesall|vec(mestmp)'; 
        endif;

    endfor;

    dgv_mese = stdc(mesall);
    dgv_mese = reshape(dgv_mese,capJ,rows(dgv_mes))';

endif;

if mtype eq 2; 
	goto finish;
endif;

opbest:


			/* ESTIMATE ORDERED PROBIT MODEL */

if capj eq 3;
	clear ordpbt,opbse,topb,maxfopb,pordp, retopb;
	clear ordpbt,maxfopb,gopb,covopb,retopb,opbse,ordpbt,topb,maxfopb,covopb,op_mes, op_mese, pordp, retopb;		
	goto finish;
endif;
	
cmlset;


if rows(startopb) eq 1;
	startmu1			= meanc(d_ij);
	startmu 			= cdfni(startmu1[2]+startmu1[1]);

	startrep			= 3;
	do while startrep le capJ-1;
		mark	= startmu1[startrep] + cdfn(startmu[rows(startmu)]);
		startmu 	= startmu|cdfni(mark);
	startrep=startrep+1;
	endo;

	startmu = cdfn(startmu); 

	startopb 			= (inv(z'z)*z'(y .ge 1))|startmu;
endif;

let mu			= mu;
c				= ((eye(capJ-3) .*-1)~zeros(capJ-3,1)) + (zeros(capJ-3,1)~eye(capJ-3));
_cml_c			= zeros(capJ-3,cols(z))~c;
_cml_d			= zeros(rows(_cml_c),1);
_cml_dirTol 	= 0.0001;
_cml_maxIters 	= 100;
__output 		= 1;
_cml_parnames 	= varsz|(mu .*. ones(capJ-2,1));  
_cml_covpar 	= 3;
_cml_algorithm	= 3; 		/* 1 = BFGS, 2=DFP, 3= NEWTON, 4=BHHH */
__title 			= "Ordered Probit";

{ ordpbt,maxfopb,gopb,covopb,retopb } = cmlprt(cml(dset,0,&opb,startopb));

if ismiss(ordpbt) eq 1;
	ordpbt	= ones(rows(startopb),1).*-999;
	opbse	= ordpbt;
	topb		= ordpbt;
	maxfopb	= -999;
else;
	opbse 	= sqrt(diag(covopb));
	opbse 	= missrv(opbse,10000);
	topb 	= abs(ordpbt ./ opbse);
	maxfopb	= maxfopb*n;
endif;

if retopb eq 0;
	let retopb = normal convergence;
elseif retopb eq 1;
	let retopb = forced exit;
elseif retopb eq 2;
	let retopb = max iters exceeded;
elseif retopb eq 3;
	let retopb = function failed;
elseif retopb eq 4;
	let retopb = gradient failed;
elseif retopb eq 5;
	let retopb = Hessian failed;
elseif retopb eq 6;
	let retopb = line search failed;
elseif retopb eq 7;
	let retopb = bad starting values;
elseif retopb eq 8;
	let retopb = error with gradient;
elseif retopb eq 9;
	let retopb = error with constraints;
elseif retopb eq 10;
	let retopb = secant update failed;
elseif retopb eq 11;
	let retopb = max time exceeded;
elseif retopb eq 12;
	let retopb = error with weights;
elseif retopb eq 13;
	let retopb = quadratic program failed;
elseif retopb eq 20;
	let retopb = Hessian failed to invert;
elseif retopb eq 34;
	let retopb = couldn't open data set;
else;
	let retopb = unknown;
endif;


			/* CONSTRUCT ORDERED PROBIT PROBABLITIES */

betahat	= ordpbt[1:k];
muhat	= ordpbt[k+1:rows(ordpbt)]; 

b_me    = betahat;
mu_me   = muhat;

if rows(ordpbt) ne (rows(betahat)+rows(muhat));
	stop;
endif;


			/* DETERMINISTIC PART & PROBABILITIES */

zbetahat = z*betahat;

p_ij	= 1-(cdfn(zbetahat));
p_ij	= p_ij~(cdfn(muhat[1]-zbetahat) - cdfn(-zbetahat));
jrep2= 2;
do while jrep2 le capJ-2;
	p_ij = p_ij~(cdfn(muhat[jrep2]-zbetahat) -
          	            cdfn(muhat[jrep2-1]-zbetahat));
jrep2 = jrep2 + 1;
endo;

pordp = p_ij~(1 - cdfn(muhat[capJ-2]-zbetahat));

check = meanc(sumc(pordp'));
if check gt 1.00001 or check lt (1-0.00001);
	errorlog "ORDERED PROBIT PROBABILITIES DON'T SUM TO 1: AVERAGE PROB IS";check;stop;
endif;


        /* CONSTRUCT MARGINAL EFFECTS AT ZBAR */

op_mes  = zeros(rows(zbar),1);

j = 1;
do while j le capJ;
    g       = gradp(&ordp_p,zbar);
    op_mes  = op_mes~g';
    j=j+1;
endo;

op_mes = op_mes[.,2:cols(op_mes)];

if sumc(sumc(op_mes')) gt 0.00001;
    "ERROR: MARGINAL EFFECTS SUM TO";;sumc(sumc(op_mes'));
    stop;
endif;


    /* SIMULATED MARGINAL ERRORS FOR ME'S */

if ismiss(covopb) eq 1;

    op_mese = ones(k,capJ) .*-999;

else;


    for i (1,M,1);

        gam_sim = ((rndn(1,cols(covopb)) * chol(covopb)) + (ordpbt' .* ones(1,1)))';
        b_me	= gam_sim[1:k];
        mu_me	= gam_sim[k+1:rows(gam_sim)]; 

	    if rows(b_me) le 1;
		    "ONLY ONE PARAMETER";stop;
	    endif;

        mestmp     = zeros(rows(zbar),1);

        j = 1;
        do while j le capJ;    
            g       = gradp(&ordp_p,zbar);
            mestmp  = mestmp~g';
            j=j+1;
        endo;

        mestmp = mestmp[.,2:cols(mestmp)];

        if i eq 1;
            mesall = vec(mestmp)';
        else;
            mesall = mesall|vec(mestmp)'; 
        endif;

    endfor;

    op_mese = stdc(mesall);
    op_mese = reshape(op_mese,capJ,rows(op_mes))';

endif;

finish:

/*
startmnl;
logit;
lgtse;
tlogit;
mnl_mes;
mnl_mese;
covlgt;
maxflgt;
plgt;
retlgt;
dogit;
dgtse;
tdogit;
dgt_mes;
dgt_mese;
covdgt;
maxfdgt;
pdgt;
retdgt;
ogev;
ogevse;
togev;
ogv_mes;
ogv_mese;
covogev;
maxfogv;
pogv;
retogv;
dogev;
dgvse;
tdgev;
dgv_mes;
dgv_mese;
covdgev;
maxfdgv;
pdgv;
retdgv;
ordpbt;
opbse;
topb;
op_mes; 
op_mese;
covopb;
maxfopb;
pordp;
retopb;
*/

			/* RETURN OUTPUTS */

retp(startmnl, 
	logit~lgtse~tlogit, mnl_mes, mnl_mese, covlgt,  maxflgt, plgt,  retlgt,
	dogit~dgtse~tdogit, dgt_mes, dgt_mese, covdgt,  maxfdgt, pdgt,  retdgt,
	ogev~ogevse~togev,  ogv_mes, ogv_mese, covogev, maxfogv, pogv,  retogv,
	dogev~dgvse~tdgev,  dgv_mes, dgv_mese, covdgev, maxfdgv, pdgv,  retdgv,
	ordpbt~opbse~topb,  op_mes,  op_mese,  covopb,  maxfopb, pordp, retopb);
endp;


			/* MLE PROCEDURE FOR MNL EQUATION */

proc lgt(params,data);
local y, x, k, capJ, n, beta, theta, rho, v_ij, p_a, p_b, p_ogev;
local v_ijall, p_denom, p_dgt, p_ij, jrep, d_ij, logl;
local vnam,mean,var,std,min,max,valid,mis;
local res, betar, noxs, b_rep, row_rep, beta1, betanr;

load res, noxs, capJ;


/* DEFINE DATA */

y	= data[.,1];
x 	= data[.,2:cols(data)];
k	= cols(x);
n	= rows(data);

y 	= y - minc(y);						@ run y from 0 to J-1 @


/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

jrep = 0;
d_ij	= zeros(n,1);
do while jrep lt capJ;
	d_ij = d_ij~(y .eq jrep);
jrep=jrep+1;
endo;

d_ij = d_ij[.,2:capJ+1];


/* DELETE OVERALL CONSTANT AND INCLUDE CHOICE SPECIFIC CONSTANT MINUS ONE (FIRST ONE) */

if sumc(res') eq k;
    if stdc(x[.,1]) eq 0;
        x = x[.,2:cols(x)];
    endif;
    /* x       = d_ij[.,2:cols(d_ij)]~x;    */
endif;


/* DEFINE PARAMETERS */

if sumc(res') eq 0; 

	beta		= reshape(params,capJ-1,k)';

elseif sumc(res') ge 1 and sumc(res') lt k;

	beta1	= params[1:k];
	betanr	= reshape(params[k+1:noxs],capJ-2,k-sumc(res'))';
	beta		= betanr[1,.]; 
	
	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	beta = beta1~beta;

    elseif sumc(res') eq k;

    beta    = params;

endif;


/* CONSTRUCT OGEV-TYPE PROBABILITIES */

v_ij		= exp(x*beta); 

if sumc(res') lt k;

    p_ij		= (1/(1 + sumc(v_ij'))) ~ (v_ij ./ (1 + sumc(v_ij')));

    if fne(meanc(sumc(p_ij')),1) eq 1;
	    errorlog "ERROR IN DEFINING MNL PROBABILITIES; PROGRAM TERMINATED";
			    "mean probablitiy =";;meanc(sumc(p_ij'));stop;
    endif;

    logl = ln(sumc((d_ij .* p_ij)'));

else;

    p_ij = v_ij ./ sumc(v_ij);
    logl = ln(p_ij);
    
endif;

retp(logl);
endp;


			/* MLE PROCEDURE FOR DOGIT EQUATION */

proc dgt(params,data);
local y, x, k, capJ, n, beta2, beta, theta1, theta, rho, v_ij, p_a, p_b, p_ogev;
local v_ijall, p_denom, p_dgt, p_ij, jrep, d_ij, logl, phat;
local infile, fh, seltheta, jj, i, check, thetas;
local res, betar, noxs, b_rep, row_rep, beta1, betanr, beta_1;

load res, noxs, capJ;


/* DEFINE DATA */

y	= data[.,1];
x 	= data[.,2:cols(data)];
k	= cols(x);
/* J	= maxc(y) - minc(y) +1;	*/
y 	= y - minc(y);						@ run y from 0 to J-1 @
n	= rows(data);


/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

jrep = 0;
d_ij	= zeros(n,1);
do while jrep lt capJ;
	d_ij = d_ij~(y .eq jrep);
jrep=jrep+1;
endo;

d_ij = d_ij[.,2:capJ+1];

load thetas,phat;

seltheta 	= thetas;
phat 	= selif(phat,seltheta);


/* DEFINE PARAMETERS */
 
theta1	= (params[rows(params)-sumc(seltheta)+1:rows(params)]);  
beta1	= params[1:rows(params)-sumc(seltheta)]; 

if rows(params) ne (rows(theta1)+rows(beta1));
	stop;
endif;


/* PULL OUT BETAS - RESTRICTED IF APPROPRIATE */

if sumc(res') eq 0; 

	beta2	= reshape(beta1,capJ-1,k)';

elseif sumc(res') ge 1;

	beta_1	= params[1:k];
	betanr	= reshape(beta1[k+1:noxs],capJ-2,k-sumc(res'))';
	beta		= betanr[1,.]; 
	
	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	beta2 = beta_1~beta;

endif;

beta		= zeros(k,1)~beta2;
rho		= 1;

theta	= 0;
jj		= 1;
i		= 1;
do while i le capJ;
	if seltheta[i] eq 1;
		theta	= theta|theta1[jj];
		jj		= jj+1;
	else;
		theta	= theta|0;
	endif;
i=i+1;
endo;

theta=theta[2:rows(theta)];


/* CONSTRUCT OGEV-TYPE PROBABILITIES */

v_ij		= exp((x*beta) ./ rho); 
v_ijall	= zeros(n,1)~v_ij~zeros(n,1);

p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]) .^rho;
p_a		= (v_ij ./ sumc(p_denom'));  

p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho-1)) + 
			((v_ij + v_ijall[.,3:capJ+2]) .^ (rho-1)); 
			
p_ogev	= p_a .* p_b;


/* DOGEV-TYPE PROBABILITIES */

p_ij = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

/*
if fne(meanc(sumc(p_ij')),1) eq 1;
	errorlog "ERROR IN DEFINING DOGIT PROBABILITIES; PROGRAM TERMINATED";
			"mean probablitiy =";;meanc(sumc(p_ij'));stop;
endif;
*/

logl = ln(sumc((d_ij .* p_ij)'));

retp(logl);
endp;



			/* MLE PROCEDURE FOR OGEV EQUATION */

proc ogv(params,data);
local y, x, k, capJ, n, beta2, beta, theta, v_ij, p_a, p_b, p_ogev;
local v_ijall, p_denom, p_dgt, p_ij, jrep, d_ij, logl;
local rho;
local res, betar, noxs, b_rep, row_rep, beta1, betanr, beta_1;

load res, noxs, capJ;

/* DEFINE DATA */

y	= data[.,1];
x 	= data[.,2:cols(data)];
k	= cols(x);
/* J	= maxc(y) - minc(y) +1;	*/
y 	= y - minc(y);						@ run y from 0 to J-1 @
n	= rows(data);

/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

jrep = 0;
d_ij	= zeros(n,1);
do while jrep lt capJ;
	d_ij = d_ij~(y .eq jrep);
jrep=jrep+1;
endo;

d_ij = d_ij[.,2:capJ+1];


/* DEFINE PARAMETERS */
 
rho 		= params[rows(params)]; 
beta1	= params[1:rows(params)-1];

if rows(params) ne (rows(rho)+rows(beta1));
	stop;
endif;

/* PULL OUT BETAS - RESTRICTED IF APPROPRIATE */

if sumc(res') eq 0; 

	beta2	= reshape(beta1,capJ-1,k)';

elseif sumc(res') ge 1;

	beta_1	= params[1:k];
	betanr	= reshape(beta1[k+1:noxs],capJ-2,k-sumc(res'))';
	beta		= betanr[1,.]; 
	
	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	beta2 = beta_1~beta;

endif;

beta		= zeros(k,1)~beta2;
theta	= zeros(capJ,1);


/* CONSTRUCT OGEV PROBABILITIES */

v_ij		= exp((x*beta) ./ rho); 
v_ijall	= zeros(n,1)~v_ij~zeros(n,1);

p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]).^rho;

p_a		= (v_ij ./ sumc(p_denom'));  

p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho-1)) + 
			((v_ij + v_ijall[.,3:capJ+2]) .^ (rho-1)); 
			
p_ogev	= p_a .* p_b;


/* DOGEV-TYPE PROBABILITIES */

p_ij = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

if meanc(sumc(p_ij')) gt 1.0001 or meanc(sumc(p_ij')) lt .9999;
	errorlog "ERROR IN DEFINING OGEV PROBABILITIES; PROGRAM TERMINATED";
			"mean probablitiy =";;meanc(sumc(p_ij'));stop;
endif;

logl = ln(sumc((d_ij .* p_ij)'));


retp(logl);
endp;


			/* MLE PROCEDURE FOR DOGEV EQUATION */

proc dgv(params,data);
local y, x, k, capJ, n, beta2, beta, theta1, theta, rho, v_ij, p_a, p_b, p_ogev;
local v_ijall, p_denom, p_dgt, p_ij, jrep, d_ij, logl;
local infile, fh, seltheta, jj, i, thetas;
local res, betar, noxs, b_rep, row_rep, beta1, betanr, beta_1;

load res, noxs, capJ;
save params;

/* DEFINE DATA */

y	= data[.,1];
x 	= data[.,2:cols(data)];
k	= cols(x);
/* J	= maxc(y) - minc(y) + 1;	*/
y 	= y - minc(y);						@ run y from 0 to J-1 @

n	= rows(data);

load thetas;
seltheta 	= thetas;

/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

jrep = 0;
d_ij	= zeros(n,1);
do while jrep lt capJ;
	d_ij = d_ij~(y .eq jrep);
jrep=jrep+1;
endo;

d_ij = d_ij[.,2:capJ+1];


/* DEFINE PARAMETERS */

beta1	= params[1:rows(params)-sumc(seltheta)-1]; 
rho 		= params[rows(params)]; 
theta1	= (params[rows(params)-sumc(seltheta):rows(params)-1]); 

if rows(params) ne (rows(theta1)+rows(rho)+rows(beta1));
	stop;
endif;
 

/* PULL OUT BETAS - RESTRICTED IF APPROPRIATE */

if sumc(res') eq 0; 

	beta2	= reshape(beta1,capJ-1,k)';
	
elseif sumc(res') ge 1;

	beta_1	= params[1:k];
	betanr	= reshape(beta1[k+1:noxs],capJ-2,k-sumc(res'))';
	beta		= betanr[1,.]; 
	
	b_rep 	= 2;
	row_rep 	= 2;
	do while b_rep le k;

		if res[1,b_rep] eq 0;
			beta		= beta|betanr[row_rep,.];
			row_rep	= row_rep + 1;
		else;
			beta		= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
		endif;

	b_rep = b_rep + 1;
	endo;

	beta2 = beta_1~beta;

endif;

beta		= zeros(k,1)~beta2;

theta	= 0;
jj		=1;
i		=1;
do while i le capJ;
	if seltheta[i] eq 1;
		theta	= theta|theta1[jj];
		jj		= jj+1;
	else;
		theta	= theta|0;
	endif;
i=i+1;
endo;

theta=theta[2:rows(theta)];


/* CONSTRUCT OGEV PROBABILITIES */

v_ij		= exp((x*beta) ./ rho); 
v_ijall	= zeros(n,1)~v_ij~zeros(n,1);

p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]).^rho;

p_a		= (v_ij ./ sumc(p_denom'));  

p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho-1)) + 
			((v_ij + v_ijall[.,3:capJ+2]) .^ (rho-1)); 
			
p_ogev	= p_a .* p_b;


/* DOGEV PROBABILITIES */

p_ij = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

if meanc(sumc(p_ij')) gt 1.0001 or meanc(sumc(p_ij')) lt .9999;
	errorlog "ERROR IN DEFINING DOGEV PROBABILITIES; PROGRAM TERMINATED";
			"mean probablitiy =";;meanc(sumc(p_ij'));stop;
endif;

logl = ln(sumc((d_ij .* p_ij)'));

retp(logl);
endp;


			/* ORDERED PROBIT MODEL */

proc opb(params,data);
local y, x, k, capJ, n, jrep, jrep2, d_ij, beta, mu, xbeta, p_i1, p_icapj, p_irest, p_ij, logl, bands;
load capJ;

/* DEFINE DATA */

y	= data[.,1];
x 	= data[.,2:cols(data)];
k	= cols(x);
/* J	= maxc(y)  - minc(y) +1;	*/
y 	= y - minc(y);						@ run y from 0 to J-1 @
n	= rows(data);

/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

jrep = 0;
d_ij	= zeros(n,1);
do while jrep lt capJ;
	d_ij = d_ij~(y .eq jrep);
jrep=jrep+1;
endo;

d_ij = d_ij[.,2:capJ+1];


/* DEFINE PARAMETERS */
 
beta		= params[1:k];
mu		= params[k+1:rows(params)]; 

if rows(params) ne (rows(beta)+rows(mu));
	stop;
endif;

/* DETERMINISTIC PART & PROBABILITIES */

xbeta = x*beta;

p_ij	= 1-(cdfn(xbeta));
p_ij	= p_ij~(cdfn(mu[1]-xbeta) - cdfn(-xbeta));
jrep2= 2;
do while jrep2 le capJ-2;
	p_ij = p_ij~(cdfn(mu[jrep2]-xbeta) -
          	            cdfn(mu[jrep2-1]-xbeta));
jrep2 = jrep2 + 1;
endo;

p_ij = p_ij~(1 - cdfn(mu[capJ-2]-xbeta));


if meanc(sumc(p_ij')) gt 1.0001 or meanc(sumc(p_ij')) lt .9999;
	errorlog "ERROR IN DEFINING ORDERED PROBIT PROBABILITIES; PROGRAM TERMINATED";
			"mean probablitiy =";;meanc(sumc(p_ij'));stop;
endif;

logl = ln(sumc((d_ij .* p_ij)'));


retp(logl);
endp;


			/* PROCEDURE FOR UPPER BOUNDS FOR THETA'S - DOGIT MODEL */

proc ineqp(b);
local infileb, fhb, sel, theta, thetas, phat, noxs, mark, mark2;
	load thetas, phat, noxs;
	sel		= thetas;
	phat		= selif(phat,sel);
	theta	= b[noxs+1:noxs+sumc(sel)];
	mark 	= theta ./ (1+sumc(theta));
	mark2	= sumc(mark .le phat);
	/* retp( mark2 - sumc(sel) ); */
	retp( (-mark)-(-phat)); 
endp;


			/* PROCEDURE FOR UPPER BOUNDS FOR THETA'S & RHO - DOGEV MODEL: NOT USED */

proc ineqp1(b);
local infileb, fhb, sel, thetas, theta, rho, phat, noxs, mark, mark2;
	load thetas, phat, noxs;
	sel		= thetas;
	phat	= (selif(phat,sel))|1;
	theta	= b[noxs+1:noxs+sumc(sel)];
	rho		= b[rows(b)];
	mark 	= (theta ./ (1+sumc(theta))) | rho;
	retp( (-mark)-(-phat)); 
endp;




                        /* PROCEEDURES FOR ME'S */


			/* CONSTRUCT GEV-TYPE PROBABLITIES */

proc 1 = dogev_p(zbar);
local beta, beta1, beta_1, betanr, betatmp, b_rep, row_rep, theta, v_ij, v_ijall, p_ogev, p_ogv, pdgv, check, dogev, thetatmp;
local p_denom, rho, jj, i, p_a, p_b;
local thetas, capJ, k, res, noxs;
load thetas, capJ, k, res, noxs;


    beta = b_me[1:noxs];

    if sumc(res') eq 0; 

	    beta	= reshape(beta,capJ-1,k)';
	
    elseif sumc(res') ge 1;

	    beta_1	= beta[1:k];

	    betanr	= reshape(beta[k+1:rows(beta)],capJ-2,k-sumc(res'))';
	    beta	= betanr[1,.]; 

	    b_rep 	= 2;
	    row_rep 	= 2;
	    do while b_rep le k;

		    if res[1,b_rep] eq 0;
			    beta	= beta|betanr[row_rep,.];
			    row_rep	= row_rep + 1;
		    else;
			    beta	= beta|(zeros(1,capJ-2) + beta_1[b_rep]);
		    endif;

	    b_rep = b_rep + 1;
	    endo;

	    beta 	= beta_1~beta;

	    dogev	= vec(beta)|theta_me|rho_me;

    endif;

    beta	= zeros(k,1)~beta;

    theta	=0;
    jj		=1;
    i		=1;
    do while i le capJ;
	    if thetas[i] eq 1;
		    theta	= theta|theta_me[jj];
		    jj		= jj+1;
	    else;
		    theta	= theta|0;
	    endif;
    i=i+1;
    endo;

    theta = theta[2:rows(theta)];

    v_ij	= exp((zbar'*beta) ./ rho_me); 
    v_ijall	= zeros(rows(v_ij),1)~v_ij~zeros(rows(v_ij),1);

    p_denom	= (v_ijall[.,1:capJ+1] + v_ijall[.,2:capJ+2]).^rho_me;

    p_a		= (v_ij ./ sumc(p_denom'));  

    p_b		= ((v_ijall[.,1:capJ] + v_ij) .^ (rho_me-1)) + 
	    		((v_ij + v_ijall[.,3:capJ+2]) .^ (rho_me-1)); 
			
    p_ogev	= p_a .* p_b;

    pdgv = (theta' ./ (sumc(theta) + 1)) + ((1./ (sumc(theta) + 1)) .* p_ogev);

    check = meanc(sumc(pdgv'));
    if check gt 1.00001 or check lt (1-0.00001);
	    errorlog "DOGEV PROBABILITIES DON'T SUM TO 1: AVERAGE PROB IS";check;stop;
    endif;

retp(pdgv[j]);
endp;


                        /* CONSTRUCT ORDERED PROBIT PROBABLITIES */

proc 1 = ordp_p(zbar);
local ordp, zbeta, mu, z, p_ij, jrep2, cdfn, pordp, gama, k, capJ;
load k, capJ;

    if cols(zbar) eq 1;
        zbar = zbar';
    endif;

	zbeta   = zbar*b_me;

	p_ij    = 1-(cdfn(zbeta));
	p_ij    = p_ij~(cdfn(mu_me[1]-zbeta) - cdfn(-zbeta));
	jrep2= 2;

	for jrep2 (2,capJ-2,1);
        	p_ij = p_ij~(cdfn(mu_me[jrep2]-zbeta) -
                            cdfn(mu_me[jrep2-1]-zbeta));
	endfor;

	pordp = p_ij~(1 - cdfn(mu_me[capJ-2]-zbeta));

retp(pordp[j]);
endp;

