//////////////////////
// Miguel Portela  	//
// Feb 2017			//
// Panel Simulation	//
// 2 fixed effects	//
//////////////////////

	clear all
	set more off
	set rmsg on
	
	timer on 1
	
// 	quiet {
	

capture log close
log using panel_simulation_2fe_cycle.txt, text replace

	mat def mat_seed = .
	mat def mat_size = .
	mat def mat_beta_educ = .

// 1st seed

	set seed 444555666	// specify initial value of random-number seed

// start here ///

// sample size

	local nnn = "50000 100000"		// we will work with two sample sizes
	tokenize `nnn'

		foreach nobs of local nnn {

		// number of simulations for each sample size
		foreach v of numlist 1/10 {

		clear

		// renew 'seed'

				set obs 1
				gen a = int(uniform()*1000000000)
				format %16.0f a
				l a
				sum a
				drop a
				ret li
					local seed2 = r(mean)

					di _new _new "/////////////////////////////////////////////////////////////////////////////" _new
					display _new _new _new "PARAMETERS:	SIZE: `nobs' SEED: `seed2'" _new _new
					di _new "/////////////////////////////////////////////////////////////////////////////" _new

		// generate the variables relevant for the estimations


		set obs `nobs'		// set the number of observation in the sample

		// define worker's id and year of observation
		gen workerid = int(_n/10.1)+1
		bysort workerid: gen ano = _n

		sort workerid ano
		tsset workerid ano

		set seed `seed2'

		gen a = uniform()
			sort a
			drop a

		gen firmid = int(_n/70.1)+1

		// generate individual's specific effect, ui
		bysort workerid (ano): gen double ui = runiform() if _n == 1
			 bysort workerid (ano): replace ui = ui[_n-1] if _n > 1

		// build the instrument, q1
		bysort workerid (ano): gen quarter= int(4*uniform()+1) if _n == 1
			 bysort workerid (ano): replace quarter = quarter[_n-1] if _n > 1
			 
			gen q1=(quarter==1)

		// build education and experience; education is a function of ui, the worker specific effect, as well as the instrument, q1
		bysort workerid (ano): gen educ= int(16*uniform())*ui if _n == 1
			bysort workerid (ano): replace educ= educ - 3*q1 if _n == 1
			
			bysort workerid (ano): replace educ = educ[_n-1] + round(uniform()) if _n > 1
			egen sd = sd(educ),by(workerid)
				sum sd, detail
				drop sd

		bysort workerid (ano): gen exper= int(20*uniform()) if _n == 1
			bysort workerid (ano): replace exper = exper[_n-1] + 1 if _n > 1
			
			gen exper2 = exper^2
			
		// firm fixed effect
		bysort firmid (ano): gen double tmp = 3*runiform() if _n == 1

			egen fi = max(tmp),by(firmid)
				drop tmp
			
			egen tmp = mean(ui),by(firmid)
			
			replace fi = fi + 9.3*tmp
				drop tmp
			
		// generate the true log wage equation, where the error term is defined by 'uniform()/ln(485)'; 485 is the minimum wage
		gen double lnwage = ln(485) + 0.06*educ + 0.007*exper - 0.00001*exper2 + 0.02*ano + ui + 0.4*fi + uniform()/ln(485)

		// # 1. OLS

		reg lnwage educ exper exper2
		reg lnwage educ exper exper2 i.ano ui
		reg lnwage educ exper exper2 i.ano ui fi

		// # 2. REGHDFE

		reghdfe lnwage educ exper exper2 i.ano,absorb(ui fi)


			mat mat_size = mat_size\(`nobs')
			mat mat_seed = mat_seed\(`seed2')
			mat mat_beta_educ = mat_beta_educ\_b[educ]
			
		}
		}

clear

// make the graphs with the empirical distribution of the parameters
	// transform matrices in data/variables

	svmat mat_seed, n(mat_seed)
	svmat mat_size, n(mat_size)
	svmat mat_beta_educ, n(mat_beta_educ)

	format %16.0f mat_seed mat_size
	format %5.4f mat_beta_educ
	drop if mat_size == .
	sum
	save data_beta_educ, replace

use data_beta_educ, clear

// graph the distribution of betas
kdensity mat_beta_educ, xline(.06)
	twoway kdensity mat_beta_educ1 if mat_size1 == `1',yaxis(1) legend(on label(1 "Sample Size: `1'"))/*
	*/ ytitle("Returns to E ducation") || kdensity mat_beta_educ1 if mat_size1 == `2', yaxis(2) legend(on label(2 "Sample Size: `2'")) /*
	*/ xline(.06) ytitle("Returns to E ducation",axis(2))

	graph export beta_educ1.png, replace

// end the `quiet' command

// }
	
timer off 1
timer list 1

log close