/*
ENDOGENOUS SWITCHING REGRESSION

FINALSWT.SRC
Written by Parke Wilde 6/22/98


CONTENTS:
PART 1.  Procedures for Probit estimation
PART 2.  Procedures for Heckman-Lee Two-Step estimation
PART 3.  Procedures for Maximum Likelihood Switching Regression estimation
PART 4.  Utilities for reading in data or making simulated data
PART 5.  Demonstration program using simulated data


FORMALITIES:
This code is written and submitted for public, non-commericial use.  Because code must often be modified for particular tasks, the writing format favors transparency over efficiency.  There are no performance guarantees.  I am to blame for any errors, but I accept no liability.  Use of this code is the user's responsibility alone.  I would love feedback.  My e-mail address is pwilde@econ.ag.gov.  If this e-mail address does not work, please write ckr2@cornell.edu to ask for a current address.  Please use this code freely.


PURPOSE:
This program estimates an endogenous switching regression model.  See Maddala, Limited-Dependent and Qualitative Variables in Econometrics (1983) and Greene, Econometric Analysis, 2nd Edition (1993).


OVERVIEW:
Estimates all parameters in a system of three equations:
  (1) y0=x0*beta0+u0,
  (2) y1=x1*beta1+u1,
  (3) y=d*y1+(1-d)*y0,
      where d=1 if z*gam+e>0, d=0 otherwise.
The stochastic specification is:
  (u0,u1,e)~N(0,SIGMA), where
  SIGMA= { sig0     0  sig0e,
              0  sig1  sig1e,
          sig0e sig1e      1}.

Estimation uses Maximum Liklihood Estimation Application Module (Aptech Systems Inc., 1996).  First, the procedure Maxswt calls the subroutine Hek to estimate beta0, beta1, sig0, sig0e, sig1, sig1e by the Heckman-Lee two-step procedure; second, it uses the returns from Hek as starting values to estimate the full model simultaneously by maximum likelihood.

The documentation here is for the use of Maxswt, but the subroutines Maxprb and Hek are included below and may be used independently.  After the procedures, the code below includes a sample program that creates a simulated data set and demonstrates Maxswt.  The sample program also indicates how to use this program with real data sets, either in the form of a GAUSS matrix or an ascii file.  If the maxlik library is installed on your computer, this whole text file should be ready to execute the demonstration (just hit "run").


FORMAT:
{betswt,f,g,vbetswt,retcode}=maxswt(ind);


GLOBALS:
Dat: a data matrix.  This matrix is not passed to subroutines.  Rather, to save memory, an index matrix is passed to subroutines, telling them where to find variables in the columns of Dat.

All globals used by Maxlik are also available.

Global variables "ynam," "xnam," etc., containing the names of the various dependent and independent variables, are used by the utility Makind below, although they are not needed for Maxswt.


INPUTS:
Ind: a matrix with 5 columns that must have exactly the following format.  Each column contains index numbers identifying the column(s) of dat in which each of the following variable types may be found: y x0 x1 z d.  If the indices in a column of Ind do not take up the whole column, the ramaining rows of that column should contain zeros.  The number of rows in Ind equals the number of variables in whichever variable type has the largest number of variables (usually z).  The procedure Makind below conveniently creates the matrix Ind, so most users will not need to worry about the structure of this matrix.


OUTPUTS:
Outputs are in the same format as output from Maxlik.

betswt: a vector of parameter estimates in the following order:
b0 | b1 | gam | sig0 | sig0e | sig1 | sig1e.

f: value of the log-likelihood function.

g: gradient.

vbetswt: the asymptotic estimated variance-covariance matrix for betswt.

retcode: return code (0 for normal convergence).
*/


/* Refer to Maximum Likelihood Routines and Initialize Dat */
new;
library maxlik;
#include maxlik.ext;
maxset;
output file=tstswt.out reset;
dat=0;


