/*
Probit and Logit (Rizzo 1995)
This source code is provided free of any charge for public,
 non-commercial use. Please provide appropriate acknowledgment
 if it supports published research or if you incorporate
 portions of this code into your own.
This code has not been subject to formal test criteria and is
 provided without any performance guarantees.
*/

/** INSTRUCTIONS: THIS FILE MUST BE BROKEN INTO 5 FILES AT THE
    ``BREAKPOINTS'' BELOW: PROBLOG.PGM (AN EXAMPLE PROGRAM)
                           PROBLOG.DTA (DATA FOR THE EXAMPLE)
                           PROBLOG.LCG (THE LIBRARY FILE)
                           PROBLOG.DEC (DECLARATIONS FILE)
                           PROBLOG.SRC (THE PROGRAM FILE)
    PLACE THESE 5 FILES INTO C:\GAUSS\PROBLOG, A NEW SUBDIRECTORY.
    COPY PROBLOG.LCG INTO YOUR C:\GAUSS\LIB SUBDIRECTORY.
    RUN THE EXAMPLE PROGRAM. ENJOY!
**/

/******** BREAKPOINT: PROBLOG.PGM FOLLOWS ******************/
new;
/** PROBLOG.PGM: AN EXAMPLE USE OF PROBLOG.SRC **/
library problog;  globset;
proc (0) = crosstab(str,names);retp;endp;@DELETE IF YOU HAVE BSTAT MODULE@
/***********************************************/
/* ADJUST THIS SECTION ACCORDING TO YOUR NEEDS */
/***********************************************/

/** SET FILE NAMES (optional) **/
_fname  = "c:\\gauss\\problog\\problog.dta";      @ data file    @
_pname  = "c:\\gauss\\problog\\prloblog.pgm";    @ program file @
_ouname = "c:\\gauss\\problog\\problog.out";     @ output file  @
output file = ^_ouname, reset;
/**********************************************/

/** LOAD DATA **/
load dat[32,5]=^_fname;
one   = ones(rows(dat),1);                   @ ind. vars    @
gpa   = dat[.,2];                            @   ::::::     @
tuce  = dat[.,3];                            @   ::::::     @
psi   = dat[.,4];                            @   ::::::     @
x     = one~gpa~tuce~psi;                    @   ::::::     @
grade = dat[.,5];                            @ dep. var.    @

/** NAMES OF THE VARIABLES (STOF TO GET LOWER CASE) **/
let __altnam=grade constant gpa tuce psi;

/******************************************************
/** RESTRICTIONS FOR THE WALD TEST **/
     /* the last 2 coeffs = 0 */

let _rwald = {0 0 1 0,
              0 0 0 1};
let _qwald =  0 0;
/** EVALUATION DATA POINT FOR THE MARGINAL ANALISYS **/
/**    OTHER THAN THE DEFAULT (MEAN OF THE DATA)    **/
_xstar=x[4,.]'; /*in vector form*/
******************************************************/

/** PROBIT AND LOGIT ESTIMATION **/
{bp,covbp,lkp}=problog(grade,x,"probit");
{bl,covbl,lkl}=problog(grade,x,"logit");


/* marginal effect analysis at each data point (probit/logit) *

mdl=0;
do until mdl gt 1;
     i=1;
     do until i gt rows(x);
          _xstar=x[i,.]';
          "DATA POINT " I;
          if mdl ==0;
               call derivati(bp,x,covbp,mdl);
          elseif mdl==1;
               call derivati(bl,x,covbl,mdl);
          endif;
          i=i+1;
     endo;
mdl=mdl+1;
endo;

*/


