*! v.1.0.0 M.Bottai, N.Orsini 10aug2012

capture program drop laplace
program laplace, eclass properties(mi) byable(onecall)
version 10 

if _by() {
		local BY `"by `_byvars'`_byrc0':"'
}
	
`version' `BY' _vce_parserun laplace, mark(OFFset CLuster) : `0'

if "`s(exit)'" != "" {
		version 10: ereturn local cmdline `"laplace `0'"'
		exit
	}

if replay() {
		if ("`e(cmd)'"!="laplace")  error 301  
		Replay `0'
	}
else `version' `BY' Estimate `0'
end

capture program drop Estimate 
version 10, missing

program Estimate, eclass byable(recall)
	syntax varlist  [if] [in]  ///
	[, Level(integer $S_level) ///
	   Quantiles(numlist) ///
	   Failure(string) ///
	   Reps(string) ///
	   TOLerance(real 1e-10) ///
	   MAXiter(integer 2000) /// 
	   SEED(string) ///
	   sigma(varlist)   ]
	   
    local cmdline : copy local 0
	marksample touse 

// Check the failure variable
		
		local nofailure = 0
		if "`failure'" == "" {
					tempname failure
					qui gen `failure' = 1 if `touse'
					local nofailure = 1
		}
		else {
				confirm variable `failure'
		 
				confirm numeric variable `failure'
			 
				qui replace `failure' = (`failure' != 0) if `touse'
 
		}

// set the seed

			if "`seed'" != "" {
				set seed `seed'
			}
			local seed `c(seed)'
 
// Check specified quantiles 

		if "`quantiles'" == "" local quantiles = .5
 		
 		foreach i of local quantiles   {
 			local tempq "`i'"
  			if `tempq' >= 1 & `tempq' < 100 { 
						if `tempq' >= 1 & `tempq' < 10 { 
 						local tempq "0`tempq'"
 						}
 		    local tempq : subinstr local tempq "." "", all   
 		    local tempq ".`tempq'" 
 			}
 			
 			local newquantiles "`newquantiles' `tempq'"
			if (`tempq'<=0  | `tempq'>=1) {
						di as err "specified quantiles out of the range (0,1)"
						exit 198
			}

 		}
 		
 		 local quantiles "`newquantiles'"
 		  local quantiles : list uniq  quantiles
 		 local quantiles:   list sort quantiles
   		 
		tempname quantlist
		local s = 1
 		local ml = 2
		foreach i of local quantiles   {
 		 
 			if (length("`i'")-1) > `ml' local ml = length("`i'")-1
 			if `s++' == 1 mat `quantlist' = `i'
			else mat `quantlist' = (`quantlist' , `i')
 		}
 				
		local newquantiles ""

 		foreach i of local quantiles   {
				local newquantiles "`newquantiles' `: di %`=`ml'+1'.`ml'f `i''" 
  		}
  		
        local quantiles `newquantiles' 
  				  					 			
// Check number of bootstrap replications 
// Set asymptotic variance estimator as default
 
 		if "`reps'" == "" {
							local nreps = 0
		}
		else {
				if `reps'<2  {
						di  as err "specify more than 1 replication"
						exit  
				}
				local nreps = `reps'
		}
			
// Remove collinear variables

		gettoken depv indepv : varlist		
		_rmcoll `indepv' if `touse' 
		local indepv `r(varlist)'  
		
		local nx : word count `indepv'
 		local eqnames : subinstr local quantiles "0." "q", all   
 
 // Define labels for the regression coefficients for each quantile
 
	local nq : word count `quantiles'
	forv i = 1/`nq' {
		local pq : word `i' of `eqnames' 
		local conams "`conams' `indepv' _cons"
 		forv j = 1/`=`nx'+1' {
			local eqnams "`eqnams' `pq'"
		}
 	}
    
	
// Define the maximization algorithm 
	qui count if  `touse'		
	local nobs = r(N)
	qui count if `failure' == 1 & `touse' == 1
	local notcensored = r(N)
	
	if (`notcensored' == `nobs') {
		local lfname "laplace_gs0"
	}
	else {
		local lfname "laplace_gs"
	}

// Check variable list for sigma. It applies only if failure is specified.
		
	if "`sigma'" == "" {
			local varsigma ""
	}
	else {
			_rmcoll `sigma'  if `touse'
			local varsigma  `r(varlist)' 
	}
	if (`nofailure' == 1 & "`sigma'" != "") {
		* the option sigma() is ignored if the option failure() is not specified  
	    local varsigma "" 
	}
	
