Title

xtspj: Split-panel jackknife estimation for fixed-effect panel data models

Description

xtspj implements the first-order split-panel jackknife (also termed the half-panel jackknife) of Dhaene and Jochmans (2015)
for possibly nonlinear panel data models with fixed effects. These models produce a potentially inconsistent maximum likelihood (ML)
estimate of the common parameters theta (the parameters common to all cross-sectional units and time periods).
The data, (Y_it,X_it), are indexed with i referring to the cross-sectional unit (or individual) and t to the time period.
For each individual there is an individual-specific but time-invariant parameter a_i.

xtspj accepts both balanced and unbalanced data sets, but there must be no gaps in the time series. Missing data are excluded
automatically. The data must be xtset with both panelvar and timevar. If the regressand is a lagged
variable, it has to be generated manually prior to estimation. When the regressors contain lagged values of the regressand
up to p lags, the model becomes autoregressive of order p and is estimated conditional on the first p observations.
The maximizaton of the likelihood function is carried out by the Newton-Raphson algorithm.

xtspj may also be used as a fast (but less flexible) alternative to ml for maximum likelihood
estimation of panel data models with fixed effects.

Syntax

Split-panel jackknife estimation of probit, logit, linear, negative binomial, Poisson, exponential, gamma,
and Weibull regression models with fixed effects:

xtspj depvar [indepvars] [if] [in] , required_options
[other_options]

Split-panel jackknife estimation of linear models (using closed-form expressions):

xtspj depvar [indepvars] [if] [in] , model(regress) required_options
[other_options] [regress_options]

Split-panel jackknife estimation of user-written fixed-effect models:

xtspj eq [eq ...] [if] [in] , required_options
[other_options]

where eq is similar to eq of ml, i.e.,

( [eqname:] [varlist_y] = [indepvars] [, nocons] )

Note that "=" is required even if an equation has an empty right-hand side, i.e., if there is only a constant term.
In the first eq, which is always without constant term due to the presence of fixed effects, nocons is
implicitly assumed, without specifying it. When a factor variable appears as in more than one ,
it must appear with a different name in each . Hence one has to generate a duplicate of the underlying variable
for each additional where the factor variable appears.

indepvars, but not depvar or varlist_y, may contain time-series operators and factor variables;
see tsvarlist and fvvarlist.

regress_options are options of the regress command. They do not affect the estimation
of the model parameters but do affect the estimated covariance matrix of the estimator. See	Remarks below.

required_options

model(string)
specifies the type of regression model to be estimated: probit (model(probit)), logit (model(logit)),
linear (model(linear)), negative binomial (model(negbin)),
Poisson (model(poisson)), exponential (model(exponential)),
gamma (model(gamma)), Weibull (model(weibull)),
or some other, user-written model.
For user-written models, see User Written Models.
The linear model can also be estimated by specifying model(regress); see Remarks below.

method(string)
takes none, like, and parm for, respectively, no correction
(the usual ML estimate theta_hat), the split-panel jackknife estimate based on the jackknifed
log-likelihood (theta_dot), and the split-panel jackknife ML estimate (theta_tilde).

other_options