/******************* THE ABOVE PROGRAM GENERATES THE FOLLOWING OUTPUT *********
================================================================================
 PROBLOG.SRC  -  PROBIT MODEL :  VERSION 1.02 (R1)            6/29/95   4:25 PM
================================================================================
 ANALYZING FILE: C:\GAUSS�PROBLOG�PROBLOG.DTA
 PROGRAM   FILE: C:\GAUSS�PROBLOG�PROBLOG.PGM
 OUTPUT    FILE: C:\GAUSS�PROBLOG�PROBBLOG.OUT

 _stval = 0 --> Starting values are OLS (default)

>>>>> _diagn = 0   No diagnostic output option selected <<<<<

Binary Probit Variable     :  grade
Number of observations     :  32   grade = 0 --> 21      grade = 1 --> 11
Percent of positive obs    :  34.37500
Convergence achieved in    :  6 iterations (max =  25 )
Convergence criteria       :  maxc(abs(db)) where db = -inv(H)*g
   value at convergence    :  1.1526E-09
Convergence tolerance      :  1.0000E-06
--------------------------------------------------------------------------------
Error sum of squares       :   4.174660       Total sum of squares  :  7.218750
Residual sum of squares    :   3.044090       R                     :  0.421692
Unrestricted log-lik       : -12.818804       R  - Laitila '93      :  0.570418
Restr.log-lik(slopes=0)    : -20.591730       R  - McFadden         :  0.377478
Likelihood Ratio=2(LU-LR)  :  15.545851       Wald test ( =0)       : 11.254100
   p-value                 :   0.001405          p-value            :  0.023852
   d.f.                    :   3.000000          d.f.               :  4.000000
--------------------------------------------------------------------------------
variable       _ols      _probit     std.err.     t-val    prob|t| x   gradient
--------------------------------------------------------------------------------
constant     -1.4980    -7.45232     2.54247    -2.93113    0.0034   -7.903e-10
gpa           0.4639     1.62581     0.69388     2.34306    0.0191   -2.164e-09
tuce          0.0105     0.05173     0.08389     0.61663    0.5375   -1.426e-08
psi           0.3786     1.42633     0.59504     2.39704    0.0165    2.986e-11
--------------------------------------------------------------------------------
================================================================================
 PROBIT MARGINAL EFFECT ANALYSIS (at the mean of the data)    6/29/95   4:25 pm
================================================================================
                             scale factor  =  0.32805
                             utility index = -0.62554
--------------------------------------------------------------------------------
variable    data point    slope      std.err.     t-val.     prob|t| x   exp(b)
--------------------------------------------------------------------------------
constant       1.000    -2.44473      0.7589     -3.2216      0.0013    0.00058
gpa            3.117     0.53335      0.2325      2.2943      0.0218    5.08253
tuce          21.938     0.01697      0.0271      0.6257      0.5315    1.05309
psi            0.438     0.46791      0.1876      2.4936      0.0126    4.16340
--------------------------------------------------------------------------------
********************************************************************************
>>>>>>>>>>>>>>>>>>>>>>>>>>> end of probit estimation <<<<<<<<<<<<<<<<<<<<<<<<<<********************************************************************************
================================================================================
 PROBLOG.ARC  -  logit model :  Version 1.02 (R1)             6/29/95   4:25 pm
================================================================================
 ANALYZING FILE: C:\GAUSS�PROBLOG�PROBLOG.DTA
 PROGRAM   FILE: C:\GAUSS�PROBLOG�PROBLOG.PGM
 OUTPUT    FILE: C:\GAUSS�PROBLOG�PROBLOG.OUT

 _stval = 0 --> Starting values are OLS (default)

>>>>> _diagn = 0   No diagnostic output option selected <<<<<
Binary Logit Variable      :  grade
Number of observations     :  32   grade = 0 --> 21      grade = 1 --> 11
Percent of positive obs    :  34.37500
Convergence achieved in    :  6 iterations (max =  25 )
Convergence criteria       :  maxc(abs(db)) where db = -inv(H)*g
   value at convergence    :  3.8029E-07
Convergence tolerance      :  1.0000E-06
--------------------------------------------------------------------------------
Error sum of squares       :   4.144171       Total sum of squares  :  7.218750
Residual sum of squares    :   3.074579       R                     :  0.425916
Unrestricted log-lik       : -12.889634       R  - Laitila '93      :  0.543687
Restr.log-lik(slopes=0)    : -20.591730       R  - McFadden         :  0.374038
Likelihood Ratio=2(LU-LR)  :  15.404191       Wald test ( =0)       :  8.873129
   p-value                 :   0.001502          p-value            :  0.064350
   d.f.                    :   3.000000          d.f.               :  4.000000
--------------------------------------------------------------------------------
variable       _ols      _logit      std.err.     t-val    prob|t| x   gradient
--------------------------------------------------------------------------------
constant     -1.4980   -13.02135     4.93132    -2.64054    0.0083   -6.258e-08
gpa           0.4639     2.82611     1.26294     2.23772    0.0252   -1.619e-07
tuce          0.0105     0.09516     0.14155     0.67223    0.5014   -1.148e-06
psi           0.3786     2.37869     1.06456     2.23442    0.0255   -1.545e-09
--------------------------------------------------------------------------------
================================================================================
 LOGIT MARGINAL EFFECT ANALYSIS (at the mean of the data)     6/29/95   4:25 pm
================================================================================
                             scale factor  =  0.18890
                             utility index = -1.08363
--------------------------------------------------------------------------------
variable    data point    slope      std.err.     t-val.     prob|t| x   exp(b)
--------------------------------------------------------------------------------
constant       1.000    -2.45976      0.8177     -3.0081      0.0026    0.00000
gpa            3.117     0.53386      0.2370      2.2522      0.0243   16.87971
tuce          21.938     0.01798      0.0262      0.6851      0.4933    1.09983
psi            0.438     0.44934      0.1968      2.2837      0.0224   10.79073
--------------------------------------------------------------------------------
********************************************************************************
>>>>>>>>>>>>>>>>>>>>>>>>>>> end of logit  estimation <<<<<<<<<<<<<<<<<<<<<<<<<<********************************************************************************

***********************/

/*************** BREAK POINT: PROBLOG.DTA FOLLOWS *****************/

@PROBLOG.DTA: LISTS THE DATA FROM
table 21.1  Greene 1993 - page 639 -
====================================
OBS  GPA    TUCE  PSI  GRADE  (32x5)
@
  1   2.66     20    0      0
  2   2.89     22    0      0
  3   3.28     24    0      0
  4   2.92     12    0      0
  5   4.00     21    0      1
  6   2.86     17    0      0
  7   2.76     17    0      0
  8   2.87     21    0      0
  9   3.03     25    0      0
 10   3.92     29    0      1
 11   2.63     20    0      0
 12   3.32     23    0      0
 13   3.57     23    0      0
 14   3.26     25    0      1
 15   3.53     26    0      0
 16   2.74     19    0      0
 17   2.75     25    0      0
 18   2.83     19    0      0
 19   3.12     23    1      0
 20   3.16     25    1      1
 21   2.06     22    1      0
 22   3.62     28    1      1
 23   2.89     14    1      0
 24   3.51     26    1      0
 25   3.54     24    1      1
 26   2.83     27    1      1
 27   3.39     17    1      1
 28   2.67     24    1      0
 29   3.65     21    1      1
 30   4.00     23    1      1
 31   3.10     21    1      0
 32   2.39     19    1      1

