*! 1.0 MB & NO 29 Dec 2003

capture program drop xtvc
program define xtvc, eclass

	syntax  [, Level(int $S_level)  H0(numlist min=1 max=1)  ] 
 	
	version 6 
 
if replay() {
		if `"`e(cmd)'"'==`"xtreg"' {

                if "`e(model)'" != "ml"   {
                        noi di in red "xtreg, mle not found"
                        exit 301 
					}
				}
			else {
                        noi di in red "xtreg, mle not found"
                        exit 301 
				}						
		}
    
  preserve

    /* CHECK THE OPTION H0 FOR SIGMA_U */
   
			if "`h0'" != ""  {	 

					if `h0' < 0 {
							di as err "option h0(#): null value must be greater than or equal to zero"
							exit 198
							}
 
                              scalar PT = `h0'
                                     }
   
    /* CHECK SIGNIFICANCE LEVEL */

	 if `level' <10 | `level'>99 { 
						di in red "level() invalid"
						exit 198
						}   
                                
   	local levelci = `level'/100	
	scalar chi2 = invchi2(1,`levelci')
   
    /* GET RESULTS FROM XTREG, MLE */
                              
	local ivar "`e(ivar)'"
      local dep  "`e(depvar)'"
      local m = e(N_g)
      scalar sigma_e = e(sigma_e)
	scalar sigma_u = e(sigma_u)
	
 	tempvar xb pr res n  qres  sqres   Sqres 
	
 	qui predict `xb'
      qui gen double `res' = `dep' - `xb'

    /*  DROP MISSING VALUES */
    
    qui keep if `res' != .

    quietly  {
	sort  `ivar'
		by  `ivar': gen `n' = cond(_n == _N,_N,.)
		by  `ivar': gen `qres' = cond(_n==_N, sum(`res'/_N),.)		
		gen `sqres' = .
		gen `Sqres' = .
              }
            
/* SCORE TEST AT SIGMA_U == # */

 if  e(sigma_u)  == 0 {                    

                     scalar S0 = 0
                     scalar pval = 1
                     scalar PTEST = 0
                
            if "`h0'" != ""  {
                    scalar PTEST = PT
                    qui replace `sqres' =  `qres'^2 / ( PTEST ^2 + sigma_e^2/`n')
			  qui replace `Sqres' = sum(`sqres') if  `n' != .

			  scalar SUMSQRES = `Sqres' in l
			  scalar S0 = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) )  
                    scalar pval = chi2tail(1, S0)
                         	}
                     }
                else {
                	  scalar PTEST = 0
                
                      if "`h0'" != ""  {
                              	    scalar PTEST = PT
                                   	    }
 
					qui replace `sqres' =  `qres'^2 / ( PTEST ^2 + sigma_e^2/`n')
					qui replace `Sqres' = sum(`sqres') if  `n' != .

					scalar SUMSQRES = `Sqres' in l
					scalar S0 = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) )  
                    		scalar pval = chi2tail(1, S0)
 					}
 					
est scalar score = S0
est scalar pval = pval
est local pcmd "xtvc"

/* EQUATION SOLVING BY BISECTION METHOD - CONFIDENCE INTERVAL FOR SIGMA_U */
   
* SEEK THE UPPER BOUND FOR SIGMA_U
 
        scalar IMP = 0
        
        scalar NOT = 0
       	    	    
	  scalar p1  = sigma_u
		
	  qui replace `sqres' =  `qres'^2 / (p1^2 + sigma_e^2/`n')
  	  qui replace `Sqres' = sum(`sqres') if  `n' != .
  
  	  scalar SUMSQRES = `Sqres' in l
    	  scalar SUP = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) ) - chi2
		           
			if SUP > 0 {
				      scalar NOT = 1
			            scalar IMP = 1
			            scalar uppb = .
			            }
		           
      scalar p2  = sigma_u + 2  
					  			    	  		        	  		    
    	scalar STEP = p2
 
	while SUP < 0 & IMP != 1  & NOT != 1 {  
		    
		scalar p2 = p2 + STEP
		
					qui replace `sqres' =  `qres'^2 / (p2^2 + sigma_e^2/`n')
  		  			qui replace `Sqres' = sum(`sqres') if  `n' != .
  
  		  		    scalar SUMSQRES = `Sqres' in l
    	  		    	    scalar SUP = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) ) - chi2	
                     
    	 if sigma_u == 0 | sigma_u < 1 {
    	            		    
         		if p2 > (sigma_u + 100 )   { 
                   					       scalar IMP = 1
                   					       scalar uppb = .
                     					     }
                       }
           
          if sigma_u > 1 {
                     
                       if p2 > (sigma_u * 100 )   {   
                   					       scalar IMP = 1
                   					       scalar uppb = .
                     					     }
                     	}
                     					             					     
                     }   
 

   if IMP != 1  & NOT != 1  {
 
   while  abs(p1-p2) >  1e-9  {
					
                 local t = (p1+p2)/2
          
             quietly {
          			replace `sqres' =  `qres'^2 / (`t'^2 + sigma_e^2/`n')
  		  		replace `Sqres' = sum(`sqres') if  `n' != .
  
  		  		scalar SUMSQRES = `Sqres' in l
    	  		      scalar St = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) ) - chi2
                                          
          			replace `sqres' =  `qres'^2 / (p1^2 + sigma_e^2/`n')
  		  	      replace `Sqres' = sum(`sqres') if  `n' != .
  
  		  		scalar SUMSQRES = `Sqres' in l
    	  		      scalar Sp1 = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) ) - chi2
        		 	}      
    		 	
	    	   		   
    		 if (St*Sp1) > 0 {
    		 		 	scalar p1 = `t' 
 					}
 							
 		  if (St*Sp1) < 0 {
    	   		             scalar p2 =  `t' 
    		 		       }
    		 	   }	
   		
	           scalar uppb =  p1
        }
        	 
