/* 

Barbara Sianesi

University College London and Institute for Fiscal Studies
barbara_s@ifs.org.uk

last modified 28 May 01

*/



program define psmatch, rclass

	version 7.0

	syntax varname , ON(varlist numeric min=1 max=2) CALiper(real) [ID(varname numeric)] [OUTcome(varname numeric)] [SMooth(varname)] [EPan] [BOth] [noCOMMON] [noCOUnt]

	quietly  {

      if `"`id'"' != `""'  { 
	      sort `id'
	      count if `id'==`id'[_n-1]
		      if r(N)!=0  {
			      n di in red "more than one obs per `id'"
			      exit
		      }
		 }

   }

	set more off
   cap log off

   tempvar tempid   
   gen `tempid'=_n  


   tokenize `on'  

   if `"`2'"' != `""'  { 


      qui {

         preserve 
      
         drop if `1'==. | `2'==.

         tempname N1 N0 sigma1 sigma2 sigma12

         mat ac A0=`1' `2' if `varlist'==0
         mat XX0=A0[1..2, 1..2]
         mat ac A1=`1' `2' if `varlist'==1
         mat XX1=A1[1..2, 1..2]
         qui tabstat `1' `2', by(`varlist') save
         mat xx0=r(Stat1)'*r(Stat1)
         mat xx1=r(Stat2)'*r(Stat2)
         qui count if `varlist'==1
         scalar `N1'=r(N)
         qui count if `varlist'==0
         scalar `N0'=r(N)
         mat Sinv=(`N1'+`N0'-2)*syminv(XX1-`N1'*xx1+XX0-`N0'*xx0)
         scalar `sigma1'=Sinv[1,1]
         scalar `sigma2'=Sinv[2,2]
         scalar `sigma12'=Sinv[2,1]

         restore

      }

  }


	if `"`smooth'"'== `""' { 

      		quietly {
		
            cap drop _times
            cap drop _matchdif
            cap drop _matchedid
            
			   preserve
	
				tempvar  newid indic dif 
				tempname maxid		
				tempfile matchid
		      local identifier `"`id'"'

				keep if `1'!=.
            cap keep if `2'!=.
				keep `varlist' `tempid' `1' `2' `id'
				gsort -`varlist'
				gen `newid'=_n if `varlist'==1
				sum `newid' if `varlist'==1
				scalar `maxid'=r(max)
				gen byte _times=0
				gen _matchdif=0 if `varlist'==1

            if `"`id'"' != `""'  { 
				   gen _matchedid=0 if `varlist'==1
            }

				local i = 1
				while `i'<= `maxid'   {

					cap drop `indic'
					gen byte `indic' = 1 if `newid'==`i' 

					sort `indic' `varlist' 
					cap drop `dif'

               if `"`2'"' != `""'  {
                  gen `dif'=`sigma1'*(`1'[1]-`1')^2+`sigma2'*(`2'[1]-`2')^2+2*`sigma12'*(`1'[1]-`1')*(`2'[1]-`2') if `varlist'==0
               }

               else {
                  gen `dif'=abs(`1'-`1'[1]) if `varlist'==0
               }

					sort `dif' `id'
					replace _times=_times+1 if `dif'[1]<`caliper' & (`indic'==1 | _n==1)
					replace _matchdif=`dif'[1] if _times==1 & `indic'==1 
					
					if `"`id'"' != `""'  { 
					   replace _matchedid=`id'[1] if _times==1 & `indic'==1 
               }

               if `"`count'"' == `""'  { 
					   n di _skip(6) `maxid'-`i'
               }
					local i = `i'+1
				
				} 

				drop if _times==0
				replace _times=1 if `varlist'==1
				keep `tempid' _times _match*
				sort `tempid'

				compress
				save `matchid'

				restore
            sort `tempid'
				merge `tempid' using `matchid'
				drop _merge
				lab var _times "no. of times used"
				lab var _matchdif "difference in score"
				cap lab var _matchedid "`identifier' of matched control"

            if `"`count'"' == `""'  { 
				   n di _newline(3) _skip(3) in w "Now use _times~=. to identify the matched treated and matched controls"
				   n di _skip(3) in w "Use [fw=_times] in your analyses, or else type: expand _times"
			   }

     if `"`outcome'"' != `""'  {      
		
               tempname mean1 mean0 number weight var1 var0 stderr
               local out `outcome'

               sum `outcome' if `varlist' == 1 [fw=_times]  
               scalar `mean1'  = r(mean)
               scalar `number' = r(N)
               scalar `var1'   = r(Var)
               sum `outcome' if `varlist' == 0 [fw=_times]  
               scalar `mean0' =  r(mean)
               sum `outcome' if `varlist' == 0 & _times!=.
               scalar `var0'   = r(Var)

               tempvar we
               gen `we'=_times^2 if `varlist'==0
               sum `we'
               scalar `weight' = r(sum)
               scalar `stderr' = sqrt( (`var1')/`number' + (`weight'/(`number'^2))*`var0' )

               cap log on

               n di _newline(3) _skip(3) in g "Mean " in w "`outcome'" in g " of matched treated  = " in y `mean1'
               n di _newline(1) _skip(3) in g "Mean " in w "`outcome'" in g " of matched controls = " in y `mean0'
               n di _newline(1) _skip(3) in g "Effect  = " in y `mean1'-`mean0'
               n di _newline(1) _skip(3) in g "Std err = " in y `stderr' 
               n di _newline(1) _skip(3) in w "Note: takes account of possibly repeated use of control observations" 
               n di _skip(3) in w "      but NOT of estimation of propensity score."
               n di _newline(1) _skip(3) in g "T-statistics for " in w "H0: effect=0" in g " is " in y (`mean1'-`mean0')/`stderr'

               return scalar effect= `mean1'-`mean0'
            }


			}
	}

	else if `"`smooth'"' != `""' { 

			quietly  {

            local outcome `"`smooth'"'
			   cap drop _m`outcome'

            if `"`both'"' != `""'  {               
                cap drop _s`outcome'
            }                                      

				preserve
	
				tempvar  newid indic dif weight woutc denom numer denom1 numer1
				tempname min max min2 max2 maxid	
				tempfile matchid

				keep `tempid' `varlist' `1' `2' `smooth'

            keep if `1'!=.
            cap keep if `2'!=.

   			gsort -`varlist'
				gen `newid'=_n if `varlist'==1
				qui sum `newid' if `varlist'==1
				scalar `maxid'=r(max)

				gen _m`outcome'=.
            gen _s`outcome'=.

	         local i=1
				while `i'<=`maxid'   {

					cap drop `indic'
					gen byte `indic' = 1 if `newid'==`i' 
					sort `indic' `varlist'

					cap drop `dif'

               if `"`2'"' != `""'  {
                  gen `dif'=`sigma1'*(`1'[1]-`1')^2+`sigma2'*(`2'[1]-`2')^2+2*`sigma12'*(`1'[1]-`1')*(`2'[1]-`2') 
               }

               else {
                  gen `dif'=abs(`1'-`1'[1]) 
               }


					if `"`epan'"'!= `""'  {
						gen `weight' = 1-(`dif'/ `caliper')^2 if abs(`dif'/`caliper')<=1 
					}
					else  {
						gen `weight' = normden(`dif'/ `caliper') 
					}

					gen `woutc' = `smooth'*`weight' 

					sum `weight' if `varlist'==0, meanonly
					scalar `denom'=r(mean)
					sum `woutc' if `varlist'==0, meanonly
					scalar `numer'=r(mean)
					replace _m`outcome' = `numer'/`denom' if  `indic'==1 

               if `"`both'"' != `""'  {          
                  sum `weight' if `varlist'==1, meanonly
					   scalar `denom1'=r(mean)
					   sum `woutc' if `varlist'==1, meanonly
					   scalar `numer1'=r(mean)
					   replace _s`outcome' = `numer1'/`denom1' if  `indic'==1 
                  scalar drop `denom1' `numer1'
                }

					drop `indic' `weight' `woutc' 
					scalar drop `denom' `numer' 
					
					if `"`count'"' == `""'  {       
					   n di _skip(6) `maxid'-`i'
					}
					local i = `i'+1

					} 

            
			  if `"`common'"' == `""'  {

				      **** imposing COMMON SUPPORT -- on the treated only

				      sum `1' if `varlist'==0
				      scalar `min'=r(min)
				      sum `1' if `varlist'==0
				      scalar `max'=r(max)

				      drop if (`1'>`max' | `1'<`min') & `varlist'==1

                  if `"`2'"' != `""'  { 
   			         sum `2' if `varlist'==0
				         scalar `min'=r(min)
				         sum `2' if `varlist'==0
				         scalar `max'=r(max)

				         drop if (`2'>`max' | `2'<`min') & `varlist'==1

				      }      
   
                  count if `varlist'==1

                  if r(N)==0  {                                
	                  n di in red "none of the treated lies within the common support"  
	                  exit                                      
                  }                                            
                                             

            }



				keep if `varlist'==1
            if `"`both'"' != `""'  {          
                keep `tempid' _s`outcome' _m`outcome'
	         }                              
				else  {
				keep `tempid' _m`outcome'
				}

				sort `tempid'

				compress
				save `matchid'

				restore
            sort `tempid'
				merge `tempid' using `matchid'
				drop _merge
				lab var _m`outcome' "matched smoothed `outcome'"
				cap lab var _s`outcome' "treated smoothed `outcome'"

			}

   if `"`count'"' == `""'  { 
	   n di _newline(3) _skip(3) in g "_m`outcome'" in w " is defined only for treated (and within the common support);"
	   n di _skip(3) in w "It should be compared to the " in g `"`smooth'"' in w " (or" in g "_s`outcome'" in w ") for those treated only," 
	   n di _skip(3) in w "i.e. with _m`outcome'~=.
   }


   tempname mean1 mean0

   if `"`both'"' != `""'  {          
	   qui sum _s`outcome'
      scalar `mean1'  = r(mean)
   }                              
   else  {
	   qui sum `smooth'  if `varlist'==1 & _m`outcome'!=.
      scalar `mean1'  = r(mean)
   }

	qui sum _m`outcome' 
   scalar `mean0'  = r(mean)

   cap log on

   n di _newline(3) _skip(3) in g "Mean " in w "`outcome'" in g " of matched treated  = " in y `mean1'                        
   n di _newline(1) _skip(3) in g "Mean " in w "`outcome'" in g " of matched controls = " in y `mean0'                        
   n di _newline(1) _skip(3) in g "Effect  = " in y `mean1'-`mean0'                                                           

   return scalar effect= `mean1'-`mean0'

	}


set more on
cap log on

end