/****************** BREAK POINT: PROBLOG.LCG FOLLOWS *******************/
c:\gauss\problog\PROBLOG.SRC
    problog  : proc
    llikfun  : proc
    derivati : proc
    globset  : proc
    globshow : proc
    prlogr   : proc

c:\gauss\problog\PROBLOG.DEC
    _corr    : matrix
    _cross   : matrix
    _inter   : matrix
    _marg    : matrix
    _headr   : matrix
    _tst     : matrix
    _toler   : matrix
    _maxit   : matrix
    _fname   : string
    _rwald   : matrix
    _qwald   : matrix
    _diagn   : matrix
    _xstar   : matrix
    _pname   : string
    _ouname  : string
    _stval   : matrix

c:\gauss\src\MATWRT.SRC
    matwrt   :proc

@
This final section allows use of the crosstable
option in the problog.arc program.

Uncomment if you have the "bstat" module
-----------------------------------------

c:\gauss\src\BSTATSET.ARC
    bstatset : proc

c:\gauss\src\BSTAT.DEC
    _bs_ver  : matrix
    _ctstat  : matrix
    _ctpause : matrix
    _ctrow   : matrix
    _ctcol   : matrix
    _ctrowp  : matrix
    _ctcolp  : matrix
    _cttotp  : matrix
    _cttype  : matrix
    _ct_tp   : matrix
    _ctcase  : matrix
    _ctnomis : matrix
    _frpause : matrix
    _frstat  : matrix
    _frfreq  : matrix
    _frtype  : matrix
    _fr_tp   : matrix
    _frcase  : matrix
    _frmiss  : matrix
    _sqndep  : matrix
    _regcor  : matrix
    _regcol  : matrix
    _regvc   : matrix
    _regmmt  : matrix
    _regcort : matrix
    _regdesc : matrix
    _regres  : matrix
    _regrlbl : matrix
    _regrnam : string
    _regln   : matrix
    _reghc   : matrix
    _ttgrpnm : matrix
    _ttcut   : matrix
    _mbcval  : matrix
    _mbfrq   : matrix
    _mbnvars : matrix
    _corcor  : matrix
    _corvc   : matrix
    _cormmt  : matrix
    _cort    : matrix
    _cordesc : matrix
    _rgmcor  : matrix
    _rgmvc   : matrix
    _rgmmmt  : matrix
    _rgmcort : matrix
    _rgmdesc : matrix

c:\gauss\src\TTEST.ARC
    ttest    : proc

c:\gauss\src\REGM.ARC
    regm     : proc

c:\gauss\src\MEANS.ARC
    means    : proc

c:\gauss\src\MEANSBY.ARC
    _celpos  : proc
    meansby  : proc

c:\gauss\src\CORREG.ARC
    reg      : proc
    regres   : proc
    corr     : proc

/* CROSSTAB */

c:\gauss\src\CROSSTAB.ARC
    crosstab : proc
    _ctprt   : proc
    _cttable : proc
    _ctinit  : proc
    _ctmkcat : proc
    sbscatch : proc

/* FREQ */

c:\gauss\src\FREQ.ARC
    freq     : proc
    _frinit  : proc
    _frmkcat : proc
    _frtable : proc
    countch  : proc
    cntchwt  : proc

c:\gauss\src\FREQSTAT.ARC
    freqstat : proc

c:\gauss\src\GETFREQ.ARC
    getfreq  : proc

c:\gauss\src\TBLSTAT.ARC
    tblstat  : proc

/* GAUSS globals */

c:\gauss\src\GAUSS.SRC
    getnr    : proc
    header   : proc
    printfmt : proc
    formatnv : proc
    formatcv : proc
    gausset  : proc

c:\gauss\src\GAUSS.DEC
    __altnam : matrix
    __con    : matrix
    __ff     : matrix
    __fmtcv  : matrix
    __fmtnv  : matrix
    __header : string
    __miss   : matrix
    __output : matrix
    __row    : matrix
    __rowfac : matrix
    __sort   : matrix
    __title  : string
    __tol    : matrix
    __vpad   : matrix
    __vtype  : matrix
    __weight : matrix

/* DTRAN */

c:\gauss\src\DTRAN.ARC
    _lt      : proc
    _le      : proc
    _gt      : proc
    _ge      : proc
    _clt     : proc
    _cle     : proc
    _cgt     : proc
    _cge     : proc
    _ceq     : proc
    _cne     : proc
    _cilt    : proc
    _cile    : proc
    _cigt    : proc
    _cige    : proc
    _cieq    : proc
    _cine    : proc
    _slctobs : proc
    dtran    : proc
    dtranchk : proc
    dtranobs : proc
    dtranprt : proc
    dtranset : proc

c:\gauss\src\DTRAN.DEC
    _dt_ver  : matrix
    _dtsel   : matrix
    _dtproc  : matrix
    _dtstyp  : matrix
    _dtprt   : matrix
    _dtname  : string
    _dtcs    : matrix
    _dt_in   : matrix
    _dt_mnm  : matrix
    _dt_snm  : matrix

/* DUTIL */

c:\gauss\src\NUMLST.ARC
    numlst   : proc

c:\gauss\src\DUTIL.DEC
    _dl_ver  : matrix
    _dlsort  : matrix
    _dlstrow : matrix
    _dlstpau : matrix
    _dlstdgt : matrix
    _dlstst  : matrix
    _dlstvar : matrix
    _dlmkey  : matrix
    _dlmsort : matrix
    _dlmdot  : matrix
    _dlsordr : matrix
    _dlgetmk : matrix

@

