*xtserialpml V1.0.1   1-Apr-2019
*Koen Jochmans University of Cambdridge and Vincenzo Verardi, FNRS-UNamur

program define xtserialpm, rclass

    version 10.1

    if replay()& "`e(cmd)'"=="xtserialpm" {
        ereturn  display
        exit
    }

    syntax varlist [if] [in] , [Center]

    tempname touse gid2

    mark `touse' `if' `in'
    markout `touse' `varlist'

    capture tsset

    capture local ivar "`r(panelvar)'"
    if "`ivar'"=="" {
        di as err "must tsset data and specify panelvar"
        exit 459
    }
    capture local tvar "`r(timevar)'"
    if "`tvar'" == "" {
        di as err "must tsset data and specify timevar"
        exit 459
    }

    preserve
    qui keep if `touse'
    tsfill, full
    qui xtdescribe
    if r(min)!=r(max) {
        di "Panel is not balanced"
    }

    local dv: word 1 of `varlist'
    local exp: list varlist -dv
    egen `gid2'=group(`ivar')

    tempvar center0

    if "`center'"=="" {
		gen `center0'=0
    }

    else {
		gen `center0'=1
    }

    mata: test("`varlist'","`gid2'","`touse'","`center0'")
    local q=res_MAT_xxx[1,3]
    local stat=res_MAT_xxx[1,1]
    local p=res_MAT_xxx[1,2]

    di as txt ""
    di as txt "Jochmans portmanteau test for within-group correlation in panel data. "
    di as txt "H0: no within-group correlation"
    di as txt _col(5) "Chi-sq(" %3.0f `q'  ")   =  "     ///
        as res %10.3f `stat'
    di as txt _col(5) "Prob > Chi-sq =  " as res %11.4f `p'

    ret scalar stat = `stat'
    ret scalar df   = `q'
    ret scalar p    = `p'


end

mata

    void test(string scalar varlist, string scalar gid2, string scalar touse, string scalar center0)
    {
        st_view(X0=.,.,tokens(varlist))
        st_view(gid=.,.,tokens(gid2))
        st_view(center=.,.,tokens(center0))

        X0=(X0,J(rows(X0),1,1))


        v=J(1,cols(X0),1)
        v[1,1]=0
        Y=X0[,1]
        dX=select(X0,v)
        Y=rowshape(Y,max(gid))
        K=cols(dX)
        X=asarray_create("real",1)
        for(k=1;k<=K;k++) {
            asarray(X,k,rowshape(dX[,k],max(gid)))
        }

        G=rows(Y)
        N=cols(Y)
        df = N*(N-1)/2-1
        Dy=(Y :!= .)

        Dx=J(rows(Y),cols(Y),1)
        for(k=1;k<=K;k++) {
			Dx=Dx:*(asarray(X,k) :!= .)
        }

        D=Dy:*Dx
        Ng= rowsum(D)
        Ig = (Ng:>1)
        Y=select(Y,Ig)
        D=select(D,Ig)
        Ng=select(Ng,Ig)

        for(k=1;k<=K;k++) {
            Xk=asarray(X,k)
            Xk=select(Xk,Ig)
            Xk=mm_cond(D:==0,-999, Xk)
            asarray(X,k,Xk)
        }

        G=rows(Y)
        N=cols(Y)

        Y=mm_cond(D:==0,-999, Y)
        dY=D:*Y:-(rowsum(D:*Y):/Ng)*J(1,N,1)

        dX=asarray_create("real",1)
        for(k=1;k<=K;k++) {
			Xk=D:*asarray(X,k)-(rowsum(D:*asarray(X,k)):/Ng)*J(1,N,1)
            Xk=mm_cond(D:==0,-999, Xk)
            asarray(dX,k,Xk)
        }

        dY=mm_cond(D:==0,-999, dY)
        XY=J(K,1,0)
        XX=J(K,K,0)

        for(k1=1;k1<=K;k1++) {
            XY[k1, 1] = mean(rowsum(D:*(asarray(dX,k1):*dY)))
            for(k2=1;k2<=K;k2++) {
                XX[k1,k2] = mean(rowsum(D:*(asarray(dX,k1):*asarray(dX,k2))))
            }
        }

        b=invsym(XX)*XY

	E=dY

for(k=1;k<=K;k++) {
E=E-asarray(dX,k)*b[k]
}

e=E
e=mm_cond(D:==0,0,e)
A=J(G,K,0)

for(k=1;k<=K;k++) {
A[,k]=rowsum(D:*(asarray(dX,k):*E))
omega=(invsym(XX)*A')'
}

asyvar = (omega'*omega)/G^2
se = sqrt(diag(asyvar))

E=Y
for(k=1;k<=K;k++) {
E=E-asarray(X,k)*b[k]
}

DE=E[,2..cols(E)]-E[,1..cols(E)-1]
DD =D[,1..cols(D)-1]:*D[,2..cols(D)]


s=J(1,0,0)
for(t=1;t<=N-2;t++) {
s=(s,range(1,t,1)')
}

ss=J(0,1,0)
for(t=2;t<=N-1;t++) {
ss=(ss\t*J(t-1,1,1))
}

ZDE = (((D[,s]:*E[,s]):*(DD[,ss]:*DE[,ss])), (D[ ,3..N]:*E[,3..N]:*(DD[,1..cols(DD)-1]:*DE[,1..cols(DE)-1])))'
ZDX = J(df,K,0)

for(k=1;k<=K;k++) {
DX = asarray(X,k)[,2..cols(asarray(X,k))]-asarray(X,k)[,1..cols(asarray(X,k))-1]
ZDX[,k] = mean((((D[ ,s]:*E[,s]):*(DD[,ss]:*DX[,ss])) , (D[,3..N]:*E[,3..N]:*(DD[,1..cols(DD)-1]:*DX[,1..cols(DD)-1]))))'
}


if (cols(asarray(dX,1))>1) {
vg = ZDE - ZDX*omega'
}

else {
vg=ZDE
}

v=mean(vg')'
V1 = (vg*vg')/G
V2 = ((vg:-v)*(vg:-v)')/G

if (max(center)==1) {
V=V2
}

else {
V=V1
}

st = G*v'*invsym(V)*v
pv = 1-chi2(df,st)
res=(st,pv,df)
st_matrix("res_MAT_xxx",res)
}
end