program define ar_sim
	version 6.0
	if "`1'" == "?" {
	global S_1 "x1 x2 x3 x4 cons se_x1 se_x2 se_x3 se_x4 se_cons pv_x1 pv_x2 pv_x3 pv_x4 pv_cons rho"
	set seed 1234567
		exit
	}
	drop _all
	
	local prime 5000
	local M 1000
	local T 5
	local unbal yes
	local unbaltb yes /* may lose from both top and bottom */
			/* required s.t. nr_ij != min(n_i, n_j) */

	local unbaltb = "`unbal'"=="yes" & "`unbaltb'"=="yes"

	local obs = (("`unbal'"=="yes") + 1)*(`unbaltb'+1)*2*`M'*`T'

	/* form an AR1 u_t */

	local tobs = `prime' + `obs'
	set obs `tobs'
	gen int t = _n
	tsset t , daily

	gen double e_t = 10*invnormal(uniform())
	gen double u_t = . in 1
	_byobs {

		update u_t =  /*
			   --------AR----------------------
			*/ 0.3*l.u_t + e_t, missval(0)  
		
	}

	keep if _n > `prime'

	/* form balanced panels */
	drop t
	gen long id = int((_n - 1) / `T')
	sort id
	qui by id:  gen t = _n

	/* unbalance the panels if requested */

	if "`unbal'" == "yes" {
		gen int keep_t = int(uniform()*5) + 1
		qui by id: keep if t <= keep_t[1]
		drop keep_t
	}

	if `unbaltb' {
		gen int keep_t = int(uniform()*5) + 1
		qui by id: keep if t >= keep_t[_N]
	}
		
	
	/* Keep every other panel
	*/
	
	keep if mod(id, 2) == 0 
	tsset id t
	
	
	gen double x1 = 2*uniform()
	gen double x2 = x1^2
	gen double x3 = 2*uniform() - 1
	gen double x4 = id / (3*`M') + uniform()
	
	gen double y =  1*x1 + 2*x2 + 3*x3 + 4*x4 - 1  + u_t

	xtgee y x1 x2 x3 x4, corr(ar1) 

	mat Rt=e(R)
	tempname p1 p2 p3 p4 pcons
	
	scalar `p1'=(_b[x1]-1)/_se[x1]
	scalar `p1'=1-2*(normprob(abs(`p1'))-.5)
	
	scalar `p2'=(_b[x2]-2)/_se[x2]
	scalar `p2'=1-2*(normprob(abs(`p2'))-.5)
	
	scalar `p3'=(_b[x3]-3)/_se[x3]
	scalar `p3'=1-2*(normprob(abs(`p3'))-.5)
	
	scalar `p4'=(_b[x4]-4)/_se[x4]
	scalar `p4'=1-2*(normprob(abs(`p4'))-.5)
		
	scalar `pcons'=(_b[_cons]+1)/_se[_cons]
	scalar `pcons'=1-2*(normprob(abs(`pcons'))-.5)
	
	post `1' _b[x1]  _b[x2]  _b[x3]  _b[x4]  _b[_cons] _se[x1] /* 
		*/_se[x2] _se[x3] _se[x4] _se[_cons] `p1' `p2' `p3' /*
		*/ `p4' `pcons' Rt[2,1] 
end