/************ BREAK POINT: PROBLOG.DEC FOLLOWS *******************/
/** PROBLOG.DEC  -- globals for problog.arc **/

declare string __header != "";

declare matrix _corr   ! = 0;
declare matrix _cross  ! = 0;
declare matrix _inter  ! = 1;
declare matrix _marg   ! = 1;
declare matrix _headr  ! = 1;
declare matrix _tst    ! = 1;
declare matrix _toler  ! = 1e-6;
declare matrix _maxit  ! = 25;
declare string _fname  ! = "";
declare matrix _rwald  ! = 0;
declare matrix _qwald  ! = 0;
declare matrix _diagn  ! = 0;
declare matrix _xstar  ! = 0;
declare string _pname  ! = "";
declare string _ouname ! = "";
declare matrix _stval  ! = 0;

/******************** BREAKPOINT: PROBLOG.SRC FOLLOWS **************/
/*
** PROCEDURE PROBLOG.SRC
** last modified: 06-14-1995
** (with minor alterations by Alan G. Isaac on 3-15-96)
**
** XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
** Written by:
** Fernando Rizzo
** Via Nicotera 74
** 73015 Salice Salentino
** Italy
** email: C.A.Bianchi@sussex.ac.uk  (respect case)
**        (friend's account)
** This procedure was written using GAUSS 2.2
** (but runs on GAUSS 3.2.13--Alan G. Isaac)
**
** DISCLAIMER:
** ~~~~~~~~~~~
** The code is written and submitted for public and non-commercial use.
** The code is provided as is, without any guarantee of any kind. The
** author shall have no liability or responsibility to the user for any
** occurrence.
**
** USE OF THE CODE:
**
** The code can be used and distributed as long as proper attribution
** to the author is included. If you have any comment or if you use the
** code or portion of it, I would appreciate if you inform me at the
** above address (copy of papers using the code are also appreciated).
** XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
**
** PURPOSE:
**       Compute maximum likelihood probit and logit estimation using
**       Newton-Rapson method with analytical derivatives. An additional
**       procedure is included that computes marginal effects analysis at
**       any given data point.
**
** User must set the library problog.LCG in the usual way.
** BE SURE TO CHANGE PROBLOG.LCG IF YOU HAVE THE BSTAT MODULE.
**
** FORMAT:
**       {b,covb,lkdu}=problog(y,x,"_mdl").
**
** INPUT:
**       y        -- tx1 matrix of dependent variable.
**
**       x        -- txk matrix of independent variables.
**                   User must include a column of ones.
**
**       _mdl     -- string: "probit" or "logit"
**
** SAVED OUTPUT:
**
**       b        -- vector of estimates (probit or logit).
**       covb     -- asymptotic covariance matrix.
**       lkdu     -- value of the log likelihood function at convergence.
**
** OUTPUT (printed out with default globals setting):
**
**                -- number of obs (total, y=1, y=0)
**                -- % of positive observations
**       iter     -- number of iterations required for convergence
**       _maxit   -- maximum number of iterations
**       crit     -- convergence criteria and its value at convergence
**       _toler   -- convergence tolerance
**       ess      -- error sum of squares
**       tss      -- total sum of squares
**       rss      -- residual sum of squares
**       lkdu     -- unrestricted log-lik function
**       lkdr     -- restricted (slopes=0) log-lik function
**       lr       -- likelihood ratio statistic
**       plr_v    -- p-value for lr (chi2, df=k-1)
**       r2       -- R2 = (1-ess/tss)
**       r2lai    -- pseudo-R2 measure from T. Laitila, A pseudo R2
**                   measure for limited and qualitative dependent
**                   variables, Journal of Econometrics 56 (1993).
**       r2mf     -- pseudo-R2 measure from McFadden.
**       waldu    -- wald statistic (b=0) (see globals below for testing
**                   restrictions on the coefficients using the Wald test)
**       pwaldu   -- pvalue for wald (chi2, df=k)
**                -- independent variables names
**       (y/x)    -- ols estimates.
**       b        -- probit or logit estimates
**       std      -- standard errors
**       tv       -- t-values
**       p_tv     -- t-statistics for probit or logit coefficients (using
**                   the normal cdf)
**
**
**                   MARGINAL EFFECT ANALYSIS
**                   (by default the analysis is performed at the mean
**                   of the data. The data point can be selected by the
**                   the user by changing the global _xstar -see below)
**
**       dens     -- scale factor
**       xstar'b  -- utility index
**                -- independent variables names
**       xstar    -- selected data point
**       part     -- value of the partial derivatives
**       stdpart  -- partial derivatives' standard errors.
**       tvpart   -- t-values of partial derivatives.
**       p_tvpart -- t-statistics for partial derivatives.
**       exp(b)   --
**
**********************************************************************
** GLOBAL VARIABLES (defalt value).      PROBLOG.DEC
**    variable type
**********************************************************************
**  GLOBSET is a procedure that resets the globals to their
**          respective default values.
**
**  GLOBSHOW is a procedure that shows the actual setting
**           for the globals
**
**
**       _corr(0)  -- option for correlation matrix of the coeffs.
**         scalar     if _corr=1, covb./std./std' will be printed
**
**       _cross(0) -- option for crosstable (need crosstab.arc -in
**         scalar     the bstat module)
**
**       _inter(1) -- option for intermediate results (#of iterations,
**         scalar     converg. criteria and its value at converg., and
**                    converg. tolerance).
**
**       _diagn(0) -- option for diagnostic results (iter,grad,etc.)
**         scalar     at each iteration.
**
**       _marg(1)  -- option for computation of marginal effect analysis.
**         scalar     (related: _xstar)
**
**       _xstar(0) -- data point at which the marginal effects analysis
**         matrix     is carried out.
**                    _xstar=0 --> data point is the mean of the Xs
**                    _xstar=x[j,.]' data point is the jth observation
**
**     __altnam(0) -- user defined variable names. Specify dependent
**        matrix      variable first.
**                    If __altnam is specified the routine will use the
**                    names contained in __altnam, e.g., let __altnam=
**                    grade constant gpa tuce psi. Use stof("..names..")
**                    if you want case sensitivity.
**                    If __altnam is not specified the routine will
**                    provide standard names for the dependent variables,
**                    e.g., Y X00...XK.
**
**       _headr(1) -- print header for the output
**        scalar
**
**       _fname(0) -- dataset file name (printed in the header)
**        string
**
**       _pname(0) - program file name
**        string
**
**      _ouname(0) - output file name
**       string
**
**        _tst(1)  - determine whether the output from "ess" to "pwald"
**        scalar     will be computed and printed (see above)
**
**        _toler   - tolerance for convergence - default is 1e-06
**        scalar
**
**      _maxit(25) - maximum number of iterations
**       scalar
**
**      _rwald(0)  - matrix of restricions on the coefficients to be
**       matrix      used in the Wald test. The test is of the form:
**                   (r*b-q)'inv(r*covb*r')*(r*b-q).
**                   If both _rwald and _qwald are set to zero (default),
**                   no test will be performed. Otherwise, both the
**                   more general test (b=0) and the chosen one will be
**                   performed and printed.
**                   Example of restrictions - the last 2 coeffs = 0:
**                   let _rwald = {0 0 1 0,
**                                 0 0 0 1};
**                   let _qwald =  0 0;
**
**      _qwald(0)  - vector of values implied by the restrictions
**       matrix
**
**      _stval(0)  - vector of starting values. The default is OLS
**       matrix      estimates. This can be changed by selecting
**                   another set of starting value. However, the
**                   program seems to work quite well with OLS.
**
**
**
** Examples:
**           call problog(y,x,"probit")
**                            will perform probit estimation and will not
**                            save output.
**
**           {b,covb,lkdu}=problog(y,x,"logit")
**                            will perform logit estimation and will save
**                            output.
**
*/

