#delimit;
capture program drop machado_sel;
capture program drop mach;
mata mata clear;

mata;


real rowvector mean(real matrix X, real colvector w)
{
        real rowvector  CP
        real scalar     n 

        CP = quadcross(w,0, X,1)
        n  = cols(CP)
        return(CP[|1\n-1|] :/ CP[n])
}

real matrix variance(real matrix X, real colvector w)
{
        real rowvector  CP
        real rowvector  means
        real scalar     n 

        CP = quadcross(w,0, X,1)
        n  = cols(CP)
        means = CP[|1\n-1|] :/ CP[n]
        return(crossdev(X,0,means, w, X,0,means) :/ (CP[n]-1))
}

numeric scalar quantile(numeric x, numeric scalar q)
{
	numeric matrix hulp1;
	numeric scalar n, hulp2;
	
	hulp1 = sort(x,1);
	n = rows(x);
	hulp2 = trunc(n * q);

	return (hulp1[hulp2,.]);
}

transmorphic matrix sort(transmorphic matrix x, real rowvector idx)
{
        return(x[order(x,idx), .])
}


numeric matrix bereken_var2(numeric matrix q, numeric matrix theta, numeric matrix theta1,
			numeric matrix yb, numeric matrix y1a, numeric scalar m) {;
	real matrix var2;
	real scalar n2, i, f1b;

	n2 = rows(yb);
	
	for(i=1;i<=rows(theta);i++)
	{
		f1b = kerncv2(yb,theta1[i],0);
		if (i == 1)
		{	
			var2 = (((m/n2) * q[i]*(1-q[i]))/(f1b^2));
		}
		else
		{	
			var2 = var2 \ (((m/n2) * q[i]*(1-q[i]))/(f1b^2));
		}
	}
	
	return (var2);
};

numeric scalar quanc1(numeric matrix x, numeric scalar q, numeric matrix s)
{
	real matrix s1, hulp, x1;
	real scalar n, i, hulp2, hulp3;
		   
	s1 = s :/ mean(s, 1);

	hulp = x, s1;	
	hulp = sort(hulp,1);
	
	x1 = hulp[.,1];
	s1 = hulp[.,2];
	n = rows(s1);
	i = 0;
	hulp2 = 0;
	hulp2a = 0;

	while(i < rows(s1))
	{
		i = i + 1;
		hulp2 = hulp2 + s1[i];		
		hulp2a = hulp2a \ hulp2;		
	}
	
	hulp3 = mean(hulp2a :<= J(rows(hulp2a),1,q * n),1) * rows(hulp2a);

	return (x1[hulp3,.]);
}

numeric matrix kerncv2(numeric matrix y, numeric scalar theta, numeric scalar h)
{
	real matrix V;
	real scalar h1, n, c;

	n=rows(y);
	c = sqrt(variance(y,1)) * 1.1;
	
	if (h == 0) {
		h1=(n^(-0.2))*c;
	}
	else {
		h1 = h;
	}
	
	V=(1/(h1)):*(normalden((theta * J(rows(y),1,1) - y):/h1));

	return(mean(V,1));
}

numeric matrix variance_machmat(numeric matrix q, numeric matrix gamma,
				numeric scalar b0, numeric scalar b1, numeric matrix y, numeric matrix  x1,
						numeric matrix x2, 
			numeric scalar l, numeric scalar c, numeric scalar ba,
			numeric matrix cov_step1, numeric scalar m, numeric matrix beta, numeric matrix u,
			numeric matrix y1a,
			numeric matrix x1b, numeric matrix theta,
			numeric matrix y1, numeric matrix x1a, numeric matrix x2a,
			numeric matrix s1)