/*  PART 1.  Probit estimation */
/*  pwprb.src */
/*  procedures needed for estimating probit */
/*  Parke Wilde 6-18-98 */

proc betols(y,x);
  /* ols beta vector */
retp(invpd(x'*x)*x'*y);
endp;

proc likprb(b,z);
    /* Returns the log-likelihood */
    /* The last column of matrix z is the binary dep. variable */
    local d, cdf;
    d=z[.,cols(z)];
    cdf = cdfn(z[.,1:(cols(z)-1)]*b);
    retp(d.*ln(cdf) + (1-d).*ln(1-cdf));
endp;

proc grdprb(b,z);
   /* Returns the gradient */
   /* See Greene p. 644 eq. 21-21 */
   local x, d;
    x=z[.,1:(cols(z)-1)];
    d=z[.,cols(z)];
    retp(lambda(d,x*b).*x);
endp;

proc hesprb(b,z);
   /* Returns the hessian. */
   /* See Greene p. 645 eq. 21-23 */
   local m,d,x;
    d=z[.,cols(z)];
    x=z[.,1:(cols(z)-1)];
    m = x*b;
    retp(-(x.*lambda(d,m).*(lambda(d,m)+m))'x);
endp;

proc (5)=maxprb(d,x);
local bstart ;
  bstart=betols(d,x);
  __output = 0;
  _max_GradProc = &grdprb;
  _max_HessProc=&hesprb;
  _max_Algorithm = 4;
  _max_GradCheckTol = 1e-3;
retp(maxlik(x~d,0,&likprb,bstart));
endp;


/* PART 2.  Heckman-Lee Two-Step estimation */
/* mills */
/* calculates inverse mills ratio for d=1 */
proc mills(zgam);
retp(-pdfn(zgam)./cdfnc(zgam));
endp;

/* millsc*/
/* calculates inverse mills ratio for d=0 */
proc millsc(zgam);
retp(pdfn(zgam)./cdfn(zgam));
endp;

proc lambda(d,m);
    /* gives correct lambda under each regime */
    retp(d.*millsc(m)+(1-d).*mills(m));
endp;

/* Two-Step Heckman Lee Switching Regression */
proc (7)=hek(y,x0,x1,z,d);
  local gam, f, g, cov, retcode,
        alpha, b0, b1,ystar0,xstar0,ystar1,xstar1,var0,var1,
        dbar0,dbar1,sig0,sig0e,sig1,sig1e;

  /* call probit first */
  {gam,f,g,cov,retcode}=maxprb(d,z);
  call maxprt(gam,f,g,cov,retcode);

  /* call betols for regression using mills ratios */
  alpha=z*gam;
  ystar0=delif(y,d);
  xstar0=delif(x0~millsc(alpha),d);
  b0=betols(ystar0,xstar0);
  var0=meanc((ystar0-xstar0*b0)^2);
  ystar1=selif(y,d);
  xstar1=selif(x1~mills(alpha),d);
  b1=betols(ystar1,xstar1);
  var1=meanc((ystar1-xstar1*b1)^2);

  /* Standard deviations for selection. */
  /* Not corrected for heterosk'ty, but adequate for starting values. */
  /* See Greene 2nd edition p. 711-713 */
  dbar0=meanc( -(millsc(alpha).*(alpha+millsc(alpha)))   );
  dbar1=meanc( -(mills(alpha) .*(alpha+mills(alpha )))   );
  sig0=sqrt(var0-(b0[rows(b0)]^2)*dbar0);
  sig0e=b0[rows(b0)];
  sig1=sqrt(var1-(b1[rows(b1)]^2)*dbar1);
  sig1e=b1[rows(b1)];

  print "Heckman regression results (b0, b1): " b0~b1;
  print "Standard errors (sig0, sig1): " sig0~sig1;
  print "Heckman covariances (sig0e, sig1e): " sig0e~sig1e;
retp(b0[1:rows(b0)-1],b1[1:rows(b0)-1],gam,sig0,sig0e,sig1,sig1e);
endp;


/*  PART 3.  Maximum Likelihood Switching Regression estimation */
proc nonzero(invec);
/* retrieves nonzero elements of a column of the index matrix Ind */
retp(invec[1:sumc(invec.>0)]);
endp;

proc likswt(b,ind);
    /* Calculates Log-likelihood function.  */
    /* See Appendix A of PW dissertation. */
    /* See Maddala p. 224 and p. 284. */
    /* b is structured: b0 b1 gam sig0 sig0e sig1 sig1e */

    local mark,b0,b1,gam,sig0,sig0e,sig1,sig1e,
          y,x0,x1,z,d,u0,u1,q0,q1,r0,r1;

    /* make local variables for elements of beta parameter */
    mark=0;
    b0=   b[mark+1:mark+sumc(ind[.,2].>0)];  mark=mark+sumc(ind[.,2].>0);
    b1=   b[mark+1:mark+sumc(ind[.,3].>0)];  mark=mark+sumc(ind[.,3].>0);
    gam=  b[mark+1:mark+sumc(ind[.,4].>0)];  mark=mark+sumc(ind[.,4].>0);
    sig0= b[mark+1]; mark=mark+1;
    sig0e=b[mark+1]; mark=mark+1;
    sig1= b[mark+1]; mark=mark+1;
    sig1e=b[mark+1]; mark=mark+1;

    /* data types are y x0 x1 z d */
    /* make local variables for indices for columns of matrix dat */
    y =ind[1:sumc(ind[.,1].>0),1];
    x0=ind[1:sumc(ind[.,2].>0),2];
    x1=ind[1:sumc(ind[.,3].>0),3];
    z =ind[1:sumc(ind[.,4].>0),4];
    d =ind[1:sumc(ind[.,5].>0),5];

    /* create variables used in calculating likelihood */
    u0=(dat[.,y]-dat[.,x0]*b0) ;
    u1=(dat[.,y]-dat[.,x1]*b1) ;
    r0=sqrt(1-sig0e^2/sig0^2) ;
    r1=sqrt(1-sig1e^2/sig1^2) ;
    q0=(( dat[.,z]*gam+(sig0e/sig0^2)*u0)/r0);
    q1=(( dat[.,z]*gam+(sig1e/sig1^2)*u1)/r1);
retp( (1-dat[.,d]).*(-ln(sig0)+lnpdfn(u0/sig0)+ln(cdfnc(q0))) +
         dat[.,d] .*(-ln(sig1)+lnpdfn(u1/sig1)+ln(cdfn(q1))) );
endp;

proc grdswt(b,ind);
    /* Calculates Gradient */
    local mark,b0,b1,gam,sig0,sig0e,sig1,sig1e,
          y,x0,x1,z,d,u0,u1,q0,q1,r0,r1,
          lam0,lam1,dq0db0,dq1db1,dq0dgam,dq1dgam,dq0ds0e,dq1ds1e,
          dq0ds0,dq1ds1;

    /* b is structured: b0 b1 gam sig0 sig0e sig1 sig1e */
    /* make local variables for elements of beta parameter */
    mark=0;
    b0=   b[mark+1:mark+sumc(ind[.,2].>0)];  mark=mark+sumc(ind[.,2].>0);
    b1=   b[mark+1:mark+sumc(ind[.,3].>0)];  mark=mark+sumc(ind[.,3].>0);
    gam=  b[mark+1:mark+sumc(ind[.,4].>0)];  mark=mark+sumc(ind[.,4].>0);
    sig0= b[mark+1]; mark=mark+1;
    sig0e=b[mark+1]; mark=mark+1;
    sig1= b[mark+1]; mark=mark+1;
    sig1e=b[mark+1]; mark=mark+1;

    /* data types are y x0 x1 z d */
    /* make local variables for indices for columns of matrix dat */
    y =ind[1:sumc(ind[.,1].>0),1];
    x0=ind[1:sumc(ind[.,2].>0),2];
    x1=ind[1:sumc(ind[.,3].>0),3];
    z =ind[1:sumc(ind[.,4].>0),4];
    d =ind[1:sumc(ind[.,5].>0),5];
    u0=(dat[.,y]-dat[.,x0]*b0) ;
    u1=(dat[.,y]-dat[.,x1]*b1) ;

    r0=sqrt(1-sig0e^2/sig0^2) ;
    r1=sqrt(1-sig1e^2/sig1^2) ;
    q0=(( dat[.,z]*gam+(sig0e/sig0^2)*u0)/r0);
    q1=(( dat[.,z]*gam+(sig1e/sig1^2)*u1)/r1);
    lam0=(-pdfn(q0))./cdfnc(q0);
    lam1=( pdfn(q1))./ cdfn(q1);

    /* Derivatives of Q0 and Q1 */
    d=dat[.,d];
    dq0db0 =(1-d).*(-(sig0e*dat[.,x0])/(r0*sig0^2));
    dq1db1 =    d.*(-(sig1e*dat[.,x1])/(r1*sig1^2));
    dq0dgam=(1-d).*(( dat[.,z])/r0);
    dq1dgam=    d.*(( dat[.,z])/r1);
    dq0ds0 =(1-d).*(((-2*sig0e*(u0.*r0))-(q0*sig0e^2))/((r0^2)*(sig0^3)));
    dq0ds0e=(1-d).*(( u0.*r0+q0*sig0e)/((r0^2).*(sig0^2)));
    dq1ds1 =    d.*(((-2*sig1e*(u1.*r1))-(q1*sig1e^2))/((r1^2)*(sig1^3)));
    dq1ds1e=    d.*(( u1.*r1+q1*sig1e)/((r1^2).*(sig1^2)));

    /* Gradients: b0,b1,gam,sig0,sig0e,sig1,sig1e   */
    retp((1-d).*((u0.*dat[.,x0]/(sig0^2)) + (lam0.*dq0db0))  ~
             d.*((u1.*dat[.,x0]/(sig1^2)) + (lam1.*dq1db1))  ~
         (1-d).*(lam0.*dq0dgam) + d.*(lam1.*dq1dgam) ~
         (1-d).*((-1/sig0)+((u0^2)/sig0^3)+lam0.*dq0ds0) ~
         (1-d).*(lam0.*dq0ds0e) ~
             d.*((-1/sig1)+((u1^2)/sig1^3)+lam1.*dq1ds1) ~
             d.*(lam1.*dq1ds1e)  );
endp;

proc (5) = maxswt(ind);
  local startvec,b0,b1,gam,yind,x0ind,x1ind,zind,dind,sig0,sig0e,sig1,sig1e;
  yind=nonzero(ind[.,1]);
  x0ind=nonzero(ind[.,2]);
  x1ind=nonzero(ind[.,3]);
  zind=nonzero(ind[.,4]);
  dind=nonzero(ind[.,5]);

  /* Get starting values */
  {b0,b1,gam,sig0,sig0e,sig1,sig1e}=
  hek(dat[.,yind],dat[.,x0ind],dat[.,x1ind],dat[.,zind],dat[.,dind]);
  startvec=b0|b1|gam|sig0|sig0e|sig1|sig1e;
  print "rows startvec" rows(startvec);
  /* print "startvec" startvec; */

  /* Call Maxlik */
  maxset;
  /* __output = 20; */
  _max_GradCheckTol=1e-3;
  print "Maxlik checking analytic gradient against numerical.";
  _max_GradProc=&grdswt;
  print "Maxlik using gradients in grdswt.";
  _max_Algorithm=4;
  print "max algorithm" _max_Algorithm;
  /* _max_ParNames= */
retp(maxlik(ind,0,&likswt,startvec));
endp;


/* PART 4.  Utilies for reading in data or making simulated data */
/* readasc reads in an ascii file with names in the first row */
/* readmat reads in a GAUSS matrix with names in the first row */
/* makind makes the index matrix "Ind" */
/* swtsim makes simulated data for switching regression estimation */

proc(2)=readin(dat,datrows,datcols);
  local datnams;
  print;
  "                    READ IN DATA";
  if rows(dat)*cols(dat)==(datrows+1)*datcols;
    print "Successfully read in data set of the expected size.";
  else;
    print "CAUTION: Data set is not of the expected size.";
    print "expected number of elements" (datrows+1)*datcols;
    print "total number of elements" rows(dat);
  endif;
  if ismiss(dat);
    print "CAUTION: Missing values were encountered.";
  else;
    print "No missing values were encountered.";
  endif;
  dat=reshape(dat,datrows+1,datcols);
  datnams=dat[1,.]';
  dat=dat[2:rows(dat),.];
  print;
  "observations (n):       " datrows;
  "variables (k):         " datcols;
retp(dat,datnams);
endp;

proc (2)=readasc(datname,datrows,datcols);
  local dat;
  load dat[]=^datname;
  retp(readin(dat,datrows,datcols));
endp;

proc (2)=readmat(datname,datrows,datcols);
  local dat;
  load dat=^datname;
  retp(readin(dat,datrows,datcols));
endp;

proc makind(dat,datnams,types);
  /* uses globals xnam ynam, etc., if types include x and y */
  local ind, count, thisnam, thisstr, thisindx;
  ind=zeros(cols(dat),rows(types));
  count=0;
  do until count==rows(types);
    count=count+1;
    thisnam=types[count];
    thisstr=varget(thisnam$+"nam");
    thisindx=indcv(thisstr,datnams);
    ind[1:rows(thisindx),count]=thisindx;
  endo;
retp(ind);
endp;

proc (2) = swtsim(trueb0,trueb1,truegam,truesig,nobs);
  local err, x, ystar0, ystar1, ystar, d, y, alldat, allname, intcept;
  err=rndn(nobs,3)*chol(truesig);
  x=rndu(nobs,2);
  intcept=ones(nobs,1);
  ystar0=(intcept~x[.,1])*trueb0+err[.,1];
  ystar1=(intcept~x[.,1])*trueb1+err[.,2];
  ystar=(intcept~x)*truegam+err[.,3];
  print "cov" (err'*err)/rows(err);
  print "corr" corrx(err);
  d=ystar .> 0;
  y=(1-d).*ystar0+d.*ystar1;
  alldat=intcept~x~y~d;
  let allname=intcept x1 x2 Y D;
retp(alldat,allname);
endp;


/* PART 5.  Demonstration program */

/* Make Simulated Data */
trueb0={10,20};
trueb1={15,25};
truegam={5,-5,-5};
truesig = {1 0 .3, 0 1 .1, .3 .1 1};
nobs=1000;
{alldat,allname}=
  swtsim(trueb0,trueb1,truegam,truesig,nobs);
tstswt=allname'|alldat;
save tstswt;

/* Read Data from GAUSS Matrix */
let datname=tstswt;
datrows=1000;
datcols=5;
{dat,datnams}=readmat(datname,datrows,datcols);

/* Call makind to make index matrix "Ind" */
let types=y x0 x1 z d;
let ynam= y ;
let x0nam= intcept x1;
let x1nam= intcept x1;
let znam= intcept x1 x2;
let dnam= d ;
ind=makind(dat,datnams,types);

/* Call maxswt and print out maximum likelihood estimation results */
call maxprt(maxswt(ind));

/* Print "true" parameters for comparison */
print;
print "True parameters, for comparison to maximum likelihood estimates";
"True b0:   " trueb0;
"True b1:   " trueb1;
"True bgam: " truegam;
"True sig0 and sig0e: " truesig[1,1]|truesig[1,3];
"True sig1 and sig1e: " truesig[2,2]|truesig[2,3];
end;