proc (3) = problog(y,x,_modl);

local _version, iter, crit, t, k, b, pdf, cdf, g, d, h, db, bn, depnam,
      _v, pn, s, ofn1, ofn2, covb, std, tv, p_tv, lkdu, ratio, lkdr, lr,
      plr_v, waldu, waldr, pwaldu, pwaldr, mm, sig, ess, tss, rss, r2,
      r2lai, r2mf, varstr, mask, fmat, final, prfit, pp, names, xstar,
      dens, part, scalfac, f, covpart, stdpart, tvpart, p_tvpart,
      partfin, _model, oldfmt;


_modl=lower(_modl);
if _modl $== "probit";
  _model=0;
elseif _modl $== "logit";
  _model=1;
elseif _modl $/= "probit" .or "logit";
  errorlog
  "===>> ERROR: --->_modl="; chrs(34) upper(_modl) chrs(34);;
               " _MODL MUST BE EITHER \34PROBIT\34 OR \34LOGIT\34";
                 end;
endif;


/** Print header **/
if _headr==1;
     __header="tldvf";
     _version=1.02|1;
     if _model==0;
          header("PROBLOG.ARC  -  probit model ",_fname,_version);
     elseif _model==1;
          header("PROBLOG.ARC  -  logit model ",_fname,_version);;
     endif;
endif;

/** Print program file and output file **/
if _pname $/= "";
     print " PROGRAM   FILE: ";; upper(_pname);
endif;
if _ouname $/= "";
     print " OUTPUT    FILE: ";; upper(_ouname);
endif;


/** Define constants and initial values (ols default) **/
t    = rows(x);
k    = cols(x);
iter = 1;
crit = 1;

if _stval == 0 and rows(_stval) == 1;
  b=y/x; "";
  " _stval = 0 --> Starting values are OLS (default)";
else; "";
  " Starting values selected by the user: ";_stval';
  if rows(_stval) /= k;
    errorlog "_stval is not conformable with the regressors matrix
             STARTING VALUES RESET TO OLS VALUE";
    pause(2);
    _stval=y/x;
  endif;
    b=_stval;
endif;


/** BEGIN LOOP **/

do until (iter gt _maxit) or (crit le _toler);

     if _model==0;  /*probit model*/

        pdf=pdfn(x*b);                  /** Normal pdf **/
        cdf=cdfn(x*b);                  /** Normal cdf **/

        /** Gradient vector  (g is the same for probit and logit) **/
        g = y.*(pdf./cdf).*x-(1-y).*(pdf./(1-cdf)).*x;
        g = sumc(g);
        /** Hessian matrix -eq. 19.2.21- H=-X'DX where D is diagonal **/
        D = pdf.*((y.*(pdf+(x*b).*cdf)./cdf^2)
            +((1-y).*(pdf-(x*b).*(1-cdf))./(1-cdf)^2));
        H = -(x.*D)'x;

     elseif _model==1;  /*logit model*/

        pdf = exp(-x*b)./(1+exp(-x*b))^2;   /** Logit pdf **/
        cdf = 1./(1+exp(-x*b));             /** Logit cdf **/

        /** Gradient vector **/
        g   = y.*(pdf./cdf).*x-(1-y).*(pdf./(1-cdf)).*x;
        g   = sumc(g);
        /** Hessian matrix -eq.19.2.22- **/
        H   = -(x.*pdf)'x;

     endif;

     db   = - inv(H)*g;    /** compute Newton step **/
          gosub  step;
     bn   = b + s*db;
     crit = maxc(maxc(abs(db)));


     if _diagn==1;      /** diagnostic results option **/
          gosub prnta;
     elseif _diagn==0;
     endif;

     b    = bn;           /** update b **/
     iter =iter+1;        /** update iteration **/

