name bivord 'Bivariate Ordered Probit Example in ML' ; ? ? TSP 4.2B (November 1994) ? Author: Jan Hanousek, CERGE-Prague and U of Pennsylvania ? -------------------------------- tsp program ----------------------------- ? ? ordered bivariate probit (6x6 states) ? The following specification is going to be estimated ? Y1 = f(X1, X2, X3, X4, X13, X14) ? and ? Y2 = g(Y1, X1, X2, X3, X4, X15, X16, X17, X18, X19) ? date: November 6, 1994, Jan Hanousek freq n; smpl 1,453; in r453, w453; ? input of TSP databank ? ? Q1 = X1B1 + U1 (where Q1 is unobserved) ? Y1 = 0 for Q1 <= 0 ( U1 <= - X1B1 ) ? 1 for 0 < Q1 <= A1 ( - X1B1 < U1 <= A1 - X1B1 ) ? 2 for A1 < Q1 <= A2 ( A1 - X1B1 < U1 <= A2 - X1B1 ) ? 3 for A2 < Q1 <= A3 ( A2 - X1B1 < U1 <= A3 - X1B1 ) ? 4 for A3 < Q1 <= A4 ( A3 - X1B1 < U1 <= A4 - X1B1 ) ? 5 for A4 < Q1 ( A4 - X1B1 < U1 ) ? ? Q2 = X2B2 + U2 (where Q2 is unobserved) ? Y2 = 0 for Q2 <= 0 ( U2 <= - X2B2 ) ? 1 for 0 < Q2 <= C1 ( - X2B2 < U2 <= C1 - X2B2 ) ? 2 for C1 < Q2 <= C2 ( C1 - X2B2 < U2 <= C2 - X2B2 ) ? 3 for C2 < Q2 <= C3 ( C2 - X2B2 < U2 <= C3 - X2B2 ) ? 4 for C3 < Q2 <= C4 ( C3 - X2B2 < U2 <= C4 - X2B2 ) ? 5 for C4 < Q2 ( C4 - X2B2 < U2 ) ? ? As standard, to insure the ordering, we one can use method of squaring, i.e. ? to parametrize: ? A1 = D1**2 ? A2 = A1 + D2**2 ? A3 = A2 + D3**2 ? A4 = A3 + D4**2 ? and ? C1 = G1**2 ? C2 = C1 + G2**2 ? C3 = C2 + G3**2 ? C4 = C3 + G4**2 ? ? Note, it is not used here. ? ? Bivariate cumulative normal distribution: ? Abramovitz and Stegun, Handbook of Mathematical Functions, formula 26.3.29 ? (Taylor expansions for small rho, average abs. error rho=.9 is 0.00887, ? and is much smaller for smaller rho. ? ? FRML biv pr = cnorm(X1B1)*cnorm(X2B2) + norm(X1B1)*norm(X2B2)*( ? rho + rho**2*X1B1*X2B2/2 + rho**3*(X1B1**2 - 1)*(X2B2**2 - 1)/6 ); ? ? ? For bivariate normal distribution we keep notation: ? pr00 = Pr( Y1=0, Y2=0) = Pr( e1<= -X1B1, e2<= -X2B2) ? pr01 = Pr( Y1=0, Y2=1), etc. ? ? Be careful about sign of X1B1 and X2B2 (see formulas for them, at the end) ? ? pr00 frml eqp00 p00 = cnorm(X1B1)*cnorm(X2B2) + norm(X1B1)*norm(X2B2)*( rho + rho**2*X1B1*X2B2/2 + rho**3*(X1B1**2 - 1)*(X2B2**2 - 1)/6 ); ? For likelihood equation we need to create equations ? pij = Pr(y1=i,y2=j) = ? = Pr{ a[i-1]-X1B1 < Q1<= a[i]-X1B1, c[j-1]-X2B2 < Q2 <= c[j]-X2B2 } ? where ? i = 0,1,...,5 and j = 0, 1,...,5 (states) ? a[-1] = c[-1] = - infinity ? a[0] = c[0] = 0 ? a[1],...,a[4] and c[1],..c[4] are current treasholds usually denited by mu(.) ? and a[5] = c[5] = + infinity ? as standard. ? ? To write likelihood equation correctly, we first specify ? qij = Pr{ Q1<= a[i]-X1B1, Q2 <= c[j]-X2B2 } ? where i=0,...4 and j=0,..,4 ? by using above mentioned approximation (see equation biv) ? ? ? To express pij using qij is straightforward: ? p00 = q00 ? p0j = q0j - q0(j-1) for j>0 ? ? symmetrically pi0 = qi0 - q(i-1)0 for i>0 ? pij = qij - q(i-1)j - qi(j-1) + q(i-1)(j-1) for i,j >0 and i<>5, j<>5 ? p5j = Pr(y2=j) - p0j - pj1 - p2j - p3j - p4j ? analogously ? pi5 = Pr(y1=i) - pi0 - pi1 - pi2 - pi3 - pi4 ? and for p55 = Pr{ a[4]-X1B1 < Q1, c[4]-X2B2 < Q2 } = ? = Pr{ a[4]-X1B1 < Q1 } - Pr{ Q1 <= a[4]-X1B1, c[4]-X2B2 < Q2 } = ? = 1 - Pr{ Q1 <= a[4]-X1B1 } - Pr{ Q1 <= a[4]-X1B1 } + ? Pr{ Q1 <= a[4]-X1B1, Q2 <= c[4]-X2B2 < Q2 } = ? = 1 - Pr{ Q1 <= a[4]-X1B1 } - Pr{ Q1 <= a[4]-X1B1 } + q44 ? ? q00 frml eq00 q00 = cnorm(X1B1)*cnorm(X2B2) + norm(X1B1)*norm(X2B2)*( rho + rho**2*X1B1*X2B2/2 + rho**3*(X1B1**2 - 1)*(X2B2**2 - 1)/6 ); ? === 0 === DOT 1-4; FRML eq0. q0. = cnorm(X1B1)*cnorm(X2B2+C.) + norm(X1B1)*norm(X2B2+C.)*( rho + rho**2*X1B1*(X2B2+C.)/2 + rho**3*(X1B1**2 - 1)*( (X2B2+C.)**2 - 1)/6 ); ? FRML eq.0 q.0 = cnorm(X1B1+A.)*cnorm(X2B2) + norm(X1B1+A.)*norm(X2B2)*( rho + rho**2*(X1B1+A.)*X2B2/2 + rho**3*((X1B1+A.)**2 - 1)*( X2B2**2 - 1)/6 ); ENDDOT; ? === 1 === DOT 1-4; FRML eq1. q1. = cnorm(X1B1+A1)*cnorm(X2B2+C.)+norm(X1B1+A1)*norm(X2B2+C.)*( rho + rho**2*(X1B1+A1)*(X2B2+C.)/2 + rho**3*((X1B1+A1)**2 - 1)*((X2B2+C.)**2 - 1)/6 ); ENDDOT; ? === 2 ==== DOT 1-4; FRML eq2. q2. = cnorm(X1B1+A2)*cnorm(X2B2+C.)+norm(X1B1+A2)*norm(X2B2+C.)*( rho + rho**2*(X1B1+A2)*(X2B2+C.)/2 + rho**3*((X1B1+A2)**2 - 1)*((X2B2+C.)**2 - 1)/6 ); ENDDOT; ? === 3 === DOT 1-4; FRML eq3. q3. = cnorm(X1B1+A3)*cnorm(X2B2+C.)+norm(X1B1+A3)*norm(X2B2+C.)*( rho + rho**2*(X1B1+A3)*(X2B2+C.)/2 + rho**3*((X1B1+A3)**2 - 1)*((X2B2+C.)**2 - 1)/6 ); ENDDOT; ? === 4 === DOT 1-4; FRML eq4. q4. = cnorm(X1B1+A4)*cnorm(X2B2+C.)+norm(X1B1+A4)*norm(X2B2+C.)*( rho + rho**2*(X1B1+A4)*(X2B2+C.)/2 + rho**3*((X1B1+A4)**2 - 1)*((X2B2+C.)**2 - 1)/6 ); ENDDOT; ? ? ? Let's go to create equations for pij ? eqp00 is above ... eqp5. and eqp.5 are already included in likelihood ? (it is necessary) FRML eqp01 p01 = q01 - p00; FRML eqp02 p02 = q02 - q01; FRML eqp03 p03 = q03 - q02; FRML eqp04 p04 = q04 - q03; FRML eqp11 p11 = q11 - q01 - q10 + p00; FRML eqp12 p12 = q12 - q02 - q11 + q01; FRML eqp13 p13 = q13 - q03 - q12 + q02; FRML eqp14 p14 = q14 - q04 - q13 + q03; FRML eqp21 p21 = q21 - q11 - q20 + q10; FRML eqp22 p22 = q22 - q12 - q21 + q11; FRML eqp23 p23 = q23 - q13 - q22 + q12; FRML eqp24 p24 = q24 - q14 - q23 + q13; FRML eqp31 p31 = q31 - q21 - q30 + q20; FRML eqp32 p32 = q32 - q22 - q31 + q21; FRML eqp33 p33 = q33 - q23 - q32 + q22; FRML eqp34 p34 = q34 - q24 - q33 + q23; FRML eqp41 p41 = q41 - q31 - q40 + q30; FRML eqp42 p42 = q42 - q32 - q41 + q31; FRML eqp43 p43 = q43 - q33 - q42 + q32; FRML eqp44 p44 = q44 - q34 - q43 + q33; FRML eqp10 p10 = q10 - p00; FRML eqp20 p20 = q20 - q10; FRML eqp30 p30 = q30 - q20; FRML eqp40 p40 = q40 - q30; ? ? ? Starting likelihood equation have a form: FRML eqY1Y2 logL = LOG{ Y00*p00 + Y01*p01 + Y02*p02 + Y03*p03 + Y04*p04 + Y05*(cnorm(X1B1) - p00 - p01 - p02 - p03 - p04) + Y10*p10 + Y11*p11 + Y12*p12 + Y13*p13 + Y14*p14 + Y15*(cnorm(A1+X1B1)-cnorm(+X1B1) - p10 - p11 - p12 - p13 - p14) + Y20*p20 + Y21*p21 + Y22*p22 + Y23*p23 + Y24*p24 + Y25*(cnorm(A2+X1B1)-cnorm(A1+X1B1)- p20 - p21 - p22 - p23 - p24) + Y30*p30 + Y31*p31 + Y32*p32 + Y33*p33 + Y34*p34 + Y35*(cnorm(A3+X1B1)-cnorm(A2+X1B1)- p30 - p31 - p32 - p33 - p34) + Y40*p40 + Y41*p41 + Y42*p42 + Y43*p43 + Y44*p44 + Y45*(cnorm(A4+X1B1)-cnorm(A3+X1B1)- p40 - p41 - p42 - p43 - p44) + Y50*(cnorm(+X2B2) - p00 - p10 - p20 - p30 - p40) + Y51*(cnorm(C1+X2B2)-cnorm(+X2B2) - p01 - p11 - p21 - p31 - p41) + Y52*(cnorm(C2+X2B2)-cnorm(C1+X2B2) -p02 - p12 - p22 - p32 - p42) + Y53*(cnorm(C3+X2B2)-cnorm(C2+X2B2)-p03 - p13 - p23 - p33 - p43) + Y54*(cnorm(C4+X2B2)-cnorm(C3+X2B2)-p04 - p14 - p24 - p34 - p44) + Y55*(1 - cnorm(A4+X1B1) - cnorm(C4+X2B2) + q44)}; ? ? ? Variables having an influence of Y1 are: FRML X1B1 - (b0 + b1*x1 + b2*x2 + b4*x4 + b13*x13 + b14*x14); ? Variables having an influence of Y2 are: FRML X2B2 - (e0 + ey1*y1 + e1*x1 + e2*x2 + e3*x3 + e4*x4 + e15*x15 + e16*x16 + e17*x17 + e18*x18 + e19*x19); PARAM ? results of the first ordered probit y1 = ... B0 3.08437 B1 -.615433E-02 B2 -.222003 B4 .078798 B13 .125386E-02 B14 -.031427 A1 1.41680 A2 1.47321 A3 1.50369 A4 3.34789 ? results of the second ordered probit y2 = ... E0 -.742711 EY1 .084460 E1 -.247361E-03 E2 .206275 E3 .037145 E4 .408900 E15 -.261408 E16 .187293 E17 .098448 E18 -.253650 E19 .374671 C1 1.04858 C2 1.48199 C3 1.88562 C4 3.28421 ? CORRELATION COEFFICIENT rho set to rho 0.2; ? ? ? ================= SUBSTITUTION ==== ? We have to substitute equtions, and keep proper order: ? ? === list of oll auxilary equations eq: list aux = eq01 eq02 eq03 eq04 eq10 eq11 eq12 eq13 eq14 eq20 eq21 eq22 eq23 eq24 eq30 eq31 eq32 eq33 eq34 eq40 eq41 eq42 eq43 eq44; ? === list of all equations eqp. list eqp = eqp00 eqp01 eqp02 eqp03 eqp04 eqp10 eqp11 eqp12 eqp13 eqp14 eqp20 eqp21 eqp22 eqp23 eqp24 eqp30 eqp31 eqp32 eqp33 eqp34 eqp40 eqp41 eqp42 eqp43 eqp44; ? === list of all new equations neqp list new = neqp00 neqp01 neqp02 neqp03 neqp04 neqp10 neqp11 neqp12 neqp13 neqp14 neqp20 neqp21 neqp22 neqp23 neqp24 neqp30 neqp31 neqp32 neqp33 neqp34 neqp40 neqp41 neqp42 neqp43 neqp44; ? === list of all new equations nneqp list last = nneqp00 nneqp01 nneqp02 nneqp03 nneqp04 nneqp10 nneqp11 nneqp12 nneqp13 nneqp14 nneqp20 nneqp21 nneqp22 nneqp23 nneqp24 nneqp30 nneqp31 nneqp32 nneqp33 nneqp34 nneqp40 nneqp41 nneqp42 nneqp43 nneqp44; ? We do first substitutions, removing q's: DOT eqp; EQSUB (NAME= n.) . aux; ENDDOT; EQSUB (NAME=eqY1Y2) eqY1Y2 aux; EQSUB (NAME=bivar) eqY1Y2 new eqp00 X1B1 X2B2; ? ? For practical application (to save memory), it's recommended to do the ? following: ? save everything into TSP databank and run the rest from new program ? for example: ? out ww453; ? keep all; ? stop; ? ? Start with ML estimation: ? genr bivar; select ^miss(logl); ML BIVAR; ? select 1; ML bivar; ? stop; end;