/* Tobit Model */
cls;

/* Monte Carlo Simulation */

loops = 1000;
nobs = 101;

sigma = 2;

beta1 = 50;
beta2 = 0.6;

betamat = beta1|beta2;
k = rows(betamat);

cutoff = 20; // Top percentage to censor
cut = round(nobs*cutoff/100);


/* Create fixed x's */
x1 = ones(nobs,1);
x2 = 100*rndu(nobs,k-1);
x = x1~x2;

// To Show Data:
// ystar = x*betamat + sigma*rndn(nobs,1);
// ystar~x


/* Storage matrices */
olstruebstore = zeros(loops,k);
olstruetstat = zeros(loops,k);

olsbstore = zeros(loops,k);
olststat = zeros(loops,k);

mlebstore = zeros(loops,k);


/* Actual MC loops */

for j (1,loops,1);

	// Creates latent, unobserved y's 
	
	// Homoskedastic
	ystar = x*betamat + sigma*rndn(nobs,1);
	// Heteroskedastic
	//ystar = x*betamat + sigma*rndn(nobs,1).*x[.,2]./50;

	latent = ystar~x;
	latent = sortc(latent,1);

	obsdata = latent;


	/* Makes "observed" or censcored data */
	for i (1,nobs,1);
	
		if i <= (nobs - cut + 1);
			obsdata[i,1] = obsdata[i,1];
		else;
			obsdata[i,1] = obsdata[(nobs - cut + 1),1];
		endif;
	
	endfor;
	
	ystar = latent[.,1];
	y = obsdata[.,1];
	x = obsdata[.,2:(k+1)];


	// OLS Regression on "true" or latent data 
	olstruebstore[j,.] = (inv(x'*x)*(x'*ystar))';
	resid = ystar - x*olstruebstore[j,.]';
	s2 = resid'*resid/(nobs - k);
	varb = s2*inv(x'*x);
	seb = sqrt(diag(varb));
	olstruetstat[j,.] = (olstruebstore[j,.]' ./ seb)';
	
	
	// OLS Regression on observed data 
	olsbstore[j,.] = (inv(x'*x)*(x'*y))';
	resid = ystar - x*olsbstore[j,.]';
	s2 = resid'*resid/(nobs - k);
	varb = s2*inv(x'*x);
	seb = sqrt(diag(varb));
	olststat[j,.] = (olsbstore[j,.]' ./ seb)';
	
	library cml;
	#include cml.ext;

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



	/* Call the CML Module */
	cmlset;
	_cml_DirTol = 1e-4;
	_cml_Algorithm = 3;
	//_cml_GradProc = &gradient; 
	//_cml_GradCheckTol = 0.01;
	cml_Bounds = { 0 100 } ; // constrain parameters to be positive
	startv = 1|1|1;
	__output = 0; //Only final results are shown; not by iteration
	{mltheta,f0,g,cov,ret} = CMLPrt(CML(obsdata,0,&tobit,startv));
	
	mlebstore[j,.] = mltheta[1:k,1]';
	
	


endfor;


print "OLS (True) Beta Mean" "     OLS (True) t-Stat Mean"; 
print meanc(olstruebstore)~meanc(olstruetstat);
print "OLS Mean            " "     OLS t-Stat Mean";
print meanc(olsbstore)~meanc(olststat);
print "MLE Mean" meanc(mlebstore);

print "OLS 'Shortcut'";
print meanc(olsbstore)/(1 - cutoff/100);


/* Need:
 variance for ols, both true and "false"
 variance for mle 
t-stats for mlebstore
graphs for all regressions
truncated regression
censored from below regression option


*/

/* Graphs */
//library pgraph;
//graphset;
//title("OLS (True) Slope Distribution");
//xtics(0,1.2,0.2,1);
//plothist(olstruebstore[.,2],10);
//plothist(olsbstore[.,2],10);  
//plothist(mlebstore[.,2],10);




yhat = x*mlebstore[loops,.]';
presid = ystar - yhat;
//plotXY(x[.,2],presid);





/* Change Points

beta1
sigma
cutoff
heteroskedastic

*/