##	PROGRAM TO IMPLEMENT RANDOM EFFECT MODEL, W/COVARIATES
##	MICHAEL BAILEY
## 	March 1, 1998

##  MANUALLY SET PARAMETERS     ##
H_100   ##  Set counter for maximum number of iterations
        ##  Also:  set covariate list (xcov) and vote list (votemat)
asig_1/4	## Manually set sigma for alpha distribution (Mode=1; Mean=1.28)
		## (not same as variance!)
CONVERGE_0.01	## Note: using new convergence - parameter by parameter
FIT_"yes"		## Fitted values for contributions:  "yes" "no"

for(y in c("1991S"))	{
#"1993S", "1991S", "1987S", "1985S", "1993H", "1991H", "1987H", "1983H"
Year_y
print(Year)

if(Year=="1993S"){
data_read.table("C:\\Mydocs\\Research\\ResearchData\\Senate\\Sen93EM.txt", header=T)
N_length(data$dem)
avex_scale( data$Ex92 + data$Ex93 + data$Ex94 )
avpic_scale( data$pImCh92 +data$pImCh93+data$pImCh94 )
avi_scale( data$Im92 + data$Im93+ data$Im94 )
avic_scale( data$ImCh92 +data$ImCh93+data$ImCh94 )
vuln_as.numeric(data$PrevMarg<.55)
Skill_scale((data$College +data$Grad)/
	(data$NoHS+ data$HS+data$SomeColl+data$AssDeg+
	data$College+data$Grad))
if(FIT=="yes"){
	FittedMNC_scale(lm(scale(data$MNC9094PCT)~Skill+ avex+avic+
		scale(data$union93)+data$dem*scale(data$union93)+ 					scale(data$Labor9094PCT)+ scale(data$Perot92)+ 
		data$dem + scale((data$CleanCC94)/2)+data$Fin+ 	data$YrElected
		+data$DemLead +data$RepLead, 	na.action=na.omit)$fitted.values)
	FittedLabor_scale(lm(scale(data$Labor9094PCT)~Skill+ avex+avic+
		scale(data$union93)+data$dem*scale(data$union93)+ 					scale(data$MNC9094PCT)	+ scale(data$Perot92)+ data$dem + 				scale((data$CleanCC94)/2) + data$Fin + data$YrElected + 
		data$DemLead +data$RepLead, na.action=na.omit)$fitted.values)}
xcov_cbind(Skill, avex, FittedMNC, avic, 
	scale(data$union93), data$dem*scale(data$union93), 
FittedLabor, scale(data$Perot92), 
	data$dem, scale((data$CleanCC94)/2), rep(1, N))
cov.names_c(" Skill", "ex", " MNC9094PCT", 
	" imch", " Union", " Dem*Union", " Labor9094PCT", " Perot",
	"dem", " CleanCCCUS",  " const")
yobs_ cbind(data$Nafta, data$NSA, data$GATTFT, data$Gatt94, 
	data$GattBW )
yobs_replace(yobs, yobs==9, NA)
vote.names_ c("Nafta", "NSA", " GattFT", " Gatt94", " GattBW" )
}

if(Year=="1991S"){	## READ IN DATA        ## 
data_read.table("C:\\Mydocs\\Research\\ResearchData\\Senate\\Sen91EM.txt", header=T)
NewCC_scale ( data$newCC91 + data$newCC92 )
avex_scale ( data$Ex90+data$Ex91+data$Ex92 )
avpic_scale ( data$pImch90+data$pImch91+data$pImch92 )
avi_scale ( data$Im90+data$Im91+data$Im92 )
avic_scale ( data$Imch90+data$Imch91+data$Imch92 )
N_length(data$dem)
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)

Skill_scale((data$College +data$Grad)/
	(data$NoHS+ data$HS+data$SomeColl+data$AssDeg+data$College+data$Grad))
if(FIT=="yes"){
	FittedMNC_scale(lm(scale(data$MNC8892PCT)~Skill+ avex+avic+
		scale(data$ManUnion)+data$dem*scale(data$ManUnion)+ 					scale(data$Labor8892PCT)+ 
	data$dem + scale(NewCC)+data$Fin91 + data$PrevElec + 
	data$DemLead +data$RepLead, 	na.action=na.omit)$fitted.values)
	FittedLabor_scale(lm(scale(data$Labor8892PCT)~Skill+ avex+avic+
		scale(data$ManUnion)+data$dem*scale(data$ManUnion)+ 					scale(data$MNC8892PCT)	+ data$dem + 							scale(NewCC) + data$Fin91 + data$PrevElec + 
		data$DemLead +data$RepLead, na.action=na.omit)$fitted.values)}
xcov_cbind(Skill, avex, FittedMNC, 	
	avic, scale(data$ManUnion), (data$dem*scale(data$ManUnion)),
	FittedLabor, 	
	data$dem, NewCC, scale(data$PrevMarg), 
	rep(1, N))
cov.names_c(" Skillwage",  "ex", " MNC8892PCT",
	" imch", "Union", " Dem*Union", " Labor8892PCT ",
	"dem", " NewCC", " PrevMargin ", 
	" const")
yobs_ cbind(data$Minivan, data$FT91)
vote.names_c(" Minvan", " FT91")
}				##	End 1991 Senate Data

