/* 
** 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.                                                                    **
*/


new;

          /* PROGRAM TO ESTIMATE GEV AND ORDERED PROBIT MODELS */
library cml;
#include cml.ext;
#include c:\mark\gauss\procs\gevs.prc;
#include c:\mark\gauss\procs\hitmiss.prc; 
cmlset;

/* SET DIRECTORIES */

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


/* 1 LOAD RAW DATA AND NAME VARIABELS */

infile1 = "c:\\mark\\data\\gauss\\Eritrea";
open fh = ^infile1 varindxi;
vnames 	= getname(infile1);
xx 		= readr(fh,rowsf(fh)); 
age     = xx[.,iHHAGE_A] ./ 10;
age2    = (age .^2) ./ 10;
trans1 =xx[.,itrans3_a]./1000;
trans2=xx[.,itrans4_a]./1000;
house1  = xx[.,ihouse]./10;
sewage1=xx[.,isewage]./10;
/*
moon = (xx[.,ichilda_1] + xx[.,ichild2_a]);@
deprate = xx[.,imoon]./xx[.,ihlab_a];
*/



/* DEFINE DEPENDENT & EXPLANATORY VARIABLES */

y	= xx[.,iINDEX]; 
y   =y - minc(y);
/*
x 	= ones(rows(xx),1)~age~age2~trans1~trans2~ deprate~xx[.,iHSIZE_A iFIGHTER iHEMP_A iunrate iEMPSTAT1 
       iEMPSTAT2 iEMPSTAT5 ihouse isewage iEDUC2 iEDUC3 iEDUC4 iEDUC5  imarital1 
       imarital2 itown1 itown2 itown3 itown4 itown5 itown6 itown7 itown8 itown9 ireligio1 iethnic iRETURNE1];
*/

x 	= ones(rows(xx),1)~age~age2~trans1~trans2~xx[.,iHSIZE_A iFIGHTER iHEMP_A iunrate iEMPSTAT1 
       iEMPSTAT2 iEMPSTAT5 ihouse isewage iEDUC2 iEDUC3 iEDUC4 iEDUC5  imarital1 
       imarital2 itown1 itown2 itown3 itown4 itown5 itown6 itown7 itown8 itown9 ireligio1 iethnic iRETURNE1];

x   = ones(rows(xx),1)~age~age2~trans1~trans2~xx[.,iHSIZE_A iFIGHTER iHEMP_A iunrate ihouse isewage iEDUC2 iEDUC3 iEDUC4 iEDUC5 imarital1
       imarital2 iethnic];

yx 	= y~x;
@ yx 	= delif(y~x,del); @


/* DELETE ANY MISSING VALUES */

yx	= packr(yx);    @ deletes any rows where there are any missing values @

y 	= yx[.,1];
x 	= yx[.,2:cols(yx)];
capj	= maxc(y) + 1;


			/* NAME X'S & SAVE */
/*
let yxvars = INDEX CONSTANT AGE AGE2 TRANS1 TRANS2 deprate HSIZE_A FIGHTER HEMP unrate EMPSTAT1 
       EMPSTAT2 EMPSTAT5 house sewage EDUC2 EDUC3 EDUC4 EDUC5 marital1 
  marital2 town1 town2 town3 town4 town5 town6 town7 town8 town9 religion ethnic RETURNE1;
*/

let yxvars = INDEX CONSTANT AGE AGE2 TRANS1 TRANS2 HSIZE_A FIGHTER HEMP unrate EMPSTAT1 
       EMPSTAT2 EMPSTAT5 house sewage EDUC2 EDUC3 EDUC4 EDUC5 marital1 
  marital2 town1 town2 town3 town4 town5 town6 town7 town8 town9 religion ethnic RETURNE1;

let yxvars = INDEX CONSTANT age age2 trans1 trans2 HSIZE_A FIGHTER HEMP_A unrate ihouse isewage EDUC2 EDUC3 EDUC4 EDUC5 marital1 
                marital2 iethnic;

            if not saved(y~x,"ind2_dat",yxvars);
                errorlog "Disk full";stop;
                end;
            endif;

       
/* DESCRIPTIVE STATS */

{ vnam,mean,var,std,min,max,valid,mis } = dstat("ind2_dat",yxvars); 

/* CONSTRUCT INDICATOR VARIABLE FOR j=J */

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

d_ij = d_ij[.,2:maxc(y)+2];

phat = meanc(d_ij);


	/* ENTER ROW VECTOR WITH 1 IF COEFFECIENT IS TO BE EQUAL ACROSS j (IN STRUCTURAL MODEL)*/

res = zeros(1,cols(x));

    /* ENTER M THE NUMBER OF BOOTSRAPPED ITERATIONS FOR THE STANDRAD ERRORS OF MARGINAL EFFECTS */

M   = 1000;

    /* ENTER logit, dogit, ogev, dogev */

/* 
logit 	- MNL only;
dogit 	- dogit and MNL
ogev 	- ogev and MNL
dogev 	- all models including ordered Probit
ordered	- ordered probit and MNL 
*/

let model		= dogev;

/* ENTER A ROW VECTOR WITH ELEMENTS EQUAL TO 1 IF THETA IS TO BE ESTIMATED AND ZERO OTHERWISE */

thetas		= ones(capj,1);
thetas[1]=0;


/* ENTER START VALUES - OPTIONAL, SET TO 0 FOR DEFAULT START VALUES */