endo;

/** END LOOP **/


@ Print final results @  "";
if _diagn==1;
  "Probit Estimation (programmed analytical deriv.) - Final Results: ";
elseif _diagn==0;
  ">>>>> _diagn = 0   No diagnostic output option selected <<<<<";
endif;
if _marg==0;
  ">>>>> _marg  = 0   No marginal effects  option selected <<<<<";
endif;

if __altnam $/=0;             /** get name dependent variable **/
     depnam=__altnam[1];
else;
     depnam="Y";
endif;


if _model==0;
     "Binary Probit Variable     :  ";;format /ld 8,0; $depnam;
elseif _model==1;
     "Binary Logit Variable      :  ";;format /ld 8,0; $depnam;
endif;


gosub results;

format /rd 2,0;
"Number of observations     :  ";;t;;"  ";;format /ld strlen(depnam),0;
$depnam "= 0 -->" pn[1] "   ";; $depnam "= 1 -->" pn[2]; format 2,5;
"Percent of positive obs    : ";; 100*pn[2]/t;
if _inter==1;
 format 2,0;
 "Convergence achieved in    : ";;iter-1;; "iterations (max = " _maxit ")";
 "Convergence criteria       :  maxc(abs(db)) where db = -inv(H)*g";
 "   value at convergence    :  ";;format /m1 /ros 5,5;crit;
 "Convergence tolerance      :  ";;_toler;
 format 2,5;
elseif _inter==0;
endif;


if _tst==1;
  print chrs(45*ones(80,1));;

  if _rwald == 0 and _qwald == 0;

  print ftos(ess,"Error sum of squares       : %lf",10,6);;
  print ftos(tss,"       Total sum of squares  :%lf",10,6);
  print ftos(rss,"Residual sum of squares    : %lf",10,6);;
  print ftos(r2,"       R\253                    :%lf",10,6);
  print ftos(lkdu,"Unrestricted log-lik       : %lf",10,6);;
  print ftos(r2lai,"       R\253 - Laitila '93      :%lf",10,6);
  print ftos(lkdr,"Restr.log-lik(slopes=0)    : %lf",10,6);;
  print ftos(r2mf,"       R\253 - McFadden         :%lf",10,6);
  print ftos(lr,  "Likelihood Ratio=2(LU-LR)  : %lf",10,6);;
  print ftos(waldu,"       Wald test (\225=0)       :%lf",10,6);
  print ftos(plr_v,"   p-value                 : %lf",10,6);;
  print ftos(pwaldu,"          p-value            :%lf",10,6);
  print ftos(k-1 ,"   d.f.                    : %lf",10,6);;
  print ftos(k   ,"          d.f.               :%lf",10,6);
  print chrs(45*ones(80,1));;


  else;

  print ftos(lkdu,"Unrestricted log-lik       : %lf",10,6);;
  print ftos(waldr,"       Wald test (w/restrs.) :%lf",10,6);
  print ftos(lkdr,"Restr.log-lik(slopes=0)    : %lf",10,6);;
  print ftos(pwaldr,"          p-value            :%lf",10,6);
  print ftos(lr,  "Likelihood Ratio=2(LU-LR)  : %lf",10,6);;
  print ftos(rows(_qwald)   ,"          d.f.               :%lf",10,6);
  print ftos(plr_v,"   p-value                 : %lf",10,6);;
  print ftos(waldu,"       Wald test (\225=0)       :%lf",10,6);
  print ftos(k-1 ,"   d.f.                    : %lf",10,6);;
  print ftos(pwaldu,"          p-value            :%lf",10,6);
  print ftos(r2lai,"R\253 - Laitila '93           : %lf",10,6);;
  print ftos(k   ,"          d.f.               :%lf",10,6);
  print ftos(r2mf,"R\253 - McFadden              : %lf",10,6);
  print ftos(r2,"R\253                         : %lf",10,6);;
  print ftos(tss,"       Total sum of squares  :%lf",10,6);
  print ftos(rss,"Residual sum of squares    : %lf",10,6);;
  print ftos(ess,"       Error sum of squares  :%lf",10,6);

  print chrs(45*ones(80,1));;

  endif;

elseif _tst==0;

  print chrs(45*ones(80,1));;

endif;

if _model==0;
     "variable      \225_ols     \225_probit     std.err.     t-val"\
     "    prob|t|\242x   gradient";
elseif _model==1;
     "variable      \225_ols     \225_logit      std.err.     t-val"\
     "    prob|t|\242x   gradient";
endif;

print chrs(45*ones(80,1));;

gosub prntb;"";     /* print final results */

if _marg==1;        /* marginal effects analysis option */
     call derivati(b,x,covb,_model);
elseif _marg==0;
endif;

