// Simulation modeled on Blundell-Bond (1998) simulation A

cap program drop simBB
program define simBB, rclass
	syntax [, alpha(real 0.8) sig_u(real 1) sig_mu(real 1) rho(real 0)]

	gen mu = `sig_mu' * invnormal(uniform()) if t == 0
	by id: replace mu = mu[1]
	gen u = `sig_u' * (sqrt(1-(`rho')^2) * invnormal(uniform()) + `rho' * mu/`sig_mu') if t==0
	gen v = invnormal(uniform()) 
	gen y = mu/(1-`alpha') + u if t == 0
	replace y = `alpha' * L.y + mu + v if t > 0

	foreach laglim in . 1 {
		foreach c in "" c {
			xtabond2 y L.y, gmm(l.y, `c' lag(1 `laglim')) h(2) two nodiff artest(0)
			return scalar alpha`c'`limname' = _b[L.y]
			return scalar j`c'`limname' = e(j)
			return scalar hansenp`c'`limname' = e(hansenp)
		}
		local limname 1	
	}
	drop mu v y u
end

cap program drop runBB
program define runBB
	clear
	syntax [, ITERations(integer 100) N(integer 100) T(integer 7) alpha(real 0.8) sig_u(real 1) sig_mu(real 1) rho(real 0)]

	set obs `=`t' * `n''
	gen t = mod(_n - 1, `t')
	gen id = int((_n - 1) / `t')
	tsset id t
	
	simulate alpha=r(alpha) hansenp=r(hansenp) j=r(j) alphac=r(alphac) hansenpc=r(hansenpc) jc=r(jc) ///
		alpha1=r(alpha1) hansenp1=r(hansenp1) j1=r(j1) alphac1=r(alphac1) hansenpc1=r(hansenpc1) jc1=r(jc1), reps(`iterations'): ///
		simBB, sig_u(`sig_u') sig_mu(`sig_mu') alpha(`alpha') rho(`rho')
	sum, sep(0)
end

mata: mata set matafavor speed
set seed 0
set more off

cap mat drop results
forvalues rho=0(0.3)0.9 {
	forvalues t=3/20 {
		runBB, alpha(0.8) t(`t') iter(500) rho(`rho') sig_u(2)
		collapse * (sd) alphasd=alpha alphasdc=alphac alphasd1=alpha1 alphasdc1=alphac1
		mkmat *, mat(a)
		mat results = nullmat(results) \ (`rho', `t', a)
	}
}
mat list results