program define HN, eclass
version 13.0


***************************************************************
parse "`e(cmdline)'", parse(",")
tokenize "`1'", parse(" ")
macro shift 2
local rest `*'

if "`e(cmd)'" != "heckman"  {
		di in red "last estimates not heckman "
		exit 301
	}

if "`e(method)'" !=  "two-step" {
		di in red "No two-step estimator"
		exit 301
	}		
***************************************************************

***************************************************************
preserve  
tempvar y1 y2 y2p pp up 
tempname N  pp

capture drop one
qui gen one=1
qui gen  `y2'= `e(depvar)'
qui gen  `y1'=`y2'~=.

qui predict `y2p' if `y1'==1, ycond
qui predict lp if `y1'==1, mills
qui predict zg if `y1'==1, xbsel  
qui predict p, psel  

sca N=r(N)
sca t=e(lambda)
sca s=e(sigma)
qui sum `p'
sca `pp'=r(mean)

qui drop if `y1'==0
qui gen `up'=`y2'-`y2p' if `y1'==1

***************************************************************
/*
sum lp zg  `y2' `y2p' `up' one if `y1'==1  
qui summ `up', det
local JB = (_result(1)/6)*((_result(14)^2)+[(1/4)*(_result(15)-3)^2]) 
local JBsig = chiprob(2, `JB')
noi di in gr "Jarque-Bera normality test: " /*
*/in ye %6.5g `JB', in gr /*
*/ "Chi(" in ye "2" in gr ")" , in ye %6.4f (`JBsig')
*/
**************************************************************

**************************************************************
qui mata: lmp("`up'")
sca drop s t N

local LMsig = chiprob(2, lmtest1)
noi di in gr "Pseudo-score LM normality test:" /*
*/in ye %6.5g lmtest1 in gr /*
*/", p-value=",in ye %6.4f (`LMsig') 
noi di in gr "(Pfaffermayr, M. (2014), Econometrics 2, 151-168)"
*, Chi(" in ye "2" in gr ")

ereturn scalar  lmtest=lmtest1
**************************************************************
restore
end

**************************************************************
version 13.1
        mata:
        mata set matastrict on
        void lmp(string scalar v1)
        {
                real colvector  u
                st_view(u, ., v1)
				lp=st_data(., "lp")
				zg=st_data(., "zg")
				W = st_data(., "`rest' one")
				
				s=st_numscalar("s")			
				t=st_numscalar("t")
				no=st_numscalar("e(N)")-st_numscalar("e(N_cens)")
				n=st_numscalar("e(N)")
	            
				/*moments of epsilon*/
				me=J(8,1,0)
				me[1,1]=0
				me[2,1]=s^2-t^2
				me[3,1]=0
				me[4,1]=3*me[2,1]^2
				me[5,1]=0
				me[6,1]=15*me[2,1]^3
				me[7,1]=0
				me[8,1]=105*me[2,1]^4
							
				/*Moments of u1i|di=1 */
				ml=J(no,8,0)
				ml[.,1]=lp
				ml[.,2]=J(no,1,1)-zg:*lp
				for (i=3; i<=8; i++) {
					ml[.,i]=(i-1)*ml[.,i-2]+((-zg):^(i-1)):*lp
				}
				/*printf("done ml\n")*/   
				
				/*moments of vi=u1i-lpi*/
				phi=(-lp:^1,lp:^2,-lp:^3,lp:^4,-lp:^5,lp:^6,-lp:^7,lp:^8)
				for (k=1; k<=8; k++) {
					phi[.,k]=phi[.,k]+ml[.,k]
					for (j=1; j<=k-1; j++) {
					cnk=comb(k,j)
					phi[.,k]=phi[.,k] + ((-1)^j)*cnk*ml[.,k-j]:*(lp:^j)
					}
				/*printf("done j \n")*/
				}
				/*printf("done phi \n")*/
				
				/*moments of epsilon+tau*vi*/
				psi=J(no,1,1)#me'
				for (k=1; k<=8; k++) {
						psi[.,k]=psi[.,k]+(t^k)*(phi[.,k])
						/*printf("k1=%g\n", k)*/
						for (l=1; l<=k-1; l++) {
						ckl=comb(k,l)
						psi[.,k]=psi[.,k]+(t^l)*ckl*(me[k-l])*(phi[.,l])
						}
						/*printf("done k \n")*/
				}
				/*printf("done psi \n")*/
				
				/*moment h*/
				h2=(colsum( (u:^3)-psi[.,3]), colsum( (u:^4)-psi[.,4]) )'/n
				
				/*variance of h */
				F22=	( 
						colsum( (psi[.,6]-psi[.,3]:^2)/n       ),
						colsum( (psi[.,7]-psi[.,3]:*psi[.,4])/n )  \
						colsum( (psi[.,7]-psi[.,3]:*psi[.,4])/n), 
						colsum( (psi[.,8]-psi[.,4]:^2)/n        )
						)
				
				F12=	(
						W'*psi[.,4]/n,
						W'*psi[.,5] /n  \ 
						colsum( (psi[.,5]-psi[.,2]:*psi[.,3])/n ), 
						colsum(  (psi[.,6]-psi[.,2]:*psi[.,4])/n )  
						)
				
				F11=	(
						(W'*diag(psi[.,2])*W/n),
						W'*psi[.,3]/n \
						(W'*psi[.,3]/n)',
						colsum( (psi[.,4]-psi[.,2]:^2)/n)
						)
				
				lmtest1=n*h2'*invsym(F22-F12'*invsym(F11)*F12)*h2 

				printf("LM test1 \n")
				lmtest1
				
				st_numscalar("lmtest1",lmtest1 )			
				
		}
        end
**************************************************************