pareto.src

没有说明

/*
** pareto.src - Pareto Regression Model for Duration Data (Censoring Optional)
**
**
** (C) Copyright 1988-1995  Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**-------------------**------------------**-------------------**-----------**
**-------------------**------------------**-------------------**-----------**
**
**  FORMAT:     { bg,vc,llik } = pareto(dataset,dep,ind);
**
**  INPUT:
**      dataset = name of Gauss dataset or name of matrix in memory
**      dep     = dependent variable (duration) name or column number
**      ind     = vector of independent variable names or column numbers
**
**  OUTPUT:
**    _cn_Inference    = MAXLIK for maximum likelihood estimates
**                    = BOOT for bootstrapped estimates
**                    = PROFILE for likelihood profile and profile t traces
**
**      bg    = vector of effect parameters that maximize the likelihood
**             on top of parameter(s) corresponding to vind.
**             PARAMETERIZATION: bg=b|g;
**                 E(Y) = exp(ind*b)
**                 V(Y) = E(Y)*(1+exp(g))
**      vc   = variance-covariance matrix of b
**      llik = value of the log-likelihood at the maximum
**
**  GLOBALS:
**
**       _cn_Censor     0 = pareto model, no censoring (default).
**               symbol = use this variable from dataset to censor.  Code
**                        each row as 0 if censored or 1 if not.
**              integer = use this column of input matrix to censor. Code
**                        each row as 0 if censored or 1 if not.
**
**      _cn_Start   choose method of calculating starting values.
**                     0 = LS (default),
**                     1 = set to vector stored in _cn_StartValues,
**                     2 = rndu-0.5,
**                     3 = zeros, or set to vector
**
**      _cn_Dispersion   starting value for scalar dispersion parameter.
**                 Default = 3.
**
**      __output    1  =  print output to screen (default),
**                  0  =  do not print to screen
**
**  OTHER GLOBALS:
**      see MAXLIK.
**
**  EXAMPLE:
**      let dep=time;
**      let ind=age party unem;
**      dataset="\\gauss\\prg\\sample";
**      call pareto(dataset,dep,ind,0);
**
*/
#include count.ext
#include gauss.ext
#include maxlik.ext

proc _cn_svpar(dataset,dep,ind);
    local res,b0,b1,b,pars;
    if _cn_Dispersion == 3;
        _cn_Dispersion = .5;
    endif;
    pars = 2;
    if ind/=0;
        pars = pars+rows(ind);
    endif;
    if _cn_Start==0;
        if ind==0;
            b0 = 0;
        else;
            b0 = -lols(dataset,dep,ind);
        endif;
        res = b0|_cn_Dispersion;
    elseif _cn_Start==1;
        b = _cn_StartValues;
        res = b;
        if rows(b)/=pars;
            "b is the wrong size for _cn_Start\g";
            end;
        endif;
    elseif _cn_Start==2;
        res = rndu(pars,1)-0.5;
    elseif _cn_Start==3;
        res = zeros(pars,1);
    else;
        res = _cn_Start;
        if rows(res)/=pars;
            errorlog "\nERROR:  Wrong number of rows in _cn_Start.\n";
            end;
        endif;
    endif;
    retp(res);
endp;

proc _cn_lipar(b,dta);
    local y,n,x,ezg,xb,ps,res,resc,aa,exb,dd,z,g,bb;
    y = dta[.,1];
    n = rows(y);
    x = ones(n,1);
    if ((_cn_Censor/=0) and (cols(dta)>2));
        x = x~dta[.,2:cols(dta)-1];
    elseif (_cn_Censor==0) and (cols(dta)>1);
        x = x~dta[.,2:cols(dta)];
    endif;
    dd = ones(n,1);
    if _cn_Censor/=0;
        dd = dta[.,cols(dta)];
    endif;
    z = 1;          /* use for variance function */

    g = trimr(b,cols(x),0);
    b = trimr(b,0,1);
    xb = x*b;
    exb = exp(xb);
    ezg = exp(g);
    ps = exb./ezg;
    aa = 1+(y./exb)./(1+2*ps);
    if _cn_Censor/=0;
        resc = ln(1-(aa^(-2-2*ps)));
    else;
        resc = 1;
    endif;
    res = ln(1+ps)-xb-ln(1+2*ps)-(3+2*ps).*ln(aa);
    res = dd.*res+(1-dd).*resc;
    retp(res);
endp;

proc 3 = pareto(dataset,dep,ind);
    local vars,b,logl,g,vc,se,parnames,st,ret;
    _max_CovPar = 3;
    _cn_fn = dataset;
    if dep$==0;
        errorlog "\nERROR:  DEP must be variable name or number.\n";
        end;
    endif;
    if ((type(dataset)/=13) and ((maxc(ind)>cols(dataset)) or
        (dep>cols(dataset))) );
        errorlog "\nERROR:  If dataset is a matrix, DEP and IND "\
              "must be column numbers of the input matrix.\n";
        end;
    endif;
    if (type(dataset)==13) and (type(_cn_Censor)==13);
        _cn_Censor = indcv(_cn_Censor,getname(dataset));
    endif;
    vars = dep;
    if ind/=0;
        vars = vars|ind;
    endif;

    st = _cn_svpar(dataset,dep,ind);
    if __title $== "";
      if _cn_Censor/=0;
          __title = "Censored ";
      endif;
      __title = __title$+"Pareto Regression Model of Duration Data";
    endif;
    if _cn_Censor/=0;
        vars = vars|_cn_Censor;
    endif;
    local infm,inf0,lcInf;
    infm = { MAXLIK, BOOT };
    inf0 = { 1, 2 };
    LcInf = _ml_check(_cn_Inference,1,infm,inf0,1);
    if LcInf == 1;
        { b,logl,g,vc,ret } = maxlik(dataset,vars,&_cn_lipar,st);
    elseif LcInf == 2;
        { b,logl,g,vc,ret } = maxboot(dataset,vars,&_cn_lipar,st);
    endif;
    if ret /= 0;
        errorlog "ERROR: Model estimation failed.";
        end;
    endif;
    if type(dataset)==13;
        vars = "beta0";
        if ind/=0;
            vars = vars|ind;
        endif;
        vars = vars|"gamma";
    else;
        vars = "beta0";
        if ind/=0;
            vars = vars|
            ((0 $+ "Col." $+ zeros(rows(ind),1))$+_cn_ftosm(ind,2));
        endif;
        vars = vars|"gamma";
    endif;
    _cn_vr = vars;
    _cn_dp = dep;
    ndpclex;
    retp(b,vc,logl*_max_NumObs);
endp;