****************************************************************************************
* Program generating wage data for males and females and estimating the access function
*
* Log-wages of males and females are supposed to follow normal laws with mean and 
* variance equal to those in our real data 
*
* Two graphs (one with confidence intervals and one without confidence intervals) are 
* also provided
*
* For an application, one just has to replace the simulated data by real data 
****************************************************************************************

#delimit ;

clear matrix;
clear mata;
clear;

set mem 500m ;

set more off;

local path  "C:\dossiers\wage_discri\DADS\FilesToTransfer\ProgramsToTransfer\";		/* Path containing the procedures for computation of access function */

do "`path'kernelhu.do"; 	/* File containing all the procedures */

local date="131026";			/* Date used to index generated files */

local nbp=100;					/* Number of ranks at which the access function is computed */

local nbsim=100;				/* Number of simulations for bootstrap */

local lsal ="lsal";				/* Log-wage variable */
local sexe ="sexe";				/* Sex variable: 1: male, 2: female */


************************************************************;
*                 Data selection / data simulation         *;
************************************************************;

* Next lines generate wage data and can be replace by the selection of real data;

local Nm=10000; /* Number of males   */
local Nf=2240;	/* Number of females; le proportion of females in our sample is 22.4% */

local tot=`Nf'+`Nm';

set obs `tot';

gen     sexe=1 if _n<=`Nm';
replace sexe=2 if _n>`Nm';

gen		  lsal=rnormal(4.79,.386) if _n<=`Nm';
replace   lsal=rnormal(4.57,.304) if _n>`Nm';


* To be kept in any case;

keep if missing(`sexe')==0 & missing(`lsal')==0;
keep `sexe' `lsal';

save "`path'data_select_`date'tot", replace;



**************************************************;
*                 Simulations                    *;
**************************************************;

foreach s of numlist 0/`nbsim' {;	/* Number of bootstrap simulations; if nbsim=0 no bootstrap */

* Initial sample used when s=0;

if `s'==0 {;
	use "`path'data_select_`date'tot", clear;
};

* Simulation of a sample for each simulation s>0 by drawing n observations with replacement in the initial sample;

if `s'>0 {; 
	use "`path'data_select_`date'tot", clear;
	gen numobs=int(uniform()*_N)+1;
	local nbvar=wordcount("`lsal' `sexe'");
	foreach i of numlist 1/`nbvar' {;
    local nvar=word("`lsal' `sexe'",`i');
	gen `nvar'b=`nvar'[numobs[_n]];
	drop `nvar';
	rename `nvar'b `nvar';
	};
};

egen nbh=sum((`sexe'==1));
egen nbf=sum((`sexe'==2));
gen ntot=nbh+nbf;

sort `lsal';
gen rangt=_n/ntot;

keep rangt `sexe' `lsal';

sort `lsal';
local N=min(_N,`nbp');

* Use of procedures stored in kernelh120326.do;

mata: calcquant("`lsal'","`sexe'","epan",0.5,`N',"lsaljourh","lsaljourf","rang");
mata: calckh("`lsal'","`sexe'","epan",2,`N');

keep if _n<=`N';
sort rang;

gen s=`s';

keep hu rang s;

if `s'>0 {; append using "`path'hu_boot_`nvarh'_`date'.dta"; };
saveold "`path'hu_boot_`nvarh'_`date'.dta", replace; 
};

sort s rang;
saveold "`path'hu_boot_`nvarh'_`date'.dta", replace; 

erase "`path'data_select_`date'tot.dta";



/* Two graphs : one with confidence interval; one without confidence interval */

use "`path'hu_boot_`nvarh'_`date'.dta", replace;

keep if s==0;

sort rang;
saveold "hu_boot_`nvarh'_`date'_shaped", replace;

use "`path'hu_boot_`nvarh'_`date'.dta", clear;
sort rang hu;

collapse (p5) hu_p5=hu (p95) hu_p95=hu, by(rang);

sort rang;

merge rang using "hu_boot_`nvarh'_`date'_shaped";
drop _merge;

sort rang;
saveold "hu_boot_`nvarh'_`date'_shaped", replace;

set scheme s2mono;

gen ligne1=1;

twoway line hu rang if rang>0.01 & hu<=1.6,
   xtitle("Rank") ytitle("Access function") legend(off) xscale(range(0 1)) xmtick(0(.2)1) xlabel(0(.2)1) 
   yscale(range(0 1.6)) ymtick(0(.2)1.6) ylabel(0(.2)1.6) legend(off)
   lpattern(solid dash) lcolor(black black) graphregion(color(white));

graph export "`path'\hu.ps", as(ps) pagesize(letter) orientation(landscape) mag(180) logo(off) tmargin(0.5) lmargin(0.5) replace;

twoway line hu hu_p5 hu_p95 ligne1 rang if rang>0.01 & hu<=1.6 & hu_p5<=1.6 & hu_p95<=1.6,
   xtitle("Rank") ytitle("Access function") legend(off) xscale(range(0 1)) xmtick(0(.2)1) xlabel(0(.2)1) 
                               yscale(range(0 1.6)) ymtick(0(.2)1.6) ylabel(0(.2)1.6) legend(off)
			 	   	           lpattern(solid shortdash shortdash dash) lcolor(black black black black) graphregion(color(white));

graph export "`path'\hu_CI.ps", as(ps) pagesize(letter) orientation(landscape) mag(180) logo(off) tmargin(0.5) lmargin(0.5) replace;

erase "hu_boot_`nvarh'_`date'_shaped.dta";