/* 'xtpsse' runs a conditional fixed-effects Poisson Panel Regression, computes sandwich and spatial SE, and tests for time-invariant spatial dependence according to Bertanha and Moser (2016) For help, see xtpsse.sthlp */ /* Marinho Bertanha, Version 05/01/2017 I'd appreciate your feedback: mbertanha@nd.edu Some of the structure here was inspired by Tim Simcoe's xtpqml.ado command */ program xtpsse, eclass version 11 syntax varlist(numeric) [if] [in], COORDinates(varlist) CUToffs(numlist >0) [ I(varname num) T(varname num) TEst(integer -1)] ** Make sure user specifies cross-section ID variable AND time index variable if length("`i'") == 0 { local i = "`r(panelvar)'" if length("`i'") == 0 { di as error "You must specify a cross-section ID variable. Use i(varname) option or 'xtset' command." exit } } if length("`t'") == 0 { local t = "`r(timevar)'" if length("`t'") == 0 { di as error "You must specify a time index variable. Use t(varname) option or 'xtset' command." exit } } ** Dimension of both coordinates and cutoffs must be equal if wordcount("`coordinates'") != wordcount("`cutoffs'") { di as error "Number of coordinate variables is different than the number of cutoffs." exit } scalar K2 = `test' if K2==-1 { scalar K2=wordcount("`varlist'")-1 } xtpoisson `varlist' `if' `in', fe i(`i') di as text "Computing Standard Errors..." tempvar n_i xb_hat_it mu_hat_it sum_mu_i p_it u_it crossec crossec2 idsamp tempname b quietly{ matrix b = e(b) local varlist : colnames b bysort `i' : egen `n_i' = sum(`e(depvar)') if e(sample) predict `xb_hat_it', xb gen `mu_hat_it' = exp(`xb_hat_it') if e(sample) by `i' : egen `sum_mu_i' = sum(`mu_hat_it') if e(sample) gen `p_it' = `mu_hat_it' / `sum_mu_i' if e(sample) gen `u_it' = `e(depvar)' - `p_it'*`n_i' if e(sample) by `i' : egen `crossec' = min(`t') if e(sample) foreach v of var `coordinates' { tempvar cd_`v' gen `cd_`v'' = `v' replace `cd_`v'' = . if `t'!=`crossec' local cdnames `cdnames' `cd_`v'' } foreach v of var `varlist' { tempvar sum_mu_`v' score_i_`v' by `i' : egen `sum_mu_`v'' = sum(`v'*`mu_hat_it') if e(sample) by `i': egen `score_i_`v'' = sum((`v' - (`sum_mu_`v''/`sum_mu_i')) * `u_it') if e(sample) replace `score_i_`v'' = . if `t'!=`crossec' local scorename `scorename' `score_i_`v'' } ** units that are not dropped in estimation by `i' : egen `crossec2' = min(`t') gen `idsamp' = `i' if e(sample)==0 & `t'==`crossec2' } mata: varcovar() ereturn repost V = V_SW di as result "Robust Standard Errors (Sandwich Formula)" ereturn display ereturn repost V = V_SP di as result "Spatial Standard Errors" ereturn display if K2>0 { di as text "`msg1'" di as text "`msg2'" di as result "Null hypothesis of time-invariant spatial dependence" di as text "test-statistic: " as result T_hat di as text "p-value (Chi2 - " Kstar " df): " as result pval } end mata: void varcovar() { K2=st_numscalar("K2") scorename=st_local("scorename") nameid=st_local("idsamp") A_hat = st_matrix("e(V)") idsamp = st_data(.,nameid,0) scores=st_data(.,scorename,0) nobs=rows(scores) nvars=cols(scores) cutoffs=strtoreal(tokens(st_local("cutoffs"))) C_hat = J(nvars,nvars,0) coordinates = st_local("cdnames") coord=st_data(.,coordinates,0) nc=2 //***********************fits spatial coord into regular lattice coordinates dist=J(nobs*(nobs-1)/2,1,0) l=0 for (i=1; i<=(nobs-1); i++) { for (j=(i+1); j<=nobs; j++) { temp=((coord[i,1]-coord[j,1])^2+(coord[i,2]-coord[j,2])^2)^.5 if (temp>0) { l=l+1 dist[l,1]=temp } } } dist=dist[1::l,.] dmin=0.95*min(dist) dstar=(sqrt(2)/2)*dmin coord2=J(nobs,2,0) spmin=(min(coord[.,1]),min(coord[.,2])) for (i=1; i<=nobs; i++) { coord2[i,.]=(floor( diag((1/dstar,1/dstar))*(coord[i,.]-spmin)')+(1\1))' } cutoffs2=round( diag((1/dstar,1/dstar))* cutoffs' )' //*****************************************computes covariance matrix C_hat for (i=1; i<=nobs; i++) { for (j=1; j<=nobs; j++) { temp=abs(coord2[i,.]-coord2[j,.]) if (temp<=cutoffs2) { prod=1 if (min(cutoffs2)>0) { for (k=1; k<=nc; k++) { prod=prod*(1-temp[1,k]/cutoffs2[1,k]) } } C_hat = C_hat + prod*scores[i,.]'*scores[j,.] } } } V_SP=A_hat*C_hat*A_hat st_matrix("V_SP",V_SP) B_hat = J(nvars,nvars,0) for (i=1; i<=nobs; i++) { B_hat = B_hat + scores[i,.]'*scores[i,.] } V_SW=A_hat*B_hat*A_hat st_matrix("V_SW",V_SW) //********************************************computes test-statistic if (K2>0) { st_local("msg1"," ") st_local("msg2"," ") Kstar=K2*(K2+1)/2 spmomg=J(nobs,K2,0) for (i=1; i<=nobs; i++) { Nl=(2*cutoffs2[1]+1)*(2*cutoffs2[2]+1)-1 for (j=1; j<=nobs; j++) { temp=abs(coord2[i,.]-coord2[j,.]) if (i~=j && temp[1]<=cutoffs2[1] && temp[2]<=cutoffs2[2]) { //Ni=Ni+1 spmomg[i,.]=spmomg[i,.]+scores[j,1::K2] } } spmomg[i,.]=(1/Nl)*spmomg[i,.] } Z=J(nobs,1,NULL) A=J(Kstar,Kstar,0) for (i=1; i<=nobs; i++) { Z[i]=&J(Kstar,K2,0) l=1 for (k=1; k<=K2; k++) { (*Z[i])[l::(l+k-1),k]=spmomg[i,1::k]' l=l+k } A=A+(*Z[i])*(*Z[i])' } Gamma=(-1/nobs)*A Omega = J(Kstar,Kstar,0) for (i=1; i<=nobs; i++) { for (j=1; j<=nobs; j++) { temp=abs(coord2[i,.]-coord2[j,.]) if (temp<=2*cutoffs2) { prod=1 if (min(cutoffs2)>0) { for (k=1; k<=nc; k++) { prod=prod*(1-temp[1,k]/(2*cutoffs2[1,k])) } } Omega=Omega+(1/nobs)*prod*( (*Z[i]) * scores[i,1::K2]' )*( scores[j,1::K2]*(*Z[j])') } } } if (min((rank(Gamma),rank(Omega)))0) { for (k=1; k<=nc; k++) { prod=prod*(1-temp[1,k]/(2*cutoffs2[1,k])) } } Omega=Omega+(1/nobs)*prod*( (*Z[i]) * scores[i,1::K2]' )*( scores[j,1::K2]*(*Z[j])') } } } } } W=luinv(Gamma)*Omega*luinv(Gamma) B=J(Kstar,1,0); for (i=1; i<=nobs; i++) { B=B+(*Z[i])*scores[i,1::K2]' } Theta=luinv(A)*B T_hat = nobs*Theta'*luinv(W)*Theta pval=1-chi2(Kstar,T_hat) st_numscalar("T_hat",T_hat) st_numscalar("Kstar",Kstar) st_numscalar("pval",pval) st_numscalar("K2",K2) st_numscalar("nvars",nvars) if (K2