@ compute and print correlation matrix of the coefficients @
if _corr==1;
     if _model==0;
          matwrt("Correlation matrix of the probit coefficients",
          (covb./std./std'),varstr,varstr,3);
     elseif _model==1;
          matwrt("Correlation matrix of the logit coefficients",
          (covb./std./std'),varstr,varstr,3);
     endif;
elseif _corr==0;
endif;


@ computes and print crosstable @
if _cross==1;
     gosub crtab;
endif;


print chrs(42*ones(80,1));;
if _model==0;
">>>>>>>>>>>>>>>>>>>>>>>>>>> end of probit estimation <<<<<<<<<<<<<<<<"\
"<<<<<<<<<<<";;
elseif _model==1;
">>>>>>>>>>>>>>>>>>>>>>>>>>> end of logit  estimation <<<<<<<<<<<<<<<<"\
"<<<<<<<<<<<";;
endif;
print chrs(42*ones(80,1));;

retp(b,covb,lkdu); /** Values available for further computation **/
"";

/***********************************************************************/
/**************************** SUBROUTINES ******************************/
/***********************************************************************/

@=================================================@
@================== step length ==================@
@=================================================@
@ begin subroutine step @
step:

s = 2;
ofn1 = 0;
ofn2 = 1;
do until ofn1 >= ofn2;
     s = s/2;
     ofn1=llikfun(b+s*db,y,x,_model);
     ofn2=llikfun(b+s*db/2,y,x,_model);
endo;

return;
@ end subroutine step @

@=================================================@
@========= print values at each iteration ========@
@=================================================@
@ begin subroutine prnta @
prnta:

format 2,0; "ITER. NR. " iter;;
format 4,2; " Step length = " s;;
format 5,15; " Likelihood = " llikfun(bn,y,x,_model);
format 3,10; "gradient ";;g';
format /m1 /ros 5,5;"crit.    ";;crit;
format /ld 5,5; "Estimates = " bn';

return;
@ end subroutine prnta @

@=================================================@
@================= final results =================@
@=================================================@
@ begin subroutine results @
results:

let _v=0 1;
pn=counts(y,_v);           /** compute # obs.               **/

covb = -inv(H);            /** asymptotic covariance matrix **/
std  = sqrt(diag(covb));   /** standard errors              **/
tv   = b./std;             /** t-values                     **/
p_tv = 2*(cdfnc(abs(tv))); /** p-values (use normal distr.) **/

lkdu=llikfun(b,y,x,_model);/** unrestriced log-lik function **/

if _model==0;              /** fitted probabilities         **/
     prfit=cdfn(x*b);              /*probit*/
elseif _model==1;
     prfit = 1./(1+exp(-x*b));     /*logit*/
endif;

if _tst==1;

ratio=pn[2]/t;        /** restricted log-lik. and lr test   **/
lkdr = t*(ratio*ln(ratio)+(1-ratio)*ln(1-ratio));   /*::::*/
lr   = 2*(lkdu-lkdr);          /* likelihood ratio           */
plr_v= cdfchic(LR,k-1);        /* p-value                    */

if _rwald ==0 and _qwald == 0;
     waldu = -b'H*b;  /** overall significance of the mod. **/
     pwaldu= cdfchic(waldu,k);
else;
     waldu = -b'H*b;  /** overall significance of the mod. **/
     pwaldu= cdfchic(waldu,k);
     waldr=(_rwald*b-_qwald)'inv(_rwald*covb*_rwald')*(_rwald*b-_qwald);
     print "Matrices of restrictions for the Wald test:
     _rwald~_qwald";
     oldfmt=formatnv("*.*lf"~6~2);
     call printfmt(_rwald~_qwald,1);"";
     pwaldr= cdfchic(waldr,rows(_qwald));
endif;

/**** R2 measures ****/          /* prfit defined above  */
ess=(y-prfit)'(y-prfit);         /* error sum of squares */
tss=(y-meanc(y))'(y-meanc(y));   /* total sum of squares */
rss=tss-ess;                     /* resid sum of squares */
r2=1-(ess/tss);                  /** "classical" R2      **/

r2mf = 1-(lkdu/lkdr);            /** McFadden pseudo-R2  **/

if _model==0;        @probit@    /** Laitila '93 - R2    **/
     sig=1;
elseif _model==1;    @logit @
     sig=(pi^2)/3;
endif;
mm=vcx(x);
r2lai = (b'mm*b)/(sig+b'mm*b);

endif;   /* (_tst==1) */

return;
@ end subroutine results @


@=================================================@
@========= format and print final results ========@
@=================================================@
@ begin subroutine prntb @
prntb:

if __altnam $/=0;        /*look for user defined var.names    */
     varstr=__altnam[2:k+1];
else;                    /*create a column of X00..X0K        */
     varstr=0 $+ "X" $+ ftocv(seqa(0,1,(rows(b))),2,0);
endif;

mask=0~ones(1,6);
let fmat[7,3] = "-*.*s" 8 8  "*.*lf" 12 4 "*.*lf" 12 5 "*.*lf" 12 5
                "*.*lf" 12 5 "*.*lf" 10 4 "*.*le" 13 3;

final=varstr~(y/x)~b~std~tv~p_tv~g;
call printfm(final,mask,fmat);
print chrs(45*ones(80,1));;

return;
@ end subroutine prntb @

@=================================================@
@========= compute and print a crosstable ========@
@=================================================@
@ begin subroutine crtab @
crtab:

call varput(0,"_ctstat");  /* set globals for the crosstab proc */
call varput(0,"_ctpause");

pp = prfit .> 0.5;       /*prfit defined above*/
let names=actual predicted;
call saved(y~pp,0,names);
"==========================================";
"FREQUENCIES OF ACTUAL & PREDICTED OUTCOMES";
"PREDICTED OUTCOME HAS MAXIMUM PROBABILITY.";
"==========================================";
call crosstab("temp",names);

return;
@ end subroutine crtab @


         /**************************************************/
endp;    /*********** END OF PROCEDURE PROBLOG.ARC *********/
         /**************************************************/

/***********************************************************************/
/******************** DECLARATION OF PROCEDURES ************************/
/***********************************************************************/

@=================================================@
@========= Probit/Logit log-lik. function ========@
@=================================================@
proc llikfun(b,y,x,_model);

  local cdf,li;

  if _model==0;                     /* probit cdf       */
       cdf=cdfn(x*b);

  elseif _model==1;                 /* logit cdf        */
       cdf=1./(1+exp(-x*b));

  endif;

  li=y.*ln(cdf)+(1-y).*ln(1-cdf);   /* log-lik function */

  retp(sumc(li));

endp;



@=================================================@
@============ marginal effect analysis ===========@
@=================================================@
proc derivati(b,x,covb,_model);

/**************************************
** Reference for  this  procedure is **
** Greene,  _Econometric  Analysis_, **
** 1993, New York, MacMillan, ch. 21 **
**************************************/

local k, varstr, xstar, dens, pp, part, scalfac, f, covpart,
      stdpart, tvpart, p_tvpart, mask, fmat, numres, partfin;

k = cols(x);

if __altnam $/=0;         /*look for user defined var.names   */
     varstr=__altnam[2:k+1];
else;                     /*create a column of X00..X0K       */
     varstr=0 $+ "X" $+ ftocv(seqa(0,1,(rows(b))),2,0);
endif;

if _xstar==0;             /*evaluation vector - mean is the   */
     xstar=meanc(x);      /*default data point                */
else;
     xstar=_xstar;        /* else, user supplied data point   */
endif;

/** Probit marginal analysis **/
if _model==0;
     dens    = pdfn(xstar'b);
     pp      = dens*(eye(rows(b))) - dens*((b'xstar)*(b*xstar'));
     part    = dens.*b ;     /* derivatives at the data point */
     scalfac = dens;         /* Probit scale factor           */

/** Logit marginal analysis **/
elseif _model==1;
     dens    = (exp(-xstar'b))/((1+exp(-xstar'b)));
     f       = dens*(1-dens);
     pp      = f*(eye(rows(b))) - (1-2*dens)*f*(b*xstar');
     part    = f*b;         /* derivatives at the data point  */
     scalfac = f;           /* Logit scale factor             */
endif;

/** Part in common **/
covpart  = pp*covb*pp';             /* covariance formula     */
stdpart  = sqrt(diag(covpart));     /* s.e. for partials      */
tvpart   = part./stdpart;           /* t-values for partials  */
p_tvpart = 2*(cdfnc(abs(tvpart)));  /* p-values for partials  */

mask=0~ones(1,6);
let fmat[7,3] = "-*.*s" 8,8 "*.*lf" 12,3 "*.*lf" 12,5 "*.*lf" 12,4
                "*.*lf" 12,4 "*.*lf" 12,4 "*.*lf" 11,5;

/** Print header **/

__header="ldf";

if xstar==meanc(x);

   if _model==0;
header("PROBIT MARGINAL EFFECT ANALYSIS (at the mean of the data)","",0);
   elseif _model==1;
header("LOGIT MARGINAL EFFECT ANALYSIS (at the mean of the data)","",0);
   endif;

else;    /* if xstar /= meanc(x) */

if _model==0;
header("PROBIT MARGINAL EFFECT ANALYSIS","",0);
elseif _model==1;
header("LOGIT MARGINAL EFFECT ANALYSIS","",0);
endif;

endif;

__header="tldvf";  /*reset __header*/

/** print scale factor, utility index, and final results **/
format /ld 5,5;
"                             scale factor  = ";;scalfac;
"                             utility index = ";;xstar'b;

print chrs(45*ones(80,1));;
"variable    data point    slope      std.err.     t-val. "\
"    prob|t| x   exp(b)";
print chrs(45*ones(80,1));;

numres=xstar~part~stdpart~tvpart~p_tvpart~exp(b);
partfin=varstr~numres;
call printfm(partfin,mask,fmat);
print chrs(45*ones(80,1));;

retp(numres);

endp;   /* END OF PROCEDURE DERIVATI */


@=================================================@
@=========== procedure to set globals ============@
@=================================================@
proc (0) = globset;
  __altnam  = 0;
     _corr  = 0;
     varput(0,"_cross");
     _inter = 1;
     _marg  = 1;
     _headr = 1;
     _tst   = 1;
     _toler = 1e-06;
     _maxit = 25;
     _fname = "";
     _rwald = 0;
     _qwald = 0;
     _diagn = 0;
     _xstar = 0;
     _pname = "";
     _ouname= "";
     _stval = 0;
endp;


@=================================================@
@======= procedure to show globals setting =======@
@=================================================@
proc (0) = globshow;
format /lz 8,5;
"-------------------------------------------------------";
"PROBLOG.ARC ******************** GLOBALS ACTUAL SETTING";
"-------------------------------------------------------";
"__altnam=  ";;$__altnam';
"_corr   = ";; _corr;
"_cross  = ";; _cross;
"_inter  = ";; _inter;
"_marg   = ";; _marg;
"_headr  = ";; _headr;
"_tst    = ";; _tst;
"_toler  = ";; _toler;
"_maxit  = ";; _maxit;
"_fname  =  ";; _fname;
"_rwald  = ";; _rwald;
"_qwald  = ";; _qwald;
"_diagn  = ";; _diagn;
"_xstar  = ";; _xstar';
"_pname  =  ";; _pname;
"_ouname =  ";; _ouname;
"_stval  = ";; _stval';
"-------------------------------------------------------";
endp;

