set echo on set messages on set seed 1234 # Section 4: The DPB function package install DPB.zip include DPB.gfn # Generate artifical data scalar N = 1024 scalar T = 6 scalar NT = N*T nulldata NT --preserve setobs T 1:1 --stacked-time-series genr time series id = $unit series x1 = normal() series x2 = x1 + normal() series x3 = randgen(t,2) list X = x1 x2 x3 gamma = 0.6 matrix beta = ones(nelem(X),1)*0.5 s2a = 1 Ti = seq(1,N)' ~ ones(N,1)*T A = mnormal(N,1)*s2a ndx = lincomb(X,beta) ini = 1 matrix Y = zeros(sumc(Ti[,2]),1) loop i=1..N --quiet T = Ti[i,2] if T>0 fin = ini + T - 1 ndy = ndx[ini] + A[i] p = exp(ndy)/(1+exp(ndy)) Y[ini] = (randgen1(u,0,1)
0 # hyslop scalar rho = par[3] s = seq(1, T-1)' matrix r = rho .^ abs(s .- s') r[diag] = 0 S[2:,2:] += r # printf "At [1], r = \n%6.2f", r tau = (AR1==1) ? 1 : par[4] r = rho .^ s # printf "At [2], r = \n%6.2f", r S[2:,1] += tau*r S[1,2:] += tau*r' endif if !err catch C = cholesky(S) if $error err = 2 endif endif return C end function function series gendata(list X, list Z, matrix param, \ matrix Ti, scalar AR1) N = rows(Ti) nk = nelem(X) nz = nelem(Z) scalar gamma = param[1] matrix beta = param[2:nk+1] matrix pivec = param[nk+2:nk+nz+1] scalar theta = param[nk+nz+2] scalar sig = exp(param[nk+nz+3]) matrix ndx0 = lincomb(Z,pivec) matrix ndx1 = lincomb(X,beta) TT = maxc(Ti[,2]) + 1 matrix vpar = theta | sig if AR1==1 vpar |= tanh(param[nk+nz+4]) elif AR1==2 vpar |= tanh(param[nk+nz+4]) | param[nk+nz+5] endif scalar err = 0 matrix C = C_theta(vpar, TT, AR1, &err) matrix E = mnormal(N,TT)*C' ini = 1 matrix Y = zeros(sumc(Ti[,2]),1) loop i=1..N --quiet T = Ti[i,2] if T>0 fin = ini + T - 1 Y[ini] = (ndx0[ini] + E[i,1])>0 loop t=2..T --quiet Y[ini+t-1] = (ndx1[ini+t-1] + E[i,t] + gamma*Y[ini+t-2])>0 endloop ini = fin + 1 endif endloop series ret = Y return ret end function #--------- GENERATE ARTIFICIAL DATA ---------- scalar N = 4096 scalar T = 6 scalar NT = N*T nulldata NT --preserve setobs T 1:1 --stacked-time-series genr time series id = $unit series x = normal() series z = uniform() list X = const x list Z = const x z gamma = 0.6 matrix beta = ones(nelem(X),1)*0.5 matrix pivec = ones(nelem(Z),1) theta = 1.2 s2a = 1 Ti = seq(1,N)' ~ ones(N,1)*T # DP: Heckman (1981) rho = 0 tau = 0 matrix param = gamma | beta | pivec | theta | ln(sqrt(s2a)) | \ atanh(rho) | atanh(tau) AR1 = 0 series y = gendata(X, Z, param, Ti, AR1) store DP_artdata.gdtb # ADP: Hyslop (1999) rho = 0.3 tau = 1 matrix param = gamma | beta | pivec | theta | ln(sqrt(s2a)) | \ atanh(rho) | atanh(tau) AR1 = 1 series y = gendata(X, Z, param, Ti, AR1) store ADP_artdata.gdtb # GADP: Keane and Sauer (2009) rho = 0.3 tau = 0.6 matrix param = gamma | beta | pivec | theta | ln(sqrt(s2a)) | \ atanh(rho) | atanh(tau) AR1 = 2 series y = gendata(X, Z, param, Ti, AR1) store GADP_artdata.gdtb # Section 4.1: Random-effects dynamic probit # ------------ DP model: Heckman (1981b) ----------------- open DP_artdata.gdtb list X = const x list Z = const x z bundle b b = DPB_setup("DP", y, X, Z) DPB_estimate(&b) print b # bwrite(b, "mod") # b = bread("mod") scalar npar = b.npar inipar = muniform(npar, 1) DPB_estimate(&b, &inipar) DPB_printout(&b) DPB_printape(&b) err = DPB_setoption(&b, "nrep", 32) DPB_estimate(&b) err = DPB_setoption(&b, "method", 1) err = DPB_setoption(&b, "nrep", 100) DPB_estimate(&b) DPB_printout(&b) # ------------------------------- TABLE 3 ---------------------------------- b = DPB_setup("DP", y, X, Z) err1 = DPB_setoption(&b, "verbose", 2) err2 = DPB_setoption(&b, "vcv", 1) loop foreach i 16 24 32 set stopwatch err_$i = DPB_setoption(&b, "nrep", $i) DPB_estimate(&b) DPB_printout(&b) t_$i = $stopwatch print t_$i endloop err3 = DPB_setoption(&b, "method", 1) loop foreach i 128 192 256 set stopwatch err_$i = DPB_setoption(&b, "nrep", $i) DPB_estimate(&b) DPB_printout(&b) t_$i = $stopwatch print t_$i endloop #------- ADP MODEL: HYSLOP (1999) --------- open ADP_artdata.gdtb list X = const x list Z = const x z bundle b b = DPB_setup("ADP", y, X, Z) DPB_estimate(&b) DPB_printout(&b) #------- GADP MODEL: KEANE and SAUER (2009) --------- open GADP_artdata.gdtb list X = const x list Z = const x z bundle b b = DPB_setup("GADP", y, X, Z) DPB_estimate(&b) DPB_printout(&b) # Section 4.2: Quadratic exponential model open QE_artdata.gdtb list X = x1 x2 x3 bundle mod mod = DPB_setup("QE", y, X) err = DPB_setoption(&mod, "vcv", 2) DPB_estimate(&mod) DPB_printout(&mod) genr time list TIME = dummify(time) TIME -= Dtime_2 Dtime_6 X += TIME mod = DPB_setup("QE", y, X) err = DPB_setoption(&mod, "vcv", 2) DPB_estimate(&mod) DPB_printout(&mod) # Section 5.1: the union dataset, Example 1 open union_full.gdtb summary union age grade south not_smsa year --simple series valid = 0 smpl year>=78 --restrict smpl year!=83 --restrict matrix m = aggregate(const, idcode) series nwav = replace(idcode, m[, 1], m[, 2]) smpl nwav==6 --restrict series valid = 1 store union.gdtb open union.gdtb --quiet smpl full smpl valid --dummy setobs idcode year --panel-vars list X = const age grade south list Z = X not_smsa # --------------------------- TABLE 5 -------------------------------- #------------ POOLED PROBIT------------------- probit union union(-1) X --p-values #------------ STATIC RE PROBIT------------------- probit union union(-1) X --p-values --random-effects # --------------------------- TABLE 6 -------------------------------- #------------ RANDOM-EFFECTS MODELS------------------- # DP MODEL: Heckman (1981) # Replication of Stewart (2006) Table 1, column 3 # 24 quadpoints # Standard errors based on Hessian b10 = DPB_setup("DP", union, X, Z) err10 = DPB_setoption(&b10, "vcv", 2) DPB_estimate(&b10) DPB_printout(&b10) # Replication of Stewart, DP MODEL with GHK nrep = 500 b11 = DPB_setup("DP", union, X, Z) err11 = DPB_setoption(&b11, "vcv", 2) err11 = DPB_setoption(&b11, "method", 1) err11 = DPB_setoption(&b11, "nrep", nrep) DPB_estimate(&b11) DPB_printout(&b11) # ADP MODEL: Hyslop (1999) # Replication of Stewart (2006) Table 1, column 4 # 500 replications # Standard errors based on Hessian b2 = DPB_setup("ADP", union, X, Z) err2 = DPB_setoption(&b2, "vcv", 2) err2 = DPB_setoption(&b2, "nrep", nrep) DPB_estimate(&b2) DPB_printout(&b2) # -------------------------------- TABLE 7 --------------------------- # GADP MODEL: Keane & Sauer (2009) # 500 replications # Standard errors based on Hessian b3 = DPB_setup("GADP", union, X, Z) err3 = DPB_setoption(&b3, "vcv", 2) err3 = DPB_setoption(&b3, "nrep", nrep) DPB_estimate(&b3) DPB_printout(&b3) # QUADRATIC EXPONENTIAL MODEL b4 = DPB_setup("QE", union, X) err4 = DPB_setoption(&b4, "vcv", 2) DPB_estimate(&b4) DPB_printout(&b4) # ---------------------------- TABLE 8 --------------------------- open union_full.gdtb setobs idcode year --panel-vars list X = const age grade south list Z = X not_smsa # DP MODEL b10 = DPB_setup("DP", union, X, Z) err10 = DPB_setoption(&b10, "vcv", 2) DPB_estimate(&b10) DPB_printout(&b10) nrep = 500 # ADP MODEL b2 = DPB_setup("ADP", union, X, Z) err2 = DPB_setoption(&b2, "vcv", 2) err2 = DPB_setoption(&b2, "nrep", nrep) DPB_estimate(&b2) DPB_printout(&b2) # GADP MODEL b3 = DPB_setup("GADP", union, X, Z) err3 = DPB_setoption(&b3, "vcv", 2) err3 = DPB_setoption(&b3, "nrep", nrep) DPB_estimate(&b3) DPB_printout(&b3) # QUADRATIC EXPONENTIAL MODEL b4 = DPB_setup("QE", union, X) err4 = DPB_setoption(&b4, "vcv", 2) DPB_estimate(&b4) DPB_printout(&b4) # Section 5.2: the union dataset, Example 2 open union_wooldridge.gdt #-------------------------- TABLE 9 ---------------------------- list TIME = dummify(year) TIME -= Dyear_2 list MARR = marr8* list X = const union_1 union80 married MARR TIME # CDP MODEL probit union X --random-effects --quadpoints=12 probit union X educ black --random-effects --quadpoints=12 # DP MODEL list X1 = const married MARR TIME list Z1 = const married b1 = DPB_setup("DP", union, X1, Z1) DPB_setoption(&b1, "nrep", 12) DPB_estimate(&b1) DPB_printout(&b1) list X2 = const married MARR TIME educ black list Z2 = const married educ black b2 = DPB_setup("DP", union, X2, Z2) DPB_setoption(&b2, "nrep", 12) DPB_estimate(&b2) DPB_printout(&b2) #-------------------------- TABLE 10 ---------------------------- m_marr = pmean(married) list X = const union_1 union80 married m_marr TIME # CDP MODEL probit union X --random-effects --quadpoints=12 # DP MODEL list X = const married m_marr TIME list Z = const married b = DPB_setup("DP", union, X, Z) DPB_setoption(&b, "nrep", 12) DPB_estimate(&b) DPB_printout(&b) # ----------------------------- TABLE 11 --------------------------- # ---- computation of predicted probabilities -------------------- list TIME = dummify(year) TIME -= Dyear_2 list MARR = marr8* list X = const union_1 union80 married MARR TIME probit union X --random-effects --quadpoints=12 -q m = nelem(MARR) k = nelem(X) bw = $coeff[1:k] s2a = exp($coeff[k+1]) # p11_m = P(u87=1 | u86=1, married=1) series ndx = bw[1] + bw[2] + bw[3]*union80 + bw[4] + lincomb(MARR, bw[5:m+4]) + bw[m+10] series ndx /= sqrt(1+s2a) wp11_m = cnorm(ndx) # p10_m = P(u87=1 | u86=0, married=1) series ndx -= bw[2]/sqrt(1+s2a) wp10_m = cnorm(ndx) # p11 = P(u87=1 | u86=1, married=0) series ndx = bw[1] + bw[2] + bw[3]*union80 + lincomb(MARR, bw[5:m+4]) + bw[m+10] series ndx /= sqrt(1+s2a) wp11 = cnorm(ndx) # p10 = P(u87=1 | u86=0, married=0) series ndx -= bw[2]/sqrt(1+s2a) wp10 = cnorm(ndx) summary wp11_m wp10_m wp11 wp10 --simple #---- DP MODEL: Heckman (1981b) ---------- list X = const married MARR TIME list Z = const married b = DPB_setup("DP", union, X, Z) DPB_setoption(&b, "nrep", 12) DPB_setoption(&b, "verbose", 0) DPB_estimate(&b) DPB_printout(&b) k = b.nk z = b.nz bz = b.coeff[k+2:k+z+1] bh = b.coeff[1:k+1] sig_a = b.coeff[k+z+3] s2a = sig_a^2 theta = b.coeff[k+z+2] # hp11_m = P(u87=1 | u86=1, married87=1) series ndx = bh[1] + bh[2] + bh[3] + lincomb(MARR, bh[4:m+3]) + bh[m+9] series ndx /= sqrt(1+s2a) series hp11_m = cnorm(ndx) # hp10_m = P(u87=1 | u86=0, married87=1) series ndx -= bh[1]/sqrt(1+s2a) series hp10_m = cnorm(ndx) # hp11 = P(u87=1 | u86=1, married87=0) series ndx = bh[1] + bh[2] + lincomb(MARR, bh[4:m+3]) + bh[m+9] series ndx /= sqrt(1+s2a) series hp11 = cnorm(ndx) # hp10 = P(u87=1 | u86=0, married87=0) series ndx -= bh[1]/sqrt(1+s2a) series hp10 = cnorm(ndx) # ----------------------------- TABLE 11 --------------------------- summary hp11_m hp10_m hp11 hp10 --simple