* SEEK THE LOWER BOUND FOR SIGMA_U
      
  if NOT != 1 {

        local c 0
		scalar p1  = sigma_u
		scalar p2  = 0
 
   while  abs(p1-p2) >  1e-9  {

   		 		 local c = `c' + 1

				local t  = (p1+p2)/2
                  
             quietly {
 
          			replace `sqres' =  `qres'^2 / (`t'^2 + sigma_e^2/`n')
  		  		replace `Sqres' = sum(`sqres') if  `n' != .
  
  		  		scalar SUMSQRES = `Sqres' in l
    	  		      scalar St = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) ) - chi2
                                          
          			replace `sqres' =  `qres'^2 / (p1^2 + sigma_e^2/`n')
  		  		replace `Sqres' = sum(`sqres') if  `n' != .
  
  		  		scalar SUMSQRES = `Sqres' in l
    	  		      scalar Sp1 = ( ( SUMSQRES - (`m'-1) )^2 / (2*(`m'-1)) ) - chi2
        		 	}         
    		 	
    		 if (St*Sp1) > 0 {
    		 		      scalar p1 =  `t'  
 					}
 							
 		  if (St*Sp1) < 0 {
    	    		             scalar p2 =  `t' 
    		 		       }
         }	
   		
   			if abs(p1) < 1e-5 {
   			              scalar p1 = 0
   			              }
   			              
	         scalar lowb =  p1
          est scalar suuppb = uppb
          est scalar sulowb = lowb   
      }
 else {
          scalar lowb =  0
          est scalar suuppb = uppb
          est scalar sulowb = lowb   
       }

 
/* SET FORMAT  */

local fmt = "%1.0f" 
                
  if "`h0'" != ""  {
                      local fmt = "%3.2f"
                     }
 
/* DISPLAY RESULTS */

 		noi di _n in smcl in gr "{hline 13}{c TT}{hline 47}"
		noi di in smcl in gr  _col(5) %8s abbrev("`dep'",8) _col(14) in gr "{c |}" _col(20) "ML Estimate" _col(36) in gr "[" "`level'" "% Conf. Interval]"     
 		noi di in smcl in gr "{hline 13}{c +}{hline 47}"
		noi di in smcl in gr   _col(5) in gr "/sigma_u"  _col(14) in gr "{c |}" _col(22)in y %9.0g  sigma_u  _col(35) %9.0g lowb _col(47) %9.0g uppb    
 		noi di in smcl in gr "{hline 13}{c BT}{hline 47}"
 		noi di in smcl in gr "Score test of sigma_u=" in y `fmt'  PTEST  in gr ": chi2(1)= " in y %3.2f S0  in gr " Prob>=chi2 = " in y %4.3f  pval 
end