*! 1.0.0  June 1998 Jeroen Weesie/ICS   STB-47 sg100
*! 2-stage estimation subject to linear constraints

program define linest
   version 5.0
      
   local options "Eform(str) Level(int $S_level) Constraints(str) Keep Modify"
   parse "`*'"   
   if "`constra'" == "" {
      if "$S_E_cmd" ~= "linest" {
         di in red "constraints required"
         exit 198
      }
      if "`keep'`modify'" ~= "" {
         di "keep and modify not permitted during replay"
         exit 198
      }   
      local replay 1
   }
   else {   
      if "$S_E_cmd" == "" {
         error 301
      }
                        
      if "`modify'" ~= "" {
         local keep 1
         * modify not supported for internal estimation commands
         intcmd `cmd'
         if $S_1 == 1 {
            di in re "modify is not allowed with internal estimation commands"
            error 198
         }            
      }
         
    * scratch

      tempname est bu Vu IVu T a C br Vr d X npu npr

    * process linear constraints       
     
      local cmd $S_E_cmd
      if "`keep'" ~= "" { msave }
      mat `bu' = get(_b)
      mat `Vu' = get(VCE)
      estimates hold `est'        
      mat post `bu' `Vu'
      mat makeCns `constra'
      mat `bu' = get(_b)
      mat `Vu' = get(VCE)
      matcproc `T' `a' `C'

    * GLS-estimator (br,Vr) of bu on T, weighted by Vu
    * Note: al restricted vectors br satisfy br = a + T b          
    * Note: If Vr is singular, we have be careful to use the 
    *       Moore-Penrose inverse rather than the g2-inverse!
    
      * X = Inv(T*Inv(Vu)*T')
      matginv `Vu', ginv(`IVu') rank(`npu')
      mat `X' = `IVu' * `T'
      mat `X' = `T'' * `X'
      matsymm `X'
      matginv `X', ginv(`X') rank(`npr')

      * br = a + T'*X*T*Inv(Vu)*(bu-a)      
      mat `br' = `X' * `T''
      mat `br' = `br' * `IVu'
      mat `d'  = `bu' - `a'
      mat `br' = `d' * `br''
      mat `br' = `br' * `T''
      mat `br' = `br' + `a'         /* restricted estimate */

      * Vr = T*(inv(T'*Inv(Vu)*T)*T' = T*X*T'
      mat `Vr' = `X' * `T''
      mat `Vr' = `T' * `Vr'         /* variance of restricted estimate */
      matsymm `Vr'
      
    * Wald's test for equality constraints

      * Wald-statistic = (bu-br)'*Inv(Vu)*(bu-br)
      mat `d' = `bu' - `br'         
      mat `X' = `d' * `IVu'
      mat `X' = `X' * `d''
                
    * post and display results
      if "`modify'" != "" {
		estimates unhold `est'
      		capture qui estimates post `br' `Vr' `C', noclear
		if _rc {
      			mat post `br' `Vr' `C'
		}
      }
      else	mat post `br' `Vr' `C'
      
    * set globals  
      global S_E_cmd  linest  
      global S_E_cmdu `cmd'
      global S_E_npu  `npu'
      global S_E_npr  `npr' 
      global S_E_rx2  = `X'[1,1]
  }

 * display results
 
   local nc = $S_E_npu - $S_E_npr
   local c "_n in gr _col(46)"  
   #del ;  
   di  in gr _n "Two-step constrained " in ye "$S_E_cmdu"
    `c' "Dim unrestricted model = " in ye %8.0f $S_E_npu
    `c' "Dim restricted model   = " in ye %8.0f $S_E_npr      
    `c' "# restrictions         = " in ye %8.0f `nc'
    `c' "X2 for restrictions    = " in ye %8.4f $S_E_rx2
    `c' "Prob > chi2(" in ye `nc' in gr ")" _col(69) /*
                        */ "= " in ye %8.4f chiprob(`nc',$S_E_rx2) _c ;
   #del cr   
  
   matrix dispCns       
   local level "level(`level')"             
   if "`eform'" ~= "" { 
      local eform "eform(`eform')"
   }
   mat mlout, `level' `eform'
   if "`replay'" ~= "" { exit }

 * keep/modify stored results
   
   if "`keep'" == "" {
      di _n in bl "Constrained estimation results are not stored."
      di    in bl "Post-estimation commands use the unconstrained model!" 
      estimates unhold `est'
   }
   else {
      if "`modify'" == "" {
         di 
         di in bl "linest saves results as-an-estimation command."
         di in bl "Beware (see online help!)"
         di in bl "  - chi2/ll/pr2 etc may been modified"
         di in bl "  - predict now returns linear predictors (index)" 
	 mrestore `nc' $S_E_rx2
         global S_E_cmd "linest"
      }
      else {
         * we replace the results of the original estimation command!
         di 
         di in bl "linest modifies the stored results of $S_E_cmdu`cmdu'"
         di in bl "Beware! (see online help!)"
	 capture local x_x = caller()
	 if !_rc { lin_m6r `nc' $S_E_rx2 "`cmd'" }
	 mrestore `nc' $S_E_rx2
         global S_E_cmd "`cmd'"
      }   
  } 
end

* symmetrizes a matrix: X = (X+X')/2
program define matsymm 
   local m `1'
   tempname mt
   mat `mt' = `m''
   mat `m' = `m' + `mt'
   mat `m' = `m' * 0.5 
end   

* store all S_E_X macros using SE_X  
program define msave
   global SE_tdf  $S_E_tdf
   global SE_ttl  $S_E_ttl   
   global SE_depv $S_E_depv
   global SE_nobs $S_E_nobs
   global SE_chi2 $S_E_chi2
   global SE_mdf  $S_E_mdf         
   global SE_pr2  $S_E_pr2         
   global SE_ll0  $S_E_ll0            
   global SE_ll   $S_E_ll            
   global SE_if   $S_E_if
   global SE_in   $S_E_in   
   global SE_vce  $S_E_vce
   global SE_vcn  $S_E_vcn   
   global SE_cn   $S_E_cn   
end

* restore all S_E_X macros from SE_X, updating relevant statistics  
program define mrestore
   local nc `1'         /* number of independent constraints */
   local x2 `2'         /* Wald-statistic for constraints */ 
      
   global S_E_ttl  "$SE_ttl (constrained)"  
   global S_E_depv $SE_depv
   global S_E_nobs $SE_nobs
   
   * approximate the new ll by ll+x2, and adapt degrees of freedom (!)
   
   if "$SE_chi2" ~= "" { global S_E_chi2 = $SE_chi2 - `x2' }
   if "$SE_mdf"  ~= "" { global S_E_mdf  = $SE_mdf  - `nc' }
   if "$SE_tdf"  ~= "" { global S_E_tdf  = $SE_tdf  + `nc' }            
   if "$SE_ll"   ~= "" { global S_E_ll   = $SE_ll   - 0.5*`x2' }            
   global S_E_ll0 $SE_ll0
   if "$SE_pr2" ~= ""  { 
      capt global S_E_pr2  = 1 - $S_E_ll / $S_E_l0
   }           
   
   global S_E_if  $SE_if
   global S_E_in  $SE_in   
   global S_E_vce $SE_vce
   global S_E_vcn $SE_vcn   
   global S_E_cn  $SE_cn      
   
   global drop SE_tdf SE_ttl SE_depv SE_nobs SE_chi2 SE_mdf SE_pr2 SE_ll0 /*
       */ S_E_ll SE_if SE_in SE_vce SE_vcn SE_cn
end

* returns in S_1 1 (0) if estimation cmd depends on an internal estimation
* command. 
program define intcmd 
   local c `1'
   if index(" `c' "," anova areg clogit cox logit logistic mlogit ")>0 | /*
   */ index(" `c' "," ologit oprobit poisson probit regress rreg ")>0  | /*
   */ index(" `c' "," stcox tobit ")>0 /*
   */ {
        global S_1 1
   }
   else global S_1 0
end