clear all
* 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

xtdpdml_simul, nreps(1000) n(100) t(4)