* 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