if(Year=="1987S"){
data_read.table("C:\\Mydocs\\Research\\ResearchData\\Senate\\Sen87EM.txt", header=T)
CCUS _ scale( data$cleanCC87 + data$cleanCC88 )
avex_scale( data$ExB86 +  data$ExB87 + data$ExB88)
avpic_ scale( data$pImchB86 +  data$pImchB87 + data$pImchB88)
avic_ scale( data$ImchB86 +  data$ImchB87 + data$ImchB88 )
avi_ scale( data$ImB86 +  data$Im87 + data$ImB88 )
N_length(data$dem) 
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)
Skill_scale((data$College +data$Grad)/
	(data$NoHS +data$HS+data$SomeColl+data$AssDeg+
	data$College+data$Grad))
if(FIT=="yes"){
	FittedMNC_scale(lm(scale(data$MNC8488PCT)~Skill+ avex+avic+
		scale(data$union87)+data$dem*scale(data$union87)+ 					scale(data$Labor8488PCT)+ 
	data$dem + scale(NewCC)+data$Fin87 + data$PrevElec + 
	data$DemLead +data$RepLead, 	na.action=na.omit)$fitted.values)
	FittedLabor_scale(lm(scale(data$Labor8488PCT)~Skill+ avex+avic+
		scale(data$union87)+data$dem*scale(data$union87)+ 					scale(data$MNC8488PCT)	+ data$dem + 							scale(NewCC) + data$Fin87 + data$PrevElec + 
		data$DemLead +data$RepLead, na.action=na.omit)$fitted.values)}
xcov_cbind(Skill, avex, FittedMNC,
	avic, scale(data$union87), data$dem*scale(data$union87), 	
	FittedLabor,
	data$dem, CCUS, scale(data$PrevMargin), 
	rep(1, N))
cov.names_c(" Skillwage", "ex", " MNC8488PCT",
	" imch", " Union", " Dem*Union", " Labor8488PCT ", 
	"dem", " cleanCC", " PrevMargin ", 
	" const")
yobs_cbind(data$Tex2, data$Cement16, data$Oil175, data$Uran80, data$BuyAm132, 
	data$Coal138, data$OTBAuth, data$OTB301, , data$Canada, data$TexClot, 
	data$TexPOO, 
	data$TexProf, data$TexBurd, data$TexDef, data$TexCl, data$TexC, 
	data$TexFt1, data$TexCBI, data$TexFt2, data$Tex1, 
	data$OPIC343)
vote.names_c( "Tex2",  " Cement ", "Oil", " Uranium", " BuyAm", " Coal", 
	" OTBPres", " Ag301", " Canada", " TexCloture", " TexPoint", 
	" TexProf", " TexBurd", " TexDef", "TexCl", "TexC", " TexFt1", 
	" TexCBI", " TexFt2", "Tex1", " OPIC343" )
}		## End 1987 Senate Data