// Estimation and bootstrap confidence intervals

		mata: laplace("`depv'", "`indepv'", "`failure'",  "`touse'", "`quantlist'", `nreps', `maxiter', `tolerance',  &`lfname'(), "`varsigma'")
		tempname coefs VCE  
		mat `coefs' = _beta
		mat `VCE' = VCE_beta 

		mat colnames `coefs' = `conams'
		mat coleq `coefs' = `eqnams'
		mat rownames `VCE' = `conams'
		mat roweq `VCE' = `eqnams'
		mat colnames `VCE' = `conams'
		mat coleq `VCE' = `eqnams'

		ereturn post `coefs' `VCE', obs(`nobs')  depn(`depv') 	
		ereturn local qlist "`quantiles'"
		ereturn repost, esample(`touse')
		
		ereturn local depvar "`depv'"
		if (`nofailure' != 1) ereturn local failvar "`failure'"
		if "`reps'" != "" ereturn scalar reps = `reps'
		ereturn scalar N = `nobs'

		if (`nofailure' == 0) ereturn scalar  N_fail = `notcensored'
		ereturn scalar n_q = `nq'
		ereturn local eqnames "`eqnames'"

		if "`reps'" != "" ereturn local vcetype "Bootstrap"
		else ereturn local vcetype "Robust"
		ereturn local predict "sqreg_p"
		ereturn local cmdline `"laplace `cmdline'"'
		ereturn local cmd "laplace"

		ereturn repost
 			
		di _n in gr "Laplace regression"  _col(50) "No. of subjects  = " in ye %10.0g e(N)
		
		if (`nofailure' != 1) di as text "" _col(50) in gr "No. of failures  = " in y %10.0f  e(N_fail) 
			 
		Replay, level(`level')  
end

capture program drop Replay
program Replay
	syntax [, Level(cilevel) ]
	ereturn display, level(`level')  
end


mata
mata clear
mata set matastrict off

real scalar laplace_lf(real rowvector p, real colvector y, real matrix X, real matrix XS, real colvector d, real colvector wgt, real scalar q)
{
	real scalar s , lf
	real colvector z
    s = exp(XS*p[1,(cols(X)+1)::(cols(X)+cols(XS))]')
    z = (y:-X*p[1,(1::cols(X))]'):/s
    lf = mean(wgt:*((z:<=0):*(d:*((1:-q):*z:+log(q*(1-q):/s)) :+ (1:-d):*log(abs(1:-q:*exp((1:-q):*z))) ) + (z:>0):*( d :*(-q:*z:+log(q*(1-q):/s)) :+ (1:-d):*(-q:*z:+log(1-q)) )))
	return(lf)
}
 
real scalar laplace_lf0(real rowvector p, real colvector y, real matrix X, real colvector wgt, real scalar q)
{
	real scalar lf
	real colvector z
    z = (y:-X*p[1::cols(X)]')
    lf = mean(wgt:*z:*(q:-(z:<=0)))
	return(lf)
}
 	
real matrix laplace_gs(real rowvector bi, real scalar maxiter, 
					   real scalar tol, real scalar step, real colvector y, real matrix X, real matrix XS,
					   real colvector d, real colvector wgt, real scalar q)
{  
	real rowvector b, nb, lb, ls  
	real scalar n, l, nl,  i, s 
	real colvector z, h
	real matrix tx
	b = bi	
	n = rows(y)
	tx = (wgt:*X)'	
 
	l = laplace_lf(b, y, X, XS, d, wgt, q)	
	s = exp(XS*b[1,(cols(X)+1)::(cols(X)+cols(XS))]')
	z = (y:-X*b[1,1::cols(X)]'):/s
	h = (q:-(z:<=0):-(1:-d):*(z:<=0):*(q:-1):/(1:-q:*exp((1:-q):*z)))
		
	lb = (tx*(h:/s)/n)'
	ls = (wgt:*XS)'*(z:*h:-d):/n:/s
	for (i=1; i<=maxiter; i++) {
 	 	    nb = b :+ (lb, ls'):*step
			nl = laplace_lf(nb, y, X, XS, d, wgt, q)
		      
			if(nl==.) {
                  nl=l
              }
 			if (nl>l) {			
					if  ( (abs(nl-l)<tol)  )  return(nb)		
					b = nb
					l = nl
					s = exp(XS*b[1,(cols(X)+1)::(cols(X)+cols(XS))]')
					z = (y:-X*b[1,1::cols(X)]'):/s
					h = (q:-(z:<=0):-(1:-d):*(z:<=0):*(q:-1):/(1:-q:*exp((1:-q):*z)))
					lb = (tx*(h:/s)/n)'
					ls = (wgt:*XS)'*(z:*h:-d):/n:/s
					step = step*1.25 
			}
	    	else {		
				 step = step*( (1.2 /(1+.4*exp(.005*min((n,1000))))) + .1 )
			   }		 
 		}
 	
	printf("{err:x}")
	displayflush()
 	return(b)


}
 

real matrix laplace_gs0(real rowvector bi, real scalar maxiter, 
					    real scalar tol, real scalar step, real colvector y, real matrix X, real matrix XS,
					    real colvector d, real colvector wgt, real scalar q)
{  
	real rowvector b, nb, lb
	real scalar n, l, nl, i
	real matrix tx

	b = bi[1,1::cols(X)]
	n = rows(y)
	tx = (wgt:*X)':/n
	l = laplace_lf0(b, y, X, wgt, q)		

 	lb = (tx*(q:-(y:<=X*b')))'   
		
		for (i=1; i<=maxiter; i++) {
	 	    nb = b :+ lb:*step
			nl = laplace_lf0(nb, y, X, wgt, q)	

			if (nl<l+tol) { 
					if  (  (abs(nl-l)<tol)  )  return((nb , J(1,cols(XS),0)))
					b = nb
					l = nl			
					lb = (tx*(q:-(y:<=X*b')) )'  
					step = step*1.25
			}
	    	else {	
				 step = step*.25
			   }		 
 		}
	
	printf("{err:x}")
	displayflush()
	return(b , J(1,cols(XS),0))
 }

real matrix laplace(string scalar depv, string matrix indepv, string scalar censored, 
                   string scalar touse, string scalar qlist, 
                   real scalar boot, real scalar maxi, real scalar tol, pointer(real scalar function) scalar f, 
                   string matrix vsigma)
{
	real rowvector th, binit, beta, beta2, covlist, quantiles 
	real colvector  r, s, wgt, y, d, bwgt
	real scalar 	step, i,   j, ya, yb
	real matrix 	tb, X,   XS
	
	covlist = st_varindex(tokens(indepv))
	st_view(y=., ., depv, touse)
    st_view(X=., ., covlist, touse)
    st_view(d=., ., censored, touse)
	slist = st_varindex(tokens(vsigma))
	st_view(XS=., ., slist, touse)
	yinit = y
	wgt = J(rows(y) , 1, 1)
	
	// Rescale the dependent variable to avoid problems of overflow
	 ya = min(y)
     yb = (max(y)-min(y))/10
     y = (y:-ya):/yb
	
	// Transform the covariates (standardize)
	if (cols(X)>0) {
	meanX = mean(X)
	sdX = sqrt(diagonal(variance(X)))'
	Z = (X :- meanX):/sdX
	Z = (Z,J(rows(y),1,1))
	}
	else {
	Z = J(rows(y),1,1)
	}
 
	// Transform the covariates (standardize) for Sigma
	if (cols(XS)>0) {
	meanXS = mean(XS)
	sdXS = sqrt(diagonal(variance(XS)))'
	ZS = (XS :- meanXS):/sdXS
	ZS = (ZS,J(rows(y),1,1))
	}
	else {
	ZS = J(rows(y),1,1)
	}
 
    binit = ((invsym(Z'Z)*Z'y)', J(1, cols(ZS), 0)) 
 	ostep = sqrt(variance(y, wgt))    
 	quantiles = st_matrix(qlist)
  
	th = J(1, (cols(quantiles)*(cols(Z)+cols(ZS))), .)	
	 
	// Get point estimates 
	
	for (i=1; i<=cols(quantiles); i++) { 
		step = ostep
		th[1,((i-1)*(cols(Z)+cols(ZS))+1)::(i*(cols(Z)+cols(ZS)))] = (*f)(binit, maxi,  tol, step, y, Z, ZS, d, wgt, quantiles[1,i])
 	    beta = th[1,((i-1)*(cols(Z)+cols(ZS))+1)::(i*(cols(Z)+cols(ZS)))]
		beta = beta[1,1::cols(Z)]
 
     // back transform the betas because of transformed response
     
		 beta = beta[1,1::cols(Z)]:*yb :+ (J(1,cols(Z)-1,0),1):*ya
			  
	 // back transform the betas because of transformed covariates 

		 if (cols(X)>0) {
						   beta[1,cols(Z)] = beta[1,cols(Z)] - (beta[1, 1::cols(Z)-1]*(meanX:/sdX)')
						   beta[1, 1::cols(Z)-1] = beta[1, 1::cols(Z)-1]:/sdX
		  }
		  
		  	if (i==1)  beta2 = beta
			else beta2 = beta2, beta			
			 			
		// Get the asymptotic variance/covariance matrix and create a block-diagonal matrix 

		if (boot==0) {	
				
						beta = th[1,((i-1)*(cols(Z)+cols(ZS))+1)::(i*(cols(Z)+cols(ZS)))]
						BZ = beta[1,1::cols(Z)]
						BZS = beta[(cols(Z)+1)::(cols(Z)+cols(ZS))]
                        
						s = exp(ZS*BZS') 
	                    z = (y:-Z*BZ') :/ s
					    q = quantiles[1,i]
						 
				     	h = sqrt(variance(z))/log(rows(z)) 
						S = Z:*((q:-(z:<=0):-(1:-d):*(z:<=0):*(q:-1):/(1:-q:*exp((1:-q):*z)))):/s
						
						z1 = z :+ h 
						z2 = z :- h
						
						K1 = (q:-(z1:<=0):-(1:-d):*(z1:<=0):*(q:-1):/(1:-q:*exp((1:-q):*z1))):/s
						K2 = (q:-(z2:<=0):-(1:-d):*(z2:<=0):*(q:-1):/(1:-q:*exp((1:-q):*z2))):/s
					
						H = invsym(Z'*(Z:*(K1:-K2):/2:/h:/s))
						
						VBV  =  H*S'*S*H            
						VBVN = VBV:*yb^2

						if(cols(Z)>1) {
							VBVN = VBVN:/( (sdX, 1)'* (sdX, 1) )
							for (r=1; r<=(cols(Z)-1); r++) {                                                                    
											VBVN[cols(Z), cols(Z)] = VBVN[cols(Z), cols(Z)] + VBVN[r,r]*(meanX[1,r]^2) - 2*VBVN[cols(Z),r]*meanX[1,r]
											VBVN[cols(Z), r] = VBVN[cols(Z), r] - VBVN[r, r]*meanX[1,r]
											 VBVN[r, cols(Z)] = VBVN[cols(Z), r]  
							}             
						}

						for (k=1; k<=cols(VBVN); k++) { 
								for (h=1; h<=cols(VBVN); h++) { 
									if (VBVN[k,h] == .) VBVN[k,h]=0
								}					
						}
						
						_makesymmetric(VBVN)
						
						if (i==1) AV = VBVN
						else AV = blockdiag(AV, VBVN)
					 
		}
	
	}
	 
	// Get bootstrap confidence intervals
		tb = J(cols(quantiles)*(cols(Z)),boot, .)
 
		for (i=1; i<=boot; i++) {

			s = ceil(rows(y) * uniform(rows(y), 1))
 			bwgt =  wgt 
			
			// Rescale the dependent variable to avoid problems of overflow
			yab = min(yinit[s])
			ybb = (max(yinit[s])-yab)/10
			yboot = (yinit[s]:-yab):/ybb
 
			// Transform the covariates (standardize)
			if (cols(X)>0) {
			meanXb = mean(X[s,], bwgt)
			sdXb = sqrt(diagonal(variance(X[s,], bwgt)))'
			Z = (X[s,] :- meanXb):/sdXb
			Z = (Z,J(rows(s),1,1))
			}
			else {
			Z = J(rows(s),1,1)
			}
			 
			// Transform the covariates (standardize) for Sigma
			if (cols(XS)>0) {
			meanXSb = mean(XS[s,], bwgt)
			sdXSb = sqrt(diagonal(variance(XS[s,], bwgt)))'	 
			ZS = (XS[s,] :- meanXSb):/sdXSb
			ZS = (ZS,J(rows(s),1,1))
			}
			else {
			ZS = J(rows(s),1,1)
			}
			binitb = ((invsym(Z'Z)*Z'yboot)', J(1, cols(ZS), 0)) 
  			step = sqrt(variance(yboot, bwgt))  					 

			for (j=1; j<=cols(quantiles); j++) {
					 TB = (*f)(binitb, maxi,  tol, step, yboot, Z, ZS, d[s], bwgt, quantiles[1,j])	 
					 TB = TB[1,1::cols(Z)]:*ybb :+ (J(1,cols(Z)-1,0),1):*yab
 
						if (cols(X)>0) {	 
								   TB[1, cols(Z)]  = TB[1,cols(Z)] - (TB[1, 1::(cols(Z)-1)]*(meanXb:/sdXb)')
								   TB[1, 1::(cols(Z)-1)] = TB[1, 1::(cols(Z)-1)]:/sdXb
						}
						 
					 tb[((j-1)*(cols(Z))+1)::(j*(cols(Z))),i] = TB'
				 }
 	 
	              
 
					
			
		}
 
  
    st_matrix("_beta", beta2)
    if (boot!=0) st_matrix("VCE_beta", variance(tb'))
	else st_matrix("VCE_beta",  AV)
   
}
 
end

exit

// Example
clear  
set obs 1000
local q = .5
gen x1 = invnorm(uniform())
gen x2 = uniform() < .5
gen t = 30 + 5*x1 + 2*x2 + invnorm(uniform()) - invnorm(`q')
gen c = 30 + 5*x1 + 2*x2 + invnorm(uniform()) - invnorm(`q')
gen y = min(t, c)
gen d = t <= c
laplace y x1 x2, f(d) q(`q')