经管之家App
让优质教育人人可得
立即打开
Title
[XT] xtprobit -- Random-effects and population-averaged probit models
Syntax
Random-effects (RE) model
xtprobit depvar [indepvars] [if] [in] [weight] [, re RE_options]
Population-averaged (PA) model
xtprobit depvar [indepvars] [if] [in] [weight] , pa [PA_options]
RE_options description
--------------------------------------------------------------------------------------------------------------
Model
noconstant suppress constant term
re use random-effects estimator; the default
offset(varname) include varname in model with coefficient constrained to 1
constraints(constraints) apply specified linear constraints
collinear keep collinear variables
SE
vce(vcetype) vcetype may be oim, bootstrap, or jackknife
Reporting
level(#) set confidence level; default is level(95)
noskip perform likelihood-ratio test
Int opts (RE)
intmethod(intmethod) integration method; intmethod may be mvaghermite, aghermite, or ghermite; default
is intmethod(mvaghermite)
intpoints(#) use # quadrature points; default is intpoints(12)
Max options
maximize_options control the maximization process; seldom used
--------------------------------------------------------------------------------------------------------------
PA_options description
--------------------------------------------------------------------------------------------------------------
Model
noconstant suppress constant term
pa use population-averaged estimator
offset(varname) include varname in model with coefficient constrained to 1
Correlation
corr(correlation) within-group correlation structure
force estimate even if observations unequally spaced in time
SE/Robust
vce(vcetype) vcetype may be conventional, robust, bootstrap, or jackknife
nmp use divisor N-P instead of the default N
scale(parm) override the default scale parameter; parm may be x2, dev, phi, or #
Reporting
level(#) set confidence level; default is level(95)
Opt options
optimize_options control the optimization process; seldom used
--------------------------------------------------------------------------------------------------------------
correlation description
--------------------------------------------------------------------------------------------------------------
exchangeable exchangeable
independent exchangeable
unstructured unstructured
fixed matname user-specified
ar # autoregressive of order #
stationary # stationary of order #
nonstationary # nonstationary of order #
--------------------------------------------------------------------------------------------------------------
A panel variable must be specified. For xtprobit, pa, correlation structures other than exchangeable and
independent require that a time variable also be specified. Use xtset.
depvar and indepvars may contain time-series operators; see tsvarlist.
by, statsby, and xi are allowed; see prefix.
iweights, fweights, and pweights are allowed for the population-averaged model, and iweights are allowed in
the random-effects model; see weight. Weights must be constant within panel.
See [XT] xtprobit postestimation for features available after estimation.
Description
xtprobit fits random-effects and population-averaged probit models. There is no command for a conditional
fixed-effects model, as there does not exist a sufficient statistic allowing the fixed effects to be
conditioned out of the likelihood. Unconditional fixed-effects probit models may be fitted with probit
command with indicator variables for the panels. The appropriate indicator variables can be generated using
tabulate or xi. However, unconditional fixed-effects estimates are biased.
By default, the population-averaged model is an equal-correlation model; xtprobit assumes corr(exchangeable).
See [XT] xtgee for information on how to fit other population-averaged models.
See logistic estimation commands for a list of related estimation commands.
Options for RE model
+-------+
----+ Model +-------------------------------------------------------------------------------------------------
noconstant; see [XT] estimation options.
re requests the random-effects estimator. re is the default if neither re not pa is specified.
offset(varname), constraints(constraints), collinear; see [XT] estimation options.
+----+
----+ SE +----------------------------------------------------------------------------------------------------
vce(vcetype) specifies the type of standard error reported, which includes types that are derived from
asymptotic theory and that use bootstrap or jackknife methods; see [XT] vce_options.
+-----------+
----+ Reporting +---------------------------------------------------------------------------------------------
level(#), noskip; see [XT] estimation options.
+---------------+
----+ Int opts (RE) +-----------------------------------------------------------------------------------------
intmethod(intmethod), intpoints(#); see [XT] estimation options.
+-------------+
----+ Max options +-------------------------------------------------------------------------------------------
maximize_options: difficult, technique(algorithm_spec), iterate(#), [no]log, trace, gradient, showstep,
hessian, shownrtolerance, tolerance(#), ltolerance(#), gtolerance(#), nrtolerance(#), nonrtolerance,
from(init_specs); see [R] maximize. Some of these options are not available if intmethod(ghermite) is
specified. These options are seldom used.
Options for PA model
+-------+
----+ Model +-------------------------------------------------------------------------------------------------
noconstant; see [XT] estimation options.
pa requests the population-averaged estimator.
offset(varname); see [XT] estimation options.
+-------------+
----+ Correlation +-------------------------------------------------------------------------------------------
corr(correlation), force; see [XT] estimation options.
+-----------+
----+ SE/Robust +---------------------------------------------------------------------------------------------
vce(vcetype) specifies the type of standard error reported, which includes types that are derived from
asymptotic theory, that are robust to some kinds of misspecification, and that use bootstrap or jackknife
methods; see [XT] vce_options.
vce(conventional), the default, uses the conventionally derived variance estimator for generalized
least-squares regression.
nmp, scale(x2|dev|phi|#); see [XT] vce_options.
+-----------+
----+ Reporting +---------------------------------------------------------------------------------------------
level(#); see [XT] estimation options.
+-------------+
----+ Opt options +-------------------------------------------------------------------------------------------
optimize_options control the iterative optimization process. These options are seldom used.
iterate(#) specifies the maximum number of iterations. When the number of iterations equals #, the
optimization stops and presents the current results, even if the convergence tolerance has not been
reached. The default value of iterate() is 100.
tolerance(#) specifies the tolerance for the coefficient vector. When the relative change in the
coefficient vector from one iteration to the next is less than or equal to #, the optimization process is
stopped. tolerance(1e-6) is the default.
nolog suppress the display of the iteration log.
trace specifies that the current estimates should be printed at each iteration.
Technical note
The random-effects model is calculated using quadrature, which is an approximation whose accuracy depends
partially on the number of integration points used. We can use the quadchk command to see if changing the
number of integration points affects the results. If the results change, the quadrature approximation is not
accurate given the number of integration points. Try increasing the number of integration points using the
intpoints() option and again run quadchk. Do not attempt to interpret the results of estimates when the
coefficients reported by quadchk differ substantially. See [XT] quadchk for details and [XT] xtprobit for an
example.
Because the xtprobit, re likelihood function is calculated by Gauss Hermite quadrature, on large problems, the
computations can be slow. Computation time is roughly proportional to the number of points used for the
quadrature.
Examples
Setup
. webuse union
Random-effects model
. xtprobit union age grade not_smsa south southXt
Equal-correlation population-averaged model
. xtprobit union age grade not_smsa south southXt, pa
Equal-correlation population-averaged model with robust variance
. xtprobit union age grade not_smsa south southXt, pa vce(robust)
Title
[XT] xttobit -- Random-effects tobit models
Syntax
xttobit depvar [indepvars] [if] [in] [weight] [, options]
options description
--------------------------------------------------------------------------------------------------------------
Model
noconstant suppress constant term
ll(varname|#) left-censoring variable/limit
ul(varname|#) right-censoring variable/limit
offset(varname) include varname in model with coefficient constrained to 1
constraints(constraints) apply specified linear constraints
collinear keep collinear variables
SE
vce(vcetype) vcetype may be oim, bootstrap, or jackknife
Reporting
level(#) set confidence level; default is level(95)
tobit perform likelihood-ratio test comparing against pooled tobit model
noskip perform likelihood-ratio test
Int opts (RE)
intmethod(intmethod) integration method; intmethod may be mvaghermite, aghermite, or ghermite; default
is intmethod(mvaghermite)
intpoints(#) use # quadrature points; default is intpoints(12)
Max options
maximize_options control the maximization process; seldom used
--------------------------------------------------------------------------------------------------------------
A panel variable must be specified; use xtset.
depvar and indepvars may contain time-series operators; see tsvarlist.
by, statsby, and xi are allowed; see prefix.
iweights are allowed; see weight. Weights must be constant within panel.
See [XT] xttobit postestimation for features available after estimation.
Description
xttobit fits a random-effects tobit models. There is no command for a conditional fixed-effects model, as
there does not exist a sufficient statistic allowing the fixed effects to be conditioned out of the
likelihood. Honore has developed a semiparametric estimator for fixed-effect tobit models. Unconditional
fixed-effects tobit models may be fitted with tobit command with indicator variables for the panels. The
appropriate indicator variables can be generated using tabulate or xi. However, unconditional fixed-effects
estimates are biased.
Options
+-------+
----+ Model +-------------------------------------------------------------------------------------------------
noconstant; see [XT] estimation options.
ll(varname|#) and ul(varname|#) indicate the censoring points. You may specify one or both. ll() indicates
the lower limit for left-censoring. Observations with depvar<ll() are left-censored, observations with
depvar>ul() are right-censored, and remaining observations are not censored. See [R] tobit.
offset(varname), constraints(constraints), collinear; see [XT] estimation options.
+----+
----+ SE +----------------------------------------------------------------------------------------------------
vce(vcetype) specifies the type of standard error reported, which includes types that are derived from
asymptotic theory and that use bootstrap or jackknife methods; see [XT] vce_options.
+-----------+
----+ Reporting +---------------------------------------------------------------------------------------------
level(#); see [XT] estimation options.
tobit specifies that a likelihood-ratio test comparing the random effects model with the pooled (tobit) model
be included in the output.
noskip; see [XT] estimation options.
+---------------+
----+ Int opts (RE) +-----------------------------------------------------------------------------------------
intmethod(intmethod), intpoints(#); see [XT] estimation options.
+-------------+
----+ Max options +-------------------------------------------------------------------------------------------
maximize_options: difficult, technique(algorithm_spec), iterate(#), [no]log, trace, gradient, showstep,
hessian, shownrtolerance, tolerance(#), ltolerance(#), gtolerance(#), nrtolerance(#), nonrtolerance,
from(init_specs); see [R] maximize. Some of these options are not available if intmethod(ghermite) is
specified. These options are seldom used.
Technical note
The random-effects model is calculated using quadrature, which is an approximation whose accuracy depends
partially on the number of integration points used. We can use the quadchk command to see if changing the
number of integration points affects the results. If the results change, the quadrature approximation is not
accurate given the number of integration points. Try increasing the number of integration points using the
intpoints() option and again run quadchk. Do not attempt to interpret the results of estimates when the
coefficients reported by quadchk differ substantially. See [XT] quadchk for details and [XT] xtprobit for an
example.
Because the xttobit likelihood function is calculated by Gauss Hermite quadrature, on large problems, the
computations can be slow. Computation time is roughly proportional to the number of points used for the
quadrature.
Example
Setup
. webuse nlswork3
. xtset idcode
Fit random-effects (RE) tobit model
. xttobit ln_wage union age grade not_smsa south southXt, ul(1.9)
Same as above, but increase the number of integration points from 12 to 25
. xttobit ln_wage union age grade not_smsa south southXt, ul(1.9) intpoints(25)
Same as above, but report likelihood-ratio test comparing RE model with the pooled model
. xttobit ln_wage union age grade not_smsa south southXt, ul(1.9) intpoints(25) tobit
多谢ls的两位~~~~~
不过这些都是existing build in routine.....我想要的是用mata做的自己的小programme.....因为我们老师不喜欢我们用别人编好的东西，总说那个都是black box....光会用是不会知道真正的在后台运行的是什么......
我这个周末试一下，如果没问题的话，会把code传上来给需要的朋友做参考~~~~~~
it's been a month......Sorry for those might need the mata codes...I totally forgot it.....
Here is the code we have obtained(I and my groupmate)....based on the data we had, we also tried the existing build-in routine from stata and the results showed that out mata codes should be correct.....i think the difficult part for this code is how to express the right criterion function.........
cap log c
clear
cd "D:\data\applied_econometrics\Assignment5"
log using re_probit_mata_ml,replace t
/********************************************************************************************/
/* This is an example program how to use the mata routine optimize to */
/* to perfor maximum likelihood. Please read the help function optimize */
/* example: re probit model */
/********************************************************************************************/
mata
mata clear
/* tobit_v0 routine, numerical gradient (1xkk) and numerical hessian,
lnlikelihood contributionscore vector not computed */
void tobit_v0(real scalar todo, real rowvector p, real colvector v, real matrix g, real matrix H)
{
external y, x, d, p, k,n
real colvector nxbeta
beta=p[1..k]
sigma=p[k+1]
v = d:*(-0.5*ln(2*3.1416)*J(n,1,1)-0.5*ln(sigma^2)*J(n,1,1)-(1/(2*sigma^2))*((y-x*beta'):^2))+(J(n,1,1)-d):*ln(J(n,1,1)-normal(x*beta'/sigma))
}
end
use randdata.dta, clear
xtset zper year
gen d=1 if meddol>0
replace d=0 if d==.
*Program*
mata
y = st_data(., "meddol")
d = st_data(., "d")
x1=st_data(.,"logc")
x2=st_data(.,"xage")
n=rows(y)
one=J(n,1,1)
x=one,x1,x2
k=cols(x)
/* read in Person identifier */
zper=st_data(.,"zper")
/* panelsetup is a useful command in order to recognize the panel structure of the data (see help file) */
info = panelsetup(zper, 1)
p2=-30, 5, -50, 400
S = optimize_init()
optimize_init_params(S, p2)
optimize_init_evaluator(S, &tobit_v0())
optimize_init_evaluatortype(S, "v0")
optimize_init_which(S, "max" )
p=optimize(S)
p
loglik=optimize_result_value(S)
st_matrix("lnL_pp",loglik)
st_matrix("nn_p",n)
ihess_pp=optimize_result_V_oim(S)
scores_pp=optimize_result_scores(S)
opg_cluster=J(k+1,k+1,0)
/* calculation clustered standard errors */
for (i=1; i<=rows(info); i++) {
score_i = colsum(panelsubmatrix(scores_pp, i, info))'
opg_cluster=opg_cluster+(score_i*score_i')
}
/*
for (i=1; i<=rows(info); i++) {
score_i = panelsubmatrix(scores_pp, i, info)
x_i=panelsubmatrix(x, i, info)
opg_cluster=opg_cluster+x_i'*(score_i[.,1]*score_i[.,1]')*x_i
}
*/
V_cluster=ihess_pp*opg_cluster*ihess_pp
st_matrix("beta_pp",p)
st_matrix("Vcluster_pp",V_cluster)
end
matrix colnames beta_pp=Constant lc age sigma_u
matrix colnames Vcluster_pp=Constant lc age sigma_u
matrix rownames Vcluster_pp=Constant lc age sigma_u
ereturn post beta_pp Vcluster_pp
eret display
*****************************************************************************************
*****************************************************************************************
******************************Random Effects********************************************
mata
mata clear
void ran_tobit_v0(real scalar todo, real rowvector p, real colvector v, real matrix g, real matrix H)
{
external y, x, k, d, info, nrep
real colvector nxbeta, xbeta,x_i,y_i,iota_T,nnxbeta
real rowvector random
real scalar T, i, t
rseed(10101)
nnxbeta=J(rows(info),1,0)
beta=p[1..k]
sigma_mu=p[k+1]
sigma_anpha=p[k+2]
for (i=1;i<=rows(info);i++) {
random=rnormal(1,nrep,0,1)
x_i=panelsubmatrix(x,i,info)
y_i=panelsubmatrix(y,i,info)
d_i=panelsubmatrix(d,i,info)
n=rows(y)
ones = J(n,1,1)
ones_i=panelsubmatrix(ones,i,info)
iota_T=x_i[.,1]
T=rows(iota_T)
xbeta = ( ( (1/sigma_mu)*normalden((y_i:-sigma_anpha*iota_T*random:-x_i*beta')/sigma_mu) ):^(d_i)) :* ( ( ones_i:-normal( (sigma_anpha*iota_T*random:+x_i*beta')/sigma_mu ) ):^(ones_i-d_i ) )
nnxbeta=mean( exp( colsum( ln(xbeta) ) ) ')
}
v=ln(nnxbeta)
}
end
use randdata.dta, clear
xtset zper year
gen d=1 if meddol>0
replace d=0 if d==.
*Program*
mata
y = st_data(., "meddol")
d = st_data(., "d")
x1=st_data(.,"logc")
x2=st_data(.,"xage")
n=rows(y)
one=J(n,1,1)
x=one,x1,x2
k=cols(x)
/* read in Person identifier */
zper=st_data(.,"zper")
/* panelsetup is a useful command in order to recognize the panel structure of the data (see help file) */
info = panelsetup(zper, 1)
pp2= -87, -32, 7.3, 714, 386
nrep=100
S = optimize_init()
optimize_init_params(S, pp2)
optimize_init_evaluator(S, &ran_tobit_v0())
optimize_init_evaluatortype(S, "v0")
optimize_init_which(S, "max" )
p=optimize(S)
p
ihess=optimize_result_V_oim(S)
ihess
iopg=optimize_result_V_opg(S)
iopg
irobust=optimize_result_V_robust(S)
irobust
scores=optimize_result_scores(S)
st_matrix("beta_sml",p)
st_matrix("oim_sml",ihess)
st_matrix("opg_sml",iopg)
st_matrix("robust_sml",irobust)
loglik=optimize_result_value(S)
loglik
st_matrix("loglik_sml",loglik)
st_matrix("n",n)
ng=rows(info)
st_matrix("ng",ng)
end
matrix colnames beta_sml=Constant logc age sigma_mu sigma_anpha
matrix colnames oim_sml=Constant logc age sigma_mu sigma_anpha
matrix rownames oim_sml=Constant logc age sigma_mu sigma_anpha
matrix colnames opg_sml=Constant logc sigma_mu sigma_anpha
matrix rownames opg_sml=Constant logc sigma_mu sigma_anpha
matrix colnames robust_sml=Constant lc age sigma_mu sigma_anpha
matrix rownames robust_sml=Constant lc age sigma_mu sigma_anpha
ereturn post beta_sml oim_sml,depname(Ran_Tobit)
ereturn display
scalar ll=loglik_sml[1,1]
outreg2 using ran_tobit_sml,ctitle("ran tobit sml")
xttobit meddol logc xage, ll(0)
outreg2 using ran_tobit_sml,append ctitle("ran tobit sml")
** and we are done***
It will be great if BAN ZHU could give me some LUN TAN BI for the codes.....hahhh~~~~
3x in advance~~~!
[em02] [em02] 
扫码加好友，拉您进群
全部版块
我的主页
收藏
