$offsymxref offsymlist
**Tobit-THREE and TWO Stage GME MC
**Data: k=10 and k=90 condition numbers and Bill's data
*March 23, 1994
* FOR COMPARISON with "order4n.gms"
* Based on "tobitwo.gms"

set  i     index /1*4/
     j     index /1*3/
     k     index /1*3/
     m     index /1*50/
     g     index /1*100/
;
set  k1(k)
     j1(j)
;

parameter  x(m,i)  observed aggregate outcomes (RHS vars)
           xdata(g,m,i)
           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
           SEbt3(g) Third stage SE(beta)
           avsebt3 Third stage SE(beta)
           varsebt3 Third stage V(SE(beta))
           betat3(g,i) Third stage beta
           avbetat3(i) Third stage average betas
           varbetat3(i) Third 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
           aventropb    average Entropy value for betas
           entropyb2(g)    Entropy value for second stage betas
           aventropb2    average Entropy value for second stage 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

;
gamma=0.5;
delta=0.05;
JJ=3;
KK=3;
MM=50;
II=4;
GG=100;
**Note adjust "cut" according to data
*cut = -.8;
cut = 0;
**
parameter beta0(i) correct betas /

**Old Table 1 - Censored Coeff
**1  1
**2 -1
**3  2
**4  1
1  1
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
;

x(m,"1")=1;
x(m,"2")=normal(0,1);
x(m,"3")=normal(0,1);
x(m,"4")=normal(0,1);

**Generating data for ML (Jeff)
xdata(g,m,"2")=x(m,"2");
xdata(g,m,"3")=x(m,"3");
xdata(g,m,"4")=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"      "2"  "3"
*"1"     -10     -10  -10
*"2"       0       0    0
*"3"      10      10   10
*"1"     -10     -10  -10
*"2"     -5      -5   -5
*"3"      0       0    0
*"4"      5       5    5
*"5"      10      10   10
;

parameter alphay(j) param space for y_0 /

1  -8
2   0
3   8
/;
parameter alphae1(j) param space for errors e1 /

1  -8
2   0
3   8
/;

parameter alphae2(j) param space for errors e2 /

1  0
2  0
3  0
/;

parameter alphae3(j) param space for errors e2 /

1 -.5
2  0
3  .5
/;

*Uniform priors
*qb(i,"1")=.02;
*qb(i,"2")=.32;
*qb(i,"3")=.32;
*qb(i,"4")=.32;
*qb(i,"5")=.02;
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
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
bound3(m)        Third stage 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
*objective..   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 THIRD STAGE Second eq
bound3(m)$(y(m) eq cut)..    sum(i,sum(k,x(m,i)*alphab(k,i)*pb(i,k))) + sum(j,pe2(m,j)*alphae3(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 second second-stage /objective1
                  addb
                  adde1
                  adde2
                  b1i
                  e1m1
                  e1m2
                  bound1
                  bound2a
                        /;
model third third-stage /objective1
                  addb
                  adde1
                  adde2
                  b1i
                  e1m1
                  e1m2
                  bound1
                  bound3
                        /;
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=500000;
options reslim=80000;
* conopt optimizer option is below
options solprint=off, nlp=conopt;
*options nlp=conopt;
*options seed=2;
options domlim=20;
options iterlim=50000;
*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);

***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;

**Generating data for ML (Jeff)
xdata(g,m,"1") = y(m) + 1 ;

num(m)=0;
num(m)$(y(m) eq cut) = 1;
cens(g)=sum(m,num(m));

*ME case
solve noise1 maximizing obj using nlp;
*CE case
*solve noise1 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));
entropyb2(g) = -sum(k,sum(i,pb.l(i,k)*log(pb.l(i,k)+1.e-8)))/(log(KK)*II);

***Beginning of Third Stage GME
yhat(m)$(s(m) le cut) = sum(i,x(m,i)*b1.l(i)) + e2.l(m);
*display yhat;
*ME SECOND stage 
*solve third maximizing obj using nlp;

*betat3(g,i) = b1.l(i);
options decimals=7;
*display b1.l;
*SEbt3(g)=sum(i,power((betat3(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;
*avbetat3(i) = sum(g,betat3(g,i))/GG;
aventb(i) = sum(g,entb(g,i))/GG;
aventropb = sum(g,entropyb(g))/GG;
aventropb2 = sum(g,entropyb2(g))/GG;
avseb = sum(g,SEb(g))/GG;
avsebt = sum(g,SEbt(g))/GG;
*avsebt3 = sum(g,SEbt3(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 aventropb;
display aventropb2;
display avbeta;
display avbetat;
*display avbetat3;
display beta;
display betat;
*display betat3;
display SEb;
display SEbt;
*display SEbt3;
*display SEy2;
*display SEy0;
display avseb;
display avsebt;
*display avsebt3;
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;
*varbetat3(i)=sum(g,(betat3(g,i) - avbetat3(i))*(betat3(g,i) - avbetat3(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;
*varsebt3=sum(g,(SEbt3(g) - avsebt3)*(SEbt3(g) - avsebt3)) / 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 varbetat3;
display varprob;
display varseb;
display varsebt;
*display varsebt3;
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;