"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/addons/HIP/HIP_public.inp" (25 Feb 2020, 10885 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Charmm source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. See also the last Fossies "Diffs" side-by-side code changes report for "HIP_public.inp": 2019d_vs_2020a.

    1 /***********************************************************
    2  Bundle structure. 
    3  IVPROBIT: initial values for mle are computed by
    4  a simple Rivers&Vuong - two-step by the RiversVuong 
    5  function. Cholesky decomposition of the south-east 
    6  block of the covariance matrix is used.   
    7  HETPROBIT: heteroskedasticity of known form. 
    8 
    9  y_i = Z_i'\beta + \epsilon_i
   10  Z_i = ENDOG_i EXOG_i (Y_i X_1i)
   11  Y_i = \Pi'X_i + U_i
   12  X_i = EXOG ADDIN
   13 ***********************************************************/
   14 
   15 function bundle HIP_setup(series y, list EXOG, list ENDOG[null], \
   16               list ADDIN[null], list HETVAR[null])
   17 
   18     bundle ret = null
   19 
   20     # first, perform a few sanity checks
   21     
   22     scalar id = nelem(ENDOG) - nelem(ADDIN)
   23 
   24     if (nelem(ADDIN)!=0 && nelem(ENDOG)==0) 
   25     printf "No endogenous variables specified. Are you sure this is what you want?\n"
   26     return ret
   27     elif (nelem(ADDIN)==0 && nelem(ENDOG)!=0)    
   28     scalar w = nelem(ENDOG)
   29     printf "At least %g additional instruments need to be specified.\n", w
   30     return ret
   31     elif (id>0) 
   32     printf "At least %g more instruments are needed. \n", id
   33     return ret
   34     endif
   35     
   36     # then, check that y is a dummy variable
   37     check = isdummy(y) > 0
   38 
   39     if check == 0
   40     printf "%s is not a dummy variable; aborting\n", argname(y)
   41     ret["err"] = 48
   42     return ret
   43     else
   44     ret["err"] = 0
   45     endif
   46 
   47     # next, check if HIP is really needed at all
   48 
   49     ret["simple_probit"] = nelem(ENDOG || HETVAR || ADDIN) == 0
   50 
   51     ret["het"] = (nelem(HETVAR)!=0)
   52     ret["iv"] = (nelem(ADDIN)!=0 && nelem(ENDOG)!=0)
   53     ret["vcvmeth"] = 0
   54     ret["verbose"] = 1
   55     ret["ntot"] = $nobs
   56     ret["neff"] = sum(ok(y || EXOG || ENDOG || HETVAR || ADDIN))
   57     ret["depvar"] = y
   58     ret["depvarname"] = argname(y)
   59 
   60     # Exogenous variables: X_1i
   61     ret["mk1"] = nelem(EXOG)
   62     ret["mEXOG"] = { EXOG } 
   63     ret["mEXOGnames"] = varname(EXOG)
   64 
   65     # Endogenous variables: Y_i
   66     ret["mp"] = nelem(ENDOG)
   67     ret["mENDOG"] = { ENDOG } 
   68     ret["mENDOGnames"] = varname(ENDOG)
   69      
   70     # Additional instruments: X_2i
   71     ret["mk2"] = nelem(ADDIN) 
   72     ret["mADDIN"] = { ADDIN } 
   73     ret["mADDINnames"] = varname(ADDIN)
   74 
   75     # Heteros. variables: W_i
   76     ret["mq"] = nelem(HETVAR)  
   77     ret["mHETVAR"] = { HETVAR } 
   78     ret["mHETVARnames"] = varname(HETVAR)
   79 
   80     # Total regressors in main equation
   81     list Z = EXOG ENDOG 
   82     ret["mh"] = ret["mk1"] + ret["mp"] 
   83     ret["mZ"] = { Z } 
   84     ret["mZnames"] = varname(Z)
   85 
   86     # Total instruments
   87     list X = EXOG || ADDIN    
   88     ret["mk"] = ret["mk1"] + ret["mk2"]
   89     ret["mX"] = { X } 
   90     ret["mXnames"] = varname(X)
   91 
   92     matrix Q = InitParm(&ret)
   93     ret["theta"] = Q[,1]
   94 
   95     return ret
   96 end function
   97 
   98 function scalar HIP_setoption(bundle *b, string opt, scalar value)
   99 
  100     err = 0
  101    
  102     if opt=="verbose"
  103         if (value < 0) || (value > 3)
  104             err = 1
  105         else
  106             b["verbose"] = value
  107         endif
  108     elif opt=="vcvmeth"
  109         if (value < 0) || (value > 2)
  110             err = 2
  111         else
  112             b["vcvmeth"] = value
  113         endif
  114     else
  115         err = 3
  116     endif
  117     
  118     if err>0
  119         printf "Warning in DPB_setoption\n"
  120     if err < 3
  121         printf "Illegal value %s = %d\n", opt, value
  122     else
  123         printf "%s unrecognised option\n", opt
  124     endif
  125     endif
  126 
  127     return err
  128 end function
  129 
  130 
  131 /********************************************************
  132 
  133  After mle estimation, replaces the estimated parameters
  134  and covariance matrix in the
  135  bundle structure in order to be used for print-out.     
  136 
  137  Options v and s:
  138  - verbose option v: 
  139  0 = quiet, 1 = standard output(default), 
  140  2 = first stages, 3 = verbose mle.
  141 
  142  - standard error estimation s:
  143  0 = OPG estimator(default), 1 = empirical Hessian, 2 = sandwich. 
  144 
  145 *********************************************************/
  146 
  147 function scalar HIP_estimate(bundle *b) 
  148     
  149     series y = b["depvar"]
  150     matrix EXOG   = b["mEXOG"]
  151     matrix ENDOG  = b["mENDOG"]
  152     matrix ADDIN  = b["mADDIN"]
  153     matrix HETVAR = b["mHETVAR"]
  154     scalar iv = b["iv"]
  155         
  156     scalar k1 = b["mk1"]
  157     scalar k2 = b["mk2"]
  158     scalar p =  b["mp"]
  159     scalar q =  b["mq"]
  160     
  161     scalar s = b["vcvmeth"]
  162     scalar verbose = b["verbose"] 
  163 
  164     # handle options for MLE
  165 
  166     if (verbose==3)
  167     setopt mle --verbose
  168     else
  169     setopt mle --quiet
  170     endif
  171         
  172     if (s==1)
  173     setopt mle --hessian
  174     elif(s==2)
  175     setopt mle --robust
  176     endif
  177 
  178     matrix theta = b["theta"]
  179     set warnings off
  180 
  181     matrix beta = {}
  182     matrix alpha = {}
  183     matrix Pi = {}
  184     matrix psi = {}
  185     matrix C = {}
  186 
  187     series llik = NA
  188     series sigma = NA
  189     matrix omega = {}
  190     matrix SCORE = {}
  191 
  192     # guard against weird initialisation in the conditional
  193 
  194     start_ok = 0
  195     iters = 0
  196     loop while (start_ok == 0) && (iters<10) --quiet
  197     iters++
  198     scalar err = HIP_params_shape(&theta, &beta, &alpha, &Pi, 
  199                       &psi, &C, k1, k2, p, q) 
  200 
  201     scalar err = HIP_loglik(y, &EXOG, &ENDOG, &ADDIN, &HETVAR, 
  202                 &beta, &alpha, &Pi, &psi, &C,
  203                 &llik, &omega, &sigma) 
  204     
  205     if err == 4
  206         theta[1:(k1+k2)] = theta[1:(k1+k2)] ./ 2
  207         printf "Brrr! Halving params\n"
  208     else
  209         start_ok = 1
  210     endif
  211     endloop
  212 
  213     if iters==10
  214     printf "Man, these data are weird! Giving up.\n"
  215     return 101
  216     endif
  217 
  218     SCORE = {}
  219 
  220     catch mle ll = llik
  221     scalar err = HIP_params_shape(&theta, &beta, &alpha, &Pi, 
  222                       &psi, &C, k1, k2, p, q) 
  223     scalar err = HIP_loglik(y, &EXOG, &ENDOG, &ADDIN, &HETVAR, 
  224                 &beta, &alpha, &Pi, &psi, &C,
  225                 &llik, &omega, &sigma) 
  226     SCORE = err ? zeros(1,rows(theta)) : \
  227         HIP_Score(y, &EXOG, &ENDOG, &ADDIN, &HETVAR,
  228               &beta, &alpha, &Pi, &psi, &C,
  229               &omega, &sigma)       
  230     
  231     deriv theta = SCORE
  232     end mle 
  233 
  234     errcode = $error
  235     b["errcode"] = $error
  236 
  237     if errcode == 0
  238     matrix V = $vcv
  239     b["theta"] = theta
  240     b["VCVtheta"] = V       
  241     
  242     if s==1
  243         b["vcvmeth"] = "Hessian"
  244     elif s==2
  245         b["vcvmeth"] = "Sandwich"
  246     else
  247         b["vcvmeth"] = "OPG"
  248     endif
  249 
  250     if iv
  251         matrix X = EXOG ~ ADDIN
  252         scalar J = rescale_results(&b, &theta, &V, &SCORE)
  253         scalar ll_m  = sum(loglik_m(&ENDOG, &X, &Pi, &C, &omega))
  254         b["lnl1m"] = ll_m  - J
  255         b["theta"] = theta
  256         b["VCVtheta"] = V       
  257         else 
  258         scalar J = 0
  259         ll_m = 0
  260         endif       
  261 
  262     b["T"] = $T
  263     b["t1"] = $t1
  264     b["t2"] = $t2
  265 
  266     b["llt"] = llik - J
  267     J = J*b["T"]
  268     b["lnl1"] = $lnl - J
  269     
  270     b["lnl1c"] = $lnl - ll_m
  271     b["infocrit"] = ($aic ~ $bic ~ $hqc) + 2*J
  272 
  273     b["SCORE"] = SCORE
  274 
  275     HIP_diagnostics(&b)
  276     else
  277     printf "errcode = %d\n", errcode
  278     endif
  279 
  280     return errcode
  281 end function
  282 
  283 
  284 
  285 function void HIP_printout(bundle *b)
  286 
  287     scalar k1 = b["mk1"]
  288     scalar k2 = b["mk2"]
  289     scalar p =  b["mp"]
  290     scalar q =  b["mq"]
  291     scalar v =  b["verbose"]
  292 
  293     theta = b["theta"]
  294     matrix beta = {}
  295     matrix alpha = {}
  296     matrix Pi = {}
  297     matrix psi = {}
  298     matrix C = {}
  299     err = HIP_params_shape(&theta, &beta, &alpha, &Pi, &psi,
  300                &C, k1, k2, p, q)
  301 
  302     stderr = sqrt(diag(b["VCVtheta"]))
  303 
  304     matrix sd_beta = {}
  305     matrix sd_alpha = {}
  306     matrix sd_Pi = {}
  307     matrix sd_psi = {}
  308     matrix sd_C = {}
  309     err = HIP_params_shape(&stderr, &sd_beta, &sd_alpha, 
  310                &sd_Pi, &sd_psi, 
  311                &sd_C, k1, k2, p, q)
  312 
  313     scalar iv = b["iv"]
  314     scalar het = b["het"]
  315     s = b["vcvmeth"]
  316 
  317     matrix V = alpha ~ sd_alpha
  318 
  319     if (het==1 && iv==1)
  320     printf "Heteroskedastic probit model with endogenous regressors\n"
  321     elif (het==0 && iv==1) 
  322     printf "Probit model with endogenous regressors\n"
  323     elif (het==1 && iv==0)
  324     printf "Heteroskedastic probit model \n"
  325     else
  326     printf "Probit model \n"
  327     endif
  328 
  329     printf "ML, using observations %d-%d\n", b["t1"], b["t2"]
  330 
  331     string depvar = b["depvarname"]
  332     printf "Dependent Variable: @depvar \n"
  333     if iv 
  334     string endog = b["mENDOGnames"]
  335     string inst_names = b["mXnames"]
  336     printf "Instrumented: %s\n", strsub(endog, ",", ", ")
  337     printf "Instruments: %s \n", strsub(inst_names, ",", ", ")
  338     endif
  339 
  340     printf "Parameter covariance matrix: %s\n", b["vcvmeth"]
  341 
  342     mnames = b["mZnames"]
  343     matrix M = beta ~ sd_beta
  344     modprint M mnames
  345      
  346     if het
  347     printf "Variance \n"
  348     mnames = b["mHETVARnames"]
  349         modprint V mnames
  350     endif
  351 
  352     if v>1 &&  iv
  353         p = b["mp"]
  354         printf "\"First-stage\" regressions\n"
  355         loop i=1..p -q
  356             matrix K = Pi[,$i] ~ sd_Pi[,$i]
  357             modprint K inst_names
  358         endloop
  359     endif
  360     infocrit = b["infocrit"]
  361     printf "Log-likelihood   %14.4f  Akaike criterion %12.4f\n", \
  362     b["lnl1"], infocrit[1]
  363     printf "Schwarz criterion  %12.4f  Hannan-Quinn   %14.4f\n", \
  364     infocrit[2], infocrit[3]
  365     if iv
  366     printf "Conditional ll     %12.6f  Cragg-Donald stat. %10.3f\n\n", \
  367     b["lnl1c"], b["CraggDonald"] 
  368     else
  369         printf "\n"
  370     endif
  371 
  372     WT = b["WaldAll"]
  373     printf "Overall test (Wald) = %g (%d df, p-value = %6.4f)\n", \
  374     WT[1], WT[2], WT[3] 
  375 
  376     if !(iv || het) # ordinary probit
  377     CI = b["normtest"]
  378     printf "Chesher and Irish normality test = %g (%d df, p-value = %6.4f)\n", \
  379     CI[1], CI[2], CI[3]
  380     endif
  381     
  382     if iv
  383     WT = b["WaldEnd"]
  384     printf "Endogeneity test (Wald) = %g (%d df, p-value = %6.4f)\n", \
  385     WT[1], WT[2], WT[3] 
  386         if (k2>p)
  387         LM = b["LMOverid"]
  388         printf "Test for overidentifying restrictions (LM) = %g (%d df, p-value = %6.4f)\n", \
  389         LM[1], LM[2], LM[3] 
  390     endif
  391     endif
  392 
  393     if het
  394     T = b["HETtest"]
  395     if iv
  396          printf "Heteroskedasticity test (Wald) = %g (%d df, p-value = %6.4f)\n", \
  397          T[1], T[2], T[3] 
  398      else
  399          printf "Heteroskedasticity test (LR) = %g (%d df, p-value = %6.4f)\n", \
  400          T[1], T[2], T[3] 
  401 
  402          CI = b["normtest"]
  403          printf "Chesher and Irish normality test = %g (%d df, p-value = %6.4f)\n", \
  404          CI[1], CI[2], CI[3]
  405      endif
  406 
  407      endif
  408 
  409  end function
  410 
  411 function bundle HIP(series y "Dependent variable", 
  412             list EXOG "Exogenous variables", 
  413             list ENDOG[null] "Endogenous variables",\ 
  414             list ADDIN[null] "Instruments", 
  415             list HETVAR[null] "Variance regressors", \
  416             int v[0:3:1] "Verbosity level", 
  417             int s[0:2:0] "Covariance matrix estimation" \
  418             {"OPG", "Hessian", "Sandwich"})
  419 
  420     # first thing, drop all obs with missing values anywhere
  421     list EVERYTHING = y || EXOG || ENDOG || ADDIN || HETVAR
  422     smpl EVERYTHING --no-missing
  423 
  424 
  425     # first, force the constant to be the first exogenous variable
  426     EXOG = EXOG - const
  427     EXOG = const EXOG
  428 
  429     # let the dance begin
  430     set warnings off
  431 
  432     bundle b = HIP_setup(y, EXOG, ENDOG, ADDIN, HETVAR)
  433 
  434     if exists(b)
  435     if v != 1
  436         HIP_setoption(&b, "verbose", v)
  437     endif
  438 
  439     if s != 0
  440         HIP_setoption(&b, "vcvmeth", s)
  441     endif
  442 
  443     b["depvarname"] = argname(y)
  444     err = b["err"]
  445     if ! err
  446         err = HIP_estimate(&b)
  447     endif
  448 
  449     if !err 
  450         if v > 0
  451         HIP_printout(&b)
  452         else
  453         printf "Done.\n"
  454         endif
  455     else
  456         printf "Error = %s\n", errmsg(err)
  457     endif
  458     endif
  459    
  460    return b
  461  end function   
  462