let startmnl    = 0;
let startdgt    = 0;
let startogv    = 0;
let startdgv    = 0;
let startopb	= 0; 

{reg, 
    logit, mnl_mes, mnl_mese, vclgt, logllgt, plgt, retlgt,
	dogit, dgt_mes, dgt_mese, vcdgt, logldgt, pdgt, retdgt,
	ogev, ogv_mes, ogv_mese,vcogv, loglogv, pogv, retogv,
	dogev, dgv_mes, dgv_mese, vcdgv, logldgv, pdgv, retdgv,
	ordprb, op_mes,  op_mese,  vcopb, loglopb, pordp, retopb} = mnmodels("ind2_dat",model,thetas,
										startmnl,startdgt,startogv,startdgv,startopb,res,M); 



/* PRINT RESULTS */

/* CLEAR SCREEN */

/* cls;  */

xvars 	    = yxvars[2:rows(yxvars)];
xvars_me    = xvars;
let mu	    = mu;
opbvars	    = xvars|(mu .*. ones(capj-2,1));
xvars 	    = ones(capj-1,1) .*. xvars;
let theta 	= theta;


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


dgtvars 	= xvars|(theta .*. ones(sumc(thetas),1)); 
ogvvars 	= xvars|"RHO";
dgvvars 	= xvars|(theta .*. ones(sumc(thetas),1))|"RHO";

output file = c:\mark\gauss\output\DOGEV\Eritrea_2a.out reset; 
output on;

format 15,4;

save plgt, y;
print;
print;
"MULTINOMIAL LOGIT";
print;
"Convergence code:";;$retlgt;
print;

i=1;
do while i le rows(logit);
    $xvars[i];;logit[i,.];
    i=i+1;
endo;

print;
"Max L";;logllgt;
print;
print;

"MARGINAL EFFECTS AND STANDARD ERRORS FOR AMOUNT OF PARTICIPATION: MNL";

print;
for i (1,rows(xvars_me),1);
   		$xvars_me[i];;mnl_mes[i,.];
        "                ";;mnl_mese[i.,];
endfor;
print;

if whichm eq 1;
   goto hitmiss;
endif;

if whichm eq 2 or whichm eq 4;
	print;
	print;
	"DOGIT";
	print;
	"Convergence code:";;$retdgt;
	print;
	i=1;
	do while i le rows(dogit);
		$dgtvars[i];;dogit[i,.];
		i=i+1;
	endo;
	print;
	"Max L";;logldgt;
	print;
	print;

    "MARGINAL EFFECTS AND STANDARD ERRORS FOR AMOUNT OF PARTICIPATION: DOGIT";

    print;
    for i (1,rows(xvars_me),1);
   	    	$xvars_me[i];;dgt_mes[i,.];
            "                ";;dgt_mese[i.,];
    endfor;
    print;
	
	save pdgt;

endif;

if whichm eq 3 or whichm eq 4;
	print;
	print;
	"ORDERED EXTREME VALUE";
	print;
	"Convergence code:";;$retogv;
	print;

	i=1;
	do while i le rows(ogev);
		$ogvvars[i];;ogev[i,.];
		i=i+1;
	endo;

	print;
	"Max L";;loglogv;
	print;
	print;

    "MARGINAL EFFECTS AND STANDARD ERRORS FOR AMOUNT OF PARTICIPATION: OGEV";

    print;
    for i (1,rows(xvars_me),1);
   	    	$xvars_me[i];;ogv_mes[i,.];
            "                ";;ogv_mese[i.,];
    endfor;
    print;

    save pogv;

endif;

if whichm eq 4;
	print;
	print;
	"DOGIT ORDERED EXTREME VALUE";
	print;
	"Convergence code:";;$retdgv;
	print;
	i=1;
	do while i le rows(dogev);
		$dgvvars[i];;dogev[i,.];
		i=i+1;
	endo;
	print;
	"Max L";;logldgv;
	print;
	print;

    "MARGINAL EFFECTS AND STANDARD ERRORS FOR AMOUNT OF PARTICIPATION: DOGEV";

    print;
    for i (1,rows(xvars_me),1);
   	    	$xvars_me[i];;dgv_mes[i,.];
            "                ";;dgv_mese[i.,];
    endfor;
    print;


    save pdgv;

endif;

if whichm eq 2 or whichm eq 4;
	print;
	"THETA VECTOR";
	print;
	thetas;
	print;
endif;

if whichm eq 4 or whichm eq 5;
	print;
	print;
	"ORDERED PROBIT";
	print;
	"Convergence code:";;$retopb;
	print;
	i=1;
	do while i le rows(ordprb);
		$opbvars[i];;ordprb[i,.];
		i=i+1;
	endo;
	print;
	"Max L";;loglopb;
	print;
	print;

	save pordp;

endif;

hitmiss:


/* PRINT OUT HIT AND MISS TABLE */

print;
"MNL";
{ table1} = hitmiss(plgt,y);
print;
if whichm eq 2 or whichm eq 4;
	print;
	"DOGIT";
	{ table2} = hitmiss(pdgt,y);
	print;
endif;

if whichm eq 3 or whichm eq 4;
	print;
	"OGEV";
	{ table3} = hitmiss(pogv,y);
	print;
endif;

if whichm eq 4;
	print;
	"DOGEV";
	{ table4} = hitmiss(pdgv,y);
	print;
endif;
output off;

stop;