if(Year=="1985S"){
data_read.table("C:\\Mydocs\\Research\\ResearchData\\Senate\\Sen85EM.txt", header=T)
CCUS_ scale(data$newCC86) 
unionp_data$union*data$man86
avex_ scale( data$Ex86 )
avic_scale( data$Imch85+ data$Imch86 )
avi_scale( data$Im85+ data$Im86 )
avpic_scale( data$pImch85+ data$pImch86 )
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)
#  Can't average over longer time period because don't have data
N_length(data$dem)
Skill_scale((data$College +data$Grad)/
	(data$NoHS+ data$HS+data$SomeColl+data$AssDeg+data$College
	+data$Grad))
if(FIT=="yes"){
	FittedMNC_scale(lm(scale(data$MNC8286PCT)~Skill+ avex+avic+
		scale(unionp)+data$dem*scale(unionp)+ 					scale(data$Labor8286PCT)+ 
	data$dem + CCUS+data$Fin85 + data$PrevElec + 
	data$DemLead +data$RepLead, 	na.action=na.omit)$fitted.values)
	FittedLabor_scale(lm(scale(data$Labor8286PCT)~Skill+ avex+avic+
		scale(unionp)+data$dem*scale(unionp)+ 							scale(data$MNC8286PCT)	+ data$dem + 							CCUS + data$Fin85 + data$PrevElec + 
		data$DemLead +data$RepLead, na.action=na.omit)$fitted.values)}
xcov_cbind(Skill, avex, FittedMNC,
	avic, scale(unionp), data$dem*scale(unionp), 
	FittedLabor, data$dem, CCUS, scale(data$PrevMargin), 
	rep(1, N)) 
cov.names_c(" Skillwage", "ex", " MNC8286PCT", 
	" Imch", " union", " DemUn", 	" Labor8286PCT ", 
	"dem", "CCUS", " PrevMargin ",  
	" const")
yobs_cbind(data$ft587, data$ft200, data$ft201, data$ft258, data$ft300,
      data$ft301, data$ft302, data$ft303, data$ft305, data$ft340, data$ft553 )
vote.names_c(" ft587", " ft200", " ft201", " ft258", " ft300", " ft301", 
	" ft302", " ft303", " ft305", " ft340", " ft553")
}

if(Year=="1993H"){
data_read.table("C:\\Mydocs\\Research\\ResearchData\\House\\Hou93EM.txt", header=T)
avex_ scale( data$ex91 + data$ex92 + data$ex93 )
avi_ scale( data$im91 + data$im92 + data$im93 )
avic_ scale( data$imch91 + data$imch92 + data$imch93 )
avpic_ scale( data$pImch91 + data$pImch92 + data$pImch93 )
CC_ scale( data$cleanCC93 + data$cleanCC94 )
vuln_as.numeric(data$Marg92<.55)
N_length(data$dem)
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)
 
Skill_scale((data$College +data$Grad)/
	(data$NoHS+data$SomeHS+data$HS+data$SomeColl+data$AssDeg+data$College+data	$Grad))
Home_ scale(data$Owner/(data$Owner+data$Renter))
xcov_cbind(Skill, avex, scale(data$MNCPCT), avic, 
	scale(data$disunion), data$dem*scale(data$disunion), 	scale(data$LaborPCT), scale(data$perot), 
	data$dem, CC, scale(data$Marg92), 
	rep(1, N)) 
cov.names_c(" Skill", "ex", " MNC9094PCT", 
	" imch", " Union", " Dem*Union", " Labor9094PCT", 
	" Perot", "dem", " cleanCC", " PrevMargin ", 
	" const")

yobs_cbind(data$gattft, data$nafta, data$gatt)
vote.names_c(" gattft ", " nafta", " gatt")
}		## End 1993 House Data

if(Year=="1991H"){	##  READ IN DATA        ## 
data_read.table("C:\\Mydocs\\Research\\ResearchData\\House\\Hou91EM.txt", header=T)
avex_ scale( data$ex90 + data$ex91 + data$ex92 )
avi_ scale( data$im90 + data$im91 + data$im92 )
avic_ scale( data$imch90 + data$imch91 + data$imch92 )
avpic_ scale( data$pImch90 + data$pImch91 + data$pImch92 )
clCCUS_ scale ( data$newCC91 + data$newCC92 )
vuln_as.numeric(data$Marg90<.55)
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)
N_length(data$dem)
xcov_cbind(scale(data$CollegePlus), avex , scale(data$MNCPCT),
	avic,	scale(data$manunion), (scale(data$manunion)*data$dem), 
	scale(data$LaborPCT), 	
	data$dem, clCCUS, scale(data$Marg90), 
	rep(1, N)) 
