[R] Help: Maximum likelihood estimation

[roach]

http://r.789695.n4.nabble.com/file/n3018477/JEC.dta JEC.dta 
This is the data.
http://r.789695.n4.nabble.com/file/n3018477/A_Study_of_Cartel_Stability_The_Joint_Executive_Committee.pdf
A_Study_of_Cartel_Stability_The_Joint_Executive_Committee.pdf 
This is the original paper.
The variable S is breakdown into four dummy variables dm1, dm2, dm3, dm4.
the following is my code.


library(foreign)
jec=read.dta("d:/JEC.dta")

## generate year variables
jec["year"]=1879+ceiling(jec$week/52)

## generate months variables
jec["month"]=ifelse(jec$week%%52==0,13,ceiling((jec$week%%52)/4))

## generate dummy variables
jec["dm1"]=0
jec[jec$week>=28&jec$week<=subset(jec,year==1883&week%%52==10)$week,"dm1"]=1
jec["dm2"]=0
jec[jec$year==1883&jec$week%%52>=11&jec$week%%52<=25,"dm2"]=1
jec["dm3"]=0
jec[jec$week>=subset(jec,year==1883&week%%52==26)$week&jec$week<=subset(jec,year==1886&week%%52==11)$week,"dm3"]=1
jec["dm4"]=0
jec[jec$week>=subset(jec,year==1886&week%%52==12)$week,"dm4"]=1

### E-M algorithm
library("alabama")
n=328
w=seq(0.05,0.9,length.out=n)
repeat{
logp=log(jec$gr)
logq=log(jec$tqg)
lakes=jec$lakes
dm1=jec$dm1
dm2=jec$dm2
dm3=jec$dm3
dm4=jec$dm4
lamda=mean(w)
## -log likelihood function
log.L=function(parm){
alpha0=parm[1]
alpha1=parm[2]
alpha2=parm[3]
beta0=parm[4]
beta1=parm[5]
beta21=parm[6]
beta22=parm[7]
beta23=parm[8]
beta24=parm[9]
beta3=parm[10]
sigma11=parm[11]
sigma12=parm[12]
sigma21=parm[13]
sigma22=parm[14]

u1=-alpha0-alpha1*logp-alpha2*lakes+logq
u21=-beta0-beta1*logq-beta21*dm1-beta22*dm2-beta23*dm3-beta24*dm4-beta3+logp
u22=-beta0-beta1*logq-beta21*dm1-beta22*dm2-beta23*dm3-beta24*dm4+logp
expon1=u1^2*sigma11+u1*u21*sigma12+u1*u21*sigma21+u21^2*sigma22
expon2=u1^2*sigma11+u1*u22*sigma12+u1*u22*sigma21+u22^2*sigma22
const=-log(2*pi)+.5*log(sigma11*sigma22-sigma12*sigma21)+log(abs(1-alpha1*beta1))
logf=const+log(lamda*exp(-0.5*expon1)+(1-lamda)*exp(-0.5*expon2))
log.L=-sum(logf)
return(log.L)
}


####
## estimate with nonlinear constraint
hin=function(parm){
h=rep(NA,1)
h[1]=parm[11]*parm[14]-parm[12]*parm[13]
h[2]=parm[11]
h[3]=parm[14]
h[4]=-parm[12]
h[5]=-parm[13]
h
}

heq=function(parm){
h=rep(NA,1)
h[1]=parm[12]-parm[13]
h
}

max.like=constrOptim.nl(par=c(-0.5,-0.5,-0.5,-0.5,0.02,-0.02,-0.02,-0.02,-0.03,0.02,1.9,-1.1,-1.1,1.9),fn=log.L,
hin=hin,heq=heq)

######
parm=max.like$par
alpha0=parm[1]
alpha1=parm[2]
alpha2=parm[3]
beta0=parm[4]
beta1=parm[5]
beta21=parm[6]
beta22=parm[7]
beta23=parm[8]
beta24=parm[9]
beta3=parm[10]
sigma11=parm[11]
sigma12=parm[12]
sigma21=parm[13]
sigma22=parm[14]
u1=-alpha0-alpha1*logp-alpha2*lakes+logq
u21=-beta0-beta1*logq-beta21*dm1-beta22*dm2-beta23*dm3-beta24*dm4-beta3+logp
u22=-beta0-beta1*logq-beta21*dm1-beta22*dm2-beta23*dm3-beta24*dm4+logp
expon1=u1^2*sigma11+u1*u21*sigma12+u1*u21*sigma21+u21^2*sigma22
expon2=u1^2*sigma11+u1*u22*sigma12+u1*u22*sigma21+u22^2*sigma22

h1_log=(-log(2*pi)+0.5*log(sigma11*sigma22-sigma12*sigma21))+(log(abs(1-alpha1*beta1))-0.5*expon1)
h2_log=(-log(2*pi)+0.5*log(sigma11*sigma22-sigma12*sigma21))+(log(abs(1-alpha1*beta1))-0.5*expon2)
w1=w

w=1/(1+(1-lamda)/lamda*exp(h2_log-h1_log))
if(cor(w,w1)>0.99) break
}





-- 
View this message in context: http://r.789695.n4.nabble.com/Help-Maximum-likelihood-estimation-tp3006584p3018477.html
Sent from the R help mailing list archive at Nabble.com.