[求助]用mata做ML estimation for Tobit model with panel data

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

<p>请问有没stata有高手用mata做过 Tobit ml estimator(error term homogeneous&no-autocorrelation),对panel data, random effect? </p><p>前段时间做过最简单的probit model with cross-sectional data.....现在因为是panel data,要考虑random effect...所以现在的optimizing function更复杂了，因为要integrate out the random effect variable.....可能用到gaussian quadrature...或者simulated maximum.....具体思路已经有了，但是因为这个周末要去别的城市参加工作面试~~~~~不知道自己还有没有时间programme出来了~~~~~不行就只有在火车上用笔纸试着写code了，汗~~~~</p><p>如果哪位大侠有编好的code,重金求购~~~！</p><p>或者哪位大虾在某本资料上介绍过用ml做tobit model panel data analysis,可以指点一下~~~~</p><p>主要是求pooled estimator, random effect estimator...and corresponding standard errors....</p><p>多谢！</p>

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

关键词：panel data Estimation ATION Tobit model programme function error

我可不会，我也想知道呢

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***

蓝色  金币 +1  金钱 +100  奖励 2009-6-4 21:46:41

It will be great if BAN ZHU could give me some LUN TAN BI for the codes.....hahhh~~~~

3x in advance~~~!

a~~~i saw Banzhu have already give&ntilde; me some LUNTANBI and 1 JINBI~~!!!!
KAI XIN~~~~ [em02] [em02]

Thanks BANZHU~~~~~~

I am busy with my bachelor thesis these days, which uses Matlab......I think I will use Stata for my master thesis, which I will probabily start with in weeks.......hope I could get some help with coins I have then!!

Thanks again.......

1# cosa

xttobit 不就是可以的么？