program define xtmoralb, eclass
	syntax varlist(ts) [if] [in] [,omethod(int 1) time robust]
	qui: tsset
	local pvar="`r(panelvar)'"
	gettoken depvar indepvar: varlist
	qui: xtreg `depvar' `indepvar'
	gen ok=e(sample)
	egen sumok=total(ok), by(`pvar')
	qui: su sumok, meanonly
	gen allok=(sumok==r(max))
	qui: replace `depvar'=. if allok==0
	qui: xtreg `depvar' `indepvar' if allok==1
	drop ok sumok allok
	gen ___kike=e(N_g)
	marksample touse

	mata: xtmoralb_mata("`varlist'","`touse'","`omethod'","`time'","`robust'")
	
	tempname b bpar V
	matrix `b'=r(b)'
	matrix `V'=r(V)
	local N=r(N)
	local t=r(t)
	local lo=r(lo)
	local nobs=r(nobs)
	local kt=r(kt)
	local vnames `indepvar'
	drop ___kike
	matrix rownames `V' = `vnames'
	matrix colnames `V' = `vnames'
	matrix colnames `b' = `vnames'
	display _newline "Moral-Benito (2011): sub-system LIML Dynamic Panel Data"
	display as text "Likelihood counterpart of first-diff GMM (i.e. Arellano-Bond)"
	display as text "Estimation results" _col(58) as text "Number of units = " as res `N'
	display as text "log likelihood = " as res `lo' _col(58) as text "Number of periods = " as res `t'
	display as text "number of parameters = " as res `kt' _col(58) as text "Number of obs = " as res `nobs'
	ereturn post `b' `V', depname(`depvar') obs(`N') esample(`touse')
	if "`robust'" != "" ereturn local vcetype "Robust"
	ereturn display
	display as txt "Number of parameters include all var-covar parameters"
	if "`omethod'" == "1" display as txt "Optimization method: Broyden-Fletcher-Goldfarb-Shanno"
	if "`omethod'" == "2" display as txt "Optimization method: Newton-Raphson"	
	if "`omethod'" == "3" display as txt "Optimization method: Davidon-Fletcher-Powell"
	if "`omethod'" == "4" display as txt "Optimization method: Berndt-Hall-Hall-Hausman"
	if "`time'" == "time" display as txt "Time dummies included in the model"
end

mata:
void xtmoralb_mata(string scalar varlist, 
		   string scalar touse, 
		   string scalar omethod,
		   string scalar time,
		   string scalar robust){

	/* data in stata is organized in panel format: first all observations for individual 1. */
	/* that is to say: sort individual time */

	real matrix dmean, data, dataliml, V
	real colvector bliml
	real rowvector t0, bl
	real scalar n, k, t, nobs, j, i, npar, ip
	
	data=st_data(.,st_tsrevar(tokens(varlist)),touse)
	n=st_data(1,("___kike"))
	nobs=rows(data)
	t=nobs/n
if(time=="time"){
	dmtx=I(n*t)-J(n,n,1)#((1/n)*I(t))
	data=dmtx*data
}
if (data[1,1]==data[2,2] & data[2,1]==data[3,2]){
	k=cols(data)-2
	lagd=1
}
else if (data[1,1]!=data[2,2] & data[2,1]!=data[3,2]){
	k=cols(data)-1
	lagd=0
}
if (k>0 & lagd==1 & time==""){
	dmean=J(n*t,1,1)*(mean(data[.,3..(k+2)]))
	data[.,3..(k+2)]=(data[.,3..(k+2)]):/(dmean)
}
else if (k>0 & lagd==0 & time==""){
	dmean=J(n*t,1,1)*(mean(data[.,2..(k+1)]))
	data[.,2..(k+1)]=(data[.,2..(k+1)]):/(dmean)
}

if (lagd==1){
	dataliml=J(n,t+1+k*t,.)
	j=1
	while(j<=n){
		dataliml[j,(cols(dataliml)-k)..(cols(dataliml))]=data[(j-1)*t+1,2..(k+2)]
		j++
	}
}
else if (lagd==0){
	dataliml=J(n,t+k*t,.)
	j=1
	while(j<=n){
		dataliml[j,(cols(dataliml)-k+1)..(cols(dataliml))]=data[(j-1)*t+1,2..(k+1)]
		j++
	}
}
	i=1
	while(i<=t){
		j=1
		while (j<=n){
			dataliml[j,i]=data[(j-1)*t+i,1]
			j++
		}
		i++
	}

if (k>0 & lagd==1){
	i=1
	while(i<=(t-1)){
		j=1
		while (j<=n){
			dataliml[j,((t+(i-1)*k)+1)..((t+(i-1)*k)+k)]=data[(j-1)*t+(i+1),3..(cols(data))]
			j++
		}
		i++
	}
}
else if (k>0 & lagd==0){
	i=1
	while(i<=(t-1)){
		j=1
		while (j<=n){
			dataliml[j,((t+(i-1)*k)+1)..((t+(i-1)*k)+k)]=data[(j-1)*t+(i+1),2..(cols(data))]
			j++
		}
		i++
	}
}


/*********************************************************************************************/
/* NOW WE DEFINE THE INITIAL VALUES IN t0 */
/*********************************************************************************************/


	Y1=dataliml[.,1..t]
if (k>0){
	Y2=dataliml[.,(t+1)..t+k*(t-1)]
}
if (lagd==1){
	Z=dataliml[.,(t+1+k*(t-1))..(t+k+1+k*(t-1))]
}
else if (lagd==0){
	Z=dataliml[.,(t+1+k*(t-1))..(t+k+k*(t-1))]
}
	Q=I(n)-Z*invsym(Z'*Z)*Z'
	d=J(t-1,t,0)
	dmx=1
	while (dmx<=t-1){
		d[dmx,dmx]=-1
		d[dmx,dmx+1]=1
		dmx=dmx++
	}
	a=cholesky(invsym(d*d'))*d
	datao=(I(n)#a)*data
if (lagd==1){
	bwg=invsym(datao[.,2..(2+k)]'datao[.,2..(2+k)])*datao[.,2..(2+k)]'datao[.,1]
	bols=invsym(data[.,2..(2+k)]'data[.,2..(2+k)])*data[.,2..(2+k)]'data[.,1]
}
else if (lagd==0){
	bwg=invsym(datao[.,2..(1+k)]'datao[.,2..(1+k)])*datao[.,2..(1+k)]'datao[.,1]
	bols=invsym(data[.,2..(1+k)]'data[.,2..(1+k)])*data[.,2..(1+k)]'data[.,1]
}
	mm=(1/t)*J(t,t,1)
	datam=(I(n)#mm)*data
	datam1=J(n,cols(data),0)
	im=1
	while (im<=n){
		datam1[im,.]=datam[(im-1)*t+1,.]
		im=im++
	}
if (lagd==1){
	eta=datam1[.,1]-datam1[.,2..2+k]*bols
}
else if (lagd==0){
	eta=datam1[.,1]-datam1[.,2..1+k]*bols
}
	g01=invsym(Z'*Z)*(Z'*eta)
	eta1=Q*eta
if (lagd==1){
	vv=data[.,1]-data[.,2..2+k]*bols-eta1#J(t,1,1)
}
else if (lagd==0){
	vv=data[.,1]-data[.,2..1+k]*bwg-eta1#J(t,1,1)
}
	v=J(n,t,0)
	i=1
	while(i<=t){
		j=1
		while (j<=n){
			v[j,i]=vv[(j-1)*t+i,1]
			j++
		}
		i++
	}
if (k>0){
	e=J(n,k*(t-1),0)
	it=1
	while (it<=t-1){
		e[.,(it-1)*k+1..(it-1)*k+k]=Q*Y2[.,(((it-1)*k)+1)..(((it-1)*k)+k)]
		it=it++
	}
}
	inits=bols\g01
if(lagd==1){
	inits=inits\(sqrt(variance(eta1)))
	inits=inits\(sqrt(diagonal(variance(v))))
}
else if(lagd==0){
	chv=cholesky(variance(v))
	i=1
	while (i<=t){
		j=1
		while (j<=t){
			if(i<j){
				chv[i,j]=chv[j,i]
			}
		j++
		}
	i++
	}
	vvr=vech(chv)
	inits=inits\vvr
}

if (k>0){
	meeta=(e,eta1)
	vcmeeta=variance(meeta)
	inits=inits\vcmeeta[rows(vcmeeta),1..cols(vcmeeta)-1]'
	vcmev=(v,e)
	vcmev=variance(vcmev)
	it=1
	while (it<=t-1){
		inits=inits\(vcmev[it,(t+1)+(it-1)*k..cols(vcmev)]')
		it=it++
	}
}
	t0=inits'



/*********************************************************************************************/
/*         FINALLY WE CALL THE OPTIMIZATION ROUTINE AND MAXIMIZE THE LIKELIHOOD              */
/*********************************************************************************************/

	S=optimize_init()
if (k>0 & lagd==1){
	optimize_init_evaluator(S,&xtliml_moralb())
}
else if (k>0 & lagd==0){
	optimize_init_evaluator(S,&xtliml_moralb_noy())
}
else if (k==0){
	optimize_init_evaluator(S,&xtliml_moralb_nox())
}
	optimize_init_evaluatortype(S,"v0")
	optimize_init_which(S,"max")
	optimize_init_params(S,t0)
	optimize_init_argument(S,1,dataliml)
	optimize_init_argument(S,2,n)
	optimize_init_argument(S,3,nobs)
	optimize_init_argument(S,4,k)
	optimize_init_singularHmethod(S,"hybrid")

if(omethod=="1"){
	optimize_init_technique(S,"bfgs")
}
else if(omethod=="2"){
	optimize_init_technique(S,"nr")
}
else if(omethod=="3"){
	optimize_init_technique(S,"dfp")
}
else if(omethod=="4"){
	optimize_init_technique(S,"bhhh")
}
else{
	optimize_init_technique(S,"bfgs")
}

	bl=optimize(S)
	lo=optimize_result_value(S)
if(robust==""){
	V=optimize_result_V_oim(S)
}
else{
	V=optimize_result_V_oim(S)
}

if (k>0 & lagd==1 & time==""){
	V=V[(1::(k+1)),(1::(k+1))]
	Vs=diagonal(V):/((1,dmean[1,.]):^2)'
	V=(Vs*(J(1,(k+1),1))):*I(k+1)
	bl[2..(k+1)]=bl[2..(k+1)]:/dmean[1,.]
	b=bl[1..k+1]
}
else if (k>0 & lagd==0 & time==""){
	V=V[(1::(k)),(1::(k))]
	Vs=diagonal(V):/((dmean[1,.]):^2)'
	V=(Vs*(J(1,(k),1))):*I(k)
	bl[1..(k)]=bl[1..(k)]:/dmean[1,.]
	b=bl[1..k]
}
else if (k>0 & lagd==1 & time=="time"){
	V=V[(1::(k+1)),(1::(k+1))]
	b=bl[1..k+1]
}
else if (k>0 & lagd==0 & time=="time"){
	V=V[(1::(k)),(1::(k))]
	b=bl[1..k]
}
else if (k==0){
	V=V[(1::(k+1)),(1::(k+1))]
	b=bl[1..k+1]
}
	kt=cols(bl)
	st_eclear()
	st_matrix("r(b)",b')
	st_matrix("r(V)",V)
	st_numscalar("r(N)",n)
	st_numscalar("r(t)",t)
	st_numscalar("r(lo)",lo)
	st_numscalar("r(nobs)",nobs)
	st_numscalar("r(kt)",kt)
}
end




/*********************************************************************************************/
/*		         DEFINITION OF THE FUNCTION WITH THE LIKELIHOOD               	     */
/*********************************************************************************************/


mata:
void xtliml_moralb(todo, t0, dataliml, n, nobs, k, L, g, H){
	t=nobs/n
	Y1=dataliml[.,1..t]
	Y2=dataliml[.,(t+1)..t+k*(t-1)]
	Z=dataliml[.,(t+1+k*(t-1))..(t+k+1+k*(t-1))]
	Q=I(n)-Z*invsym(Z'*Z)*Z'
	B0=I(t+k*(t-1))
	C0=J(t,(k+1),0)
	C0[1,1]=t0[1]+t0[k+2]
	C0[1,2..(k+1)]=t0[2..(k+1)]+t0[(k+3)..(2*k+2)]
	
	ii=2
	while(ii<=t){
		B0[ii,ii-1]=-t0[1]
		B0[ii,k*(ii-2)+t+1..k*(ii-2)+t+k]=-t0[2..(k+1)]
		C0[ii,1]=t0[k+2]
		C0[ii,2..(k+1)]=t0[(k+3)..(2*k+2)]
		ii++
	}

	B110=B0[1::t,1..t]
	B120=B0[1::t,(t+1)..(t+k*(t-1))]

	o110=J(t,t,0)
	ii=1
	while(ii<=t){
		o110[ii,ii]=t0[(2*k+3)+ii]^2
		ii++
	}
	
	o110=o110+(t0[(2*k+3)]^2)*(J(t,1,1)*J(1,t,1))

	cpkt0=2*k+4+t

	o120p=J(t,(t-1)*k,0)
	ii=1
	while(ii<=t){
		o120p[ii,.]=t0[cpkt0..((cpkt0-1)+(t-1)*k)]
		ii++
	}

	cpkt1=cpkt0+(t-1)*k
	
	o120ps=J(t,(t-1)*k,0)
	ii=1
	while(ii<=(t-1)){
		o120ps[ii,(ii-1)*k+1..k*(t-1)]=t0[cpkt1..(cpkt1+k*(t-ii)-1)]
		cpkt1=cpkt1+k*(t-ii)
		ii++
	}

	o120=o120p+o120ps
	

	U10=(B110*Y1'+B120*Y2'-C0*Z')'
	HH=(Y2-U10*invsym(o110)*o120)'*Q*(Y2-U10*invsym(o110)*o120)
	L=-(n/2)*ln(det(o110))-(1/2)*sum(diagonal(invsym(o110)*U10'U10))-(n/2)*ln(det(HH/n))


}
end

mata:
void xtliml_moralb_nox(todo, t0, dataliml, n, nobs, k, L, g, H){
	t=nobs/n
	Y1=dataliml[.,1..t]
	Z=dataliml[.,(t+1+k*(t-1))..(t+k+1+k*(t-1))]
	Q=I(n)-Z*invsym(Z'*Z)*Z'
	B0=I(t+k*(t-1))
	C0=J(t,(k+1),0)
	C0[1,1]=t0[1]+t0[k+2]
	
	ii=2
	while(ii<=t){
		B0[ii,ii-1]=-t0[1]
		C0[ii,1]=t0[k+2]
		ii++
	}

	B110=B0[1::t,1..t]

	o110=J(t,t,0)
	ii=1
	while(ii<=t){
		o110[ii,ii]=t0[(2*k+3)+ii]^2
		ii++
	}
	
	o110=o110+(t0[(2*k+3)]^2)*(J(t,1,1)*J(1,t,1))

	
	U10=(B110*Y1'-C0*Z')'
	L=-(n/2)*ln(det(o110))-(1/2)*sum(diagonal(invsym(o110)*U10'U10))

}
end


mata:
void xtliml_moralb_noy(todo, t0, dataliml, n, nobs, k, L, g, H){
	t=nobs/n
	Y1=dataliml[.,1..t]
	Y2=dataliml[.,(t+1)..t+k*(t-1)]
	Z=dataliml[.,(t+1+k*(t-1))..(t+k+k*(t-1))]
	Q=I(n)-Z*invsym(Z'*Z)*Z'
	B0=I(t+k*(t-1))
	C0=J(t,k,0)
	C0[1,1..k]=t0[1..k]+t0[(k+1)..(2*k)]
	
	ii=2
	while(ii<=t){
		B0[ii,k*(ii-2)+t+1..k*(ii-2)+t+k]=-t0[1..k]
		C0[ii,1..k]=t0[(k+1)..(2*k)]
		ii++
	}

	B110=B0[1::t,1..t]
	B120=B0[1::t,(t+1)..(t+k*(t-1))]

	o110=J(t,t,.)
	hp=1
	ii=1
	while(ii<=t){
		jj=1
		while(jj<=t){
			if(ii==jj){
			o110[ii,jj]=t0[(2*k)+hp]
			hp=hp+1
			}
			else if(ii!=jj){
				if(ii<jj){
					o110[ii,jj]=t0[(2*k)+hp]
					hp=hp+1
				}
				else if(ii>jj){
					o110[ii,jj]=o110[jj,ii]
				}
			}
			jj++
		}
		ii++
	}
	o110=lowertriangle(o110)*lowertriangle(o110)'
	
	cpkt0=(2*k)+hp

	o120p=J(t,(t-1)*k,0)
	ii=1
	while(ii<=t){
		o120p[ii,.]=t0[cpkt0..((cpkt0-1)+(t-1)*k)]
		ii++
	}

	cpkt1=cpkt0+(t-1)*k
	
	o120ps=J(t,(t-1)*k,0)
	ii=1
	while(ii<=(t-1)){
		o120ps[ii,(ii-1)*k+1..k*(t-1)]=t0[cpkt1..(cpkt1+k*(t-ii)-1)]
		cpkt1=cpkt1+k*(t-ii)
		ii++
	}

	o120=o120p+o120ps
	
	
	U10=(Y1'+B120*Y2'-C0*Z')'
	HH=(Y2-U10*invsym(o110)*o120)'*Q*(Y2-U10*invsym(o110)*o120)
	L=-(n/2)*ln(det(o110))-(1/2)*sum(diagonal(invsym(o110)*U10'U10))-(n/2)*ln(det(HH/n))

}
end