cov.names_c(" Skill ", "ex", " MNCPCT",
	" imch", " Union", " Union*Dem", " LaborPCT",
	"Dem", " cleanCCUS", " Marg90 ",  
	" const")
yobs_cbind( data$FT115, data$FText116, data$JAuto272, data$TB273 )
vote.names_c(" FT115", " FText116", " JAuto272", " TB273" )
}			## End data for 1991 House

if(Year=="1987H"){
data_read.table("C:\\Mydocs\\Research\\ResearchData\\House\\Hou87EM.txt", header=T)
clCCUS_ scale (data$cleanCC88 + data$cleanCC87)
avex_ scale( data$exB86 + data$exB87 + data$exB88 )
avi_ scale( data$imB86 + data$imB87 + data$imB88 )
avic_ scale( data$imchB86 + data$imchB87 + data$imchB88 )
avpic_ scale( data$pImchB86 + data$pImchB87 + data$pImchB88 )
vuln_as.numeric(data$Marg86<55)
N_length(data$dem)
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)
xcov_cbind(scale(data$CollegePlus), avex , scale(data$MNCPCT),
	avic,	scale(data$union87), (scale(data$union87)*data$dem), 	scale(data$LaborPCT),
	data$dem, clCCUS,  scale(data$Marg86), 
	rep(1, N))
cov.names_c(" Skill ", "ex", " MNCPCT",
	" imch", " Union", " Union*Dem", " LaborPCT", 
	"Dem", " cleanCCUS", " Marg86 ", 
	" const")	
yobs_cbind(data$Tex319, data$NegObj69, data$Gepa72, data$Michel77, 	data$OTBGep78, data$Gepb426, data$PbWk456, data$Can267, data$Tex341, 
	data$Tex426)
vote.names_c(" NegObj69 ", " Tex319", " Gepa72", " Michel77", " OTBGep78", 
	" Gepb426", " PbWk456", " Can267", " Tex341", 
	" Tex426")
}

if(Year=="1983H"){
data_read.table("C:\\Mydocs\\Research\\ResearchData\\House\\Hou83EM.txt", header=T)
avex_ scale( data$ex82 + data$ex83 + data$ex84 )
avi_ scale( data$im82 + data$im83 + data$im84 )
avic_ scale( data$imch82 + data$imch83 + data$imch84 )
avpic_ scale( data$pImch81 + data$pImch82 + data$pImch83 + data$pImch84 )
clCCUS_ scale( data$cleanCC84 )
N_length(data$dem)
mw_as.numeric(data$state=="IL"|data$state=="IN"|
	data$state=="IA"|data$state=="KS"|data$state=="MI"|
	data$state=="MN"|data$state=="MO"|data$state=="NE"|
	data$state=="ND"|data$state=="OH"|data$state=="SD"|data$state=="WI")
south_as.numeric(data$state=="AL"|data$state=="AR"|data$state=="FL"|
	data$state=="GA"|data$state=="KY"|data$state=="LA"|
	data$state=="MS"|data$state=="NC"|data$state=="OK"|data$state=="SC"|
	data$state=="TN"|data$state=="TX"|data$state=="VA")
west_as.numeric(data$state=="AZ"|data$state=="CO"|
	data$state=="ID"|data$state=="MT"|data$state=="NV"|
	data$state=="NM"|data$state=="UT"|data$state=="WY"|data$state=="CA"|
	data$state=="OR"|data$state=="WA"|data$state=="AK"|data$state=="HI") 
# NE = CT, ME, MA, NH, RI, VT, DE, NJ, NY, PA, MD (south in census)
xcov_cbind(scale(data$CollegePlus), avex , scale(data$MNCPCT),
	avic, scale(data$union84M), (scale(data$union84M)*data$dem), 	scale(data$LaborPCT),
	data$dem,clCCUS, scale(data$Marg82), 
	rep(1, N)) 
cov.names_c(" Skill ",  "ex", " MNCPCT ",
	" imch", " Union", " Union*Dem", " LaborPCT",
	"Dem", " CleanCCUS",  " PrevMargin ", 
" const") 
yobs_cbind(data$Auto417, data$Auto415, data$Auto416, data$Unfair781, 
	data$Unfair782, data$DelayLab824, data$Steel874, data$Wine875, 	data$Tar214, data$AsiaGSP876)
vote.names_c(" Auto417"," Auto415", " Auto416", " Unfair781", 
	" Unfair782", " DelayLab824", " Steel874", " Wine875", 
	" Tar214"," AsiaGSP876 ")
}



