Adapting 'xtbcfe' for Dynamic Panel Count Data models (Poisson/Negative Binomial)

[Sophie Janssen]

Dear Statalist members,

I am relatively new to Stata and kindly seek your guidance.

I am working on a study where my dependent variable, Incident_Count, is a count variable representing monthly incident counts across 15 countries (N=15) over 144 months (T=144). The panel dataset is balanced but contains many months with zero observations. My data exhibit challenges such as autocorrelation, cross-sectional dependence, heteroskedasticity, and significant overdispersion.

The dependent variable seems very persistent, thus I would like to extend the estimation of panel count models (using xtpoisson and xtnbreg) with dynamic panel count data models. For clarity, the xtpoisson model is specified as:

xtpoisson Incident_Count GPRC Log_GPI Log_GDP Log_GDP2 Log_SIS Log_SIS2 i.month, fe vce(cluster country)

Additional note on the above: because observable nonlinear relationships with regressors Log_GDP and Log_SIS, their squared terms are also included. I still doubt to add these squared terms however, due to multicollinearity issues...

With regard to the dynamic modeling: I already have tried implementing a GMM-based Arellano-Bond estimator for count data using the following code (found this on Statalist, and note: does not account for time effects):

program gmm_poipre
    version 11
    syntax varlist if, at(name) myrhs(varlist) mylhs(varlist) myidvar(varlist)

    quietly {
        tempvar mu mubar ybar 
        gen double `mu' = 0 `if'

        * Compute the linear predictor
        local j = 1
        foreach var of varlist `myrhs' {
            replace `mu' = `mu' + `var'*`at'[1,`j'] `if'
            local j = `j' + 1
        }

        * Apply exponential function to get predicted Poisson mean
        replace `mu' = exp(`mu')

        * Replace the original variable list with the transformed dependent variable
        replace `varlist' = L.`mylhs' - L.`mu'*`mylhs'/`mu' `if'
    }
end

local myrhs "GPRC Log_GPI Log_GDP Log_GDP2 Log_SIS Log_SIS2 L1_Incident_Count"

gmm gmm_poipre, nequations(1) onestep vce(cluster country) myrhs(`myrhs') ///
    mylhs(Incident_Count) myidvar(country) ///
    parameters(Incident_Count:GPRC Incident_Count:Log_GPI Incident_Count:Log_GDP Incident_Count:Log_SIS Incident_Count:L1_Incident_Count) ///
    instruments(L.GPRC L.Log_GPI L.Log_GDP L.Log_GDP2 L.Log_SIS L.Log_SIS2 L2_Incident_Count)

While the above implementation works for the Arellano-Bond estimator, I am exploring whether I can adapt the xtbcfe package to handle the sensitivity problem to identifying valid instruments. Specifically, I am exploring whether the bias-corrected fixed effects estimator in xtbcfe could be adapted to handle dynamic Poisson (or Negative Binomial) models. My primary goal is to avoid instrument specification while robustly addressing issues such as autocorrelation and heteroskedasticity.

The current xtbcfe package assumes a linear framework. Would it make sense to manually replicate its bias correction logic for a Poisson/NB framework, or would this contradict its underlying assumptions? Alternatively, are there better packages and/or approaches in Stata for estimating dynamic panel count data models from which I might not have heard yet?

I would very much like to hear your thoughts regarding this topic. Any advice, code examples, or references would be much appreciated as I refine my methodology.

Thank you in advance for your help.

Best regards,

Sophie