*! xtloglin 1.0.1 12sep2023
*! author david vincent
*! email dvincent@dveconometrics.co.uk (or davidwvincent@hotmail.com)

//----------------------------------------------------------------------//

program define xtloglin, sortpreserve rclass

	version 15.1

	syntax, null(string) [NEGative *]

	tempvar touse
	qui gen byte `touse'=e(sample)

	if !inlist("`e(cmd)'","regress","xtreg"){
		di as err "xtloglin allowed after {opt regress} or {opt xtreg}"
		exit 198
	}

	if inlist("`e(cmd)'","xtreg") & !inlist("`e(model)'","fe","re","be"){
		di as err "models {opt fe,re,be} allowed after {opt xtreg}"
		exit 198
	}

	if !inlist("`null'","log","linear"){
		di as err "must specify {opt model(linear)} or {opt model(log)} as the null model to be tested"
		exit 198
	}

	if "`null'"=="log" & "`negative'"!=""{
		di as err "can only specify {opt negative} with {opt model(linear)}"
		exit 198
	}

	cap tsset

	if "`options'"!="" {

	//error structure specified in xtloglin

		_parse_error_cov, `options'

		loc clustvar = "`r(clustvar)'"
		loc evar     = "`r(evar)'"

		if "`e(cmd)'"=="xtreg" & "`evar'"=="robust"{
			loc clustvar="`e(ivar)'"
			loc evar="cluster"
		}

		if ("`e(cmd)'"=="xtreg" & "`clustvar'"!=""){

			_check_nested `clustvar', ivar(`e(ivar)') esample(`touse')
		}
	}

	else{

	//specified at estimation

		local clustvar = "`e(clustvar)'"
		local evar     = "`e(vce)'"

		if ("`e(cmd)'"=="xtreg" & "`clustvar'"!="")  loc evar="cluster"
		if ("`e(cmd)'"=="regress" & "`evar'"=="ols") loc evar="conventional"
	}

	//match xtreg,be requirements

	if ("`e(cmd)'"=="xtreg" & "`e(model)'"=="be" & "`evar'"!="conventional") {
		di as error "option robust or cluster not allowed after xtreg,be"
		exit 198
	}

	//----------------------------------------------------------------------//

	local c=0
	local xvars: colname e(b)
	foreach v of local xvars{
		local c=`c'+strmatch("`v'","_cons")
	}

	if `c'{
		local cons "_cons"
		local xvars: list xvars-cons
	}

	tempvar res diff
	tempname npanels

	local yvar="`e(depvar)'"

	if "`null'"=="linear"{

		qui count if `yvar'<=0
		if r(N)>0 & "`negative'"==""{
			di as err "untransformed variable `yvar' must be strictly positive. Specify {opt negative} to test linear model with negative outcomes"
			exit 411
		}
	}

	//inputs for LM-statistic
	mata: _lm_inputs("`yvar'","`xvars'",`c',"`null'","`touse'","`res'","`diff'","`npanels'","`negative'")

	if "`evar'"=="cluster"{
		sort `clustvar'
		markout `touse' `clustvar'
	}

	tempname lm_name p_name

	//LM-statistic
	mata: _lm_stat("`res'","`diff'","`evar'","`clustvar'","`touse'","`lm_name'", "`p_name'","`npanels'")

	return scalar lm = `lm_name'
	return scalar p  = `p_name'

	return loc null     `null'
	return loc evar     `evar'
	return loc clustvar `clustvar'

	local ttype=cond("`negative'"=="","Box-Cox","MacKinnon-Magee")

	di
	di as txt "Robust LM test of linear and log-linear functional forms"
	di as txt "H0: `null' dependent variable"
	di as txt "H1: `ttype' transformation"

	di in smcl as text "{hline 57}"
	di as txt _col(5) "LM-chi2(1) = " as res %10.3f `lm_name'
	di as txt _col(5) "Prob > LM  = " as res %11.4f `p_name'
	di in smcl as text "{hline 57}"
	di as text "Error Variance: `evar'"

	if "`clustvar'"!=""{
		di as txt "Cluster Variable: `clustvar'"
	}