Estimation
level(#)
sets the confidence level. The default is level(95).

ltol(#)
sets the tolerance level for changes in the objective function value. When the difference
between the objective function values of the current and the previous iteration, divided
by the absolute value of the current objective function value, is positive and smaller
than level(#), the algorithm stops and reports that convergence has been achieved.
The default is ltol(1e-4).

ptol(#)
sets the tolerance level for changes in the parameter values. The algorithm stops and
reports convergence when the change in the parameter is smaller than ptol(#).
The change is computed as the absolute difference between the parameter values in the current and
the previous iteration. When the parameter is a vector, the maximum element of the vector
of absolute differences is taken. The default is ptol(1e-4).

maxiter(#)
sets the maximum number of iterations the algorithm is allowed to use. The default is maxiter(100).

diagnosis
specifies that a simple diagnostic algorithm be invoked when the Newton-Raphson
algorithm gives an updated parameter vector that does not improve the objective function value
(see Remarks below). The diagnostic algorithm is slow and disabled by default.

Results
verbose
specifies whether the iteration log of the maximization (i.e., the objective function values)
and extra notes (e.g., about omitted individuals) should be printed on the screen. The printed
content can be lengthy. It is disabled by default.

alpha(newvarname)
specifies the name of the variable to be created to store the estimates of the fixed effects
(the same value for all observations corresponding to i). When method(none) is requested,
newvarname contains the ML estimate of the fixed effects. When method(like) or
method(parm) is requested, newvarname is the ML estimates of the fixed effects
with theta held fixed in the likelihood function at, respectively,
theta_dot and theta_tilde.

Remarks

Covariance matrix

The covariance matrix of the estimated theta is computed as follows. When no correction
is requested, the usual oim covariance matrix based on the Hessian of the concentrated
log-likelihood is computed (see vce_option). When method(like) or
method(parm) is requested, the covariance matrix is obtained from the same Hessian,
but now evaluated at the corresponding estimate, theta_dot or theta_tilde.
This requires maximizing the uncorrected log-likelihood with respect to the fixed-effect
parameters while theta is kept fixed at theta_dot or theta_tilde.
For this maximization, the options ltol(#), ptol(#), and maxiter(#) are
also effective.

Option diagnosis

When diagnosis is turned on, a diagnostic algorithm is invoked every time the current
objective function value in the Newton-Raphson algorithm is smaller than that in the previous
iteration. Given the gradient vector and Hessian matrix obtained from the previous iteration,
the algorithm successively reduces the step size, used to update the parameter vector, from 1
to 0.5, 0.5^2, ..., possibly down to a minimum value of 0.5^10. The objective function is
evaluated at the parameter vector updated with the reduced step size. When the objective
function value improves on that from the previous iteration, the diagnostic algorithm stops
and maximization continues from the currently updated parameter vector. Otherwise it keeps
reducings the step size unless the minimum step size is already reached, in which case the
diagnostic algorithm reports that no improvement could be achieved by reducing the step size,
it stops, and maximization continues from the parameter vector obtained before the diagnostic
algorithm was called.

Check of multicollinearity

By default, the program runs a multicollinearity check eq-wise, for each half-panel; i.e., on the set of
regressors of each equation, with observations corresponding to each half-panel. The check is performed using the
command _rmcoll, and the coefficients of the omitted variables are kept in the regression,
but fixed at 0 and assigned a variance of 0.

When fvvarlist is present, the indicators corresponding to the base level
(marked with b.) and those marked with o. by fvexpand are omitted from the estimation
in the same way as above.

Noninformative individuals

For the probit and logit regression models, a check is performed for each individual i in each of
the subpanels to ensure that the regressand varies across t. When the regressand do not vary across t,
the individual is omitted from the estimation altogether, because it is not informative. A similar check is also
performed in, e.g., probit. Similarly, for the Poisson and negative binomial regression models,
individuals for which the regressand is zero for all t are uninformative and are therefore dropped.
For the exponential, gamma, and Weibull regression models, individuals for which the regressand is nonpositive are dropped;
this is similar to streg.

model(linear) vs. model(regress)

When (model(linear)) is specified, the linear model estimates are computed by maximization of a Gaussian log-likelihood
using the Newton-Raphson algorithm. When (model(regress)) is specified, the linear model estimates are computed
using least-squares-based closed-form expressions. Apart from small numerical differences, the two ways of
computing the estimates gives identical results.

When model(regress) is specified, no estimate of the error variance is reported, as is often the case for least-squares-based procedures.
When model(linear) is specified, the estimated error variance is reported, as is standard for ML-based procedures.
For the implementation, model(regress) computes within-group estimates. If method(none) is specified,
the resulting coefficient estimates are identical to those obtained from xtreg, fe.

Users are advised to use model(regress) when the error variance is not of interest. Since model(regress) uses
closed-form expressions to compute the estimates, it is faster and more accurate than model(linear).

When model(regress) is specified, the options ltol(#), ptol(#), maxiter(#), diagonsis, and verbose
are ineffective (i.e., they are ignored). The option alpha(newvarname) remains effective.

When model(regress) is specified, the estimated covariance matrix is obtained from the regress command by estimating the
model using the same dataset but after the within transformation. Here, when regress is called, the nocons option
is activated and cannot be overridden. However, other options for the regress command can be specified via
regress_options in order to control the estimation of the covariance matrix. In particular, robust can be specified
to obtain robust covariance matrix estimates.

Model parameterizations

The parameterizations for the probit, logit, and Poisson regression models are standard. For the negative binomial model, the
parameterization is the same as for nbreg.

For the gamma regression model the log-likelihood is specified as

log[f(Y_it|X_it; theta, a_i, kappa)] = - log[gamma(kappa)] - kappa * log(u_it) - (kappa - 1) * log(Y_it) - Y_it / u_it

with

u_it = exp(X'_it * theta + a_i)

where a_i is the ith fixed effect and kappa is a scalar shape
parameter. This parameterization is the same as for glm, family(gamma) link(log) for the
regression slope coefficients; we also report the shape parameter estimate.

For the exponential regression model the log-likelihood is specified as

log[f(Y_it|X_it; theta, a_i)] = log(1/u_it) - Y_it / u_it,

which is the same as for glm, family(gamma) link(log).

For the Weibull regression model the log-likelihood is specified as

log[f(Y_it|X_it; theta, a_i, k)] = log(kappa) + (kappa-1) * log(Y_it) + log(u_it) - u_it * Y_it ^ kappa

where kappa=exp(k).

Example

Data setup

. webuse chicken, clear

. gen worker=.

. by restaurant: replace worker=_n

. xtset restaurant worker

. label define sex 0 "female" 1 "male"

. label values gender sex

Uncorrected fixed-effect probit ML

. xtspj complain age i.gender, model(probit) method(none)

Jackknifed log-likelihood

. xtspj complain age i.gender, model(probit) method(like)

Jackknifed ML

. xtspj complain age i.gender, model(probit) method(parm)

Stored Results

xtspj stores the following in e():

Scalars

e(N) number of observations

e(empty) contains 1 if there are no covariates in the model

e(converged) contains 1 if the maximization converged

e(ll) log-likelihood value (only for method(none) or method(like) when model(regress) is not specified)

Matrices

e(b) coefficient vector

e(V) covariance matrix

Functions

e(sample) estimation sample

Macros

e(model) value of model(string)

e(method) value of method(string)

e(cmd) xtspj

e(cmdline) command as typed

e(vce) oim or vcetype specified in vce() (for model(regress) only)

e(vcetype) title used to label Std. Err. (for model(regress) only)

e(title) title of table of results

e(properties) b V

e(depvar) the regressands or the list of regressands (when user-written model with several eq were specified)

User-written Models

Parameterization of the density

To use xtspj with a user-written model, the density has to be of the form

f(Y_it|X_it; u_1_it, ..., u_M_it)

with

u_1_it = X'_1_it * theta_1 + a_i,

u_m_it = X'_m_it * theta_m

where a_i is the fixed effect of the i-th individual and m=2,...,M for an arbitrary
number of M>=1 equations depending on covariates X_m_it. The regressands Y_it can
be a vector, at most of dimension M. The first equation, u_1_it contains an additive
fixed effect. Nevertheless, this design is general in the sense that, when the fixed effects enter
non-additively, one simply needs to keep the first equation free of covariates. It is also possible
to set X_m_it=1 for one or several m>1; then the corresponding u_m_it=theta_m are
simply parameters entering the model.

User-written evaluator

The likelihood function has to be specified as a mata class containing two
member functions ::Evaluate() and ::Check(). The template for the class is

1 | mata
2 | class xtspjYourModel extends xtspjModel {
3 | public:
4 | void Evaluate(), Check()
5 | }
6 | void function xtspjYourModel::Evaluate(
7 | real matrix Y,
8 | real matrix XB,
9 | real colvector LogLikelihood,
10| real matrix Gradient,
11| pointer(real colvector) matrix Hessian) {
12| // compute the log-likelihood, gradient, and Hessian.
13| }
14| void function xtspjYourModel::Check(
15| real matrix Data,
16| real scalar Keep) {
17| // check data.
18| }
19| end

The user needs to supply the relevant mata code in the lines 12 and 17 as discussed below, and change
YourModel (lines 2, 6, and 14) to a desired name, e.g., MyModel. When xtspj is executed,
this name must be specified in the option model(string), e.g., model(MyModel), so that
xtspj calls the user-written specification of the log-likelihood function.

void function ::Evaluate()

This function computes the log-likelihood, gradient, and Hessian for a given set of observations (i,t).
The function is called repeatedly and separately for each i and various sets . The arguments
of ::Evaluate() are:

- real matrix Y,
the matrix of regressands. The k-th column of Y corresponds to the k-th regressand defined
by an eq, and each row corresponds to an observation (i,t) with t
belonging to a specific set .

- real matrix XB,
the matrix of linear indices (u_m_it for 1<=m<=M). The m-th column corresponds
to the m-th equation and each row corresponds to an observation (i,t).

- real colvector LogLikelihood,
the column vector of log-likelihood values. Each row contains the log-likelihood value corresponding to
an observation (i,t) evaluated at the corresponding row of XB and Y; i.e.,

log[f(Y_it|X_it; u_1_it, ..., u_M_it)]

where t belongs to a given .

- real matrix Gradient,
the matrix of scores. The m-th column of Gradient corresponds to the log-likelihood
derivative with respect to the m-th equation, evaluated at XB and Y; i.e.,

d log[f(Y_it|X_it; u_1_it, ..., u_M_it)]
----------------------------------------
d u_m_it

where t belongs to a given .

- pointer(real colvector) matrix Hessian,
the pointer matrix pointing to the Hessian. The matrix must be
made symmetric. The (m,n)-th element of Hessian is a pointer to the
column vector of second-order log-likelihood derivatives with respect to the m-th and n-th
equation, evaluated at XB and Y; i.e.,

d^2 log[f(Y_it|X_it; u_1_it, ..., u_M_it)]
------------------------------------------
d u_m_it d u_n_it

where t belongs to a given .

Here the arguments real matrix gradient and pointer(real colvector) matrix Hessian are optional.
If they are not changed inside ::Evaluate(), xtspj uses a numerical differentiation algorithm
to compute the derivatives.

Below is an example of ::Evaluate() for the probit model:

1 | void function xtspjprobit::Evaluate(
2 | real matrix Y,
3 | real matrix XB,
4 | real colvector LogLikelihood,
5 | real matrix Gradient,
6 | pointer(real colvector) matrix Hessian) {
7 | q=(Y:*2:-1)
8 | LogLikelihood=lnnormal(q:*XB)
9 | Gradient=q:*(normalden(q:*XB):/normal(q:*XB))
10| Hessian=&(-Gradient:*(XB+Gradient))
11| }

Here the first point to be noted is that YourModel is replaced with probit. The LogLikelihood,
Gradient, and Hessian are computed in lines 8 to 10. Lines 9 and 10 can be omitted for a numerical
differentiation algorithm to take over.

When a model contains more than one equation (M>1), the variable Gradient and Hessian may
need to be preallocated to expedite the computation. To do so, simply insert

1 | Gradient=J(rows(Data),cols(Data),.)
2 | Hessian=J(rows(Data),cols(Data),NULL)

between lines 6 and 7. In addition, one should avoid reusing variables in the calculation of the Hessian;
i.e., code such as

1 | H=...
2 | Hessian[1,1]=&H
3 | H=...
4 | Hessian[1,2]=&H

should be replaced by

1 | H11=...
2 | Hessian[1,1]=&H11
3 | H12=...
4 | Hessian[1,2]=&H12

void function ::Check()

If a check is not desired, one can simply replace line 17 in this code with
Keep=1. To perform a model-specific check of the data, one has to supply suitable code to ::Check().
Similar to ::Evaluate(), the function is called repeatedly and separately for each
i and various sets . The arguments are:

- real matrix Data,
the data to be checked. The columns of Data correspond to the regressands of the model, followed
by the regressors, in the same order as specified by the equations. The rows correspond to the
observations (i,t) with t belonging to a given set .

- real scalar Keep,
the main output. Keep is set to 1 if Data passes the check, and to 0 otherwise.

Below is an example of ::Check() for the probit model:

1 | void function xtspjprobit::Check(
2 | real matrix Data,
3 | real scalar Keep) {
4 | if (sum(Data[.,1])==0 || sum(Data[.,1])==rows(Data)) {
5 | Keep=0
6 | } else {
7 | Keep=1
8 | }
9 | }

References

Dhaene, G., and K. Jochmans. 2015. Split-panel jackknife estimation of fixed-effect models. The Review of Economic Studies
82(3): 991-1030.