set seed 23 set more off clear all cd "Y:\Allison_Williams_MoralBenito_DPD\Simulations" * cd "/Volumes/ENRIQUE/Allison_Williams_MoralBenito_DPD/Simulations" * First simulate the data * design 5 in table 1 of moral-benito (2013) local b2 .75 local b3 .5 local b1 .25 local b4 -0.17 local g 1 local h 1 local f .67 set matsize 10000 local reps = 500 mat rxtdpdml = J(`reps',2,.) mat rxtdpd = J(`reps',2,.) local it = 1 while `it' <= `reps'{ clear set matsize 10000 set obs 200 /* 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 y0= ( (`b1'*`f'+(1-`b3')) / ((1-`b3')*(1-`b2')-`b1'*`b4') ) * c + `g'*rnormal() gen x0= ( (`b4'+`f'*(1-`b2')) / ((1-`b3')*(1-`b2')-`b1'*`b4') ) * c + `g'*rnormal()*2.56 * T = 4 gen x1= `b3'*x0 + `b4'*y0 +`f'*c + `g'*rnormal()*2.56 gen y1= `b1'*x1 + `b2'*y0 +`h'*c + `g'*rnormal() gen x2= `b3'*x1 + `b4'*y1 +`f'*c + `g'*rnormal()*2.56 gen y2= `b1'*x2 + `b2'*y1 +`h'*c + `g'*rnormal() gen x3= `b3'*x2 + `b4'*y2 +`f'*c + `g'*rnormal()*2.56 gen y3= `b1'*x3 + `b2'*y2 +`h'*c + `g'*rnormal() /* * T = 6 gen x4= `b3'*x3 + `b4'*y3 +`f'*c + `g'*rnormal()*2.56 gen y4= `b1'*x4 + `b2'*y3 +`h'*c + `g'*rnormal() gen x5= `b3'*x4 + `b4'*y4 +`f'*c + `g'*rnormal()*2.56 gen y5= `b1'*x5 + `b2'*y4 +`h'*c + `g'*rnormal() * T = 8 gen x6= `b3'*x5 + `b4'*y5 +`f'*c + `g'*rnormal()*2.56 gen y6= `b1'*x6 + `b2'*y5 +`h'*c + `g'*rnormal() gen x7= `b3'*x6 + `b4'*y6 +`f'*c + `g'*rnormal()*2.56 gen y7= `b1'*x7 + `b2'*y6 +`h'*c + `g'*rnormal() * T = 10 gen x8= `b3'*x7 + `b4'*y7 +`f'*c + `g'*rnormal()*2.56 gen y8= `b1'*x8 + `b2'*y7 +`h'*c + `g'*rnormal() gen x9= `b3'*x8 + `b4'*y8 +`f'*c + `g'*rnormal()*2.56 gen y9= `b1'*x9 + `b2'*y8 +`h'*c + `g'*rnormal() * T = 12 gen x10= `b3'*x9 + `b4'*y9 +`f'*c + `g'*rnormal()*2.56 gen y10= `b1'*x10 + `b2'*y9 +`h'*c + `g'*rnormal() gen x11= `b3'*x10 + `b4'*y10 +`f'*c + `g'*rnormal()*2.56 gen y11= `b1'*x11 + `b2'*y10 +`h'*c + `g'*rnormal() */ reshape long y x, i(id) j(t) xtset id t /* if first observation is t=0, xtdpdml does not use it but xtdpd does!! */ replace t = t + 1 di "`it'" /* generating missing data i.e. unbalancedness */ /* missing at random */ gen p = invlogit(0.5*x+rnormal()) qui sum p, det replace x = . if p < `r(p10)' replace y = . if p < `r(p10)' /* xtdpdml with fiml option for unbalanced simulations */ xtdpdml y, pre(x) fiml mat b1 = _b[y3:y2] mat b2 = _b[y3:x3] if e(converged)==1{ mat rxtdpdml[`it',1] = b1 mat rxtdpdml[`it',2] = b2 } /* xtdpd */ xtdpd y L.y x, dgmmiv(x L.y , lag(1 )) mat b1 = _b[L.y] mat b2 = _b[x] mat rxtdpd[`it',1] = b1 mat rxtdpd[`it',2] = b2 local it = `it' + 1 } svmat rxtdpdml svmat rxtdpd gen bias1ml = rxtdpdml1 - 0.75 gen bias2ml = rxtdpdml2 - 0.25 gen bias1gmm = rxtdpd1 - 0.75 gen bias2gmm = rxtdpd2 - 0.25 set more off log using result_simulation_unbalanced.log, replace sum bias1ml bias2ml bias1gmm bias2gmm, det sum rxtdpdml1 rxtdpdml2 rxtdpd1 rxtdpd2, det log close