/*	function for the calculation of the standard erros of the Machado and Matta technique with
	sample selection as in Albrecht, Van Vuuren and Vroman, 2004.

	In :	q		- 	range of quantile-levels to be computed (in range 0-1) - k x 1 vector
			gamma	- 	estimates of the first stage, exluding normalized coefficients
			b0 		- 	normalized constant in first stage
			b1 		- 	normalized regressor of continuous coefficient in first stage
			y		-	left hand side variable. Observations that are missing
						should be NaN.
			x1		-	right hand side variable of the first stage, including instruments
						and constant (in first column).
			x2 		-	right hand side variable of the second stage, excluding instruments
						but also with observations of those individuals that do
						not have an observed lhs.
			m 		- 	number of repetitions of the MM-method
			theta	- 	function value obtained from function machado()
			beta	-	computed levels of quantile regression at u1 (size: m x 1).
			u		-	levels at which quantile regressions ar estimated (size: m x 1)
			y1a		-	estimated levels for Machado and Mata method (size: m x 1)
			x1b		-	sampled x's for Machado and Mata method (size: m x k); k = columns of x2
			
	Out	:	function value	=	k x 1 vector containing the quantiles using the MM-method
								computed at the levels of q.

	WARNING: Unfortunately I was not able to speed this procedure a bit more. It can take
				quite a while even with moderate levels of m and sample size.
*/
{
	numeric scalar n, n1, j1, j2;
	numeric matrix var_hulp;
	
	n = rows(y);
	n1 = rows(y1a);	
	
	var_hulp = J(cols(beta)-l,cols(beta)-l,0);

	hulp_x = x1a[.,1..cols(x1a)-1];
	hulp_gamma = b1 \ gamma;
	z = hulp_x * hulp_gamma;

	hulp_nd = normalden(-z);
	hulp_nd1 = (J(rows(z),1,1)-normal(-z));
	lambda = hulp_nd :/  hulp_nd1;

	vlambda = lambda;
	
	for (i=2;i<=l;i++)
	{
		i1 = J(rows(lambda),1,i);
		vlambda = vlambda, (lambda:^i1);
	}

	x3 = x2a[.,1..cols(x2a)-1], vlambda, x2a[.,cols(x2a)];
	psi1 = x3' * x3;
	n2 = J(rows(psi1), cols(psi1), 1 / rows(x1a));
	psi1 = n2 :* psi1;
	
	dlambda = lambda:^2 - lambda :* z;
	lambda1 = mean(s1:^2, 1);
	
	for(j1=1;j1<=m;j1++) {;
		j1;
		bhat1 = beta[j1,.]';
		m1 = kerncv3(x3,y1,bhat1);
		delta1 = bhat1[rows(bhat1)-cols(vlambda)..rows(bhat1)-1];
		v11 = y1 - x3 * bhat1;
		v21 = y1 - x2a[.,2..cols(x2a)] * bhat1[2..rows(bhat1)-cols(vlambda)];
		q1 = u[j1];
		mu21 = quanc1(v21, q1 ,s1);
		v31 = v21 - J(rows(v21),1,mu21);
		lambda21 = berekenlambda2(z, s1, v31);

		for(j2=j1+1;j2<=m;j2++) {;												  
			cov1 = J(cols(beta)-l, cols(beta)-l,n1) :* compute_cov1(gamma, b0, b1,
								y1,x1a,x2a,l,c,ba,beta[j1,.]', beta[j2,.]', u[j1],  u[j2],
								n, cov_step1, s1, z, lambda, vlambda, x3, psi1, m1, dlambda, delta1,
								lambda1, v11, v21, mu21, lambda21, n2);
			var_hulp = var_hulp + cov1 + cov1;
		};
		
	};
							   
	for(j=1;j<=m;j++)
	{
		bhat1 = beta[j,.]';
		m1 = kerncv3(x3,y1,bhat1);
		delta1 = bhat1[rows(bhat1)-cols(vlambda)..rows(bhat1)-1];
		v11 = y1 - x3 * bhat1;
		v21 = y1 - x2a[.,2..cols(x2a)] * bhat1[2..rows(bhat1)-cols(vlambda)];
		q1 = u[j];
		mu21 = quanc1(v21, q1 ,s1);
		v31 = v21 - J(rows(v21),1,mu21);
		lambda21 = berekenlambda2(z, s1, v31);
		cov1 = J(cols(beta)-l, cols(beta)-l,n1) :* compute_cov1(gamma, b0, b1,
								y1,x1a,x2a,l,c,ba,beta[j,.]', beta[j,.]', u[j],  u[j],
								n, cov_step1, s1, z, lambda, vlambda, x3, psi1, m1, dlambda, delta1,
								lambda1, v11, v21, mu21, lambda21, n2);
		var_hulp = var_hulp + cov1;
	}

	hulp2 = J(rows(y1a),1,1);
	for(i=1;i<=rows(theta);i++)
	{
		c = quantile(y1a,0.75)-quantile(y1a,0.25)
		h1 =(n^(-0.4))*c;
		V = (1/(h1))*(normalden((theta[i] * hulp2 - y1a)/h1));
		V = V :* J(rows(V),cols(x1b),1);
		hulp = x1b :* V;
		if (i==1){
			var10 = mean(hulp, 1);
		}
		else {
			var10 = var10 \ mean(hulp, 1);			
		}
	}
	
	var_hulp = (var_hulp) :/ (m^2);
	
	
	for(i=1;i<=rows(q);i++)
	{
		f1a = kerncv2(y1a,theta[i],0);
		hulp1 = var10[i,.] * var_hulp * var10[i,.]';
		if (i == 1)
		{
			var1 = (((m / n) * q[i]*(1-q[i]) + (m / n1) * hulp1) / (f1a^2));
		}
		else
		{
			var1 = var1 \ (((m / n) * q[i]*(1-q[i]) + (m / n1) * hulp1) / (f1a^2));
		}
	}

	return (sqrt(var1 / m));	 
}

numeric matrix berekenlambda2(numeric matrix z, numeric matrix s1, numeric matrix v)
{
	numeric matrix sigma, f;
	numeric scalar h, n;
	
	n = rows(v);
	sigma = sqrt(variance(v,1));
	h = 1.06 :* sigma * n:^(-0.2);
	arg = v :/ h;
	
	f = ((1 / (h)) :* normalden(arg));
	return (mean(f :* s1, 1));	
}

numeric matrix kerncv3(numeric matrix X, numeric matrix y, numeric matrix bhat){
	real scalar c, h;
	real matrix V;

	r= y - X * bhat;

	n=rows(X);

	c = quantile(r,0.75)-quantile(r,0.25);
	h=(n^-(0.2))*c;

	V= (normalden(r/h):*X)'*X;
	nh = J(rows(V),cols(V),1 / (n*h)); 
	V = invsym(V) :* nh;

	return(V);
}

