$offsymxref offsymlist **Tobit-Two Stage GME MC **Data: k=10 and k=90 condition numbers and Bill's data *October 19, 95 * Based on "tn1000.gms" set i index /1*4/ i1 index /1*5/ j index /1*3/ k index /1*5/ m index /1*30/ g index /1*10/ ; set k1(k) j1(j) ; parameter x(m,i) observed aggregate outcomes (RHS vars) xdata(g,m,i1) pb0(i,k) estimated prob for beta cut censored leve of y SEbt(g) Second stage SE(beta) avsebt Second stage SE(beta) varsebt Second stage V(SE(beta)) betat(g,i) Second stage beta avbetat(i) Second stage average betas varbetat(i) Second stage Var of betas pred(m) number of correctly predicted ypred1(m) ypred(m) yhat(m) predicted y2 from first stage estimation sumyp(g) number of misses avsumyp average number of misses avpred average number of correctly predicted sss(m) dist(m) prob(g,i,k) estimated prob for beta per sample avprob(i,k) average recovered prob for beta varprob(i,k) var of recovered prob for beta cens(g) number of censured per sample ee1(m) rand error vector u(m) for generating data u1(m) for generating data num(m) pe0(m,j) estimated prob for errors qb(i,k) prior prob for beta qy(m,j) prior prob for beta qe1(m,j) prior prob for errors qe2(m,j) prior prob for errors alphab(k,i) param space for betas * alphae(j,m) param space for errors beta(g,i) estimated betas avbeta(i) averags beta error(g,m) estimated errors error1(m) estimated errors 1 error2(m) estimated errors 2 derror(m) difference between correct and estimated errors meanerr mean value of estimated errors s0(m) pure LHS values s(m) Noisy LHS values y(m) Left hand side variable y0(m) Left hand side "correct" variable yy20(m) ystar(g,m) recovered LHS var e(m) estimated errors paramb(i,k) parameter (beta) space parame(m,j) parameter (errors) space entb(g,i) normalized entropy of estimated beta coefficients aventb(i) average nornmalized Ent of Betas sbar normalized ent of system (betas) ente(g,m) normalized entropy of estimated errors avente(m) average normalized Ent for Errors entropy Entropy value (value of objective fn) entropyb(g) Entropy value for betas entropye Entropy value for errors SEb(g) square error of betas SEb1(g) square error of betas except intercept SEy2a(g) sq error of y*_2 the minus error case SEy2b(g) sq error of y*_2 the plus error case SEy(g) sq error of y*_2 SEy0(g) sq error of PURE y*_2 avseb av sq errors of betas avseb1 av sq errors of betas except intercept varseb1 var of seb except intercept avsey av sq errors of y_2 avsey2a av sq errors of y_2 (minus case) avsey2b av sq errors of y_2 (plus case) avsey0 av sq errors of PURE y*_2 avseg av sq errors of errors avasee average sq errors of errors ave1 average val of recovered e1 SEeg(g) square error of errors ASEe(g) Average square error of errors var variance varb(i) variance of bata i vare(m) variance of error m n(m) index tot JJ number of parameters in error param space (size of j) KK number of parameters in beta param space (size of k) MM number of observations (constraints) size of m II number of coefficients beta to be estimated (size of i) GG number of MC experiments weight weight of regularization param delta delta weight factor for objective fn gamma weight factor for objective fn avy average value of y vary var of y sdy sd of y miny maxy ; gamma=0.5; delta=0.05; JJ=3; KK=5; MM=30; II=4; GG=10; **Note adjust "cut" according to data *cut = -.8; cut = 0; ** parameter beta0(i) correct betas / *1 -.75 *2 1 *3 1 1 2 2 1 3 -3 4 2 /; parameter SEe1 sum of error from y-Xb SEe1q sum of error sq from y-Xb e01(m) errors from y-Xb V(g) recovered variance avV average variance ; table x(m,i) data **NOTE X matrix is similar for 3 cases ***Bill's data set * 1 2 3 *1 1 0.0400 0.0783 *2 1 0.0800 -0.3607 *3 1 0.1200 0.9041 *4 1 0.1600 0.4323 *5 1 0.2000 0.0300 *6 1 0.2400 -0.3321 *7 1 0.2800 -0.9800 *8 1 0.3200 -0.3309 *9 1 0.3600 -0.6978 *10 1 0.4000 -0.7574 *11 1 0.4400 0.0245 *12 1 0.4800 0.0454 *13 1 0.5200 -0.0273 *14 1 0.5600 -0.9851 *15 1 0.6000 0.1877 *16 1 0.6400 -0.9579 *17 1 0.6800 -0.0500 *18 1 0.7200 0.3117 *19 1 0.7600 0.8185 *20 1 0.8000 -0.6174 *21 1 0.8400 0.3884 *22 1 0.8800 -0.9531 *23 1 0.9200 -0.7092 *24 1 0.9600 -0.2145 *25 1 1.0000 0.7504 *26 1 1.0400 -0.3051 *27 1 1.0800 0.6412 *28 1 1.1200 0.1519 *29 1 1.1600 0.8101 *30 1 1.2000 0.9129 *31 1 1.2400 0.6208 *32 1 1.2800 0.2246 *33 1 1.3200 -0.2321 *34 1 1.3600 -0.2599 *35 1 1.4000 -0.2336 *36 1 1.4400 -0.6832 *37 1 1.4800 -0.5938 *38 1 1.5200 -0.7668 *39 1 1.5600 -0.7978 *40 1 1.6000 0.3110 *41 1 1.6400 0.7477 *42 1 1.6800 0.9893 *43 1 1.7200 0.7076 *44 1 1.7600 -0.6539 *45 1 1.8000 -0.1962 *46 1 1.8400 0.0435 *47 1 1.8800 -0.8027 *48 1 1.9200 0.6889 *49 1 1.9600 -0.2986 *50 1 2.0000 0.3355 *51 1 2.0400 -0.4065 *52 1 2.0800 -0.8544 *53 1 2.1200 -0.2713 *54 1 2.1600 0.3093 *55 1 2.2000 0.0431 *56 1 2.2400 -0.9433 *57 1 2.2800 -0.0973 *58 1 2.3200 0.3408 *59 1 2.3600 0.9349 *60 1 2.4000 0.6199 *61 1 2.4400 -0.3389 *62 1 2.4800 0.6220 *63 1 2.5200 0.2959 *64 1 2.5600 0.1617 *65 1 2.6000 -0.5786 *66 1 2.6400 -0.5323 *67 1 2.6800 0.8043 *68 1 2.7200 -0.8508 *69 1 2.7600 0.7873 *70 1 2.8000 0.2097 *71 1 2.8400 -0.7478 *72 1 2.8800 0.5661 *73 1 2.9200 0.2068 *74 1 2.9600 0.4763 *75 1 3.0000 0.2391 *76 1 3.0400 -0.2794 *77 1 3.0800 -0.2303 *78 1 3.1200 -0.1173 *79 1 3.1600 0.5545 *80 1 3.2000 0.6799 *81 1 3.2400 0.8093 *82 1 3.2800 0.2088 *83 1 3.3200 0.6892 *84 1 3.3600 -0.5798 *85 1 3.4000 0.7631 *86 1 3.4400 0.8722 *87 1 3.4800 0.3443 *88 1 3.5200 -0.2180 *89 1 3.5600 -0.6593 *90 1 3.6000 -0.1709 *91 1 3.6400 -0.9701 *92 1 3.6800 0.1404 *93 1 3.7200 -0.9288 *94 1 3.7600 0.0663 *95 1 3.8000 -0.4511 *96 1 3.8400 0.0917 *97 1 3.8800 0.9582 *98 1 3.9200 -0.6056 *99 1 3.9600 -0.5759 *100 1 4.0000 -0.0377 **k=10 case (same as in multinomial paper) * 1 2 3 4 * *1 -0.052359 -0.264344 -0.058610 0.111364 *2 -0.098536 0.069810 0.060271 0.054685 *3 -0.210793 -0.270722 0.036207 -0.210084 *4 -0.134768 0.065137 -0.085801 -0.140940 *5 -0.076803 -0.146102 -0.155144 0.063655 *6 -0.154308 -0.081523 0.074593 -0.155242 *7 0.083309 -0.037870 0.072019 0.205073 *8 0.062334 0.323160 -0.097266 -0.014514 *9 -0.105366 -0.034904 0.015750 -0.531353 *10 -0.151496 -0.186632 -0.006663 -0.241347 *11 -0.287023 -0.407090 -0.030673 0.090842 *12 0.276012 0.246759 0.177114 -0.135687 *13 -0.285272 -0.043978 -0.257213 0.236365 *14 0.091886 -0.098755 0.264811 0.009687 *15 -0.279293 0.010607 -0.242569 0.047538 *16 -0.388233 -0.039411 0.014311 0.035271 *17 -0.163587 0.118327 0.015172 -0.301582 *18 -0.144417 -0.339407 0.074788 0.111695 *19 0.145756 0.185476 0.005867 -0.269005 *20 -0.039196 -0.040142 -0.311574 0.377680 *21 -0.108334 0.173895 -0.218945 -0.054025 *22 -0.161187 0.033520 -0.161315 0.129206 *23 0.258251 0.071409 0.104663 0.069091 *24 -0.174804 0.031342 -0.316047 0.108814 *25 0.210636 -0.295945 0.223393 0.013047 *26 -0.064620 -0.220928 0.148018 -0.047185 *27 0.066566 -0.048299 0.182580 0.243102 *28 -0.000621 0.078313 -0.101813 0.019408 *29 0.005704 0.108200 -0.035409 -0.118853 *30 0.444332 0.204394 0.494187 -0.102119 **k=90 condition number (same as multinomial paper) *1 -0.092828 -0.243467 -0.028175 0.128663 *2 -0.043228 0.030561 0.002750 0.031111 *3 -0.222735 -0.268274 0.039670 -0.204956 *4 -0.129042 0.058481 -0.095608 -0.143364 *5 -0.120547 -0.123416 -0.122069 0.082353 *6 -0.129731 -0.100305 0.047039 -0.165708 *7 0.096314 -0.046086 0.059999 0.199523 *8 0.079360 0.315392 -0.108560 -0.021799 *9 -0.134270 -0.011897 0.049519 -0.519050 *10 -0.175767 -0.173573 0.012391 -0.230976 *11 -0.295540 -0.415305 -0.043006 0.094564 *12 0.292501 0.250470 0.182870 -0.142813 *13 -0.274977 -0.066494 -0.290523 0.232074 *14 0.124225 -0.114996 0.241149 -0.004139 *15 -0.274563 -0.004848 -0.265468 0.045599 *16 -0.302050 -0.109422 -0.088476 -0.001406 *17 -0.130650 0.096553 -0.016707 -0.315631 *18 -0.142852 -0.348282 0.061622 0.111078 *19 0.122691 0.210565 0.042816 -0.259230 *20 -0.083727 -0.021097 -0.283920 0.396741 *21 -0.115755 0.174845 -0.217642 -0.050835 *22 -0.149544 0.017258 -0.185312 0.124294 *23 0.247686 0.087450 0.128347 0.073540 *24 -0.203038 0.039477 -0.304369 0.120925 *25 0.170849 -0.262710 0.272203 0.029973 *26 -0.055388 -0.228384 0.137072 -0.051114 *27 0.114012 -0.079140 0.137439 0.222861 *28 -0.011363 0.084531 -0.092729 0.023996 *29 0.003259 0.112042 -0.029735 -0.117825 *30 0.502239 0.189007 0.472166 -0.126964 *T = 30 K = 4 K(X'X) = 1.0 *New experimints with orthonormal case for revised cens paper, Oct 95 1 2 3 4 1 0.0872 0.4224 0.0340 0.1616 2 0.2305 -0.0517 0.2540 -0.2943 3 0.0465 -0.2464 -0.2995 0.0887 4 0.1857 -0.0447 -0.6032 0.1893 5 -0.0020 -0.2116 0.0030 0.0509 6 -0.1548 0.1933 -0.1359 0.2666 7 0.4572 -0.1171 -0.1067 0.0108 8 -0.1160 0.2803 -0.1657 0.1204 9 -0.0857 -0.1544 -0.1024 0.3533 10 -0.1834 0.2025 0.0782 0.0365 11 -0.0599 -0.0706 0.1314 -0.0777 12 0.2350 0.1311 0.1567 -0.0078 13 0.0245 -0.0016 -0.0658 -0.1874 14 -0.0881 0.1665 -0.0561 -0.1503 15 0.2378 -0.0158 0.1364 0.1757 16 -0.1673 0.0063 -0.2379 -0.3473 17 -0.3065 -0.1596 -0.0747 0.1185 18 0.1677 -0.1085 0.0522 0.0886 19 -0.1960 0.0918 0.0420 -0.1254 20 0.0162 -0.0975 0.0527 0.3841 21 -0.1958 0.2969 -0.0934 0.1065 22 0.0887 -0.0319 -0.2095 -0.1879 23 -0.0889 0.0338 -0.2134 -0.3736 24 0.2888 0.5100 -0.0604 -0.0228 25 -0.1992 -0.0387 -0.1887 -0.0837 26 0.2625 0.1372 -0.1446 -0.0247 27 -0.0956 0.0300 0.1501 -0.0206 28 0.1794 0.0401 -0.2317 -0.0841 29 -0.0073 -0.0298 -0.1852 -0.1561 30 0.1236 -0.1935 0.0273 -0.0497 ; *x(m,"1")=0; *x(m,"3")=0; *x(m,"1")=normal(0,1); *x(m,"3")=normal(0,1); **Data for ML-Jeff *xdata(g,m,i1) =0; *xdata(g,m,"2") = x(m,"1"); *xdata(g,m,"3") = x(m,"2"); *xdata(g,m,"4") = x(m,"3"); *xdata(g,m,"5") = x(m,"4"); table alphab(k,i) parameter space for betas "1" "2" "3" "4" *"1" -10 -10 -10 -10 *"2" 0 0 0 0 *"3" 10 10 10 10 "1" -10 -10 -10 -10 "2" -5 -5 -5 -5 "3" 0 0 0 0 "4" 5 5 5 5 "5" 10 10 10 10 *"1" -100 -100 -100 -100 *"2" 0 0 0 0 *"3" 100 100 100 100 ; parameter alphay(j) param space for y_0 / 1 -5 2 0 3 5 *1 -3 *2 0 *3 3 /; parameter alphae1(j) param space for errors e1 / 1 -5 2 0 3 5 *1 -3 *2 0 *3 3 /; parameter alphae2(j) param space for errors e2 / 1 -1 2 0 3 1 /; *Uniform priors *qb(i,"1")=.1; *qb(i,"2")=.2; *qb(i,"3")=.4; *qb(i,"4")=.2; *qb(i,"5")=.1; *qb(i,"1")=.02; *qb(i,"2")=.32; *qb(i,"3")=.32; *qb(i,"4")=.32; *qb(i,"5")=.02; **qb(i,"1")=.02; **qb(i,"2")=.08; **qb(i,"3")=.3; **qb(i,"4")=.58; **qb(i,"5")=.02; **qb("3","1")=.02; **qb("3","2")=.58; **qb("3","3")=.3; **qb("3","4")=.08; **qb("3","5")=.02; *qb(i,"1")=0; *qb(i,"2")=0.5; *qb(i,"3")=0.5; *qb("3","1")=0.5; *qb("3","2")=0.5; *qb("3","3")=0; qb(i,k)=1/KK; qe1(m,j)=1/JJ; qe2(m,j)=1/JJ; qy(m,j)=1/JJ; *qe(m,"6")=1; variables theta(i) coeff for OLS b1(i) coefficients e1(m) errors e1 e2(m) errors e2 yy0(m) estimated y0's theta(i) xx(m) xx1(m) xx2(m) ; positive variables pb(i,k) probabilities py(m,j) probabilities pe1(m,j) probabilities pe2(m,j) probabilities ; variables obj objective function ; equations objective objective function objectce objective function objective1 objective function objective2 objective function obj1 objective for OLS obj2 objective for OLS obj3 objective for OLS addb(i) adding up betas addy(m) adding up y0 adde1(m) adding up errors adde2(m) adding up errors *add1 adding up b1i(i) b1k coefficients eqy0(m) y0 coeff e1m1(m) e1m coefficients e1m2(m) e1m coefficients noise adding noise to zero bound1(m) stochastic conservation bound2(m) stochastic conservation bound2a(m) stochastic conservation boundg(m) stochastic conservation moment1(m) part of moment const moment2(m) part of moment const mom1(i) moment consistency constraint mom2(i) moment consistency constraint ; *The CE objective for GENERAL problem with NOISE and WEIGHTS *objective.. obj =e= (1-delta) * sum(k,sum(i,pb(i,k)*log((pb(i,k)+1.e-8)/(qb(i,k)+1.e-8)))) + delta * sum(j,sum(m,pe(m,j)*log((pe(m,j)+1.e-8)/(qe(m,j)+1.e-8)))); *The Tobit ME objective for GENERAL problem with NOISE **Four criterion case objective.. obj =e= -sum(k,sum(i,pb(i,k)*log(pb(i,k)+1.e-8))) -sum(j,sum(m,py(m,j)*log(py(m,j)+1.e-8))) - sum(j,sum(m,pe1(m,j)*log(pe1(m,j)+1.e-8))) - sum(j,sum(m,pe2(m,j)*log(pe2(m,j)+1.e-8))); *The Tobit ME objective for SECOND STAGE objective1.. obj =e= -sum(k,sum(i,pb(i,k)*log(pb(i,k)+1.e-8))) - sum(j,sum(m,pe1(m,j)*log(pe1(m,j)+1.e-8))) - sum(j,sum(m,pe2(m,j)*log(pe2(m,j)+1.e-8))); *The Tobit ME objective for REGULAR (continuous) NOISE objective2.. obj =e= -sum(k,sum(i,pb(i,k)*log(pb(i,k)+1.e-8))) - sum(j,sum(m,pe1(m,j)*log(pe1(m,j)+1.e-8))); *The Tobit CE objective for GENERAL problem with NOISE **Four criterion case objectce.. obj =e= sum(k,sum(i,pb(i,k)*log((pb(i,k)+1.e-8)/(qb(i,k)+1.e-8))))+sum(j,sum(m,py(m,j)*log((py(m,j)+1.e-8)/(qy(m,j)+1.e-6)))) + sum(j,sum(m,pe1(m,j) *log((pe1(m,j)+1.e-8)/(qe1(m,j)+1.e-6)))) + sum(j,sum(m,pe2(m,j)*log((pe2(m,j)+1.e-8)/(qe2(m,j)+1.e-6)))); ****End CE Note: CE objective is composed of TWO lines addb(i).. sum(k,pb(i,k)) =e= 1; adde1(m)$(y(m) gt cut).. sum(j,pe1(m,j)) =e= 1; adde2(m)$(y(m) eq cut).. sum(j,pe2(m,j)) =e= 1; addy(m)$(y(m) eq cut).. sum(j,py(m,j)) =e= 1; b1i(i).. sum(k,alphab(k,i)*pb(i,k)) =e= b1(i); e1m1(m)$(y(m) gt cut).. sum(j,alphae1(j)*pe1(m,j)) =e= e1(m); e1m2(m)$(y(m) eq cut).. sum(j,alphae2(j)*pe2(m,j)) =e= e2(m); eqy0(m)$(y(m) eq cut).. sum(j,alphay(j)*py(m,j)) =e= yy0(m); bound1(m)$(y(m) gt cut).. sum(i,sum(k,x(m,i)*alphab(k,i)*pb(i,k))) + sum(j,pe1(m,j)*alphae1(j)) =e= y(m); ****************************** ****Second eq. for Tobit First Stage Estimation bound2(m)$(y(m) eq cut).. sum(i,sum(k,x(m,i)*alphab(k,i)*pb(i,k))) =l= yy0(m); ****************************** ****************************** ****Consistency for SECOND STAGE Second eq bound2a(m)$(y(m) eq cut).. sum(i,sum(k,x(m,i)*alphab(k,i)*pb(i,k))) + sum(j,pe2(m,j)*alphae2(j)) =l= yhat(m); ****************************** ****************************** ****Consistency for GENERAL NOISE CASE boundg(m).. sum(i,sum(k,x(m,i)*alphab(k,i)*pb(i,k))) + sum(j,pe1(m,j)*alphae1(j)) =e= y(m); ****************************** ***OLS objective obj1.. obj =e= - sum(m, (y(m) - sum(i,x(m,i)*theta(i)))*(y(m) - sum(i,x(m,i)*theta(i)))); ***OLS objective for SECOND stage Tobit obj2.. obj =e= - sum(m, (y(m) - sum(i,x(m,i)*theta(i)))*(y(m) - sum(i,x(m,i)*theta(i)))); ***OLS objective for continuous data obj3.. obj =e= - sum(m, (s(m) - sum(i,x(m,i)*theta(i)))*(s(m) - sum(i,x(m,i)*theta(i)))); model noise1 noise-regular /objective addb adde1 adde2 addy eqy0 b1i e1m1 e1m2 bound1 bound2 /; model noisece noise-GCE /objectce addb adde1 adde2 addy eqy0 b1i e1m1 e1m2 bound1 bound2 /; model second second-stage /objective1 addb adde1 adde2 b1i e1m1 e1m2 bound1 bound2a /; model regnoise regular noise / objective2 addb adde1 b1i e1m1 boundg /; model ols least squares /obj1/; model ols1 least squares /obj2/; model ols2 least squares /obj3/; * conopt and minos are two non-linear optimization packages * minos is used by default. options limrow=0,limcol=0; options work=400000; * conopt optimizer option is below options solprint=off, nlp=conopt; *options nlp=conopt; *options seed=5; options iterlim=5000; options domlim=20; *options solprint=off; *initialization pb.l(i,k)=1/KK; pe1.l(m,j)=1/JJ; pe2.l(m,j)=1/JJ; py.l(m,j)=1/JJ; e1.l(m)=0; e2.l(m)=0; *lambda1.l=1; *lambda2.l=1; **Starts MC loop loop(g, s0(m)=sum(i,x(m,i)*beta0(i)); *display s0; **normal generator ee1(m) = normal(0,1); s(m) = s0(m) + ee1(m); options decimals=7; *display ee1; **t dist generator *u(m)=power(normal(0,1),2)+power(normal(0,1),2)+power(normal(0,1),2); *u1(m)=ee1(m)/sqrt(u(m)); *s(m) = s0(m) + u1(m); **Chi sq generator *u(m)=power(normal(0,1),2)+power(normal(0,1),2)+power(normal(0,1),2)+power(normal(0,1),2); *u1(m)=(u(m)-4)/sqrt(8); *s(m) = s0(m) + u1(m); *display s; y(m)=s(m); **Creating endogenous empirical var of the y **and create the error bounds accordingly avy = sum(m,y(m))/MM; vary=sum(m,(y(m)-avy)*(y(m)-avy))/(MM-1); sdy=sqrt(vary); *display s0; *display y; *display avy; *display vary; *display sdy; ***the ZERO case *y(m)$(s(m) gt 0) = s(m); **For 30% k=10 case *cut = -.8; cut = 0; **For 30% k=90 case *cut = -.6; **For 50% k=10 and 90 cases *cut = -.1; y(m)$(s(m) le cut) = cut; num(m)=0; num(m)$(y(m) eq cut) = 1; cens(g)=sum(m,num(m)); **Calculating the sd for observed and not observed y maxy = max( y("1"),y("2"),y("3"),y("4"),y("5"),y("6"),y("7"),y("8"),y("9"),y("10"), y("11"),y("12"),y("13"),y("14"),y("15"),y("16"),y("17"),y("18"),y("19"),y("20"), y("21"),y("22"),y("23"),y("24"),y("25"),y("26"),y("27"),y("28"),y("29"),y("30") ); miny = -((cens(g)/MM) * maxy); vary = ((maxy - miny)**2) / 12; sdy=sqrt(vary); avy = (maxy + miny)/2; *display maxy; *display miny; *display avy; *display vary; *display sdy; **Endogenously determined errors alphae1("1")=-3*sdy; alphae1("3")=3*sdy; alphay(j)=alphae1(j); *** **Data for ML-Jeff *xdata(g,m,i1) =0; *xdata(g,m,"2") = x(m,"1"); *xdata(g,m,"3") = x(m,"2"); *xdata(g,m,"4") = x(m,"3"); *xdata(g,m,"5") = x(m,"4"); *xdata(g,m,"1")=y(m)+1; *ME case solve noise1 maximizing obj using nlp; *CE case *solve noisece minimizing obj using nlp; *solve ols maximizing obj using nlp; *options decimals=7; *display theta.l; options decimals=7; *display x; **display y0; *display y; prob(g,i,k)=pb.l(i,k); options decimals=7; *display alphab; beta(g,i) = sum(k, alphab(k,i)*pb.l(i,k)); options decimals=7; **display beta; yy20(m)$(s(m) le cut) = sum(i,x(m,i)*b1.l(i)) -yy0.l(m); sss(m)=s(m); yy20(m)$(s(m) gt cut) = 0; sss(m)$(s(m) gt cut) = 0; dist(m)=yy20(m)-sss(m); SEy2a(g) = sum(m,power(dist(m),2)); *Predicting the recovered negative y's num(m)=0; pred(m)=0; ***pred(m)$(yy20(m) le cut) = 1; ***pred(m)$(yy20(m) gt cut) = 0; ***num(m)$(y(m) eq cut) = 1; ***num(m)$(y(m) gt cut) = 0; ***ypred1(m)=0; ***ypred1(m)=pred(m)-num(m); ***display ypred1; ***ypred(m)=0; ***ypred(m)$(ypred1(m) ne 0) = 1; ***sumyp(g)=sum(m,ypred(m)); pred(m)$(yy20(m) le cut) = 0; pred(m)$(yy20(m) gt cut) = 1; display pred; sumyp(g)=sum(m,pred(m)); options decimals=7; *display sumyp; *display y; *display s; *display yy20; *display dist; **End predicting *SEy2a(g) = sum(m,power((yy20(m)-sss(m)),2)); *display yy20; *display dist; *display SEy2a; **Least Sq case **yy20(m)=sum(i,x(m,i)*theta.l(i)); **SEy2(g) = sum(m,power((-yy20(m)-s(m)),2)); **display SEy2; *SEy2(g) = sum(m,power((yy20(m)-s(m)),2)); **End LS case options decimals=7; *display SEy2; ******************* ***JUST for OLS *beta(g,i)=theta.l(i); *** end ols special ******************* options decimals=7; *display gamma; *display delta; *SEb(g)=sum(i,(beta(g,i)-beta0(i))*(beta(g,i)-beta0(i))); SEb(g)=sum(i,power((beta(g,i)-beta0(i)),2)); SEb1(g)=power((beta(g,"2")-beta0("2")),2) + power((beta(g,"3")-beta0("3")),2); **SEb1(g)=power((beta(g,"2")-beta0("2")),2) + power((beta(g,"3")-beta0("3")),2)+power((beta(g,"4")-beta0("4")),2); ** Below is SEb for OLS case *SEb=sum(i,(theta.l(i)-beta0(i))*(theta.l(i)-beta0(i))); options decimals=7; *display SEb; **paramb(i,k)=pb.l(i,k)*alphab(k,i) **options decimals=7; **display paramb; **pe0(m,j)=pe.l(m,j); options decimals=7; **display pe0; options decimals=7; *display alphae1; *display alphae2; *MATRIX case *error(m) = sum(j, alphae(j,m)*pe.l(m,j)); *VECTOR case **error1(m) = sum(j, alphae1(j)*pe.l(m,j)); options decimals=7; *display e1.l; *display e2.l; *display pe1.l; *display pe2.l; n(m)$(y(m) ne 0) = 1; n(m)$(y(m) eq 0) = 0; tot = sum(m,n(m)); ave1 = sum(m,e1.l(m))/tot; *error(m) = y(m) - sum(i,x(m,i)*beta(g,i)); **options decimals=7; **display error; *e0(m)=y(m)-y0(m); options decimals=7; **display e0; **SEe=sum(m,(error(m)-e0(m))*(error(m)-e0(m))); options decimals=7; **display SEe; **ASEe=SEe/MM; options decimals=7; **display ASEe; ****Calculate sum of errors square of the esimated errors **The OLS case *SEe1=sum(m,(y(m)-sum(i,x(m,i)*theta.l(i)))); *SEe1q=sum(m,(y(m)-sum(i,x(m,i)*theta.l(i)))*(y(m)-sum(i,x(m,i)*theta.l(i)))); *display SEe1; *display SEe1q; *e01(m)=y(m)-sum(i,x(m,i)*theta.l(i)); *SEe1=sum(m,e01(m)); *SEe1q=sum(m,e01(m)*e01(m)); *options decimals=7; *display e01; *display SEe1; *display SEe1q; **End ols case ***For ME cases **SEe1=sum(m,(y(m)-sum(i,x(m,i)*beta(i)))); **SEe1q=sum(m,(y(m)-sum(i,x(m,i)*beta(i)))*(y(m)-sum(i,x(m,i)*beta(i)))); **display SEe1; **display SEe1q; ystar(g,m)$(y(m) ne 0) = y(m); *ystar(m)$(y(m) eq 0) = e2.l(m); **Regular case ystar(g,m)$(y(m) eq 0) = sum(i,x(m,i)*beta(g,i)) + e2.l(m); ************************* ***With ave1 for calculation of recovered y*_2 *ystar(g,m)$(y(m) eq 0) = sum(i,x(m,i)*beta(g,i)) - ave1; ************************* *display ave1; SEy(g)=sum(m,(ystar(g,m)-s(m))*(ystar(g,m)-s(m))); SEy0(g)=sum(m,(ystar(g,m)-s0(m))*(ystar(g,m)-s0(m))); e01(m)=ystar(g,m)-sum(i,x(m,i)*beta(g,i)); SEe1(g)=sum(m,e01(m)); *SEe1q(g)=sum(m,e01(m)*e01(m)); error(g,m) = y(m) - sum(i,x(m,i)*beta(g,i)); SEe1q(g)=sum(m,error(g,m)*error(g,m)); options decimals=7; *display ystar; *display n; *display tot; *display s; *display s0; *display e01; *display SEe1; *display SEe1q; V(g) = SEe1q(g)/(MM); options decimals=7; *display V; ******* **derror(m) = e0(m) - error(m); options decimals=7; **display derror; *var = sum(m,derror(m)*derror(m)) / (MM-II); *Calculating mean value of errors **meanerr = sum(m,error(m))/MM; options decimals=7; **display meanerr; **var = sum(m,(error(m)-meanerr)*(error(m)-meanerr)) / (MM-1); options decimals=7; **display var; *MATRIX case *parame(m,j)=pe.l(m,j)*alphae(j,m) *VECTOR case **parame(m,j)=pe.l(m,j)*alphae(j) options decimals=7; **display parame; entb(g,i) = -sum(k,pb.l(i,k)*log(1.e-7+pb.l(i,k)))/log(KK); options decimals=7; *display entb; **ente(m) = -sum(j,pe.l(m,j)*log(1.e-7+pe.l(m,j)))/log(JJ); options decimals=7; **display ente; **entropy of ME level entropyb(g) = -sum(k,sum(i,pb.l(i,k)*log(pb.l(i,k)+1.e-8)))/(log(KK)*II); options decimals=7; *display entropyb; **entropye = -sum(j,sum(m,pe.l(m,j)*log(pe.l(m,j)+1.e-8)))/(log(JJ)*MM); options decimals=7; **display entropye; entropy = obj.l; options decimals=7; *display entropy; **entropy = entropyb + entropye; options decimals=7; **display entropy; *varb(i) = sum(k, (paramb(i,k) - beta(i))*(paramb(i,k) - beta(i))*pb.l(i,k)); ***varb(i) = sum(k, (alphab(k,i) - beta(i))*(alphab(k,i) - beta(i))*pb.l(i,k)); options decimals=7; ***display varb; *varb(i) = sum(k, (paramb(i,k) - beta(i))*(paramb(i,k) - beta(i))) / (KK-1); *options decimals=7; *display varb; *vare(m) = sum(j, (alphae(j,m) - error(m))*(alphae(j,m) - error(m))*pe.l(m,j)); **vare(m) = sum(j, (alphae(j) - error(m))*(alphae(j) - error(m))*pe.l(m,j)); options decimals=7; **display vare; *vare(m) = sum(j, (parame(m,j) - error(m))*(parame(m,j) - error(m))) / (JJ-1); *options decimals=7; *display vare; *var(k) = sum(j, ((param(k,j) - beta(k))*(param(k,j) - beta(k)))*p.l(k,j)); *display var; *Note: below is the cross-entropy for each beta *ent(k) = (sum(j,p.l(k,j)*log((1.e-7+p.l(k,j))/q(k))))/log(JJ); *options decimals=7; *display ent; *entropy = sum(k,sum(j,p.l(k,j)*log((1.e-7+p.l(k,j))/1.e-7+q(k,j)))); *entropy = -sum(k,sum(j,p.l(k,j)*log(1.e-7+p.l(k,j)))); *options decimals=7; *display entropy; *options decimals=7; *display bound.m; *checking the partition functions, etc *p0(k,j) = q(k,j)*exp(sum(i,bound.m(i)*x(i,k)*alpha(j,k)) + b1.m(k)*alpha(j,k)); *omega(k) = sum(j,p0(k,j)); *p0(k,j) = p0(k,j)/omega(k); *beta(k) = sum(j,alpha(j,k)*p0(k,j)); *options decimals=7; *display p0; *options decimals=7; *display beta; *Calculating priors for next expermints *entb(i) = beta(i)/KK; *alphab("1",i) = beta(i) - 4*abs(entb(i)); *alphab("2",i) = beta(i) - 3*abs(entb(i)); *alphab("3",i) = beta(i) - 2*abs(entb(i)); *alphab("4",i) = beta(i) - 1*abs(entb(i)); *alphab("5",i) = beta(i); *alphab("6",i) = beta(i) + 1*abs(entb(i)); *alphab("7",i) = beta(i) + 2*abs(entb(i)); *alphab("8",i) = beta(i) + 3*abs(entb(i)); *alphab("9",i) = beta(i) + 4*abs(entb(i)); options decimals=7; *display entlim.m; *****Checking lambda - Sbar relatopnship **SEe=sum(m,error(m)); options decimals=7; **display SEe; **display addb.m; **entb("2") = (sum(m,2*x(m,"2")*(y(m) - sum(i,x(m,i)*beta(i)))) + addb.m("2") +b1i.m("2"))/(+1+log(beta("2"))); **entb("3") = (sum(m,2*x(m,"3")*(y(m) - sum(i,x(m,i)*beta(i)))) + addb.m("3") +b1i.m("3"))/(+1+log(beta("3"))); ***Beginning of Second Stage GME yhat(m)=yy20(m); *display yhat; *ME SECOND stage solve second maximizing obj using nlp; betat(g,i) = b1.l(i); options decimals=7; *display b1.l; SEbt(g)=sum(i,power((betat(g,i)-beta0(i)),2)); ****Regular OLS on censured data *solve ols maximizing obj using nlp; *options decimals=7; *display theta.l; *SEbsec=sum(i,power((theta.l(i)-beta0(i)),2)); *display SEbsec; *OLS SECOND stage *solve ols1 maximizing obj using nlp; *y(m)=s(m); *yhat(m)$(s(m) le cut) = yhat(m); *options decimals=7; *display theta.l; *SEbsec=sum(i,power((theta.l(i)-beta0(i)),2)); *display SEbsec; ***Regular GME case *solve regnoise maximizing obj using nlp; *options decimals=7; *display b1.l; *SEbsec=sum(i,power((b1.l(i)-beta0(i)),2)); *display SEbsec; ****Regular OLS on "continueous" y's *solve ols2 maximizing obj using nlp; *options decimals=7; *display theta.l; *SEbsec=sum(i,power((theta.l(i)-beta0(i)),2)); *display SEbsec; **End MC loop ); avbeta(i) = sum(g,beta(g,i))/GG; avbetat(i) = sum(g,betat(g,i))/GG; aventb(i) = sum(g,entb(g,i))/GG; avseb = sum(g,SEb(g))/GG; avsebt = sum(g,SEbt(g))/GG; avsumyp = sum(g,sumyp(g))/GG; sbar = sum(g,entropyb(g))/GG; avprob(i,k) = sum(g,prob(g,i,k))/GG; avseb1 = sum(g,SEb1(g))/GG; avsey2a = sum(g,SEy2a(g))/GG; avV = sum(g,V(g))/GG; options decimals=7; *display avprob; *display sbar; display aventb; display avbeta; display avbetat; display beta; display betat; *display SEb; *display SEbt; *display SEy2; *display SEy0; display avseb; display avsebt; *display avseb1; *display avsey2a; *display avsey; *display avsey0; *display V; *display avV; parameter varbeta(i) Variance of Betas varseb Variance of SEbeta varsumyp Variance of SEbeta varsey2a Var of SEy2a varsey2b Var of SEy2b varsey Var of SEy varsey0 Var of SEy0 varv Variance of estimated var avcens average number of zero cases ; varbeta(i)=sum(g,(beta(g,i) - avbeta(i))*(beta(g,i) - avbeta(i))) / GG; varbetat(i)=sum(g,(betat(g,i) - avbetat(i))*(betat(g,i) - avbetat(i))) / GG; varprob(i,k)=sum(g,(prob(g,i,k) - avprob(i,k))*(prob(g,i,k) - avprob(i,k))) / GG; varv=sum(g,(V(g) - avV)*(V(g) - avV)) / GG; varseb=sum(g,(SEb(g) - avseb)*(SEb(g) - avseb)) / GG; varsebt=sum(g,(SEbt(g) - avsebt)*(SEbt(g) - avsebt)) / GG; varsumyp=sum(g,(sumyp(g) - avsumyp)*(sumyp(g) - avsumyp)) / GG; varseb1=sum(g,(SEb1(g) - avseb1)*(SEb1(g) - avseb1)) / GG; varsey2a=sum(g,(SEy2a(g) - avsey2a)*(SEy2a(g) - avsey2a)) / GG; **varsey=sum(g,(SEy(g) - avsey)*(SEy(g) - avsey)) / GG; **varsey0=sum(g,(SEy0(g) - avsey0)*(SEy0(g) - avsey0)) / GG; options decimals=7; display varbeta; display varbetat; display varprob; display varseb; display varsebt; *display varsumyp; *display varseb1; *display varsey2a; *display varsey; *display varsey0; *display varv; avcens = sum(g,cens(g))/GG; display avcens; display avsumyp; *display sumyp; *display cens; *display xdata;