end

 //----------------------------------------------------------------------//

 program _check_nested, sortpreserve

	syntax varname, ivar(varname) esample(varname)

	local cvar `varlist'

	if "`cvar'" == "`ivar'" exit

	marksample touse
	qui replace `touse'=0 if !`esample'

	tempvar uniq tot

	sort `ivar' `cvar'
	qui egen `uniq'=tag(`ivar' `cvar') if `touse'
	qui egen `tot'=total(`uniq'), by(`ivar')

	qui count if `tot'>1 & `touse'

	if r(N)>0{
		di as err "panels are not nested within clusters"
		exit 498
	}

 end

 //----------------------------------------------------------------------//

 program _parse_error_cov, rclass

	syntax, [NOTRobust ROBust CLUSter(varname)]

	local cluster_opt=cond("`cluster'"!="","cluster(varname)","")

	opts_exclusive "`notrobust' `robust' `cluster_opt'"

	if "`cluster'"!="" {

		confirm var `cluster'
		return loc evar =  "cluster"
		return loc clustvar  = "`cluster'"
	}

	if "`robust'"!="" {

		return loc evar =  "robust"
		return loc clustvar  = ""
	}

	if "`notrobust'"!="" {

		return loc evar =  "conventional"
		return loc clustvar  = ""
	}

end

//--------------------------------------------------------------------------//

mata:

struct data
{

 string scalar yvar,xvars
 string scalar touse
 string scalar null
 string scalar negative
 string scalar idvar

 real   matrix X,y
 real   matrix T_lam
 real   matrix id_rows

 real   matrix beta
 real   matrix yhat_n
 real   matrix res
 real   matrix diff

 real   matrix MZZ
 real   matrix MXX
 real   matrix MZT_lam
 real   matrix MXT_lam

 real   scalar K
 real   scalar N
 real   scalar NT
 real   scalar cons
 real   scalar sigma_e,sigma_u,sigma_b,rmse

}

//--------------------------------------------------------------------------//

void setup(struct data scalar z,
           string scalar yvar,
		   string scalar xvars,
		   real   scalar cons,
		   string scalar null,
		   string scalar negative,
		   string scalar touse)
{

 real scalar sf
 real vector id,y_T

 st_view(z.y,.,yvar,touse)			//y-transformed = y_T, or log(y_T)
 st_view(z.X,.,tokens(xvars),touse)
 z.beta=st_matrix("e(b)")

 z.negative=negative
 z.touse=touse
 z.null=null
 z.cons=cons
 z.NT=rows(z.X)
 z.K=cols(z.X)
 z.N=z.NT


 if(st_global("e(cmd)")=="xtreg"){

	z.idvar=st_global("e(ivar)")
	st_view(id,.,z.idvar,touse)
	z.id_rows=panelsetup(id,1)
	z.N=rows(z.id_rows)

	if(st_global("e(model)")=="fe" | st_global("e(model)")=="re"){

		z.sigma_e=st_numscalar("e(sigma_e)")
		z.sigma_u=st_numscalar("e(sigma_u)")
	}
 }

 //untransformed y-variable
 y_T=(z.null=="linear" ? z.y : exp(z.y))

 if(z.negative=="negative" & z.null=="linear"){

	//derivative of MacKinnon-Magee transformation for testing linearity (theta=0)
	z.T_lam=y_T:^2

 }
 else{

	//scaling factor (Box Cox only)
	sf=exp(mean(log(abs(y_T))))

	//derivative of scaled Box-Cox transformation for testing linear & log-linear (lambda=1 | 0)
	z.T_lam=(z.null=="linear" ? (y_T:*(ln(y_T:/sf):-1):+1) : ln(y_T):*(ln(sqrt(y_T):/sf)))

 }

}

//--------------------------------------------------------------------------//

void _lm_stat(string scalar resvar,
			  string scalar diffvar,
			  string scalar evar,
			  string scalar clustvar,
			  string scalar touse,
		      string scalar lm_name,
			  string scalar p_name,
			  string scalar npanels)
{

 real scalar N,M,lm
 real vector res,diff,clust

 if(st_global("e(model)")!="be"){

	st_view(res,.,resvar,touse)
	st_view(diff,.,diffvar,touse)
 }
 else{

	//number of panels (first N-rows in dataset)
	N=st_numscalar(npanels)
 	st_view(res,(1::N),resvar)
	st_view(diff,(1::N),diffvar)

 }

 if(evar=="conventional"){

	N=rows(res)
	M=quadcross(diff,diff)*(quadcross(res,res)/N)

 }

 if(evar=="robust"){

	M=quadcross(diff,res:^2,diff)
 }

 if(evar=="cluster"){

	st_view(clust,.,clustvar,touse)
	M=_cluster_cross(clust,res,diff)

 }

 lm=(quadcross(res,diff)^2)/M

 st_numscalar(lm_name,lm)
 st_numscalar(p_name,1-chi2(1,lm))

}

//--------------------------------------------------------------------------//

real scalar _cluster_cross(real vector clust,
				           real vector res,
						   real vector diff)
{

 real scalar N,M,i
 real vector res_i,diff_i,clust_rows

 clust_rows=panelsetup(clust,1)
 N=rows(clust_rows)
 M=0

 for(i=1;i<=N ;i++){

		res_i=panelsubmatrix(res,i,clust_rows)
		diff_i=panelsubmatrix(diff,i,clust_rows)

		M=M+diff_i'(res_i*res_i')*diff_i

 }

 return(M)

}