numeric matrix compute_cov1(numeric matrix gamma, numeric scalar b0,
				numeric scalar b1, numeric matrix w, numeric matrix x,
				numeric matrix x1, numeric scalar l, numeric scalar c, numeric scalar ba,
					numeric matrix bhat1, numeric matrix bhat2,
					numeric scalar q1, numeric scalar q2, numeric scalar n,
					numeric matrix cov, numeric matrix s1, numeric matrix z, numeric matrix lambda, numeric matrix vlambda,
					numeric matrix x3, numeric matrix psi1, numeric matrix m1, numeric matrix dlambda,
					numeric matrix delta1, numeric matrix lambda1, numeric matrix v11, numeric matrix v21,
					numeric scalar mu21, numeric matrix lambda21, numeric matrix n2)
{
/* Computes variance as in Albrecht, Van Vuuren and Vroman, 2004 for off-diagonal elements
	In:
	q1		- 	level of quantile for row of the off-diagonal element
	q2		-	level of quantile for column of the off-diagonal element (q1 < q2)
	bhat1	-	estimate as returned from compute_regression
	Other inputs are the same as in  compute_regression
	Out:
	function value =	covariance matrix of quantile regression. Thes size is the same as
						the number of columns contained in x1

	*/
	numeric matrix y, delta2, v12, v22, hulp_gamma,
							hulp_x, hulp_nd, hulp_nd1, i1, hulp11, hulp21,
							m2, j, psi2, psi;
	numeric scalar i, mu22, smu2, q3, k;


	y = w;

	delta2 = bhat2[rows(bhat2)-cols(vlambda)..rows(bhat2)-1];


	v12 = y - x3 * bhat2;
	v22 = y - x1[.,2..cols(x1)] * bhat2[2..rows(bhat1)-cols(vlambda)];
	
	n1 = rows(x);
	
	mu22 = quanc1(v22, q2 ,s1);
	
	v32 = v22 - J(rows(v22),1,mu22);

	lambda22 = berekenlambda2(z, s1, v32);

	smu2 = sqrt((lambda1 / (lambda21 * lambda22)) * q1 * (1-q2) / n);

	hulp11 = J(rows(lambda),1,0);
	hulp21 = hulp11;
	
	for(i=1;i<=cols(vlambda);i++)
	{
		i1 = J(rows(lambda),1,i-1);
		i2 = J(rows(lambda),1,i);
		hulp11 = hulp11 + J(rows(lambda),1,delta1[i]) :* dlambda :* (lambda:^(i1)) :* i2;
		hulp21 = hulp21 + J(rows(lambda),1,delta2[i]) :* dlambda :* (lambda:^(i1)) :* i2;
	}
	
	
	m2 = kerncv3(x3,y,bhat2);

	q3 = J(rows(psi1), cols(psi1), q1 * (1-q2));
	psi1 = q3 :* psi1;
	psi2 = bereken_psi2b(v11, v12, x[.,2..cols(x)-1], x3, hulp11, hulp21, cov);
	  
	psi = psi1 + (n1 / n) * psi2;
	
	cov1 = (m1 * psi * m2) :* n2;
	k = cols(x1);
	
	return (cov1[1..k, 1..k]);
}

