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;