@ Loading the data set @
load vp4 presinc dsouth64 south4 south64b home4x triheat zecon41 adaaca 
neweng64 dev3 dev2 west76 ncent72 legis julyec2 a elecvote jimcrow govrace 
pcin pcapreal educat period rain outstate votpop DV DP TieBreak pcblck;
votpop = 0.1*votpop;

@ Loading the Campbell variables and preparing their scaling @
campvar=triheat~julyec2~presinc~vp4~home4x~adaaca~dev3~dev2~legis~zecon41~dsouth64~south4~south64b~west76~ncent72~neweng64;
campdiv = {100 10 10 10 10 1000 100 100 1000 100 1 10 10 10 10 10};
campdiv=campdiv';

@ Marking the number of observations @
N = rows(DP);

@ Parameters names @
Parnames ={"sigma_d ","beta_0  ","Gallup  ","GNP Grow","Incumben","VP H S  ","Pres H S",
                     "ADA-ACA ","PrevVote","Prv8Vote","StateLeg","StateGrw","South64 ","SouthDem",
                     "South 64","West7680","NorthCen","NewEng  ","sigma_p ","betad_0 ","betar_0 ",
                     "beta_gov","beta_pre","beta_jim","beta_inc","beta_blk","beta_mov","beta_edu",
                     "dummy88 ","dummy84 ","dummy80 ","dummy76 ","dummy72 ","dummy68 ",
                     "dummy64 ","dummy60 ","dummy56 ","dummy52 ","ExtState","rho     ","eta     ","S       "};


@ Max Likelihood Procedure   @
library maxlik;  #include maxlik.ext;  maxset;
Active = ones(42,1);
_max_Active = Active;
_max_GradTol = 0.000001;
_max_ParNames = Parnames;


Proc Lik(q, DP);

local sed, d0, cm,sep, betad0, betar0, bgov, bpre, bjim, bvinc, bvblk, bin, 
bedu, tm, outst, rho, eta, S, dbar, gama, alpha, d_star, eps_star, f1, f2, 
Effort, xb, prob, L;

@  Parameters  @
sed = .1*q[1];
d0  = q[2];
cm  = q[3:18]./campdiv;
sep = q[19]; 
betad0 = q[20];
betar0 = q[21];
bgov = .1*q[22];
bpre = .1*q[23];
bjim = q[24];
bvinc = q[25];
bvblk = q[26];
bin = q[27]; 
bedu = q[28];
tm = q[29:38]; 
outst = q[39];

rho = q[40];
eta = q[41];
S=q[42];

dbar = d0 + campvar*cm;
gama = exp(betad0 - betar0);  
alpha = 1/(1+gama);
d_star = 1/(1+(1/DV-1)*gama); 

dbar = dbar.*(1-a) + d_star.*a; @ a=1 for the observations with missing campvar data @

Effort = (TieBreak+S).*rho*alpha*(1-alpha).*pdfn((alpha-dbar)/sed)/sed  - eta*votpop;

Xb =  betad0 + bgov*govrace + bpre*rain + bvinc*pcapreal + bedu*educat + 
bjim*jimcrow + bvblk*pcblck + bin*pcin + period*tm + OutState*outst;


f1=(1/sed)*pdfn( (d_star-dbar)/sed ).*((d_star./DV).^2)*gama;

eps_star = ln( DP./d_star ) - rho*Effort - Xb;

f2=(1/sep)*pdfn( eps_star/sep ).*(1/DP);

prob = (f1.*(1-a)+a).*f2;

L = ln(prob);

retp(L); endp;

load start;

output file = Table1.out reset;

{ w,f,foc,cov,s} = maxlik(DP,0,&lik,start);
{ w,f,foc,cov,s } = MAXPrt(w,f,foc,cov,s);

Format 12,12;
"The value of the likelihood is ";;n*f;
output off;