numeric matrix bereken_psi2b(numeric matrix z1, numeric matrix z2,
				numeric matrix x1, numeric matrix x2, numeric matrix hulp1,
				numeric matrix hulp2, numeric matrix cov)
{
	numeric matrix psi2, arg1, arg2, f1, f2, hulpa, hulpd;
	numeric scalar i, h1, h2, c1, c2, n;
	
	r1 = z1;
	r2 = z2;

	n = rows(r1);

	c1 = sqrt(variance(r1,1)) * 1.1;
	c2 = sqrt(variance(r2,1)) * 1.1;

	h1=(n^-(0.4))*c1;
	h2=(n^-(0.4))*c2;
	
	arg1 = z1 / h1;
	arg2 = z2 / h2;

	f1 = (1 / h1) * normalden(arg1);
	f2 = (1 / h1) * normalden(arg2);

	
	hulp3 = f1 :* x2;
	hulp4 = f2 :* x2;
	
	xb1 = x1 :* hulp1;
	xb2 = x1 :* hulp2;
	
	k1 = cols(x2);
	
	for(j=1;j<=k1;j++)
	{
		if (j == 1)											   
		{
			hulp_1 = (xb2 :* (hulp4[.,j]))';
			hulp_2 = mean((cov * hulp_1)',1)';
			hulpa = hulp_2';
		}
		else
		{
			hulp_1 = (xb2 :* (hulp4[.,j]))';
			hulp_2 = mean((cov * hulp_1)',1)';
			hulpa = hulpa \ hulp_2';
		}
	}
	
	
	for(i=1;i<=k1;i++)
	{
		hulpd = mean(((hulp3[.,i]) :* xb1), 1);
		
		if (i == 1)
		{
			psi2 = (hulpd * hulpa');
		}
		else
		{
			psi2 = psi2 \ (hulpd * hulpa');
		} 
	}

	return (psi2);
}

end;

program qregb, eclass;
	version 9.0;
	syntax varlist(numeric) [if] , gender(varlist min=1 max=1)  select(varlist) [, m(int 100)] [, q(int 9)]
					[, l(int 2)] [, method(int 1)]
					[, variance(int 1)] [, ba(real 0.8)] [, c(real 0.2)]
					[, display(int 1)];


	tempfile hulp_file_overall;
	quietly save `hulp_file_overall', replace;
	
	local n = `m';
	
	if (`method' > 6) {;
		local method = 1;
	};

	if (`method' < 0) {;
		local method = 1;
	};

	mata uniformseed("Xc644806e61911463d1cc0e3d7be41beb01ba");
	
	mata hulp_u = uniform(`n',1) * 0.9998 + J(`n', 1, 0.0001);
	mata hulp_u = sort(hulp_u,1);
								 
	local vars3_hulp = "`varlist'";
	local i = 0;
	local k_hulp : list sizeof varlist;
	
	while (`i' < `k_hulp') {;
		local i = `i'+ 1;

		tempvar hulp_x;
		local name_x = word("`vars3_hulp'",1);
		
		quietly g `hulp_x' = `name_x';
		
		local vars3_hulp1 : list local(vars3_hulp) - local(name_x);
		local vars3_hulp "`vars3_hulp1'";

		if (`i' > 1) {;
			quietly drop if missing(`hulp_x');
		};
	};
								 
	local k_hulp : list sizeof select;
	local vars3_hulp = "`select'";

	local i = 0;
	
	while (`i' < `k_hulp') {;
		local i = `i'+ 1;

		tempvar hulp_x;
		local name_x = word("`vars3_hulp'",1);		
		g `hulp_x' = `name_x';
		
		local vars3_hulp1 : list local(vars3_hulp) - local(name_x);
		local vars3_hulp "`vars3_hulp1'";

		quietly drop if missing(`hulp_x');
	};							   
	
	local y = "`1'";
	tempvar fulltime;
	g `fulltime' = !missing(`y');
	local i = 0;
	local k : list sizeof varlist;
	local k = `k' - 1;

	local hulp1_varlist = word("`varlist'", 1);
	
	global vars : list local(varlist) - local(hulp1_varlist);

	global vars1 = "";
	global vars2 : list local(select) - global(vars);
	local k_iv : list sizeof select;

	local k1 = 0;

	global vars3 $vars $vars2;

	
	tempfile hulp_file4;

	quietly save `hulp_file4', replace;
	quietly keep if `gender' == 1;
	capture program drop single;
	capture program drop dsingle;
											  
	single `fulltime' $vars $vars2, h(0.2);
		
	use `hulp_file4', clear;
	matrix gamma = e(b);
	local k_totaal = colsof(gamma);
	matrix cov_step1 = e(V);
	local b1 = gamma[1,1];
	local b0 = gamma[1,colsof(gamma)];

	local i = 0;
	tempvar lambda;
	tempvar lambda1;
	tempvar z;
	
	g `z' = `b0';
	local i = 0;
	local vars3_hulp "$vars3";
	
	while (`i' < `k_totaal'-1) {;
		local i = `i'+ 1;
		local b1 = gamma[1,`i'];

		tempvar hulp_x;
		local name_x = word("`vars3_hulp'",1);
		g `hulp_x' = `name_x';
		local vars3_hulp1 : list local(vars3_hulp) - local(name_x);
		local vars3_hulp "`vars3_hulp1'";
		
		quietly replace `z' = `z' + `b1' * `hulp_x';
	};
	
	g `lambda' = normalden(-`z') / (1-normal(-`z'));
	g `lambda1' = 1;
	
	local i = 0;
	
	while (`i' < `l') {;
		local i = `i' + 1;
		tempvar vlambda_`i';
		quietly replace `lambda1' = `lambda1' * `lambda';
		quietly g `vlambda_`i'' = `lambda1';
		local vars_lambda `vars_lambda' `vlambda_`i'';
	};

	local hulp_increment = 1 / (`q'+1);
	
	if ((`method' == 2) | (`method' == 4)) {;
		local q1 = 0;

		while (`q1' < 1 - (2 * `hulp_increment')) {;
			local q1 = `q1' + `hulp_increment';
			quietly qreg `y' if gender == 1, q(`q1');
			matrix b1 = e(b);
			
			local hulp_list : list local(hulp_list) | local(q1);

			if (`q1' == `hulp_increment'){;
				matrix result3 = `q1', b1;
			};
			else {;
				matrix result3 = result3 \ (`q1', b1);
			};		
		};
/*
		quietly sqreg `y' if gender == 1, q(`hulp_list') reps(2);

		matrix result3 = result3, e(b)';*/
	};


	if ((`method' == 3) | (`method' == 5)) {;
		local q1 = 0;
		local i = 0;
		
		while (`q1' < 1 - (2 * `hulp_increment')) {;
			local q1 = `q1' + `hulp_increment';
			local i = `i' + 1;
			quietly qreg `y' if gender == 0, q(`q1');
			matrix b1 = e(b);
			if (`i' == 1){;
				matrix result3 = `q1', b1;
			};
			else {;
				matrix result3 = result3 \ (`q1', b1);
			};		
		};
	};

	if (`method' == 5) {;
		tempvar gender1;
		g `gender1' = 1 - gender;
		mata yb_gender0 = st_data(., "`y'", "`gender1'");
	};

	local i = 0;

	tempfile hulp_file2;
	quietly save `hulp_file2', replace;

	quietly keep if gender == 1;
	quietly keep if `fulltime' == 1;

	tempvar x;
	tempvar hulp;
	tempvar s;
	tempvar hulp_x;
	tempvar hulp_v;
	tempvar v2;

	tempvar y1;
	g `y1'= !missing(`y'); 
	tempvar fulltime1;
	g `fulltime1' = `fulltime' * `y1';

	
	quietly g `x' = `z' - `c';
	quietly g `hulp' = `x' / (`ba' - `x') * (`x' < `ba') + (`x' >= `ba');
	quietly g `s' = (1 - exp(-`hulp')) * ((`x' >= 0) * (`x' < `ba')) + (`x' >= `ba');

	local n = 10;
	while (`i' < `n') {;
		local i = `i' + 1;
		mata: st_matrix("hulp2", hulp_u[`i',1]);

		local hulp2 = hulp2[1,1];
		
		qreg `y' $vars `vars_lambda' if `fulltime' == 1 , q(`hulp2');
	
		matrix b = e(b);
		matrix b3 = e(b);
		matrix b = b[1,1..colsof(b)-`l'-1], b[1,colsof(b)];
		
		g `hulp_v' = 0;
		local j = 0;
		local vars3_hulp "$vars3";

		local k : list sizeof varlist;
		local k = `k' - 1;

		while (`j' < (`k' - `l')) {;
			local j = `j'+ 1;
			local b1 = b[1,`j'];

			tempvar hulp_x;
			local name_x = word("`vars3_hulp'",1);
			local vars3_hulp1 : list local(vars3_hulp) - local(name_x);
			local vars3_hulp "`vars3_hulp1'";
			quietly g `hulp_x' = `name_x';
			
			quietly replace `hulp_v' = `hulp_v' + `b1' * `hulp_x';
			drop `hulp_x';
		};

		quietly g `v2' = `y' - `hulp_v';				
		mata v2 = st_data(., "`v2'", "`fulltime1'");
		mata s = st_data(., "`s'", "`fulltime1'");
		mata q = hulp_u[`i',1];
		mata mu2 = quanc1(v2, q ,s);
		mata st_matrix("mu2", mu2);

		matrix b = b';
		matrix b3 = b3';
		
		matrix b = b[1..rowsof(b)-1,1] \ mu2;
		matrix b3 = b3[1..rowsof(b3)-1,1] \ mu2;
		  
		mata: b = st_matrix("b");
		mata: b3 = st_matrix("b3");
		
		if (`i'  == 1) {;
			mata: result = b';
			mata: result_b3 = b3';
		};
		else {;
			mata: result = result \ b';
			mata: result_b3 = result_b3 \ b3';
		};

		drop `hulp_v';
		drop `v2';
	};

	use `hulp_file2', clear;

	if (`method' < 5){;
		quietly keep if `gender' == 1;
	};
	else {;
		quietly keep if `gender' == 0;
	};

	set seed 339487731;
	
	quietly {;
	if (`n'< _N){;
		quietly bsample `n';
	};
	else
	{;
		tempfile hulp_file1;
		quietly save `hulp_file1', replace;
		local n1 = int(`n' / _N);
		local n2 = `n' - `n1' * _N;
		
		local i = 0;
		while (`i' < `n1') {;
			local i = `i' + 1;
			use `hulp_file1', clear;
			quietly bsample _N;
			tempfile hulp_file_`i';
			quietly save `hulp_file_`i'', replace;			
		};

		use `hulp_file1', clear;
		quietly bsample `n2';
		
		local i = 0;
		while (`i' < `n1') {;
			local i = `i' + 1;
			quietly append using `hulp_file_`i'';
		};	
	};
	};

	quietly g cons = 1;

	local i = 0;
	local vars3_hulp = "$vars3";

	while (`i' < `k') {;
		local i = `i' + 1;
		
		local x = word("`vars3_hulp'",1);
		local vars3_hulp1 : list local(vars3_hulp) - local(x);
		local vars3_hulp "`vars3_hulp1'";
			
		if (`i' == 1) {;
			mata X = st_data(., "`x'");
		};
		else{;
			mata X = X, st_data(., "`x'");
		};		
	};

	mata y = st_data(.,"`y'");
	
	local nobs = _N;
	mata: hulp_mat = J(`nobs',1,1);
	mata: X = X, hulp_mat;
	
	local i = 0;
	
	while (`i' < `n') {;
		local i = `i' + 1;
		mata: hulp = X[`i',.] * result[`i',.]';
		if (`i' == 1){;
			mata: result1 = hulp;
		};
		else {;
			mata: result1 = result1 \ hulp;
		};	   
	};

	drop _all;

	local nn = `n'/* * `n'*/;
	set obs `nn';
	tempvar hulp_var;
	quietly g `hulp_var' = .;
	
	local i = 0;

	while (`i' < (`n' /** `n'*/)) {;
		local i = `i' + 1;
		mata: st_store(`i',1,result1[`i',1]);
	};

	local q1 = 0;
	local i = 0;

	while (`q1' < 1 - (2 * `hulp_increment')) {;
		local q1 = `q1' + `hulp_increment';
		local i = `i' + 1;
		quietly qreg `hulp_var', q(`q1');
		matrix b = e(b);
		if (`i' == 1){;
			matrix result2 = `q1', b;
		};
		else {;
			matrix result2 = result2 \ (`q1', b);
		};		
	};


	if (`method' > 2) {;
		use `hulp_file2', clear;

		quietly keep if `gender' == 1;
		tempvar y1;
		g `y1'= !missing(`y'); 
		tempvar fulltime1;
		g `fulltime1' = `fulltime' * `y1';
	  	quietly keep if `fulltime1' == 1;

		if (`n'< _N){;
			quietly bsample `n';
		};
		else
		{;
			tempfile hulp_file1;
			quietly save `hulp_file1', replace;
			local n1 = int(`n' / _N);
			local n2 = `n' - `n1' * _N;
		
			local i = 0;
			while (`i' < `n1') {;
				local i = `i' + 1;
				use `hulp_file1', clear;
				quietly bsample _N;
				tempfile hulp_file_`i';
				quietly save `hulp_file_`i'', replace;			
			};

			use `hulp_file1', clear;
			quietly bsample `n2';
		
			local i = 0;
			while (`i' < `n1') {;
				local i = `i' + 1;
				append using `hulp_file_`i'';
			};	
		};

		tempvar cons;
		g `cons' = 1;

		local i = 0;
		local vars3_hulp = "$vars3";

		while (`i' < `k') {;
			local i = `i' + 1;
			local x = word("`vars3_hulp'",1);
			local vars3_hulp1 : list local(vars3_hulp) - local(x);
			local vars3_hulp "`vars3_hulp1'";
			if (`i' == 1) {;
				mata X1 = st_data(., "`x'");
			};
			else{;
				mata X1 = X1, st_data(., "`x'");
			};		
		};

		mata y1 = st_data(.,"`y'");

		local nobs = _N;
		mata: hulp_mat = J(`nobs',1,1);
		mata: X1 = X1, hulp_mat;
	
		local i = 0;
	
		while (`i' < `n') {;
			local i = `i' + 1;
			mata: hulp = X1[`i',.] * result[`i',.]';
			if (`i' == 1){;
				mata: result1_full = hulp;
			};
			else {;
				mata: result1_full = result1_full \ hulp;
			};
		};

		drop _all;

		set obs `n';
		tempvar hulp_var;
		g `hulp_var' = .;
	
		local i = 0;

		while (`i' < `n') {;
			local i = `i' + 1;
			mata: st_store(`i',1,result1_full[`i',1]);
		};

		local q1 = 0;
		local i = 0;

		while (`q1' < 1 - (2 * `hulp_increment')) {;
			local q1 = `q1' + `hulp_increment';
			local i = `i' + 1;
			quietly qreg `hulp_var', q(`q1');
			matrix b = e(b);
			if (`i' == 1){;
				matrix result2_full = `q1', b;
			};
			else {;
				matrix result2_full = result2_full \ (`q1', b);
			};		
		};
		matrix result2_1_full = result2_full[1..rowsof(result2_full),1];
		matrix result2_2_full = result2_full[1..rowsof(result2_full),2];
	};

	use `hulp_file2', clear;
	quietly keep if gender == 1;
	ereturn clear;

	matrix result2_1 = result2[1..rowsof(result2),1];
	matrix result2_2 = result2[1..rowsof(result2),2];
	mata q =  st_matrix("result2_1");
	mata theta1 = st_matrix("result3");


	if (`variance' == 1) {;
		mata n = `n';
		mata theta = st_matrix("result2_2");
		matrix b1 = gamma[1,1];
		matrix b0 = gamma[1,colsof(gamma)];
		matrix gamma = gamma[1, 2..colsof(gamma)-1]';
		mata b0 = st_matrix("b0");
		mata b1 = st_matrix("b1");
		mata gamma = st_matrix("gamma");

		tempvar y1;
		g `y1'= !missing(`y'); 
		mata beta = result_b3;
		mata y = st_data(.,"`y'");
		mata y1 = st_data(.,"`y'", "`y1'");
		mata yb = st_data(.,"`y'", "`gender'");
		mata cov_step1 = st_matrix("cov_step1");
		mata cov_step1 = cov_step1[2..cols(cov_step1)-1,2..cols(cov_step1)-1];

		tempvar fulltime1;
		g `fulltime1' = `fulltime' * `y1'; 

		local i = 0;

		local vars3_hulp = "$vars3";

		while (`i' < `k_totaal'-1) {;
			local i = `i' + 1;
			local x = word("`vars3_hulp'",1);
			local vars3_hulp1 : list local(vars3_hulp) - local(x);
			local vars3_hulp "`vars3_hulp1'";
			
			if (`i' == 1) {;
				mata x1 = st_data(., "`x'", "`fulltime'");
				mata x1a = st_data(., "`x'", "`fulltime1'");
			};
			else{;
				mata x1 = x1, st_data(., "`x'", "`fulltime'");
				mata x1a = x1a, st_data(., "`x'", "`fulltime1'");
			};		
		};

		mata: hulp_mat1 = J(rows(x1),1,1);
		mata: hulp_mat2 = J(rows(x1a),1,1);	
		mata x1 = x1, hulp_mat1;
		mata x1a = x1a, hulp_mat2;

		local i = 0;

		local vars_hulp = "$vars";

		local k_hulp : list sizeof vars_hulp;

		while (`i' < `k_hulp') {;
			local i = `i' + 1;
			local x = word("`vars_hulp'",1);
			local vars_hulp1 : list local(vars_hulp) - local(x);
			local vars_hulp "`vars_hulp1'";
			
			if (`i' == 1) {;
				mata x2 = st_data(., "`x'");
				mata x2a = st_data(., "`x'", "`y1'");
			};
			else{;
				mata x2 = x2, st_data(., "`x'");
				mata x2a = x2a, st_data(., "`x'", "`y1'");
			};
		};

		
		mata: hulp_mat3 = J(rows(x2),1,1);
		mata: hulp_mat4 = J(rows(x2a),1,1);	
		mata x2 = x2, hulp_mat3;
		mata x2a = x2a, hulp_mat4;
		tempvar x;
		tempvar hulp;
		tempvar s;
		quietly g `x' = `z' - `c';
		quietly g `hulp' = `x' / (`ba' - `x');
		quietly g `s' = (1 - exp(-`hulp')) * ((`x' >= 0) * (`x' < `ba')) + (`x' >= `ba');
		
		mata s = st_data(., "`s'", "`fulltime1'");

		if ((`method' ==  2) | (`method' ==  3) | (`method' == 5)) {;
			mata st_matrix("theta1", theta1[.,2]);
			matrix dtheta = theta1 - result2_2;
		};
							
		mata sd = variance_machmat(q, gamma, b0, b1, y, x1, x2, `l', `c', `ba', cov_step1, `n',
						beta, hulp_u, result1, X, theta, y1, x1a, x2a, s);

		
		if ((`method' == 3) | (`method' == 4)) {;
			mata theta_full = st_matrix("result2_2_full");		
			mata sd2 = variance_machmat(q, gamma, b0, b1, y, x1, x2, `l', `c', `ba', cov_step1, `n',
						beta, hulp_u, result1_full, X1, theta_full, y1, x1a, x2a, s);						
		};
		
		if (`method' == 2) {;
			mata var2 = bereken_var2(q, theta, theta1[.,2], yb, result1, n);
			mata sd1 = sqrt(abs(sd :* sd + var2 :/ n));
			mata st_matrix("sd", sd1);
		};

		if (`method' == 3) {;
			use `hulp_file2', clear;
			tempvar gender1;
			g `gender1' = 1 - `gender';
			mata yb = st_data(.,"`y'", "`gender1'");
			mata var2 = bereken_var2(q, theta, theta1[.,2], yb, result1, n);
			mata sd1 = sqrt(abs(sd :* sd + var2 :/ n));
			mata st_matrix("sd", sd1);
		};

		if (`method' == 1) {;
			mata st_matrix("sd", sd);
			ereturn matrix q = result2_1;
			ereturn matrix theta = result2_2;
			ereturn matrix sd = sd;
		};
								   
		if (`method' == 2) {;
			ereturn matrix q = result2_1;
			ereturn matrix theta = result2_2;
			ereturn matrix dtheta = dtheta;
			ereturn matrix sd = sd;
		};

		
		if (`method' == 3) {;
			ereturn matrix q = result2_1;
			ereturn matrix dtheta = dtheta;
			ereturn matrix sd = sd;
		};

		if (`method' == 5) {;
			mata var2 = bereken_var2(q, theta, theta1[.,2], yb_gender0, result1, n);
			mata sd1 = sqrt(abs(sd :* sd + var2 :/ n));
			ereturn matrix q = result2_1;
			ereturn matrix dtheta = dtheta;
			mata st_matrix("sd", sd1);
			ereturn matrix sd = sd;
		};


		if (`method'== 4) {;
			mata st_matrix("theta1", theta1[.,2]);
			matrix dtheta = theta1 - result2_2;
			matrix dtheta1 = theta1 - result2_2_full;
			matrix dtheta2 = result2_2_full - result2_2;
			
			mata var2 = bereken_var2(q, theta, theta1[.,2], yb, result1, n);
			
			mata sd = sqrt(abs(sd :* sd + var2 :/ n));
			mata sd_1 = sqrt(abs(sd2 :* sd2 + var2 :/ n));
			mata sd_2 = sqrt(abs(sd2 :* sd2 + sd :* sd));
			mata st_matrix("sd", sd);
			mata st_matrix("sd_1", sd_1);
			mata st_matrix("sd_2", sd_2);
 
			ereturn matrix q = result2_1;
			ereturn matrix dtheta = dtheta;
			ereturn matrix dtheta_obs = dtheta2;
			ereturn matrix dtheta_unobs = dtheta1;
			ereturn matrix sd = sd;
			ereturn matrix sd_obs = sd_1;
			ereturn matrix sd_unobs = sd_2;
		};

		if (`display' == 1) {;
			matrix q = e(q);
			if (!(`method' == 4)) {;
				display "";
				display "Estimation of the Machado-Mata method with sample selection correction";
				display "";
				display "------------------------------------------------------------------";
				display "";
				display "   quantile       point        standard       lower         upper";
				display "                estimate       deviation      bound         bound";
				display "";
				display "------------------------------------------------------------------";
				display "";

				local i = 0;
				local n_q = rowsof(q);
				while (`i' < `n_q') {;
					local i = `i' + 1;
					if (`method' == 1) {;
						matrix theta = e(theta);
						matrix sd = e(sd);
						local q_dis = q[`i',1];
						local theta_dis = theta[`i',1];
						local sd_dis = sd[`i',1];
						local lower = `theta_dis' - 2 * `sd_dis';
						local upper = `theta_dis' + 2 * `sd_dis';
						display _skip(2) %10.4f `q_dis' "   " %10.4f `theta_dis' "   " %10.4f `sd_dis'
										"    " %10.4f `lower' "    " %10.4f `upper';
					};
					else {;
						matrix dtheta = e(dtheta);
						matrix sd = e(sd);
						local q_dis = q[`i',1];
						local dtheta_dis = dtheta[`i',1];
						local sd_dis = sd[`i',1];
						local lower = `dtheta_dis' - 2 * `sd_dis';
						local upper = `dtheta_dis' + 2 * `sd_dis';
						display _skip(2) %10.4f `q_dis' "   " %10.4f `dtheta_dis' "   " %10.4f `sd_dis'
									"    " %10.4f `lower' "    " %10.4f `upper';
					};
				};
			};
			else {;
				display "";
				display "Estimation of the Machado-Mata method with sample selection correction";
				display "                    Selection on observables    ";
				display "";
				display "------------------------------------------------------------------";
				display "";
				display "   quantile       point        standard       lower         upper";
				display "                estimate       deviation      bound         bound";
				display "";
				display "------------------------------------------------------------------";
				display "";

				local i = 0;
				local n_q = rowsof(q);
				while (`i' < `n_q') {;
					local i = `i' + 1;
					matrix dtheta2 = e(dtheta_obs);
					matrix sd2 = e(sd_obs);
					local q_dis = q[`i',1];
					local dtheta_dis2 = dtheta2[`i',1];
					local sd_dis2 = sd2[`i',1];
					local lower2 = `dtheta_dis2' - 2 * `sd_dis2';
					local upper2 = `dtheta_dis2' + 2 * `sd_dis2';
					display _skip(2) %10.4f `q_dis' "   " %10.4f `dtheta_dis2' "   " %10.4f `sd_dis2'
									"    " %10.4f `lower2' "    " %10.4f `upper2';
				};

				display "";
				display "                    Selection on unobservables    ";
				display "";
				display "------------------------------------------------------------------";
				display "";
				display "   quantile       point        standard       lower         upper";
				display "                estimate       deviation      bound         bound";
				display "";
				display "------------------------------------------------------------------";
				display "";

				local i = 0;
				local n_q = rowsof(q);
				while (`i' < `n_q') {;
					local i = `i' + 1;
					matrix dtheta1 = e(dtheta_unobs);
					matrix sd1 = e(sd_unobs);
					local q_dis = q[`i',1];
					local dtheta_dis1 = dtheta1[`i',1];
					local sd_dis1 = sd1[`i',1];
					local lower1 = `dtheta_dis1' - 2 * `sd_dis1';
					local upper1 = `dtheta_dis1' + 2 * `sd_dis1';
					display _skip(2) %10.4f `q_dis' "   " %10.4f `dtheta_dis1' "   " %10.4f `sd_dis1'
									"    " %10.4f `lower1' "    " %10.4f `upper1';
				};
			};
		};	
	};
	else{;
		if (`method' == 1) {;
		
			ereturn matrix q = result2_1;
			ereturn matrix theta = result2_2;
		};
						
		if ((`method' ==  2) | (`method' ==  3)) {;
			mata st_matrix("theta1", theta1[.,2]);
			matrix dtheta = theta1 - result2_2;
			ereturn matrix theta1 = theta1;
			ereturn matrix theta = result2_2;
			ereturn matrix q = result2_1;
			ereturn matrix dtheta = dtheta;
		};

		if (`method' ==  4) {;
			mata st_matrix("theta1", theta1[.,2]);
			matrix dtheta = theta1 - result2_2;
			matrix dtheta1 = theta1 - result2_2_full;
			matrix dtheta2 = result2_2_full - result2_2;
			ereturn matrix q = result2_1;
			ereturn matrix dtheta = dtheta;
			ereturn matrix dtheta_obs = dtheta2;
			ereturn matrix dtheta_unobs = dtheta1;
		};
		
		if (`display' == 1) {;
			matrix q = e(q);
			display "";
			display "Estimation of the Machado-Mata method    ";
			display "";
			display "----------------------------";
			display "";
			display "   quantile       point        ";
			display "                 estimate      ";
			display "";
			display "----------------------------";
			display "";
			local i = 0;
			local n_q = rowsof(q);
			while (`i' < `n_q') {;
				local i = `i' + 1;
				if (`method' == 1) {;
					matrix theta = e(theta);
					matrix sd = e(sd);
					local q_dis = q[`i',1];
					local theta_dis = theta[`i',1];
					display _skip(2) %10.4f `q_dis' "   " %10.4f `theta_dis';
				};
				else {;
					matrix dtheta = e(dtheta);
					matrix sd = e(sd);
					local q_dis = q[`i',1];
					local dtheta_dis = dtheta[`i',1];
					display _skip(2) %10.4f `q_dis' "   " %10.4f `dtheta_dis';
				};
			};
		};			
	};
	
	quietly use `hulp_file_overall', clear;
end;