/*	MARYANN.PRG
**
** 	PROGRAM FOR TOBIT REGRESSION OF PANEL DATA
**	USING SIMULATED MAXIMUM LIKELIHOOD ESTIMATION
**
**	BY MARY ANN HOFMANN
**
*/

/*	NOTE ABOUT DATA SET CONFIGURATION
**	THE DATA SET MUST BE ARRANGED SO THAT THERE IS ONE ROW
**	FOR EACH INDIVIDUAL, AND THE DEPENDENT AND INDEPENDENT
**	VARIABLES ARE ARRANGED IN COLUMNS BY YEAR. 
**	THE DATA SET MUST INCLUDE A GROUP OF DUMMY VARIABLES,
**	WHERE THE VALUE IS ONE IF THE DEPENDENT VARIABLE IS NONZERO,
**	ZERO OTHERWISE.  FOR EXAMPLE:
**	
**	OBS #	DEP_T1	DEP_T2	X1_T1	X1_T2	X2_T1	X2_T2	DV_T1	DV_T2
**
**	1        100     101     5000   5025      3       3       1       1
**
**	2         0       0      4250   4300      4       5       0       0
**
**
*/

	new;

@ open output file and format it @
	output file = maryann.out reset;
	format  /rd 12, 4;

@ define n and t @
	n = 103;    @ # of individuals  @
	t = 4;      @ # of time periods @

@ cc is added to the mean function value before taking the log @
	cc = 0;

@ read data @
	load dat[n,17] = tobit.db;

@ define variables @
	id 	= dat[.,1];
	ch87	= dat[.,2];	
	ch88	= dat[.,3]; @ logs of charitable contributions ($)  @
	ch89	= dat[.,4];	
	ch90	= dat[.,5];
	pr87	= dat[.,6]; @ after-tax price of contributions: log(1-tax rate)    @
	pr88	= dat[.,7];
	pr89	= dat[.,8];
	pr90	= dat[.,9];
	inc87	= dat[.,10]; @ log of disposable y ($1000) before contribution@
	inc88	= dat[.,11];
	inc89	= dat[.,12];
	inc90	= dat[.,13];
	dv87	= dat[.,14]; @ dummy variable = 1 if ch > 0 @
	dv88	= dat[.,15];
	dv89	= dat[.,16];
	dv90	= dat[.,17];

@ define matrices @
	yy  = ch87~ch88~ch89~ch90;	@ dependent variable @
	pr  = pr87~pr88~pr89~pr90;	@ first X variable   @
	inc = inc87~inc88~inc89~inc90;	@ second X variable  @
	dv  = dv87~dv88~dv89~dv90;	@ dummy variables    @

@ create matrix of random numbers for simulations  @
	draws = 20; @ defines H (number of simulations)     @
	seed1 = 1;
	rnd   = rndns(n,draws,seed1);

@ set start values (use best guess or results from Limdep)   @
@ NOTE:  B[5] cannot be zero  @
	b0 = {1,-1,1,1,1};


@ define f-tilda function @
proc tilda(y,x1,x2,dum,rh,p);
local j, ff, ffj;
	j = 1;
	ff = 1;
	do while j <= t;
	ffj = 	(((1/abs(p[5]))*pdfn((y[1,j] - p[1] -p[2]*x1[1,j] - p[3]*x2[1,j]
		- p[4]*rh)/abs(p[5])))^dum[1,j]) *
		((cdfn((- p[1] -p[2]*x1[1,j] - p[3]*x2[1,j] - p[4]*rh)
		/ abs(p[5])))^(1-dum[1,j]));
	ff = ff*ffj;
	j = j + 1;
	endo;
	retp(ff);
	endp;

@set up the Maximum Likelihood function @

library optmum;
#include optmum.ext;
OPTSET;

@ define likelihood function @
 	proc smle(b);
	local i, h, lnf, fst, fsth, mnfst, lnfst ;


@ set counter for individuals @ 
@ set up lnf to accumulate sum of ln f(.) over individuals @
	i=1;
	lnf = 0;
	do while i <= n;

@ set counter for draws @ 
@ set up fst to accumulate sum of simulations @
		h    =  1;
		fst  =  0;
		do while h <= draws;
		fsth = 	tilda(yy[i,.],pr[i,.],inc[i,.],dv[i,.],rnd[i,h], b);
		fst  =  fst + fsth;
		h    =  h + 1;
		endo;

@ compute mean of fst and take its log @
	mnfst = fst/draws;
	lnfst = ln(mnfst+cc);

@ increase sum of ln f(.) @

	lnf = lnf + lnfst;
	i = i + 1;
	endo;
	
	retp(-lnf/n) ;
	endp;


__title = "TOBIT by SML";
_opgtol = 0.00001;
_opstmth = "BFGS HALF";

{b, func, grad, retcode} = optprt(optmum(&smle,b0));

@ estimate standard errors @

@ create initial matrix that will be used to compute covariance matrix @
sum = zeros(5,5);
i = 1;
do while i <= n;
	@ compute numerator (num) and denominator (den) for s @
	num = zeros(1,5);
	den = 0;
	h = 1;
	do while h <= draws;
		proc tildb(p);
		local til;
		til = tilda(yy[i,.],pr[i,.],inc[i,.],dv[i,.],rnd[i,h],p);
		retp(til);
		endp;
		clo = tildb(b);
		num = num + gradp(&tildb,b);
		den = den + clo;
	h = h + 1;
	endo;

s = num/den;
ss = (s's);
sum = sum + ss;
i = i + 1;
endo;

@ compute the covariance matrix as the inverse of sum @
cov = invpd(sum);	

@ obtain the standard errors  and t-statistics @
se = sqrt(diag(cov));
tst = b./se;

@ print the results @
"";
"SMLE TOBIT RESULTS";
"";
" coefficient	std. error	t-statistic";
b~se~tst; 


output off;