//--------------------------------------------------------------------------//

void _lm_inputs(string scalar yvar,
				string scalar xvars,
				real   scalar cons,
			    string scalar null,
				string scalar touse,
				string scalar resvar,
				string scalar diffvar,
				string scalar npanels,
				string scalar negative)
{

 struct data scalar z

 setup(z,yvar,xvars,cons,null,negative,touse)

 if(st_global("e(model)")=="ols"){
	_ols_inputs(z)
 }

 if(st_global("e(model)")=="re"){
	_re_inputs(z)

 }

 if(st_global("e(model)")=="fe"){
	_fe_inputs(z)
 }

 if(st_global("e(model)")=="be"){
	_be_inputs(z)
 }


 if(st_global("e(model)")!="be"){

	//row specific, order matters
	st_store(.,st_addvar("double",(resvar,diffvar)),touse,(z.res,z.diff))
 }
 else{

	//add means to first N-rows order irrelevant as clustering not allowed
	st_store((1::z.N),st_addvar("double",(resvar,diffvar)),(z.res,z.diff))
 }

 st_numscalar(npanels,z.N)

}

//--------------------------------------------------------------------------//

void _ols_inputs(struct data scalar z)
{

 real scalar min,max
 real vector beta,res,diff
 real matrix X,Z,yhat_n,MZZ_inv,MXX_inv

 beta=z.beta

 res=J(z.NT,1,.)
 diff=J(z.NT,1,.)

 X=(z.cons ? (z.X,J(rows(z.X),1,1)) : z.X)

 //use normalised_yhat(it) = a*yhat(it)+b to avoid huge yhat^r values (no impact on LM-stat)
 yhat_n=X*beta'
 min=min(yhat_n)
 max=max(yhat_n)
 yhat_n = (z.null=="linear" & z.negative!="negative" ? yhat_n:^-1 : (yhat_n:-min):/(max:-min))

//instruments
 Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
 _editmissing(Z,0)

 MZZ_inv=invsym(quadcross(Z,Z))
 MXX_inv=invsym(quadcross(X,X))

 //predict T_lam using Z then partial out X
 z.diff=(Z*MZZ_inv*quadcross(Z,z.T_lam))-(X*MXX_inv*quadcross(X,z.T_lam))

 //residuals
 z.res=(z.y-X*beta')

}

//--------------------------------------------------------------------------//

void _fe_inputs(struct data scalar z)
{

 real scalar i,max,min,T,st,fi
 real vector beta,beta1,res,diff,y,e
 real matrix X,Z,yhat_n,MZZ_inv,MXX_inv,Q

 beta=z.beta
 beta1=select(z.beta,!e(z.K+1,z.K+1))

 res=J(z.NT,1,.)
 diff=J(z.NT,1,.)

 //normalised fitted values (before FE-transformation)
 X=(z.X,J(rows(z.X),1,1))
 z.yhat_n=X*beta'
 min=min(z.yhat_n)
 max=max(z.yhat_n)
 z.yhat_n = (z.null=="linear" & z.negative!="negative" ? z.yhat_n:^-1 : (z.yhat_n:-min):/(max:-min))

 _fe_ZZ_XX_components(z)

 MZZ_inv=invsym(z.MZZ)
 MXX_inv=invsym(z.MXX)

 for(i=1;i<=z.N;i++){

	y=panelsubmatrix(z.y,i,z.id_rows)
	X=panelsubmatrix(z.X,i,z.id_rows)
	yhat_n=panelsubmatrix(z.yhat_n,i,z.id_rows)

	T=rows(X)
	e=J(T,1,1)
	st=z.id_rows[i,1]
	fi=z.id_rows[i,2]
	Q=I(T)-e*e'/T

	//instruments
	Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
	_editmissing(Z,0)

	//transformed variables
	y=Q*y
	X=Q*X
	Z=Q*Z

	//predict T_lam_FE using Z_FE then partial out X_FE
	diff[st::fi]=Z*MZZ_inv*z.MZT_lam-X*MXX_inv*z.MXT_lam

	//residuals FE=eit
 	res[st::fi]=(y-X*beta1')

 }

 z.diff=diff
 z.res=res


}

//--------------------------------------------------------------------------//

void _fe_ZZ_XX_components(struct data scalar z)
{

 real scalar i,T
 real vector T_lam,yhat_n,e
 real matrix X,Z,Q

 z.MZZ     = J(z.K+4,z.K+4,0)
 z.MXX     = J(z.K,z.K,0)
 z.MZT_lam = J(z.K+4,1,0)
 z.MXT_lam = J(z.K,1,0)

 for(i=1;i<=z.N;i++){

	X=panelsubmatrix(z.X,i,z.id_rows)
	T_lam=panelsubmatrix(z.T_lam,i,z.id_rows)
	yhat_n=panelsubmatrix(z.yhat_n,i,z.id_rows)

	T=rows(X)
	e=J(T,1,1)
	Q=I(T)-e*e'/T

	//instruments
	Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
	_editmissing(Z,0)

	z.MZZ=z.MZZ+Z'Q*Z
	z.MXX=z.MXX+X'Q*X

	z.MZT_lam=z.MZT_lam+Z'Q*T_lam
	z.MXT_lam=z.MXT_lam+X'Q*T_lam

 }


}

//--------------------------------------------------------------------------//

void _re_inputs(struct data scalar z)
{

 real scalar i,max,min,T,st,fi,phi2
 real vector beta,res,diff,y,e
 real matrix X,Z,yhat_n,MZZ_inv,MXX_inv,Q,V

 beta=z.beta

 res=J(z.NT,1,.)
 diff=J(z.NT,1,.)

 //normalised fitted values (before RE-transformation)
 X=(z.X,J(rows(z.X),1,1))
 z.yhat_n=X*beta'
 min=min(z.yhat_n)
 max=max(z.yhat_n)
 z.yhat_n = (z.null=="linear" & z.negative!="negative" ? z.yhat_n:^-1 : (z.yhat_n:-min):/(max:-min))

 _re_ZZ_XX_components(z)

 MZZ_inv=invsym(z.MZZ)
 MXX_inv=invsym(z.MXX)

 for(i=1;i<=z.N;i++){

	y=panelsubmatrix(z.y,i,z.id_rows)
	X=panelsubmatrix(z.X,i,z.id_rows)
	yhat_n=panelsubmatrix(z.yhat_n,i,z.id_rows)
	T=rows(X)
	X=(X,J(T,1,1))

	e=J(T,1,1)
	st=z.id_rows[i,1]
	fi=z.id_rows[i,2]
	Q=I(T)-e*e'/T
	phi2=(z.sigma_e^2)/(T*z.sigma_u^2+z.sigma_e^2)
	V=(Q+sqrt(phi2)*(I(T)-Q)) 							//note sqrt transform

	//instruments
	Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
	_editmissing(Z,0)

	//transformed variables
	y=V*y
	X=V*X
	Z=V*Z

	//predict T_lam_RE using Z_RE then partial out X_RE
	diff[st::fi]=Z*MZZ_inv*z.MZT_lam-X*MXX_inv*z.MXT_lam

	//residuals_RE
 	res[st::fi]=(y-X*beta')

 }

 z.diff=diff
 z.res=res

}

//--------------------------------------------------------------------------//

void _re_ZZ_XX_components(struct data scalar z)
{

 real scalar i,T,phi2
 real vector T_lam,yhat_n,e
 real matrix X,Z,Q,V

 z.MZZ     = J(z.K+5,z.K+5,0)
 z.MXX     = J(z.K+1,z.K+1,0)
 z.MZT_lam = J(z.K+5,1,0)
 z.MXT_lam = J(z.K+1,1,0)

 for(i=1;i<=z.N;i++){

 	X=panelsubmatrix(z.X,i,z.id_rows)
	T_lam=panelsubmatrix(z.T_lam,i,z.id_rows)
	yhat_n=panelsubmatrix(z.yhat_n,i,z.id_rows)
	T=rows(X)
	X=(X,J(T,1,1))

	e=J(T,1,1)
	Q=I(T)-e*e'/T
	phi2=(z.sigma_e^2)/(T*z.sigma_u^2+z.sigma_e^2)
	V=(Q+phi2*(I(T)-Q)) 							  //note no sqrt transform

	//instruments
	Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
	_editmissing(Z,0)

	z.MZZ=z.MZZ+Z'V*Z
	z.MXX=z.MXX+X'V*X

	z.MZT_lam=z.MZT_lam+Z'V*T_lam
	z.MXT_lam=z.MXT_lam+X'V*T_lam

 }

}

//--------------------------------------------------------------------------//

void _be_inputs(struct data scalar z)
{

 real scalar i,max,min,T,st,fi,w
 real vector beta,res,diff,y,e
 real matrix X,Z,yhat_n,MZZ_inv,MXX_inv,P

 beta=z.beta

 res=J(z.N,1,.)
 diff=J(z.N,1,.)

 //normalised fitted values (before BE-transformation)
 X=(z.X,J(rows(z.X),1,1))
 z.yhat_n=X*beta'
 min=min(z.yhat_n)
 max=max(z.yhat_n)
 z.yhat_n = (z.null=="linear" & z.negative!="negative" ? z.yhat_n:^-1 : (z.yhat_n:-min):/(max:-min))

 _be_ZZ_XX_components(z)

 MZZ_inv=invsym(z.MZZ)
 MXX_inv=invsym(z.MXX)


 for(i=1;i<=z.N;i++){

	y=panelsubmatrix(z.y,i,z.id_rows)
	X=panelsubmatrix(z.X,i,z.id_rows)
	yhat_n=panelsubmatrix(z.yhat_n,i,z.id_rows)
	T=rows(X)
	X=(X,J(T,1,1))

	//weight by sqrt(Ti) for WLS
	w=(st_global("e(typ)")!="" ? sqrt(T) : 1)	//sqrt not 1
	e=J(T,1,1)
	P=(e'/T)*w

	//instruments
	Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
	_editmissing(Z,0)

	//transformed variables (panel means)
	y=P*y
	X=P*X		//i.e., 1 x K
	Z=P*Z

	//predict T_lam_BE using Z_BE then partial out X_BE
	diff[i]=Z*MZZ_inv*z.MZT_lam-X*MXX_inv*z.MXT_lam

	//residuals_RE
 	res[i]=(y-X*beta')

 }

 z.diff=diff
 z.res=res

}

//--------------------------------------------------------------------------//

void _be_ZZ_XX_components(struct data scalar z)
{

 real scalar i,T,w
 real vector T_lam,yhat_n,e
 real matrix X,Z,P

 z.MZZ     = J(z.K+5,z.K+5,0)
 z.MXX     = J(z.K+1,z.K+1,0)
 z.MZT_lam = J(z.K+5,1,0)
 z.MXT_lam = J(z.K+1,1,0)

 for(i=1;i<=z.N;i++){

 	X=panelsubmatrix(z.X,i,z.id_rows)
	T_lam=panelsubmatrix(z.T_lam,i,z.id_rows)
	yhat_n=panelsubmatrix(z.yhat_n,i,z.id_rows)
	T=rows(X)
	X=(X,J(T,1,1))

	//weight by Ti for WLS
	w=(st_global("e(typ)")!="" ? 1 : 2)
	e=J(T,1,1)
	P=(e*e')/(T^w)

	//instruments
	Z=(X,yhat_n,yhat_n:^2,yhat_n:^3,yhat_n:^4)
	_editmissing(Z,0)

	z.MZZ=z.MZZ+Z'P*Z
	z.MXX=z.MXX+X'P*X

	z.MZT_lam=z.MZT_lam+Z'P*T_lam
	z.MXT_lam=z.MXT_lam+X'P*T_lam

 }

}

end