* Code from Williams, Allison, Moral-Benito paper
* Simulations Program
* Based on design 5 in table 1 of moral-benito (2013)
program xtdpdml_simul
	*** Section 4.1, Comparisons with A-B, Simulated data
	set more off
	clear all
	version 13.1, user

	* Default values: seed = 3, bylab = .75, bx = .25, t = 4, nreps = 10, N = 100
	syntax, [SEEDnum(integer 3) bylag(real .75) bx(real .25) t(integer 4) ///
		ONLYConverged nreps(integer 10) N(integer 100)]

	quietly {
		local N `n'
		* The following are currently fixed but could be changed.
		local b3 .5
		local b4 -0.17
		local g 1
		local h 1
		local f .67

		set matsize 10000
		* reps = the number of simulations
		local reps `nreps'
		set seed `seednum'

		mat rxtdpdml = J(`reps',2,.)
		mat rxtdpd = J(`reps',2,.)
		local it = 1
		while `it' <= `reps'{

			clear
			set obs `N' /* this is N (units in the panel */
			gen id=_n
			gen c=rnormal()*1.72

			* initial observations under mean stationarity: eqs (29) and (30) page 457 Moral-Benito (2013)
			gen y1= ( (`bx'*`f'+(1-`b3')) / ((1-`b3')*(1-`bylag')-`bx'*`b4') ) * c + `g'*rnormal()
			gen x1= ( (`b4'+`f'*(1-`bylag')) / ((1-`b3')*(1-`bylag')-`bx'*`b4') ) * c + `g'*rnormal()*2.56

			forval tnum = 2/`t' {
				local tlag = `tnum' - 1
				gen x`tnum'= `b3'*x`tlag' + `b4'*y`tlag' +`f'*c + `g'*rnormal()*2.56
				gen y`tnum'= `bx'*x`tnum' + `bylag'*y`tlag' +`h'*c + `g'*rnormal()
			}

			reshape long y x, i(id) j(t)
			xtset id t

			di "`it'"

			/* xtdpdml */
			capture xtdpdml y, pre(x) skipcfatransform skipconditional constinv

			local mlconverged = e(converged)
			if e(converged)== 1{
				mat bylag = _b[y2:y1]
				mat bx = _b[y2:x2]
				mat rxtdpdml[`it',1] = bylag
				mat rxtdpdml[`it',2] = bx
			}

			/* xtdpd */
			*Could also use xtdpd y L.y x, dgmmiv(x L.y , lag(1 ))
			xtabond y, pre(x)
			if "`onlyconverged'" == "" | `mlconverged' == 1 {
				mat bylag = _b[L.y]
				mat bx = _b[x]
				mat rxtdpd[`it',1] = bylag
				mat rxtdpd[`it',2] = bx
			}

			local it = `it' + 1
		}

		svmat rxtdpdml
		svmat rxtdpd

		gen biasyml = rxtdpdml1 - `bylag'
		gen biasxml = rxtdpdml2 - `bx'

		gen biasygmm = rxtdpd1 - `bylag'
		gen biasxgmm = rxtdpd2 - `bx'

	}

	display
	display as result "N = `N', T = `t', bylag = `bylag', bx = `bx', `reps' simulations
	sum biasyml biasygmm biasxml biasxgmm
	display
	quietly sum biasyml
	display "Root Mean Squared Error of lagged y ml  = " sqrt(r(mean)^2 + r(Var))
	quietly sum biasygmm
	display "Root Mean Squared Error of lagged y gmm = " sqrt(r(mean)^2 + r(Var))
	quietly sum biasxml
	display "Root Mean Squared Error of x ml         = " sqrt(r(mean)^2 + r(Var))
	quietly sum biasxgmm
	display "Root Mean Squared Error of x gmm        = " sqrt(r(mean)^2 + r(Var))
	display

end