/* Tobit Model */
cls;

/* Load Data and Define Variables */
load path = /Users/stevetargos/Documents/GAUSS;
load data[6367,10] = tobitdata.txt;
Constant = data[2:6367,1]; //C
RateMarr = data[2:6367,2]; //v1
Age = data[2:6367,3]; //v2
YrsMarr = data[2:6367,4]; //v3
NumKids = data[2:6367,5]; //v4
Religious = data[2:6367,6]; //v5
Education = data[2:6367,7]; //v6
WifeOcc = data[2:6367,8]; //v7
HusOcc = data[2:6367,9]; //v8
Y = data[2:6367,10]; // yrb
nobs = rows(Y);

Xmat = Constant~RateMarr~Age~YrsMarr~NumKids~Religious~Education~WifeOcc~HusOcc;

k = cols(Xmat);
vnum = seqa(1,1,k);
vname1 = "constant";
vname2 = "RateMarr";
vname3 = "Age";
vname4 = "YrsMarr";
vname5 = "NumKids";
vname6 = "Religious";
vname7 = "Education";
vname8 = "WifeOcc";
vname9 = "HusOcc";

let string vnames = {"Constant", "RateMarr", "Age", "YrsMarr", "NumKids", "Religious", "Education", "WifeOcc", "HusOcc"};

bols = inv(Xmat'*Xmat)*(Xmat'*Y);
resid = y - Xmat*bols;
s2 = resid'*resid/(nobs - k);
varbols = s2*inv(Xmat'*Xmat);
sebols = sqrt(diag(varbols));
tstatols = bols./sebols;
print "Variable   	   Coefficient    	Standard Error     t Statistic ";
print "--------------------------------------------------------------------";
print vnum~bols~sebols~tstatols;

print;





/* ***************************** */
/* MAXIMUM LIKELIHOOD ESTIMATION */
/* ***************************** */

library cml;
#include cml.ext;

data = Y~Xmat; //total dataset array

/* Likelihood Procedure */
proc tobit(theta,data);
	local u, sig, y, x, b;
	sig = theta[k+1];
	y = data[.,1];
	x = data[.,2:k+1];
	b = theta[1:k,1];
	u = y[.,1] ./= 0; // 1 or 0 dummy; 1 -> y>0, 0-> y=0
	retp(u.*lnpdfmvn(y-x*b,sig^2) + (1-u).*(ln(cdfnc(x*b/sig))));
endp;


/* Call the CML Module */
cmlset;
_cml_DirTol = 1e-3;
_cml_Algorithm = 3;
//_cml_GradProc = &gradient; 
//_cml_GradCheckTol = 0.01;
startv = bols|1;
__output = 0;
{mltheta,f0,g,cov,ret} = CMLPrt(CML(data,0,&tobit,startv));