#*************** 
#the function **  

mymat<-function(par,data) {

# define the parameter matrix used in following part
vbar2<-matrix(0,n,nt)
vbar3<-matrix(0,n,nt)
v8 <-matrix(0,n,nt)
regw<-matrix(0,n,nt)
wden<-matrix(0,n,nt)
lia<-matrix(0,n,nt)
ccl<-matrix(1,n,ns)
eta<-c(0,0)

# setup the parts for loglikelihood
q1<-exp(par[1])
pr1<-q1/(1+q1)
pr2<-1-pr1

eta[2]<-par[2]

a<-par[3:10]
b<-par[11:19]
w<-par[20:npar]

for(m in 1:ns){

regw<-w[1]*acwrk+w[2]*actr+w[3]+w[4]*eta[m]
	
vbar2=a[1]+     eta[m]+a[2]*acwrk+a[3]*actr+a[4]*edu+a[5]*v_refg+a[6]*v_econ+a[7]*age+
	regw*a[8]

vbar3=b[1]+b[2]*eta[m]+b[3]*acwrk+b[4]*actr+b[5]*edu+b[6]*v_refg+b[7]*v_econ+b[8]*age+
	regw*b[9]

v8=1+exp(vbar2)+exp(vbar3)

lia<-ifelse(home==1,1/v8,
	ifelse(wrk==1,exp(vbar2)/v8,
		ifelse(tr==1,exp(vbar3)/v8,1)))

wden<-ifelse(wrk==1,dnorm((lnw-regw)/w[5])/w[5],1)


ccl[,m]<-lia[,1]*lia[,2]*lia[,3]*lia[,4]*lia[,5]*lia[,6]*lia[,7]*lia[,8]*
	wden[,1]*wden[,2]*wden[,3]*wden[,4]*wden[,5]*wden[,6]*wden[,7]*wden[,8]
}


#****************************
#cat("w=",w, "\n")
#cat("vbar2=",vbar2[1,], "\n")
#cat("vbar3=",vbar3[1,], "\n")
#cat("v8=",v8[1,], "\n")
#cat("lia=",lia[1,], "\n")
#cat("wden=",wden[1,], "\n")
#cat("lnw=",head(lnw), "\n")
#cat("regw=",regw[1,], "\n")
#cat("ccl=",ccl[1:6,], "\n")
#cat("---------------", "\n")
#****************************



func<-pr1*ccl[,1]+pr2*ccl[,2]
func<-ifelse(func<.Machine$double.xmin,0.00001,func)
f<-sum(log(func))
return(-f)
}


#*********************************
mydata<-read.table("j:/new folder/mydata9x.txt", head=F)

nt<<-8	# number of periods
ns<<-2	# number of person type
n<<-50	# number of people
npar<<-24	# number of parameters

id<-as.numeric(mydata[,1])
tr<-as.matrix(mydata[,2:(nt+1)])
wrk<-as.matrix(mydata[,(nt+2):(2*nt+1)])
home<-as.matrix(mydata[,(2*nt+2):(3*nt+1)])
actr<-as.matrix(mydata[,(3*nt+2):(4*nt+1)])
acwrk<-as.matrix(mydata[,(4*nt+2):(5*nt+1)])
lnw<-as.numeric(mydata[,5*nt+2])
edu<-as.numeric(mydata[,5*nt+3])
age<-as.numeric(mydata[,5*nt+4])
v_refg<-as.numeric(mydata[,5*nt+5])
v_econ<-as.numeric(mydata[,5*nt+6])

# initial value
guess<-rep(0.1,times=npar)
guess[npar]<-1.0

system.time(r1<-optim(guess,mymat,data=mydata, hessian=F))
guess<-r1$par

system.time(r2<-optim(guess,mymat,data=mydata, method="BFGS",hessian=T, 
		control=list(trace=T, maxit=1000)))

std.err<-sqrt(diag(solve(r2$hessian)))
res<-cbind(r2$par,std.err,r2$par/std.err)
colnames(res)<-c("parameter","std.err","t test")
#res