* * ACOUP1.RAT * dep var gchangu * * See the appendix of Londregan & Poole (1988) for the details. * * LARGE SAMPLE 1950-1982 * * bma global 500 local 50 * * CHANGE 1: On next 2 lines, change 33 to 23 for the 1960-82 sample. * calendar(org=indiv) 33 allocate 40 119:33 env errors=1 warning=1 seed 16339 * declare index rhs1(1) rhs2(1) rhsw(1) declare index bser1(1) bserw2(1) strser(1) * * read in data from rats format file * * CHANGE 2: DATA60.DAT for the 1960-82 sample. * open data data50.dat data(format=rats) / grgdp gchangu $ educ wgrowth1 grgdp1 $ exadj1 $ dstham dafrica gchangu1 * * * growth = educ * gchangu = exadj{1} gchangu{1} * * both = constant dstham dafrica wgrowth1 grgdp1 * * * * CHANGE 3: number of draws for the bootstrap * ieval ndraws = 1024 * * CHANGE 4: enter left-hand-side variables for the two equations * set lhs1 = grgdp(t) set lhs2 = gchangu(t) * * CHANGE 5: enter instrument(s) for growth * enter(varying,entries=k1) rhs1 # educ * * CHANGE 6: enter variables which affect both growth & instability * (do NOT enter CONSTANT here) * enter(varying,entries=kw) rhsw # dstham dafrica wgrowth1 grgdp1 * * CHANGE 7: enter variables which affect instability * enter(varying,entries=k2) rhs2 # exadj1 gchangu1 ieval kw = kw + 1 ieval k = k1+kw+k2 ieval kk = 2*k display 'k=' k display 'k1=' k1 display 'k2=' k2 display 'kw=' kw display 'kk=' kk declare lvector sn1(k1) sn2(k2) snw(kw) * * CHANGE 8: enter variable labels (names/titles) for growth instruments * input sn1 educ * * CHANGE 9: enter variable labels (names/titles) for instability instruments * input sn2 exadj1 gchangu1 * * CHANGE 10: enter variable labels (names/titles) for common variables * (NOTE: you MUST have constant here) * input snw constant dsouth dafrica wgrowth1 grgdp1 ************************************************ * Do not make any changes below here * ************************************************ * * get rid of missing observations * make alldata / numobs numvars #lhs1 lhs2 rhs1 rhsw rhs2 display 'numobs = ' numobs display 'numvars= ' numvars ieval colnow = 1 dofor s = lhs1 lhs2 rhs1 rhsw rhs2 set s 1 numobs = alldata(t,colnow) ieval colnow = colnow + 1; end dofor smpl 1 numobs * * estimate reduced form coefficients * linreg(define=groweqn) lhs1 / reseqn1 psi1 # rhs1 constant rhsw rhs2 declare vector pihat(kk) declare vector g1(kk) declare vector g2(kk) declare symm i2k(kk,kk) declare symm ik(k,k) eval tau = sqrt(rss/nobs) prbit lhs2 / lprobs psi2 # rhs1 constant rhsw rhs2 reseqn1 eval eta = psi2(k+1) eval varadj = 1.0 + eta*eta*tau*tau eval rhohat = (eta*tau)/sqrt(varadj) eval sigmahat = tau display 'eta = ' eta display 'sigmahat = ' tau display 'rhohat= ' rhohat matrix ik = iden(1) matrix i2k = iden(1) do i=1,k eval pihat(i) = psi1(i) eval pihat(k+i) = psi2(i)/(sqrt(varadj)) eval g1(i) = pihat(k+i) eval g1(k+i) = 0.0 eval g2(i) = 0.0 eval g2(k+i) = pihat(i) end do i set vpihat 1 kk = pihat(t) set sg1 1 kk = g1(t) set sg2 1 kk = g2(t) * * get the correct reduced form significance tests * prj(cdf=pfitted) set prus = lhs2(t) - pfitted(t) make pru / pruobs # prus display 'pruobs=' pruobs display 'nobs =' nobs eval pseesq = 0.0 do i=1,pruobs eval pseesq = pseesq + (pru(i,1)*pru(i,1))/(pruobs-k) end do i display 'psee = ' sqrt(pseesq) do i=1,k eval wi = (pseesq * xx(i,i))/varadj eval ti = psi2(i)/sqrt(wi) eval coi = pihat(k+i) display 'coefficient for variable ' i ' = ' coi cdf ttest ti (pruobs-k) end do i * * do i=1,k1 set k+i 1 kk = i2k(t,i) labels k+i # sn1(i) end do i do i=1,kw set k+k1+i 1 kk = i2k(t,k1+i) set k+k1+kw+i 1 kk = i2k(t,k+k1+i) labels k+k1+i k+k1+kw+i # snw(i) snw(i) end do i do i=1,k2 set k+k1+kw+kw+i 1 kk = i2k(t,k+k1+kw+i) labels k+k1+kw+kw+i # sn2(i) end do i ieval strser(1) = k+1 do i=k+2,(k+k1+kw+kw+k2) enter(varying) strser # strser i end do i * * * bootstrap VCOV of vpihat <-- vec(pihat) * * declare vector bpihat(kk) declare vector betamean(kk) declare symm betacmom(kk,kk) matrix betamean=const(0.0) matrix betacmom=const(0.0) eval bsigmean = 0.0 eval bsigcmom = 0.0 eval brhomean = 0.0 eval brhocmom = 0.0 ieval bser1(1) = 1 do ssnow=2,k1 enter(varying) bser1 # bser1 ssnow end do ssnow ieval bserw2(1) = k1+1 do ssnow=k1+2,k-1 enter(varying) bserw2 # bserw2 ssnow end do ssnow do jj=1,ndraws boot entry set blhs1 = lhs1(entry(t)) set blhs2 = lhs2(entry(t)) ieval sernow = 1 dofor ssnow = rhs1 rhsw rhs2 set sernow = ssnow(entry(t)) ieval sernow = sernow + 1 end do ssnow linreg(noprint) blhs1 / breseqn1 bpsi1 # bser1 constant bserw2 eval btau = sqrt(rss/nobs) prbit(noprint) blhs2 / blprobs bpsi2 # bser1 constant bserw2 breseqn1 eval bteta = bpsi2(k+1) eval brhohat = (bteta*btau)/sqrt(1+bteta*bteta*btau*btau) eval brhomean = brhomean + brhohat eval brhocmom = brhocmom + brhohat*brhohat eval bsighat = btau eval bsigmean = bsigmean + bsighat eval bsigcmom = bsigcmom + bsighat*bsighat do i=1,k eval bpihat(i) = bpsi1(i) eval bpihat(k+i) = bpsi2(i)/(sqrt(1 + bteta*bteta*btau*btau)) end do i ewise betamean(i) = betamean(i) + bpihat(i) ewise betacmom(i,j) = betacmom(i,j) + bpihat(i)*bpihat(j) end do jj matrix betamean=betamean*scale(1.0/ndraws) ewise betacmom(i,j)=betacmom(i,j)/ndraws-betamean(i)*betamean(j) eval brhomean = brhomean /ndraws eval brhocmom = brhocmom/ndraws - brhomean*brhomean eval brhostd = sqrt(brhocmom) eval bsigmean = bsigmean /ndraws eval bsigcmom = bsigcmom/ndraws - bsigmean*bsigmean eval bsigstd = sqrt(bsigcmom) display 'var(rho) = ' brhocmom ' std dev = ' brhostd display 'var(sigma) = ' bsigcmom ' std dev = ' bsigstd * * * * Recover structural parameters * smpl 1 kk linreg(define=recoveqn, noprint) vpihat # sg1 sg2 strser eval gam1star = beta(1) eval gam2star = beta(2) declare rect gamstar(2,2) matrix gamstar = iden(1) eval gamstar(1,2) = -gam1star eval gamstar(2,1) = -gam2star declare symm sigmapi(kk,kk) declare rect tmpgam(kk,kk) matrix tmpgam = kroneker(ik,gamstar) matrix sigmapi = tmpgam*betacmom*tr(tmpgam) ****** GLS declare rect tmp1(kk,kk) declare rect xy1(kk,1) declare vector sigu(kk) declare symm isigpi(kk,kk) make(equation=recoveqn) xt1 matrix isigpi = inv(sigmapi) matrix tmp1 = tr(xt1)*isigpi matrix xx = inv(tmp1*xt1) matrix xy1 = tmp1*pihat matrix beta = xx * xy1 matrix sigu = pihat - xt1 * beta set ssigus 1 kk = sigu(t) linreg(equation=recoveqn,residuals=ssigus,create,noscale) declare vector betagls(2+k1+k2+kw+kw) matrix betagls = beta *** chi-square test of over identifying restrictions prj prvpihat declare vector prpihat(kk) do i=1,kk eval prpihat(i) = prvpihat(i) end do i declare real overiden matrix overiden = tr(prpihat - pihat) * inv(betacmom) * (prpihat-pihat) ieval orestrct = k1 + k2 - 2 cdf chisquare overiden orestrct fetch qcdstat = cdstat fetch qsignif = signif display 'degrees of freedom = ' orestrct display 'stat into program = ' overiden display 'test statistic = ' qcdstat display 'significance level = ' qsignif