P_length(xcov[1,])           ##  Number of covariates
if (length(b)!= P) {b_rep(1, P)}
TT_length(yobs[1,])          	##  Number of votes
if(length(a)!=TT) a_rep(1, TT)
if(length(ak)!=TT) ak_rep(0, TT)
k_ak/a
xi_0
ximinus1_0

##  INITIALIZE VARIABLES AND FUNCTIONS  ##
##  Generated data
thetaq_c(-11:11)/2
	##	NOTE: CHANGING interval between THETAQ MEANS WEIGHTS MUST CHANGE!
Q_length(thetaq)        		##  Number of quadrature nodes 
w_c(1, rep(c(4, 2), (Q-3)/2), 4, 1)/6	## Based on Simpson's Rule
	##	NOTE: CHANGING interval between THETAQ MEANS WEIGHTS MUST CHANGE!

onesTT_rep(1, TT)
onesQ_rep(1, Q)
onesN_rep(1, N)
onesP_rep(1, P)

## Splus minimizes maxfun; so set to negative to get max
maxfun_function(xi, q, r, f, w){
	log(xi[1]) + (.5/asig^2)*(log(xi[1])-asig^2)^2 +
	sum(w*(
	r*log(1 + exp(xi[2] - xi[1]*q)) + 
	(f-r)*log(1+exp(xi[1]*q - xi[2]))		)	)
	}

## Mode of a lognormal r.v. is exp(mu  - sigma^2) 
## 	--> for mode = 1 ==>, mu = sigma^2
##	Alpha is dist'ed LN(1, var) where mu = -.5, sigma = asig (above)

Ptq_matrix(NA, TT, Q)
PTtheta_array(NA, dim=c(N, TT, Q))      ## Prob(y_it|q, a, k)
Liq_matrix(NA, N, Q)                    ##  L_i(X_q) - prob Y_i|theta=theta_q
ritq_array(0, dim=c(N, TT, Q))          ##  yti * L*W/sum of L*W over q
fitq_array(0, dim=c(N, TT, Q))          ##  P(i is thetaq) for yobs!=NA
rm(thetai)           ##  Del. thetai with names attached

