/********************************************************** ** MONTE CARLO EXPERIMENTS FOR DYNAMIC PANEL DATA MODELS ** ** WRITTEN BY S. C. AHN ** ***********************************************************/ @ CLEAN MEMMORY @ new ; @ OPEN OUTPUT FILE @ output file = dynamic2.out reset ; @ FORMAT OUTPUT FILE @ format /rd 8,3 ; /****************************************** ** Starting comments on the output file ** ******************************************/ " MODEL SPECIFICATION: " ; " y_t = d*y_{t-1} + a + e_t ;" ; " y_0 = gamma*a + e_0 " ; " Non-Stationary y: gamma = 1/(1-d)+1 " ; /**************************** ** Data Generating Process ** ****************************/ @ For the size of N @ j1 = 2 ; let jj[4,1] = 50 100 300 500 ; do while j1 <= 2 ; jj1 = jj[j1,1] ; @ For the size of T @ j2 = 1 ; let jjj[3,1] = 4 7 11 ; do while j2 <= 1 ; jj2 = jjj[J2,1] ; @ For the size of delta @ j3 = 1 ; let jjjj[3,1] = 0.9 0.8 0.5 ; do while j3 <= 3 ; jj3 = jjjj[J3,1] ; @ For siga @ j4 = 3 ; let jjjjj[3,1] = 0 0.5 1 ; do while j4 <= 3 ; jj4 =jjjjj[J4,1] ; @ Values of major parameters @ ddd = 1 ; @ SEED @ iter = 100 ; @ # OF ITERATIONS @ tt = jj2 ; @ TIME HORIZON @ nn = jj1 ; @ # OF CROSS-SECTIONAL UNITS @ dd = jj3 ; @ DELTA @ se = 1 ; @ SE(E) @ sa = jj4 ; @ SE(A) @ rr = 1/(1-dd)+1 ; @ GAMMA @ s0 = sqrt(se^2/(1-dd^2)) ; @ SE(Y0) @ @ Printing out the parameter values in the output file @ "" ; "============================" ; "VALUES OF MAJOR PARAMETERS" ; "SEED = " ddd ; "ITER = " iter ; "T = " tt ; "N = " nn ; "DELTA = " dd ; "SIGA = " sa ; "SIGE = " se ; "GAMMA(1/(1-d)+1) = " rr ; "SIG0 = " s0 ; @ Vectors to store estimators from each iterations @ @ ARELLANO-BOND GMM @ abb = zeros(iter,1) ; @ AB ESTIMATE @ abv = zeros(iter,1) ; @ AB VARIANCE @ abss= 0 ; @ SIZE OF T-TEST @ @ ARELLANO-BOVER GMM @ abob = zeros(iter,1) ; @ ABo ESTIMATE @ abov = zeros(iter,1) ; @ ABo VARIANCE @ abos = 0 ; @ SIZE OF T-TEST @ @ AHN-SCHMIDT GMM1 @ as1b = zeros(iter,1) ; @ AS1 ESTIMATE @ as1v = zeros(iter,1) ; @ AS1 VARIANCE @ as1s = 0 ; @ SIZE OF T-TEST @ as1s1= 0 ; @ SIZE OF AHN-TEST @ @ AHN-SCHMIDT GMM1 WITH AS1 MOMENT CONDITIONS ONLY @ as1ob = zeros(iter,1) ; @ AS1 ESTIMATE @ as1ov = zeros(iter,1) ; @ AS1 VARIANCE @ as1os = 0 ; @ SIZE OF T-TEST @ as1os1= 0 ; @ SIZE OF AHN-TEST @ /******************** ** BEGING THE LOOP ** ********************/ i = 1; do while i <= iter ; /*********************** ** CREATING VARIABLES ** ***********************/ @ UPDATING SEED FOR EACH ITERATION @ ss1 = (ddd+i) ; @ CREATING INDIVIDUAL EFFECTS @ aaa = sa*rndns(nn,1,ss1) ; @ CREATING DATA: YY INCLUDES Y0,Y1,...,YT @ yy = rr*aaa + s0*rndns(nn,1,ss1) ; j = 1 ; do while j <= tt ; yy = yy~(dd*yy[.,j]+aaa+se*rndns(nn,1,ss1)) ; j = j+1; endo ; @ CREATING DIFFERENCED DATA: DYY INCLUDES (Y0-Y1),...,(YT-1-YT) @ dyy = yy[.,2]-yy[.,1] ; j = 2 ; do while j <= tt ; dyy = dyy~(yy[.,j+1]-yy[.,j]) ; j = j+1; endo ; @ STACKING DATA @ depy = vec(yy[.,2:tt+1]') ; @ Y @ lagy = vec(yy[.,1:tt]') ; @ Y_{-1} @ pdepy = meanc(yy[.,2:tt+1]') .*. ones(tt,1) ; @ P_v Y @ plagy = meanc(yy[.,1:tt]' ) .*. ones(tt,1) ; @ P_v Y_{-1} @ qdepy = depy - pdepy ; @ Q_v Y @ qlagy = lagy - plagy ; @ Q_v Y_{-1} @ ddepy = vec(dyy[.,2:tt]'); dlagy = vec(dyy[.,1:tt-1]'); abwdy = ddepy ; abwly = vec(yy[.,2:tt]') ; abwdly= dlagy ; @ CREATING ARELLANO AND BOND, ARELLANO-BOVER IVS @ @ CREATING L(T,P) MATRIX: NEEDED TO CREATE AB IVS @ proc ll(t,p) ; local i,u ; u = zeros(t,p) ; i = 1 ; do while i <= p ; u[t-p+i-1,i] = - 1 ; u[t-p+i,i] = 1 ; i = i + 1 ; endo ; retp(u) ; endp ; @ CREATING AHLL(T,P) MATRIX: NEEDED TO CREATE DIAGONAL ANDERSON-HSIAO IVS @ proc ahll(t,p) ; local i,u ; u = zeros(t,1) ; u[p+2] = - 1 ; u[p+1] = 1 ; retp(u) ; endp ; @ CREATING AB INSTRUMENTS @ abiv = zeros(tt*nn,1) ; j = 1 ; do while j <= tt-1 ; abiv = abiv~( yy[.,j] .*. ll(tt,tt-j) ) ; j = j + 1 ; endo ; pp = tt*(tt-1)/2 ; abiv = abiv[.,2:pp+1] ; @ CREATING SUB-AB (Diagonal form of Andersen-Hsiao) INSTRUMENTS @ ahiv = zeros(tt*nn,1) ; j = 1 ; do while j <= tt-1 ; ahiv = ahiv~( yy[.,j] .*. ahll(tt,j-1) ) ; j = j + 1 ; endo ; ahpp=tt-1 ; ahiv = ahiv[.,2:ahpp+1] ; @ CREATING ABo INSTRUMENTS @ aboiv = zeros(tt*nn,1) ; j = 1 ; do while j <= tt-1 ; mmmm = zeros(tt,1) ; mmmm[j+1] = 1 ; aboiv = aboiv~( dyy[.,j] .*. mmmm ) ; j = j + 1 ; endo ; aboiv = aboiv[.,2:cols(aboiv)] ; @ CREATING AB-Type ABo (AABo) INSTRUMENTS @ aaboiv = zeros(tt*nn,1) ; j = 1 ; do while j <= tt-1 ; mmmm = zeros(tt,1) ; mmmm[j+1] = 1 ; aaboiv = aaboiv~( dyy[.,1:j] .*. mmmm ) ; j = j + 1 ; endo ; aaboiv = aaboiv[.,2:cols(aaboiv)] ; @ CREATING DIAGONAL LAGGED Y @ diagly = zeros(tt*nn,1) ; j = 1 ; do while j <= tt ; mmmm = zeros(tt,1) ; mmmm[j] = 1 ; diagly = diagly~(yy[.,j] .*. mmmm) ; j = j + 1 ; endo ; g1pp = cols(diagly) ; diagly = diagly[.,2:g1pp] ; @ # OF AB INSTRUMENTS @ nabiv = pp ; @ # OF AH INSTRUMENTS @ nahiv = cols(ahiv) ; /********************** ** ARELLANO-BOND GMM ** **********************/ @ TSLS ESTIMATOR @ tscov = invpd( (lagy'abiv)*invpd(abiv'abiv)*(abiv'lagy) ) ; tsb = tscov*(lagy'abiv)*invpd(abiv'abiv)*(abiv'depy) ; qres = qdepy - qlagy*tsb ; tssige2 = (qres'qres/(nn*tt-tt-1)) ; tscov = tssige2*tscov ; res = depy - lagy*tsb ; opivab = zeros(1,cols(abiv)) ; j = 0 ; do while j <= nn-1 ; zz = res[tt*j+1:tt*j+tt,.]'abiv[tt*j+1:tt*j+tt,.] ; opivab = opivab|zz ; j = j + 1 ; endo ; opivab = opivab[2:nn+1,.] ; optwab = invpd(opivab'opivab) ; VARAB1 = INVPD( (LAGY'ABIV)*(OPTWAB)*(ABIV'LAGY) ) ; ABB1 = VARAB1*(LAGY'ABIV)*(OPTWAB)*(ABIV'DEPY) ; abv[I,1] = VARAB1 ; abb[I,1] = ABB1 ; wald = (ABB1-dd)^2/VARAB1 ; pval = cdfchic(wald,1) ; @ Computing rejection rate @ if pval < 0.05 ; rej05 = 1 ; else ; rej05 = 0 ; endif ; abss = rej05 + abss ; /*********************** ** ARELLANO-BOVER GMM ** ***********************/ iv = abiv~aboiv ; @ TSLS ESTIMATOR @ tscov = invpd( (lagy'iv)*invpd(iv'iv)*(iv'lagy) ) ; tsb = tscov*(lagy'iv)*invpd(iv'iv)*(iv'depy) ; qres = qdepy - qlagy*tsb ; tssige2 = (qres'qres/(nn*tt-tt-1)) ; tscov = tssige2*tscov ; res = depy - lagy*tsb ; opivab = zeros(1,cols(iv)) ; j = 0 ; do while j <= nn-1 ; zz = res[tt*j+1:tt*j+tt,.]'iv[tt*j+1:tt*j+tt,.] ; opivab = opivab|zz ; j = j + 1 ; endo ; opivab = opivab[2:nn+1,.] ; optwab = invpd(opivab'opivab) ; abov1 = INVPD( (LAGY'IV)*(OPTWAB)*(IV'LAGY) ) ; abob1 = abov1*(LAGY'IV)*(OPTWAB)*(IV'DEPY) ; abov[I,1] = abov1 ; abob[I,1] = abob1 ; wald = (abob1-dd)^2/abov1 ; pval = cdfchic(wald,1) ; @ Computing rejection rate @ if pval < 0.05 ; rej05 = 1 ; else ; rej05 = 0 ; endif ; abos = rej05 + abos ; /********************* ** AHN-SCHMIDT GMM1 ** *********************/ /* ** First step GMM */ @ define some global variables for gmm1 @ WELL = LL(TT,TT-1) ; well = well[1:tt,1:tt-2] ; m11 = abiv'depy ; m12 = abiv'lagy ; m21 = (yy[.,tt+1].*.well)'depy ; m22 = (yy[.,tt].*.well)'depy + (yy[.,tt+1].*.well)'lagy ; m23 = (yy[.,tt].*.well)'lagy ; library optmum; #include optmum.ext; optset ; @ Finding initial value @ @ step 1 for gmm1 @ m11 = abiv'depy ; m12 = abiv'lagy ; m21 = (yy[.,tt+1].*.well)'depy ; m22 = (yy[.,tt].*.well)'depy + (yy[.,tt+1].*.well)'lagy ; m23 = (yy[.,tt].*.well)'lagy ; wei1 = invpd((abiv'abiv)/nn)~zeros(cols(abiv),tt-2) ; wei2 = zeros(tt-2,cols(abiv))~eye(tt-2) ; wei = wei1|wei2 ; /* GRADIENT */ proc grad1(b) ; local gra, m ; m = (m11-m12*b)|(m21-m22*b+m23*b^2) ; gra = m12|(m22-2*m23*b) ; retp( - 2*m'wei*gra ) ; endp ; /* GMM */ proc f1(b) ; local m ; m = (m11-m12*b)|(m21-m22*b+m23*b^2) ; retp( m'wei*m ) ; ENDP ; _opgdprc = &grad1 ; @ starting values @ b0 = -1 ; __title = "GMM11"; _opgtol = 1e-4 ; _opstmth = "bfgs brent"; __output = 0 ; {b,func,grad,retcode} = optmum(&f1,b0) ; g = m12|(m22-2*m23*b) ; v = invpd(g'g) ; s = (b-dd)^2/v ; bb1 = func~b~v~s ; @ step 2 @ /* GMM */ proc f2(b) ; local m ; m = (m11-m12*b)|(m21-m22*b+m23*b^2) ; retp( m'wei*m ) ; ENDP ; _opgdprc = &grad1 ; @ starting values @ b0 = 1 ; __title = "GMM12"; _opgtol = 1e-4 ; _opstmth = "bfgs brent"; __output = 0 ; {b,func,grad,retcode} = optmum(&f2,b0) ; g = m12|(m22-2*m23*b) ; v = invpd(g'g) ; s = (b-dd)^2/v ; bb2 = func~b~v~s ; bb = bb1|bb2 ; bb = sorthc(bb,1) ; tsbb = bb[1,2] ; @ define some global variables for gmm1 @ REST = YY[.,TT+1] - YY[.,TT]*tsbb ; IVB1 = REST .*. WELL ; RES = DEPY - LAGY*tsbb ; AS1IV = ABIV~IVB1 ; opivas1 = zeros(1,cols(as1iv)) ; j = 0 ; do while j <= nn-1 ; zz = res[tt*j+1:tt*j+tt,.]'as1iv[tt*j+1:tt*j+tt,.] ; opivas1 = opivas1|zz ; j = j + 1 ; endo ; opivas1 = opivas1[2:nn,.] ; opwas1 = invpd(opivas1'opivas1) ; @ step 1 for gmm1 @ /* GRADIENT */ proc grad2(b) ; local gra, m ; m = (m11-m12*b)|(m21-m22*b+m23*b^2) ; gra = m12|(m22-2*m23*b) ; retp( - 2*m'opwas1*gra ) ; endp ; /* GMM */ proc f3(b) ; local m ; m = (m11-m12*b)|(m21-m22*b+m23*b^2) ; retp( m'opwas1*m ) ; ENDP ; _opgdprc = &grad2 ; @ starting values @ b0 = -1 ; __title = "GMM11"; _opgtol = 1e-4 ; _opstmth = "bfgs brent"; __output = 0 ; {b,func,grad,retcode} = optmum(&f3,b0) ; g = m12|(m22-2*m23*b) ; v = invpd(g'opwas1*g) ; s = (b-dd)^2/v ; bb1 = func~b~v~s ; @ step 2 @ /* GMM */ proc f4(b) ; local m ; m = (m11-m12*b)|(m21-m22*b+m23*b^2) ; retp( m'opwas1*m ) ; ENDP ; _opgdprc = &grad2 ; @ starting values @ b0 = 1 ; __title = "GMM12"; _opgtol = 1e-4 ; _opstmth = "bfgs brent"; __output = 0 ; {b,func,grad,retcode} = optmum(&f4,b0) ; g = m12|(m22-2*m23*b) ; v = invpd(g'opwas1*g) ; s = (b-dd)^2/v ; bb2 = func~b~v~s ; bb = bb1|bb2 ; bb = sorthc(bb,1) ; as1b[i,1] = bb[1,2] ; as1v[i,1] = bb[1,3] ; wald = bb[1,4] ; pval = cdfchic(wald,1) ; @ Computing rejection rate @ if pval < 0.05 ; rej05 = 1 ; else ; rej05 = 0 ; endif ; as1s = as1s + rej05 ; wald = f3(dd) - bb[1,1] ; pval = cdfchic(wald,1) ; @ Computing rejection rate @ if pval < 0.05 ; rej05 = 1 ; else ; rej05 = 0 ; endif ; as1s1 = as1s1 + rej05 ; i = i + 1 ; endo ; i = i - 1 ; "" ; "AB GMM RESULTS:" ; " MEAN OF ESTIMATES :" meanc(abb) ; " MSE :" meanc((abb-dd)^2) ; " MEAN OF ASYM. VAR :" meanc(abv) ; " SIZE OF T-TEST :" abss/iter ; ""; "ABo GMM RESULTS:" ; " MEAN OF ESTIMATES :" meanc(abob) ; " MSE :" meanc((abob-dd)^2) ; " MEAN OF ASYM. VAR :" meanc(abov) ; " SIZE OF T-TEST :" abos/iter ; ""; ""; "AS1 GMM RESULTS:" ; " MEAN OF ESTIMATES :" meanc(as1b) ; " MSE :" meanc((as1b-dd)^2) ; " MEAN OF ASYM. VAR :" meanc(as1v) ; " SIZE OF T-TEST :" as1s/iter ; " SIZE OF AHN TEST :" as1s1/iter ; j4 = j4 + 1 ; endo ; j3 = j3 + 1 ; endo ; j2 = j2 + 1 ; endo ; j1 = j1 + 1 ; endo ; OUTPUT OFF