* Version 1.2 28 Nov 2012
* By Joao Santos Silva

program define ivppml, eclass 
version 11.1
syntax varlist(numeric min=1) [if] [in] [aweight pweight fweight iweight], INSTruments(string) [STrict TWOstep ONEstep Igmm ///
vce(string) TECHnique(string) WMATrix(string) Center WINITial from(string)   tracelevel(string)  ///
conv_maxiter(integer 5000) conv_ptol(real 1e-6) conv_vtol(real 1e-6)  conv_nrtol(real 1e-6) ///
IGMMITerate(integer 5000) igmmeps(real 1e-6) igmmweps(real 1e-6) ]
marksample touse
tempname logy y _rhs _rhss _drop zeros _res instvs
gettoken y _rhs: varlist
unab instvs :  `instruments'


 qui gen `logy'=-222  if (`touse')                                                                             // Creates regressand for first step
 qui replace `logy'=log(`y') if (`touse')&(`y'>0)                                                              // Modifies regressand for first step
 qui reg `logy' `instvs' if (`touse')&(`y'>0)  [`weight'`exp']                                                   // Performs first step
 gen `zeros'=1                                                                                                 // Initialize observations selector
 local tot=1                                                                                                   // Initialize parameter counter
 local _drop ""                                                                                                // List of regressors to exclude
 local _rhss ""                                                                                                // List of regressors to include
 qui summarize `y' if (`touse')                                                                       
 if r(max)>1000000 di as error"`y' has very large values, consider rescaling"                                  // Warn if there are large values
    foreach x of local instvs{                                                                                   // Start loop over all regressors: LOOP 1.1
    local tot = `tot'+1                                                                                       // Add to parameter counter
      if (_se[`x']==0) {                                                                                       // Try to include regressors dropped in the
                                                                                                               // first stage: LOOP 1.1.1
          
          qui summarize `x' if (`y'>0)&(`touse'), meanonly
          local _mean=r(mean)
          qui summarize `x' if (`y'==0)&(`touse')
          if (r(min)<`_mean')&(r(max)>`_mean')&("`strict'" == "") {                                            // Include regressor if conditions met and
              local _rhss "`_rhss' `x'"                                                                        // strict is off 
          }
          else{
              qui su `x' if `touse', d                                                                         // Otherwise, drop regressor
              local _mad=r(p50)
              qui inspect  `x'  if `touse'                                                                         
              qui replace `zeros'=0 if (`x'!=`_mad')&(r(N_unique)==2)&(`touse')                                // Mark observations to drop
              local _drop "`_drop' `x'"
          }
      }                                                                                                        // End of LOOP 1.1.1 
      if _se[`x']>0 {                                                                                          // Include safe regressors: LOOP 1.1.2
      local _rhss "`_rhss' `x'" 
      }                                                                                                        // End LOOP 1.1.2
    }                                                                                                          // End LOOP 1.1
 qui su `touse' if `touse'                                                                                     // Summarize touse to obtain N
 local _enne=r(sum)                                                                                               // Save N
 *qui replace `touse'=0 if (`zeros'==0)&("`keep'"=="")&(`y'==0)&(`touse')                                       // Drop observations with perfect fit
                                                                                                               // if keep is off
if ("`_drop'" != "") {
di ""
di as result "Dropped instruments: `_drop'"  
}
if "`igmm'"=="" {
if "`technique'"==""{
gmm (`y' - exp({xb: `_rhs'} + {b0})) if `touse' [`weight' `exp'], instruments(`_rhss')  `twostep' `onestep'  ///
technique(nr) vce(`vce') `wmatrix' `center' `winitial' from(`from') tracelevel(`tracelevel')  conv_maxiter(`conv_maxiter') ///
conv_ptol(`conv_ptol') conv_vtol(`conv_vtol')  conv_nrtol(`conv_nrtol') ///
deriv(/xb = -1*exp({xb:} + {b0})) deriv(/b0 = -1*exp({xb:} + {b0})) 
}
else{
gmm (`y' - exp({xb: `_rhs'} + {b0})) if `touse' [`weight' `exp'], instruments(`_rhss')  `twostep' `onestep'  ///
technique(`technique') vce(`vce') `wmatrix' `center' `winitial' from(`from') tracelevel(`tracelevel')  conv_maxiter(`conv_maxiter') ///
conv_ptol(`conv_ptol') conv_vtol(`conv_vtol')  conv_nrtol(`conv_nrtol') ///
deriv(/xb = -1*exp({xb:} + {b0})) deriv(/b0 = -1*exp({xb:} + {b0})) 
}
}
else{
if "`technique'"==""{
gmm (`y' - exp({xb: `_rhs'} + {b0})) if `touse' [`weight' `exp'], instruments(`_rhss')  igmm  ///
technique(nr) vce(`vce') `wmatrix' `center' `winitial' from(`from') tracelevel(`tracelevel')  conv_maxiter(`conv_maxiter') ///
conv_ptol(`conv_ptol') conv_vtol(`conv_vtol')  conv_nrtol(`conv_nrtol') ///
deriv(/xb = -1*exp({xb:} + {b0})) deriv(/b0 = -1*exp({xb:} + {b0})) 
}
else{
gmm (`y' - exp({xb: `_rhs'} + {b0})) if `touse' [`weight' `exp'], instruments(`_rhss')  igmm  ///
technique(`technique') vce(`vce') `wmatrix' `center' `winitial' from(`from') tracelevel(`tracelevel')  conv_maxiter(`conv_maxiter')  ///
conv_ptol(`conv_ptol') conv_vtol(`conv_vtol')  conv_nrtol(`conv_nrtol') igmmiterate(`igmmiterate') igmmeps(`igmmeps') igmmweps(`igmmweps')   ///
deriv(/xb = -1*exp({xb:} + {b0})) deriv(/b0 = -1*exp({xb:} + {b0})) 
}
}
qui predict double `_res' if `touse', equation(1)
qui replace `_res'=`y'-`_res' if `touse'
qui cor `y' `_res' if `touse'
di as text  "R-squared = " _continue
di as result r(rho)^2
end