The following is a GAMS program for multinomial/binomial case.
Note: its a very long program and has a lot of features that you
probably don't need. Before you run it, you must input your data!
Look at
1) model "base", the primal;
2) model "dual" the GME dual;
3) model "dualw" GME dual with weights.
Note in the dual model there is an additional eq called "testx". If you
activate this eq. for some var it holds the Lagrange multiplier at the
level you wish (or inequality, etc.).
Also pay attention after you choose the relevant GME model
make sure you print the correct Betas (look after the solve commands).
Finally, the program calculates marginal effects for comntinuous and discrete
variables and other tests.
**************************************************************************
$offsymxref offsymlist
**The DUAL QC case (note: no weights)
* Note, creating the data based on Bill Griffith's method
* (Nov, 1994) with 3 categories.
**GME logit for Bank 1 experiment
**August 1998
**Based on Koop9 JASA example article
**Note use index 0 for the k index if want to add a var
set i index /1*350/
j index /1*3/
j0 index /1*3/
j1 index /1*3/
j2 index /1*3/
k index /1*10/
k2 index /
edti
dti
hdti
ltvd
ltv
fscore
loannumb
bad
gfund
pub
exp
hisp
black
denied/
k5 index /
deny1
deny2
deny3
cd
inc
l999
appval
m
f
lsize/
k9 index/
INSUFND
COAPPL
REFIDUM/
h index /1*9/
k8 index for discrete partials /0*1/
h0 index /1*3/
h1 index /1*3/
l index /1*3/
w index /1*3/
d index /1*3/
n index /1*2/
d1 index /1*2/
g index /1*1/
;
set i1(n)
g1(g)
;
alias (k1,k);
alias (k3,k);
i1(n)=yes;
*i1("1")=no;
parameter x(i,k) observed aggregate outcomes (inputs outputs etc.)
basic(i,k2) bank1 basic data
basic1(i,k5) bank1 second set of basic data
basic2(i,k9) bank1 basic data
datab(i,k5) bank1 second set of basic data
b0(i) left hand side var
pseudor Pseudo R-Sq for sample
xdata(i,k) keep data y and x
ind(i) index
alpha(j,k) parameter space for each RHS var and observation
alphag(n,j,k) parameter space for each RHS var and observation
psi(i,n)
w0(i,n,h)
ei0(i,n)
alphaeg(h,n) parameter space for each errors
alphab0(j0,k) parameter space for RHS vars correspond to y=0
alphab1(j1,k) parameter space for RHS var correspond to y=1
alphab2(j2,k) parameter space for RHS var correspond to y=2
qb(n,k,j) prior probabilities for betas
qb0(k,j0) prior probabilities for beta0
qb1(k,j1) prior probabilities for beta1
qd1(n,d) prior prob. for delta1
qd2(n,d) prior prob. for delta2
qe(i,h) prior probabilities for errors
qe0(i,h0) prior probabilities for errors of beta0
qee1(i,h1) prior probabilities for errors of beta1
qw(i,w) prior probabilities for x'beta
qe1(i,l) prior probabilities for errors
error(i) estimated errors
error1(i) estimated errors for betas
a no. of parameters (length of alpha = max of j)
beta(g,n,k) estimated espected value of oefficients
omega(k) partition function for k var
omeg1(k) for partials
delta Weight factor for objective functional
sume(i) normalization for ME dist
denom(i) denominator
p1(i) temporary for data generation
p2(i) temporary for data generation
p3(i) temporary for data generation
limit(i) temporary for data generation
u(i) temporary for data generation
p0(k,j) estimated probabilities
recp(i,n) recovered p's
recw(i,n,h) recovered w's
omeg(i) partition fn for i
psi(i,n) partition fn for w
tablep(i,n,n) table for predictin results
tablepred(n,n) predicted vs. exact cases
param(k,j) parameter (beta) space
paramw(i,w) parameter (x'beta) space
entropy Entropy value (value of objective fn)
MSE(g) mean sq. error
entp(g) normalized ent for p's
entpi(g,i) normalized ent for each i
entpe(g) normalized ent for errors
aventp average ent for p's
varentp var of ent for p's
aventpi(i) average ent per i
varentpi(i) var of ent per i
aventpe average ent per errors
varentpe var ent per errors
avmse average MSE
varmse var MSE
dist(g,n,k) distance between correct and estimated betas
var(k) variance of coefficients
ent(n,k) normalized entropy of coefficients
ypred(i) predicted LHS var
pred(i) predicting var
ypred1(i) predicted minus correct y
sumyp(g) total number of misses
minmiss min number of misses
maxmiss max number of misses
predprop predicted proportion of ones
exactprop correct proportion of ones
xbeta(n,i) estimated x'beta
exbeta(n,i) estimated exp(x'beta)
xbeta1(n,i) estimated recovered probabilities
xbetam(i) max xbeta
xbetae(n,i) estimated x'beta+estimated errors
SEbeta(i) square error for each observation
SEbeta1 square error for all observations
avbeta(n,k) average values of betas
varbeta(n,k) var of betas
totvar variance of betas
avtablep(n,n) average value of prediction vs correct
vartablep(n,n) var of value of prediction vs correct
avmiss averages number of misses
varmiss var of misses
JJ number of parameters in param space (size of j)
JJ0 number of parameters in param space (size of j0)
JJ1 number of parameters in param space (size of j1)
II number of observations (size of i)
HH number of param space for errors (size of h)
HH0 number of param space for errors (size of h)
HH1 number of param space for errors (size of h)
LL number of param space for consistency errors (size of l)
WW number of points for w index
DD number of points for d index
NN total number of categories (i.e., index n)
GG number of experiments
KK
**Below params for calculating Cov and Info matrices
omega1(g,i)
IME1(g,n)
IME2(g,n)
IME3(g,n)
IME4(g,n)
IML1(g,n)
IML2(g,n)
IEE(g,n)
INFO(g,n,k,k)
INFO1(g,n,k,k)
TINFO(g,n,k,k)
avIME1(n)
avIME2(n)
avIME3(n)
avIME4(n)
avIML1(n)
avIML2(n)
xsq(i)
xsqk(k,k)
;
KK=card(k);
NN=card(n);
JJ=3;
JJ0=3;
JJ1=3;
*II=card(i);
**Adgjust for missing values etc
II = card(i);
GG=1;
HH=card(h);
HH0=3;
HH1=3;
LL=3;
WW=3;
DD=3;
a=12;
delta=0.2;
*Generating the data
*First step: generating the x's
parameter xxx(k)
rho1 determines correlation between x's
t1(i)
t2(i)
t3(i)
x(i,k) the multinomial data
xx(i,k2) the multinomial data
beta0(n,k) the correct betas
e(i,n) the errors
s(i,n) value of "continuous" LHS var
b0(i) the correct LHS variable (category)
y(i,n) the 0 1 specification
b(i) transposing b0 to creat the y's
bounde endogenous error bound
alphae(h) param space for errors
dum999(i) dummy var to identify NA data points
sumx(k) total sum of x's
avx(k) mean value of data
avomeg mean value of partition fn
omeg1a
junk1(n)
avrecp(n) mean val of recovered prob
partial(n,k) partial p-j wrt x-k
partials(k,k8) final table of partial derivatives
partial2(n,k) partial p-j wrt x-k holding EDTI fixed
partials2h(k,k8) table of partial holding EDTI fixed hispanic
partial5(n,k) partial p-j wrt x-k holding BAD fixed
partials5h(k,k8) table of partial holding BAD fixed hispanic
partials2b(k,k8) table of partial holding EDTI fixed blacks
partials5b(k,k8) table of partial holding BAD fixed black
sumpart
erstat entropy ratio stat
dumx(k) dummy for discrete partial calculations
dumb(i) count no. blacks
dumh(i) count numberHispanics
totb total number of blacks in sample
toth total number of hispanics in sample
count(i)
sumerr sum error sq for category 1
rand1(i) random number
dumd(i)
totd total number of denied per sample
dumbd(i)
dumhd(i)
totbd total number of blacks denied in sample
tothd total number of hispanics denied in sample
**For predicting hispanics table
predh
ypredh(i)
ypred1h(i)
tableph(i,n,n)
tablepredh(n,n) Prediction table Hispanics
**For predicting Blacks table
predb
ypredb(i)
ypred1b(i)
tablepb(i,n,n)
tablepredb(n,n) Prediction table Blacks
dumj(i)
S1(i,k) for calculating probabilities
partiali(i,n,k) partial for each obs
countk
countki(i)
maxbeta max value of beta
minbeta min value of beta
bound(k,d1) table of lower and upper bounds for partials
lowbound(k) lower bound partials
upbound(k) upper bound partials
pxi(i)
;
rho1 = .85;
$ontext
table beta0(n,k) the correct betas
**Just for ommitting var exp
1 2 3
1 0 0 0
2 -1 1 2
3 1 -1 -2
;
$offtext
beta0(n,k) = 0;
options decimals=7;
parameter alphae1(h) param space for beta0 errors /
**1 0.0
**2 0.0
**3 0.0
1 -1000.0
2 0.0
3 1000.0
/;
*Uniform priors
*qb(n,k,j)=1/JJ;
*qe(i,h)=1/HH;
variables
e1(i,n) error coefficients
lm(n,k) Lagrange multipliers
;
positive variables
p(i,n) probabilities of i choosing j
pe(i,n,h) probabilities for errors
;
variables
obj objective function
;
equations
objective objective function
objdual dual objective function
objdualm dual objective function pure ME Logit
objdualw dual with weights
normalize(n,k) Normalizing first category
addp(i) adding up for p(in)
adde(i,n) adding up for errors
e1k0(i,n) e1k coefficients
e1k1(i,n) e1k coefficients
e1k2(i,n) e1k coefficients
consist(n,k) consistency for n and k
testx(n) Equation for testing model without some variable
;
**********
*The ME objective for Soofgi with BANG
objective.. obj =e= -sum(i,sum(n,p(i,n)*log((p(i,n)+1.e-8)))) -
sum(i,sum(n,sum(h,pe(i,n,h)*log((pe(i,n,h)+1.e-8)))));
**********
**********
*The ME objective for Soofgi with BANG and weight factor for objective
$ontext
objective.. obj =e=(1-delta)*(-sum(i,sum(n,p(i,n)*log((p(i,n)+1.e-4)))))
-delta*(sum(i,sum(n,sum(h,pe(i,n,h)*log((pe(i,n,h)+1.e-4))))));
$offtext
**********
*********
*Dual Objective With Noise
objdual.. obj =e= -sum(i$(dum999(i) ne 1),sum(n,sum(k,y(i,n)*x(i,k)*lm(n,k))))
+ sum(i$(dum999(i) ne 1),log(sum(n,exp(sum(k,x(i,k)*lm(n,k))))))
+ sum(i$(dum999(i) ne 1),sum(n,log(sum(h,exp(sum(k,x(i,k)*alphae(h)*lm(n,k)))))));
*********
*********
**From Fargo9 prog
*Dual Objective With Noise
objdualw.. obj =e= - sum(i$(dum999(i) ne 1),sum(n,sum(k,y(i,n)*x(i,k)*lm(n,k))))
+ (1-delta)*(sum(i$(dum999(i) ne 1),log(sum(n,exp((1/(1-delta))*sum(k,x(i,k)*lm(n,k)))))))
+ delta*(sum(i$(dum999(i) ne 1),sum(n,log(sum(h,exp((1/delta)*sum(k,x(i,k)*alphae(h)*lm(n,k))))))));
*********
*********
*Dual Objective PURE ME LOGIT
objdualm.. obj =e= -sum(i$(dum999(i) ne 1),sum(n,sum(k,y(i,n)*x(i,k)*lm(n,k))))
+ sum(i$(dum999(i) ne 1),log(sum(n,exp(sum(k,x(i,k)*lm(n,k))))))
;
*********
*********
*Dual Objective (PURE)
*objdual.. obj =e= - sum(i,sum(n,sum(k,y(i,n)*x(i,k)*lm(n,k))))
* + sum(i,log(sum(n,exp(sum(k,x(i,k)*lm(n,k))))));
*********
normalize(n,k).. lm("1",k) =e=0;
testx(n).. lm(n,"9") =e= 0;
addp(i).. sum(n,p(i,n)) =e= 1;
adde(i,n).. sum(h,pe(i,n,h)) =e= 1;
e1k0(i,"1").. sum(h,alphae1(h)*pe(i,"1",h)) =e= e1(i,"1");
e1k1(i,"2").. sum(h,alphae(h)*pe(i,"2",h)) =e= e1(i,"2");
**e1k2(i,"3").. sum(h,alphae(h)*pe(i,"3",h)) =e= e1(i,"3");
*The Soofgi with BANG case
consist(n,k).. sum(i,p(i,n)*x(i,k)) + sum(i,e1(i,n)*x(i,k)) =e= sum(i,y(i,n)*x(i,k));
**Note "l" means less or equal
**Note "g" means greater or equal
*option nlp=conopt ;
*Note: remove zero when alpha space is equal for all
model base test /objective
addp
e1k0
e1k1
** e1k2
adde
consist
/;
model dual /objdual
normalize
* testx
/;
model dualw dual with weights /objdualw
normalize
* testx
/;
model dualme /objdualm
normalize
* testx
/;
*model ols leas squares /obj1/;
* conopt and minos are two non-linear optimization packages
* minos is used by default.
options limrow=0,limcol=0;
options domlim=50000;
options iterlim=8000;
options reslim=150000;
*options work=1500000;
* conopt optimizer option is below
options solprint=off;
options nlp=conopt;
*rho1 = .85;
*t1(i) = normal(0,1);
*t2(i) = normal(0,1);
*t3(i) = normal(0,1);
$ontext
table xx(i,k2) koop-poirier data
1 2 3 4
1 1 3.000000 11.00000 1.000000
;
$offtext
**Read two basic bank1 data
*************************************************
**** data input and reassignment
*************************************************
$inlinecom { }
*$include "table1.prn" {import 350 obs sample data}
**File with loan number
$include "table1n.prn" {import 350 obs sample data}
$include "table1na.prn" {import 350 obs sample data}
$include "table1b.prn" {import 350 obs sample data}
**Definning scaling factors
**Note, it can be done automatic
parameter scale(k)
;
scale(k) = 1;
scale("4") = 1000;
scale("6") = 10000;
**For the small data set of 22 obs Bank 1 case
*scale(k) = scale(k)*10000;
**The Basic Data
x(i,"1") = 1;
$ontext
x(i,"0") = 1;
***Income case
*x(i,"1")$(basic1(i,"inc") ne 999) = basic1(i,"inc")/scale("1");
*x(i,"1")$(basic1(i,"inc") eq 999) = 999;
**Female var
x(i,"1")$(basic1(i,"f") ne 999) = basic1(i,"f");
x(i,"1")$(basic1(i,"f") eq 999) = 999;
$offtext
*x(i,"2") = basic(i,"edti");
**Building var EXCESDTI to be above 42 percent
*x(i,"2")$((basic(i,"dti") ge 42) and (basic(i,"dti") ne 999) ) =1;
*x(i,"2")$((basic(i,"dti") lt 42) and (basic(i,"dti") ne 999)) =0;
x(i,"2")$((basic(i,"dti") ge 41) and (basic(i,"dti") ne 999) ) =1;
x(i,"2")$((basic(i,"dti") lt 41) and (basic(i,"dti") ne 999)) =0;
x(i,"2")$(basic(i,"dti") eq 999) = 999;
**Continuous case dti
*x(i,"2") = basic(i,"dti") ;
x(i,"3") = basic(i,"ltv");
x(i,"4")$(basic(i,"fscore") ne 999) = basic(i,"fscore")/scale("4");
x(i,"4")$(basic(i,"fscore") eq 999) = 999;
x(i,"5") = basic(i,"bad");
x(i,"6")$(basic(i,"gfund") ne 999) = basic(i,"gfund")/scale("6");
x(i,"6")$(basic(i,"gfund") eq 999) = 999;
x(i,"7") = basic(i,"pub");
x(i,"8") = basic(i,"exp");
x(i,"9") = basic(i,"hisp");
x(i,"10") = basic(i,"black");
***Creating dummy var for identification of missing or NA values in data
$ontext
dum999(i)$((x(i,"2") eq 999) or (x(i,"3") eq 999) or (x(i,"4") eq 999)
or (x(i,"5") eq 999) or (x(i,"6") eq 999) or (x(i,"7") eq 999)
or (x(i,"8") eq 999) ) = 1;
display dum999;
$offtext
* deny1 deny2 deny3 cd inc l999 appval M F lsize
dum999(i)$((x(i,"2") eq 999) or (x(i,"3") eq 999) or (x(i,"4") eq 999)
or (x(i,"5") eq 999) or (x(i,"6") eq 999) or (x(i,"7") eq 999)
or (x(i,"8") eq 999)
or (basic1(i,"deny1") eq 7) or (basic1(i,"deny1") eq 9)
or (basic1(i,"deny2") eq 7) or (basic1(i,"deny2") eq 9)
or (basic1(i,"deny3") eq 7) or (basic1(i,"deny3") eq 9)
or (basic1(i,"lsize") gt basic1(i,"appval"))
or ((basic(i,"denied") eq 0) and (basic2(i,"INSUFND") eq 1))
**All are zeroes
* or (basic2(i,"REFIDUM") eq 1)
**Making only Hispanic subsample
* or (basic(i,"black") eq 1)
**Making only Hispanic subsample with Hisp coapplicants
* or (basic(i,"black") eq 1)
* or (basic2(i,"COAPPL") eq 0)
**Making only Black subsample
* or (basic(i,"hisp") eq 1)
***Only Female subsample
* or (basic1(i,"F") eq 1)
) = 1;
*****
dumj(i)=0;
dumj(i)$((basic(i,"loannumb") eq 1000036495)
or (basic(i,"loannumb") eq 0372764027)
or (basic(i,"loannumb") eq 0215432760)
or (basic(i,"loannumb") eq 0215687611)
or (basic(i,"loannumb") eq 0215771480)
or (basic(i,"loannumb") eq 0372760445)
or (basic(i,"loannumb") eq 0215686118)
or (basic(i,"loannumb") eq 0215788146)
or (basic(i,"loannumb") eq 0215786744)
or (basic(i,"loannumb") eq 0215562228)
or (basic(i,"loannumb") eq 0215721527)
or (basic(i,"loannumb") eq 0372320254)
or (basic(i,"loannumb") eq 0215731872)
or (basic(i,"loannumb") eq 0215700570)
or (basic(i,"loannumb") eq 0215769757)
or (basic(i,"loannumb") eq 0372304539)
or (basic(i,"loannumb") eq 0215320742)
or (basic(i,"loannumb") eq 0215622766)
or (basic(i,"loannumb") eq 0215512777)
or (basic(i,"loannumb") eq 0215633128)
or (basic(i,"loannumb") eq 0215576178)
or (basic(i,"loannumb") eq 0215695168)
or (basic(i,"loannumb") eq 0372360010)
or (basic(i,"loannumb") eq 0372132277)
or (basic(i,"loannumb") eq 0372757775)
or (basic(i,"loannumb") eq 0215627179)
****Additional delitions specprog 1 and nonverif 1
or (basic(i,"loannumb") eq 1000100127)
or (basic(i,"loannumb") eq 0372763484)
or (basic(i,"loannumb") eq 1000081848)
or (basic(i,"loannumb") eq 1000094783)
or (basic(i,"loannumb") eq 1000053698)
or (basic(i,"loannumb") eq 0371712795)
or (basic(i,"loannumb") eq 0372123642)
or (basic(i,"loannumb") eq 1000083121)
or (basic(i,"loannumb") eq 1000049547)
or (basic(i,"loannumb") eq 0215717780)
or (basic(i,"loannumb") eq 0372373861)
or (basic(i,"loannumb") eq 1000079562)
or (basic(i,"loannumb") eq 0215723879)
or (basic(i,"loannumb") eq 0372758114)
or (basic(i,"loannumb") eq 0215536685)
or (basic(i,"loannumb") eq 0372215699)
or (basic(i,"loannumb") eq 0215629464)
or (basic(i,"loannumb") eq 0215554753)
or (basic(i,"loannumb") eq 0215758032)
or (basic(i,"loannumb") eq 0215701719)
or (basic(i,"loannumb") eq 0215614045)
or (basic(i,"loannumb") eq 0215771522)
or (basic(i,"loannumb") eq 1000045619)
or (basic(i,"loannumb") eq 0215620588)
)=1;
display dumj;
****dum999(i)$( (dum999(i) ne 1) and (dumj(i) eq 1)) = 1;
*dum999(i)$( (dum999(i) ne 1) and (dumj(i) ne 1)) = 1;
*dum999(i)$( (dum999(i) ne 1) and(((basic(i,"hisp") eq 1) and (dumj(i) ne 1))
* or ((basic(i,"black") eq 1) and (dumj(i) ne 1)))) = 1;
*****
display dum999;
II = card(i) - (sum(i,dum999(i)));
display II;
***Random sampling to choose half the sample
*rand1(i) = uniform(0,1);
*dum999(i)$((dum999(i) ne 1) and (rand1(i) le 0.5)) = 1;
*dum999(i)$((dum999(i) ne 1) and (rand1(i) gt 0.5)) = 1;
II = card(i) - (sum(i,dum999(i)));
display II;
dumb(i)$(dum999(i) ne 1) = x(i,"10");
dumh(i)$(dum999(i) ne 1) = x(i,"9");
totb = sum(i,dumb(i));
toth = sum(i,dumh(i));
dumd(i)$((dum999(i) ne 1) and (basic(i,"denied") eq 1)) = 1;
totd = sum(i,dumd(i));
dumbd(i)$((dumb(i) eq 1) and (basic(i,"denied") eq 1)) = 1;
totbd = sum(i,dumbd(i));
dumhd(i)$((dumh(i) eq 1) and (basic(i,"denied") eq 1)) = 1;
tothd = sum(i,dumhd(i));
$ontext
datab(i,k5) = basic1(i,k5);
datab(i,k5)$(basic1(i,k5) eq 0) = -9;
display datab;
$offtext
display totb, toth, totd, totbd, tothd;
dumd(i) =0;
$ontext
dumd(i)$((dumb(i) eq 1) and (dumh(i) eq 1)) = 1;
totd = sum(i,dumd(i));
$offtext
**Calculating errors bounds
bounde = 1/(sqrt(II));
***
*bounde = 0.15;
bounde = 0.1;
*bounde = 1/(0.5*sqrt(II));
$ontext
**giving values to alphae
1 -.1
2 -.075
3 -.05
4 -.025
5 0
6 .025
7 .05
8 .075
9 .1
$offtext
alphae(h) =0;
alphae("1") = - bounde;
alphae("2") = - (0.75*bounde);
alphae("3") = - (0.5*bounde);
alphae("4") = - (0.25*bounde);
alphae("6") = (0.25*bounde);
alphae("7") = (0.5*bounde);
alphae("8") = (0.75*bounde);
alphae("9") = bounde;
***
**For simple ME-Logit just use supports 0
*alphae(h) =0;
***
display alphae;
**The LHS var
b0(i) = basic(i,"denied");
**Note adjust the index k before experimenting with number of vars
**Omitting var 1
*x(i,"1") = xx(i,"2");
*x(i,"2") = xx(i,"3");
*x(i,"3") = xx(i,"4");
*Calculate Probabilities
*First, calculate x'beta
***s(i,n) = sum(k,x(i,k)*beta0(n,k));
*Denominator is:
***denom(i) = 1 + exp(s(i,"2")) + exp(s(i,"3"));
*probability of the categories is:
***p1(i) = 1/denom(i);
***p2(i) = exp(s(i,"2"))/denom(i);
***p3(i) = exp(s(i,"3"))/denom(i);
***limit(i) = p1(i) + p2(i);
***u(i) = uniform(0,1);
**Creating the LHS category variable b0
***b0(i)$(u(i) le p1(i)) = 0;
***b0(i)$((u(i) gt p1(i)) and (u(i) le limit(i))) = 1;
***b0(i)$(u(i) gt limit(i)) = 2;
**Starting MC loops
Loop(g,
**End of data generation
*Converting the categiries b(i) into y(i,n) with 0 1 values
Loop(i,
b(i) = b0(i);
);
b(i)$(b0(i)=0) = 8;
b(i)$(b0(i)=1) = 9;
y(i,n) = 0;
Loop(i,
y(i,n)=b(i);
);
Loop(n,
y(i,"1")$(y(i,n) EQ 8)=1;
y(i,"2")$(y(i,n) EQ 9)=1;
$ontext
y(i,"3")$(y(i,n) EQ 2)=1;
y(i,"4")$(y(i,n) EQ 3)=1;
y(i,"5")$(y(i,n) EQ 4)=1;
$offtext
y(i,n)$(y(i,n) NE 1)=0;
);
*initialization
p.l(i,n)=1/NN;
pe.l(i,n,h)=1/HH;
tablepred(n,n) =0;
lm.l(n,k) = 0;
*ME case
*solve base maximizing obj using nlp;
*Dual GME Problem
solve dual minimizing obj using nlp;
***Dual With Weights GME Problem
*solve dualw minimizing obj using nlp;
*Dual pure ME Problem
*solve dualme minimizing obj using nlp;
*CE case
*solve base minimizing obj using nlp;
*solve ols maximizing obj using nlp;
options decimals=7;
**display p.l;
***Just for Primal ME case
$ontext
beta(g,n,k) = -consist.m(n,k)/scale(k);
lm.l(n,k) = beta("1",n,k);
display beta;
display lm.l;
$offtext
**End Primal ME case
**Estimated betas for dual case with Rescaling
beta(g,n,k) = lm.l(n,k) / scale(k);
**Dual with Weights case
***beta(g,n,k) = (1/(1-delta))*lm.l(n,k) / scale(k);
***USE ALWAYS WHEN working with weights
***Just a rescaling
*lm.l(n,k) = (1/(1-delta))*lm.l(n,k);
*beta(g,n,k) = lm.l(n,k) / scale(k);
****
*beta(g,n,k) = (1/(1-delta))*consist.m(n,k);
*When normalizing Beta0=0 use the following
****beta(g,"1",k) = 0;
***Calculating P's and w's from Betas
omeg(i)$(dum999(i) ne 1) = sum(n,exp(sum(k,x(i,k)*lm.l(n,k))));
recp(i,n)$(dum999(i) ne 1) = exp(sum(k,x(i,k)*lm.l(n,k))) / omeg(i);
psi(i,n)$(dum999(i) ne 1) = sum(h,exp(sum(k,x(i,k)*lm.l(n,k)*alphae(h))));
recw(i,n,h)$(dum999(i) ne 1) = exp(sum(k,x(i,k)*lm.l(n,k)*alphae(h))) / psi(i,n);
ei0(i,n)$(dum999(i) ne 1)=sum(h,alphae(h)*recw(i,n,h));
sumerr = ((sum(i$(dum999(i) ne 1),ei0(i,"2"))) *
(sum(i$(dum999(i) ne 1),ei0(i,"2"))) )/ II;
display sumerr;
***CAlculating PARTIALS
sumx(k) = sum(i$(dum999(i) ne 1),x(i,k));
avx(k) = sum(i$(dum999(i) ne 1),x(i,k))/ II;
avomeg = sum(n,exp(sum(k,avx(k)*lm.l(n,k))));
avrecp(n) = exp(sum(k,avx(k)*lm.l(n,k))) / avomeg;
display avrecp;
****Calculating av prob differently
avrecp(n) = (sum(i$(dum999(i) ne 1),recp(i,n)) )/ II;
display avrecp;
****
**General continuous case
*omeg1(k) = sum(n,avrecp(n) * lm.l(n,k));
*partial(n,k) = avrecp(n) * (lm.l(n,k) - omeg1(k));
omeg1(k) = sum(n,avrecp(n) * beta(g,n,k));
partial(n,k) = avrecp(n) * (beta(g,n,k) - omeg1(k));
display partial;
****Calculating av partials differently
S1(i,k) = sum(n,(recp(i,n)$(dum999(i) ne 1) ) * beta(g,n,k));
partiali(i,n,k) = recp(i,n)$(dum999(i) ne 1) * (beta(g,n,k) - S1(i,k));
partial(n,k) = sum(i$(dum999(i) ne 1),partiali(i,n,k)) / II;
display partial;
****
**pxi(i) = partiali(i,"2","9")$(dum999(i) ne 1);
***display pxi;
$ontext
**Comparing calculations
omeg1a = (1 + exp(sum(k,avx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,avx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,avx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
*display partial;
$offtext
****EXP with discrete cases
partials(k,k8) = partial("2",k);
countk = sum(i$(dum999(i) ne 1),x(i,"2")$(x(i,"2") eq 1));
partials("2","0") = sum(i$((dum999(i) ne 1) and (x(i,"2") eq 0)),partiali(i,"2","2")) / (II-countk);
partials("2","1") = sum(i$((dum999(i) ne 1) and (x(i,"2") eq 1)),partiali(i,"2","2")) / countk;
countk = sum(i$(dum999(i) ne 1),x(i,"5")$(x(i,"5") eq 1));
partials("5","0") = sum(i$((dum999(i) ne 1) and (x(i,"5") eq 0)),partiali(i,"2","5")) / (II-countk);
partials("5","1") = sum(i$((dum999(i) ne 1) and (x(i,"5") eq 1)),partiali(i,"2","5")) / countk;
countk = sum(i$(dum999(i) ne 1),x(i,"7")$(x(i,"7") eq 1));
partials("7","0") = sum(i$((dum999(i) ne 1) and (x(i,"7") eq 0)),partiali(i,"2","7")) / (II-countk);
partials("7","1") = sum(i$((dum999(i) ne 1) and (x(i,"7") eq 1)),partiali(i,"2","7")) / countk;
countk = sum(i$(dum999(i) ne 1),x(i,"8")$(x(i,"8") eq 1));
partials("8","0") = sum(i$((dum999(i) ne 1) and (x(i,"8") eq 0)),partiali(i,"2","8")) / (II-countk);
partials("8","1") = sum(i$((dum999(i) ne 1) and (x(i,"8") eq 1)),partiali(i,"2","8")) / countk;
***HISP
partials("9","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)),
partiali(i,"2","9")) / (II-(totb+toth));
partials("9","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 1)),partiali(i,"2","9")) / toth;
***BLACK
partials("10","0") = sum(i$((dum999(i) ne 1) and (x(i,"10") eq 0) and (x(i,"9") eq 0)),
partiali(i,"2","10")) / (II-(toth+totb));
partials("10","1") = sum(i$((dum999(i) ne 1) and (x(i,"10") eq 1)),partiali(i,"2","10")) / totb;
display countk, partials;
****END EXP
****Calculating upper and lower bounds for partials
**See Econometrics Review 1998
maxbeta = max(beta("1","2","2"),beta("1","2","3"), beta("1","2","4"),
beta("1","2","5"), beta("1","2","6"), beta("1","2","7"),
beta("1","2","8"), beta("1","2","9"), beta("1","2","10"));
minbeta = min(beta("1","2","2"),beta("1","2","3"), beta("1","2","4"),
beta("1","2","5"), beta("1","2","6"), beta("1","2","7"),
beta("1","2","8"), beta("1","2","9"), beta("1","2","10"));
lowbound(k) = 0.25 * (beta("1","2",k) - maxbeta);
upbound(k) = 0.25 * (beta("1","2",k) - minbeta );
bound(k,"1") = lowbound(k);
bound(k,"2") = upbound(k);
display maxbeta, minbeta, bound;
$ontext
**Discrete cases PARTIALS
partials(k,k8) = partial("2",k);
dumx(k) = avx(k);
dumx("2") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
*display partial;
partials("2","0") = partial("2","2");
dumx("2") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("2","1") = partial("2","2");
dumx(k) = avx(k);
dumx("5") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("5","0") = partial("2","5");
dumx("5") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("5","1") = partial("2","5");
dumx(k) = avx(k);
dumx("7") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("7","0") = partial("2","7");
dumx("7") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("7","1") = partial("2","7");
dumx(k) = avx(k);
dumx("8") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("8","0") = partial("2","8");
dumx("8") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("8","1") = partial("2","8");
dumx(k) = avx(k);
****NOTE holding both black and HISP at 0
dumx("9") = 0;
dumx("10") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("9","0") = partial("2","9");
***NOTE hold black at 0
dumx(k) = avx(k);
dumx("9") = 1;
dumx("10") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("9","1") = partial("2","9");
dumx(k) = avx(k);
dumx("10") = 0;
dumx("9") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("10","0") = partial("2","10");
***NOTE hold HISP at 0
dumx(k) = avx(k);
dumx("10") = 1;
dumx("9") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials("10","1") = partial("2","10");
$offtext
****
$ontext
omeg1(k) = (1 + exp(sum(k,avx(k)*lm.l("2",k)))) *
(1 + exp(sum(k,avx(k)*lm.l("2",k)))) ;
partial(k) = (lm.l("2",k) * (exp(sum(k,avx(k)*lm.l("2",k))) ))/ omeg1(k);
$offtext
**w0(i,n,h)=exp(sum(k,x(i,k)*alphae(h)*lm.l(n,k)))/psi(i,n);
**ei0(i,n)=sum(h,alphae(h)*recw(i,n,h));
***display beta, recp, ei0;
display beta, avrecp, partials, avx;
****
*Calculate additional partials for discrete cases
**Hispanic case
****EXP
**HISPANICS
*Holding EDTI fixed
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and (x(i,"2") eq 0)) =1;
countk = sum(i,countki(i));
partials2h("9","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"2") eq 0)),partiali(i,"2","9")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0) and (x(i,"2") eq 0)) =1;
countk = sum(i,countki(i));
partials2h("9","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0)
and (x(i,"2") eq 0)),partiali(i,"2","9")) / countk;
display countk, partials2h;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and
(x(i,"2") eq 1)) =1;
countk = sum(i,countki(i));
partials2h("9","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"2") eq 1)),partiali(i,"2","9")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0) and
(x(i,"2") eq 1)) =1;
countk = sum(i,countki(i));
partials2h("9","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0)
and (x(i,"2") eq 1)),partiali(i,"2","9")) / countk;
display countk, partials2h;
**End Hispanics
**BLACKS
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and (x(i,"2") eq 0)) =1;
countk = sum(i,countki(i));
partials2b("10","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"2") eq 0)),partiali(i,"2","10")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1) and (x(i,"2") eq 0)) =1;
countk = sum(i,countki(i));
partials2b("10","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1)
and (x(i,"2") eq 0)),partiali(i,"2","10")) / countk;
display countk, partials2b;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and
(x(i,"2") eq 1)) =1;
countk = sum(i,countki(i));
partials2b("10","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"2") eq 1)),partiali(i,"2","10")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1) and
(x(i,"2") eq 1)) =1;
countk = sum(i,countki(i));
partials2b("10","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1)
and (x(i,"2") eq 1)),partiali(i,"2","10")) / countk;
display countk, partials2b;
**End Blacks EDTI case
*Holding Var 5 BAD fixed
**HISPANICS
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and (x(i,"5") eq 0)) =1;
countk = sum(i,countki(i));
partials5h("9","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"5") eq 0)),partiali(i,"2","9")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0) and (x(i,"5") eq 0)) =1;
countk = sum(i,countki(i));
partials5h("9","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0)
and (x(i,"5") eq 0)),partiali(i,"2","9")) / countk;
display countk, partials5h;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and
(x(i,"5") eq 1)) =1;
countk = sum(i,countki(i));
partials5h("9","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"5") eq 1)),partiali(i,"2","9")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0) and
(x(i,"5") eq 1)) =1;
countk = sum(i,countki(i));
partials2h("9","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 1) and (x(i,"10") eq 0)
and (x(i,"5") eq 1)),partiali(i,"2","9")) / countk;
display countk, partials5h;
**End Hispanics
**BLACKS
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and (x(i,"5") eq 0)) =1;
countk = sum(i,countki(i));
partials5b("10","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"5") eq 0)),partiali(i,"2","10")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1) and (x(i,"2") eq 0)) =1;
countk = sum(i,countki(i));
partials5b("10","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1)
and (x(i,"5") eq 0)),partiali(i,"2","10")) / countk;
display countk, partials5b;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0) and
(x(i,"5") eq 1)) =1;
countk = sum(i,countki(i));
partials2b("10","0") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 0)
and (x(i,"5") eq 1)),partiali(i,"2","10")) / countk;
countki(i)=0;
countki(i)$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1) and
(x(i,"5") eq 1)) =1;
countk = sum(i,countki(i));
partials2b("10","1") = sum(i$((dum999(i) ne 1) and (x(i,"9") eq 0) and (x(i,"10") eq 1)
and (x(i,"5") eq 1)),partiali(i,"2","10")) / countk;
display countk, partials5b;
**End Blacks Var 5 BAD case
***END EXP
**Holding Var 2 edti at 0
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("2") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2h("9","0") = partial2("2","9");
dumx(k) = avx(k);
dumx("9") = 1;
dumx("10") = 0;
dumx("2") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2h("9","1") = partial2("2","9");
display partials2h;
***EDTI eq 1 case
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("2") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2h("9","0") = partial2("2","9");
dumx(k) = avx(k);
dumx("9") = 1;
dumx("10") = 0;
dumx("2") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2h("9","1") = partial2("2","9");
display partials2h;
**********************
**Holding Var 5 BAD at 0
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("5") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5h("9","0") = partial5("2","9");
dumx(k) = avx(k);
dumx("9") = 1;
dumx("10") = 0;
dumx("5") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5h("9","1") = partial5("2","9");
display partials5h;
***BAD var eq 1 case
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("5") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5h("9","0") = partial5("2","9");
dumx(k) = avx(k);
dumx("9") = 1;
dumx("10") = 0;
dumx("5") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5h("9","1") = partial5("2","9");
display partials5h;
******
**The Black case
*****
**Holding Var 2 edti at 0
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("2") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2b("10","0") = partial2("2","10");
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 1;
dumx("2") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2b("10","1") = partial2("2","10");
display partials2b;
***EDTI eq 1 case
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("2") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2b("10","0") = partial2("2","10");
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 1;
dumx("2") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial2(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials2b("10","1") = partial2("2","10");
display partials2b;
**********************
**Holding Var 5 BAD at 0
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("5") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5b("10","0") = partial5("2","10");
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 1;
dumx("5") = 0;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5b("10","1") = partial5("2","10");
display partials5b;
***BAD var eq 1 case
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 0;
dumx("5") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5b("10","0") = partial5("2","10");
dumx(k) = avx(k);
dumx("9") = 0;
dumx("10") = 1;
dumx("5") = 1;
omeg1a = (1 + exp(sum(k,dumx(k) * lm.l("2",k))) )*
(1 + exp(sum(k,dumx(k) * lm.l("2",k))) );
junk1(n) = exp(sum(k,dumx(k) *lm.l(n,k)));
partial5(n,k) = (beta(g,n,k) * junk1(n) )/ omeg1a;
partials5b("10","1") = partial5("2","10");
display partials5b;
****
*End Special partials
****
$ontext
***Calculating INF Matrix ME and ML
IME1(g,n) = II + sum(i,exp(sum(k,beta(g,n,k)*x(i,k))));
omega1(g,i) = sum(n,exp(sum(k,beta(g,n,k)*x(i,k))));
IME3(g,n) = sum(i,(exp(sum(k,beta(g,n,k)*x(i,k)))) /
(omega1(g,i)*omega1(g,i)));
*IME3(g,n) = sum(i,1/p.l(i,n));
IME2(g,n) =
sum(i,(1+exp(-sum(k,beta(g,n,k)*x(i,k))))/(exp(-sum(k,beta(g,n,k)*x(
i,k)))));
IEE(g,n) = sum(i,sum(h,1/pe.l(i,n,h)));
**Note 1 and 2 are similar
xsq(i) = sum(k,x(i,k)*x(i,k));
xsqk(k,k1) = sum(i,x(i,k)*x(i,k1));
*INFO(g,n,k,k1) = sum(i,(exp(sum(k3,beta(g,n,k3)*x(i,k3)))) /
* (omega1(g,i)*omega1(g,i)))*xsqk(k1,k);
INFO(g,n,k,k1) = sum(i,((exp(sum(k3,beta(g,n,k3)*x(i,k3))) /
(omega1(g,i)*omega1(g,i))*(x(i,k)*x(i,k1)))));
INFO1(g,n,k,k1) = sum(i,sum(h,((exp(sum(k3,beta(g,n,k3)*x(i,k3)*alphae(h))) /
(psi(i,n)*psi(i,n))*((x(i,k)*alphae(h))*(x(i,k1)*alphae(h)))))));
TINFO(g,n,k,k1) = INFO(g,n,k,k1) + INFO1(g,n,k,k1);
IML1(g,n) = sum(i,((exp(sum(k,beta(g,n,k)*x(i,k)))) /
(omega1(g,i)*omega1(g,i)))*xsq(i));
**original ML1 below
*IML1(g,n) = sum(i, ((exp(sum(k,beta(g,n,k)*x(i,k))))/
* ((1+exp(sum(k,beta(g,n,k)*
* x(i,k))))* (1+exp(sum(k,beta(g,n,k)*x(i,k)))))) * xsq(i));
IML2(g,n) = sum(i, ((exp(-sum(k,beta(g,n,k)*x(i,k))))/
((1+exp(-sum(k,beta(g,n,k
)*x(i,k))))* (1+exp(-sum(k,beta(g,n,k)*x(i,k)))))) );
***End
$offtext
**PREDICTING with the recovered Betas
xbeta(n,i)$(dum999(i) ne 1) = sum(k,x(i,k)*beta(g,n,k) * scale(k));
***Predicting according to Xbeta
$ontext
**The Five categories case
pred(i) = max(xbeta("1",i),xbeta("2",i),xbeta("3",i),xbeta("4",i),xbeta("5",i));
ypred(i)$(pred(i) eq xbeta("1",i)) = 0;
ypred(i)$(pred(i) eq xbeta("2",i)) = 1;
ypred(i)$(pred(i) eq xbeta("3",i)) = 2;
ypred(i)$(pred(i) eq xbeta("4",i)) = 3;
ypred(i)$(pred(i) eq xbeta("5",i)) = 4;
ypred1(i)=ypred(i)-b0(i);
**options decimals=7;
tablep(i,"1","1")$(ypred(i) eq 0 and b0(i) eq 0) = 1;
tablep(i,"1","2")$(ypred(i) eq 0 and b0(i) eq 1) = 1;
tablep(i,"1","3")$(ypred(i) eq 0 and b0(i) eq 2) = 1;
tablep(i,"1","4")$(ypred(i) eq 0 and b0(i) eq 3) = 1;
tablep(i,"1","5")$(ypred(i) eq 0 and b0(i) eq 4) = 1;
tablep(i,"2","1")$(ypred(i) eq 1 and b0(i) eq 0) = 1;
tablep(i,"2","2")$(ypred(i) eq 1 and b0(i) eq 1) = 1;
tablep(i,"2","3")$(ypred(i) eq 1 and b0(i) eq 2) = 1;
tablep(i,"2","4")$(ypred(i) eq 1 and b0(i) eq 3) = 1;
tablep(i,"2","5")$(ypred(i) eq 1 and b0(i) eq 4) = 1;
tablep(i,"3","1")$(ypred(i) eq 2 and b0(i) eq 0) = 1;
tablep(i,"3","2")$(ypred(i) eq 2 and b0(i) eq 1) = 1;
tablep(i,"3","3")$(ypred(i) eq 2 and b0(i) eq 2) = 1;
tablep(i,"3","4")$(ypred(i) eq 2 and b0(i) eq 3) = 1;
tablep(i,"3","5")$(ypred(i) eq 2 and b0(i) eq 4) = 1;
tablep(i,"4","1")$(ypred(i) eq 3 and b0(i) eq 0) = 1;
tablep(i,"4","2")$(ypred(i) eq 3 and b0(i) eq 1) = 1;
tablep(i,"4","3")$(ypred(i) eq 3 and b0(i) eq 2) = 1;
tablep(i,"4","4")$(ypred(i) eq 3 and b0(i) eq 3) = 1;
tablep(i,"4","5")$(ypred(i) eq 3 and b0(i) eq 4) = 1;
tablep(i,"5","1")$(ypred(i) eq 4 and b0(i) eq 0) = 1;
tablep(i,"5","2")$(ypred(i) eq 4 and b0(i) eq 1) = 1;
tablep(i,"5","3")$(ypred(i) eq 4 and b0(i) eq 2) = 1;
tablep(i,"5","4")$(ypred(i) eq 4 and b0(i) eq 3) = 1;
tablep(i,"5","5")$(ypred(i) eq 4 and b0(i) eq 4) = 1;
**End Five categories case
$offtext
**The Two categories case OCC
pred(i)$(dum999(i) ne 1) = max(xbeta("1",i),xbeta("2",i));
ypred(i)$((dum999(i) ne 1) and (pred(i) eq xbeta("1",i))) = 0;
ypred(i)$((dum999(i) ne 1) and (pred(i) eq xbeta("2",i))) = 1;
ypred1(i)$(dum999(i) ne 1)=ypred(i)-b0(i);
tablep(i,"1","1")$((dum999(i) ne 1) and (ypred(i) eq 0 and b0(i) eq 0)) = 1;
tablep(i,"1","2")$((dum999(i) ne 1) and (ypred(i) eq 0 and b0(i) eq 1)) = 1;
tablep(i,"2","1")$((dum999(i) ne 1) and (ypred(i) eq 1 and b0(i) eq 0)) = 1;
tablep(i,"2","2")$((dum999(i) ne 1) and (ypred(i) eq 1 and b0(i) eq 1)) = 1;
**End Two categories case
**display tablep;
loop(i1,
tablepred(n,i1) = sum(i,tablep(i,n,i1));
);
display tablepred;
ypred1(i)$(dum999(i) ne 1)=ypred(i)-b0(i);
ypred1(i)$((dum999(i) ne 1) and (ypred1(i) ne 0)) = 1;
sumyp(g)=sum(i,ypred1(i));
display sumyp;
dist(g,n,k) = beta(g,n,k) - (1*beta0(n,k));
MSE(g) = sum(n,sum(k,dist(g,n,k)*dist(g,n,k))) / ((NN-1)*KK);
*display dist;
*display MSE;
***Calculating Norm entropies and Other Statistics
entp(g) = -sum(i$(dum999(i) ne 1),sum(n,recp(i,n)*log(1.e-6+recp(i,n))))/(log(NN)*II);
**entpi(g,i) = -sum(n,p.l(i,n)*log(1.e-6+p.l(i,n)))/log(NN);
entpe(g) = -sum(i$(dum999(i) ne 1),sum(n,sum(h,recw(i,n,h)*log(1.e-6+recw(i,n,h)))))/(log(HH)*II*NN);
pseudor = 1-entp(g);
display entp, entpe, pseudor;
***Calculating Predicted table just for Hispanics
**The Two categories case OCC
predh(i)$((dum999(i) ne 1) and (basic(i,"hisp") eq 1)) = max(xbeta("1",i),xbeta("2",i));
ypredh(i)$((dum999(i) ne 1) and (pred(i) eq xbeta("1",i)) and (basic(i,"hisp") eq 1)) = 0;
ypredh(i)$((dum999(i) ne 1) and (pred(i) eq xbeta("2",i)) and (basic(i,"hisp") eq 1)) = 1;
ypred1h(i)$(dum999(i) ne 1)=ypredh(i)-b0(i);
tableph(i,"1","1")$((dum999(i) ne 1) and (basic(i,"hisp") eq 1) and (ypred(i) eq 0 and b0(i) eq 0)) = 1;
tableph(i,"1","2")$((dum999(i) ne 1) and (basic(i,"hisp") eq 1) and (ypred(i) eq 0 and b0(i) eq 1)) = 1;
tableph(i,"2","1")$((dum999(i) ne 1) and (basic(i,"hisp") eq 1) and (ypred(i) eq 1 and b0(i) eq 0)) = 1;
tableph(i,"2","2")$((dum999(i) ne 1) and (basic(i,"hisp") eq 1) and (ypred(i) eq 1 and b0(i) eq 1)) = 1;
loop(i1,
tablepredh(n,i1) = sum(i,tableph(i,n,i1));
);
display tablepredh;
***Calculating Predicted table just for Blacks
**The Two categories case OCC
predb(i)$((dum999(i) ne 1) and (basic(i,"black") eq 1)) = max(xbeta("1",i),xbeta("2",i));
ypredb(i)$((dum999(i) ne 1) and (pred(i) eq xbeta("1",i)) and (basic(i,"black") eq 1)) = 0;
ypredb(i)$((dum999(i) ne 1) and (pred(i) eq xbeta("2",i)) and (basic(i,"black") eq 1)) = 1;
ypred1b(i)$(dum999(i) ne 1)=ypredb(i)-b0(i);
tablepb(i,"1","1")$((dum999(i) ne 1) and (basic(i,"black") eq 1) and (ypred(i) eq 0 and b0(i) eq 0)) = 1;
tablepb(i,"1","2")$((dum999(i) ne 1) and (basic(i,"black") eq 1) and (ypred(i) eq 0 and b0(i) eq 1)) = 1;
tablepb(i,"2","1")$((dum999(i) ne 1) and (basic(i,"black") eq 1) and (ypred(i) eq 1 and b0(i) eq 0)) = 1;
tablepb(i,"2","2")$((dum999(i) ne 1) and (basic(i,"black") eq 1) and (ypred(i) eq 1 and b0(i) eq 1)) = 1;
loop(i1,
tablepredb(n,i1) = sum(i,tablepb(i,n,i1));
);
display tablepredb;
********
);
**End Monte Carlo Loop
*Loop(g,
*maxmiss$(sumyp(g) gt sumyp(g-1))=sumyp(g);
*minmiss$(sumyp(g) lt sumyp(g-1))=sumyp(g);
*);
*display sumyp;
*display maxmiss;
*display minmiss;
*minmiss = min(sumyp(g));
*maxmiss = max(sumyp(g));
display sumyp;
*display MSE;
avbeta(n,k) = sum(g,beta(g,n,k))/GG;
$ontext
avIME1(n) = sum(g,IME1(g,n))/GG;
avIME2(n) = sum(g,IME2(g,n))/GG;
avIME3(n) = sum(g,IME3(g,n))/GG;
avIML1(n) = sum(g,IML1(g,n))/GG;
avIML2(n) = sum(g,IML2(g,n))/GG;
$offtext
varbeta(n,k) = sum(g,(beta(g,n,k) - avbeta(n,k))*(beta(g,n,k) - avbeta(n,k)))/GG;
totvar = sum(n,sum(k,varbeta(n,k)));
*loop(i1,
*tablepred(g,n,i1) = sum(i,tablep(i,n,i1));
*);
*avtablep(n,n) = sum(g,tablepred(g,n,n))/GG;
*vartablep(n,n) = sum(g,(tablepred(g,n,n) - avtablep(n,n))*(tablepred(g,n,n) - avtablep(n,n)));
*minmiss = min(sumyp(g));
*maxmiss = max(sumyp(g));
avmiss = sum(g,sumyp(g))/GG;
varmiss = sum(g,(sumyp(g)-avmiss)*(sumyp(g)-avmiss)) / GG;
avmse = sum(g,mse(g))/GG;
varmse = sum(g,(mse(g)-avmse)*(mse(g)-avmse)) / GG;
aventp = sum(g,entp(g))/GG;
varentp = sum(g,(entp(g)-aventp)*(entp(g)-aventp)) / GG;
**aventpi(i) = sum(g,entpi(g,i))/GG;
**varentpi(i) = sum(g,(entpi(g,i)-aventpi(i))*(entpi(g,i)-aventpi(i))) / GG;
aventpe = sum(g,entpe(g))/GG;
varentpe = sum(g,(entpe(g)-aventpe)*(entpe(g)-aventpe)) / GG;
$ontext
display avbeta;
display varbeta;
display totvar;
display avmiss;
display IME1;
display IME2;
display IME3;
display IEE;
display IML1;
display IML2;
display xsq;
display xsqk;
display INFO;
display INFO1;
display TINFO;
display avIME1;
display avIME2;
display avIME3;
display avIML1;
display avIML2;
*display minmiss;
*display maxmiss;
display varmiss;
display alphae;
display avmse;
display varmse;
display aventp;
display varentp;
**display aventpi;
**display varentpi;
display aventpe;
display varentpe;
$offtext
*Note: below is the entropy for each beta
*ent(n,k) = -sum(j,pb.l(n,k,j)*log(1.e-7+pb.l(n,k,j)))/log(JJ);
*options decimals=7;
*display ent;
*entropy = sum(k,sum(j,pb.l(k,j)*log((1.e-7+pb.l(k,j))/1.e-7+qb(k,j))));
*entropy = -sum(n,sum(k,sum(j,pb.l(n,k,j)*log(1.e-7+pb.l(n,k,j)))));
*options decimals=7;
*display entropy;
erstat = (II*log(NN) + II*NN*log(HH) ) - obj.l;
display erstat;
xxx(k) = sqrt(sum(i$(dum999(i) ne 1),x(i,k)*x(i,k)));
display xxx;