##	ITERATIONS              ##
for(h in 1:H){
print ("Iteration Number ");    print(h)

##	EXPECTATION STEP        ##
for(q in 1:Q) {                 ##  Start first q loop
Ptq[,q]_(1+exp(-a*thetaq[q] +ak))^(-1)     ## TTx1 matrix for each q
PTtheta[,,q]_(outer(rep(1, N), Ptq[,q])^yobs)*(
		outer(rep(1, N), (1-Ptq[,q]))^(1-yobs))                                    
      		##  Prob y_it given q and vt parameters
PTtheta[,,q]_replace(PTtheta[,,q], PTtheta[,,q] =="NA", 1) 	
			## Cum product, so ok to replace with 1's
Liq[, q]_exp(log(PTtheta[,,q])%*%rep(1, TT))  ##  prob of i's t votes given q
}                               ##  End  first q loop
ProbThetaq_(1/sqrt(2*pi))*exp((-.5)*(outer(rep(1, N), thetaq)-
	crossprod(t(xcov), b)%*%rep(1, Q))^2)
                 	##  Prob Thetaq given X'beta (numerator (n x Q)

Ptheta_(Liq*ProbThetaq)/(((Liq*ProbThetaq)%*%onesQ)%*%t(rep(1, Q)))
	## Because bin size cancels out, do not use w
	## (not having bin size increases flexibility)

for (t in 1:TT){                ##  Start first t loop
ritq[,t,]_outer(yobs[,t],rep(1, Q))*Ptheta  ##  prob(i is q)*y_it
fitq[,t,]_Ptheta*((1+outer(yobs[,t],rep(1, Q)))/(1+outer(yobs[,t],rep(1, Q))))
                                ##  jerry-rigged way to del. missing obs from
                                ##  fitq matrix (which is length(beta)(i is 						  ##	q)|y_it!=NA)
}                               ##  End  first t loop
ritq_replace(ritq, ritq=="NA", 0)       ##  Am summing, so ok to replace w/ 
fitq_replace(fitq, fitq=="NA", 0)       ##  zeroes
rtq_apply(ritq, c(2, 3), sum)   ##  Sum y_it*E(type q) over i for each t
ftq_apply(fitq, c(2, 3), sum)   ##  Sum of fiqt over i's for each t & q


thetai_Ptheta%*%(thetaq)		## E(theta_i)
thetai_replace(thetai, Liq[,1]==1, NA)    ##  Liq[,1]=1 is for
gthetai_thetai[!is.na(thetai)]
gid_data$id[1:N][!is.na(thetai)]
gxcov_matrix(xcov[!is.na(thetai)], ,length(b))
b_solve(crossprod(gxcov)) %*%crossprod(gxcov, gthetai)
print(b)

##  MAXIMIZATION STEP       ##
ximinus1_xi 

##	IDENTIFY WITH A1_1, K1_0
a[1]_1
ak[1]_0

for(t in 2:TT){
max_nlminb(start=c(a[t], ak[t]), objective = maxfun, lower=c(0, -Inf), 
	r=rtq[t,], q=thetaq, f=ftq[t,], w=w)
ak[t]_max$parameters[2]
a[t]_max$parameters[1]
}                               ## End  Max step 

k_ak/a
print(a)
print(k)
xi_c(a, k, b)   ##  For convergence test


##      Convergence Criterion -- only set final estimates upon convergence
print(" convergence criterion ");	print(max(abs((xi-ximinus1))))

##  BEGIN conversion loop
if((max(abs((xi-ximinus1))))<CONVERGE) {
print("Convergence at"); print(h)

##  STANDARD ERRORS OF ESTIMATES        ##
dlogLda_matrix(NA, N, TT)
dlogLdk_matrix(NA, N, TT)
dlogLdb_matrix(NA, N, P)

## Create terms for Standard Errors of Ideal Point Estimates (delta method)
	# GAMMA is first derivative matrix of E(theta)
GAMMA.alpha_matrix(NA, nrow=N, ncol=TT)
GAMMA.kappa_matrix(NA, nrow=N, ncol=TT)
GAMMA.beta_matrix(NA, nrow=N, ncol=P)

## N x 1 vector of \int P(Yi|.)P(theta|.)dtheta
InfoIntegral_((Liq*ProbThetaq)%*%w)

## Fill up individual dlogLikelihood/dparameters and GAMMA matrix
for(tt in 2:TT) {
thetakappa_rep(1, N)%*%t(k[tt]- thetaq)
voteprob_yobs[,tt]%*%t(rep(1, Q))- rep(1, N)%*%t(Ptq[tt,])

dlogLda[,tt]_ (log(a[tt])-2*asig)/(a[tt]*asig)+
		(1/InfoIntegral)*(
		ProbThetaq*Liq*thetakappa*voteprob)%*%w
dlogLdk[,tt]_(1/InfoIntegral)*(ProbThetaq*Liq*a[tt]*voteprob)%*%w
GAMMA.alpha[,tt]_ (1/InfoIntegral)*((Liq*ProbThetaq*thetakappa*voteprob*0.5)%*%
	(w*thetaq)) - 
	(1/InfoIntegral^2)* ((Liq*ProbThetaq*thetakappa*voteprob)%*%w)*
	((Liq*ProbThetaq*0.5)%*%(w*thetaq))
GAMMA.kappa[,tt]_(1/InfoIntegral)*((Liq*ProbThetaq*a[tt]*voteprob*0.5)%*%
	(w*thetaq)) - 
	(1/InfoIntegral^2)* ((Liq*ProbThetaq* a[tt]*voteprob)%*%w)*
	((Liq*ProbThetaq*0.5)%*%(w*thetaq))
}                       #       End tt loop

## fill up dlogL/dbetas and GAMMA.beta
for(p in 1:P){	
dlogLdb[,p]_ (1/InfoIntegral)*(
	Liq*ProbThetaq*(
	(rep(1, N)%*%t(thetaq))-(xcov%*%b)%*%t(rep(1, Q)))*
	(xcov[,p]%*%t(rep(1, Q)))
	)%*%w
	## See notes
GAMMA.beta[,p]_(1/InfoIntegral)*
	((Liq*ProbThetaq*0.5*
	((rep(1, N)%*%t(thetaq))-(xcov%*%b)%*%t(rep(1, Q)))*
	(xcov[,p]%*%t(rep(1, Q))))%*%(w*thetaq)) 
	- (1/InfoIntegral)*(dlogLdb[,p])*
	((Liq*ProbThetaq*0.5)%*%(w*thetaq))
## Last part originally (1/InfoIntegral^2) but cancels with InfoIntegral 
## that we have to multiply by dlogLdb
}				#	End of p loop

dlogLda_replace(dlogLda, dlogLda=="NA", 0)
dlogLdk_replace(dlogLdk, dlogLdk=="NA", 0)
dlogLdb_replace(dlogLdb, dlogLdb=="NA", 0)

#  Create first derivative matrix (N x (2*TT -2 +P))
dlogLdxi_cbind(dlogLda[,2:TT], dlogLdk[,2:TT], dlogLdb)   

infoa_matrix(0, 2*TT-2+P, 2*TT -2 + P)

for(i in 1:N){          #       Start i loop
infoa_infoa + outer(dlogLdxi[i,], dlogLdxi[i,])
                        # Cumulative sum
}                       #       End i loop
xihat_c(a[2:TT], k[2:TT], b)
xi.se_sqrt(diag(solve(infoa)))
xi.t_xihat/xi.se

#  Create GAMMA matrix (for delta method) - is really GAMMA' (derivative)
GAMMA_cbind(GAMMA.alpha[,2:TT], GAMMA.kappa[,2:TT], GAMMA.beta)
GAMMA_replace(GAMMA, GAMMA=="NA", 0)
ETheta.se_sqrt(diag(GAMMA%*%(solve(infoa))%*%(t(GAMMA))))


## Calculate Percent Correctly Predicted
PITX_1/(1+exp((-1)*onesN%*%t(a)*((xcov%*%b[,1])%*%onesTT - onesN%*%t(k))))
temp_matrix(cut(PITX, breaks = c(0, .5, 1)), N, TT)
temp_replace(temp, temp==1, 0)
temp_replace(temp, temp==2, 1)
correct_(1-(yobs-temp)^2)
CorPercentX_mean(correct, na.rm=T)

PITXY_1/(1+exp((-1)*onesN%*%t(a)*(outer(thetai[,1],onesTT) - onesN%*%t(k))))
temp_matrix(cut(PITXY, breaks = c(0, .5, 1)), N, TT)
temp_replace(temp, temp==1, 0)
temp_replace(temp, temp==2, 1)
correct_(1-(yobs-temp)^2)
CorPercentXY_mean(correct, na.rm=T)

##  DISPLAY MODULE
print(Year)
print(" Convergence criteria:  ");	print(CONVERGE)

atemp_a[2:TT]
row.names(atemp)_vote.names[2:TT]
temp2_cbind(atemp, xi.t[1:(TT-1)])
names(temp2)_list("Coef.", "t-stat")
print("Discrimination Parameter Estimates");	print(temp2)

ktemp_k[2:TT]
row.names(ktemp)_vote.names[2:TT]
temp3_cbind(ktemp, xi.t[(TT):(2*TT-2)])
names(temp3)_list("Coef.", "t-stat")
print("Cutpoint Estimates");	print(temp3)

print("Number of observations ");       print(length(gid))
print ("Percent Correctly Predicted|X "); print(CorPercentX)
print ("Percent Correctly Predicted|XY "); print(CorPercentXY)
print(" asig "); print(asig)

btemp_b
row.names(btemp)_cov.names
temp1_cbind(btemp, xi.t[(2*TT-1):(2*TT-2+P)])
names(temp1)_list("Coef.", "t-stat")
print("Effect Parameter Estimates");	print(temp1)

break}                          ##  End conversion loop
}                               ##  End  h loop

}