*! qregpd2 
*! version 1.0 2026-09-28
*! This command is a modified version of qregpd by David Powell, Travis Smith, and Matthew Baker.
*! It has been modified to support multi-chain MCMC draw retention.
*! The command is intended to be called internally by qrpdmcmc.
*! authors Kamilla Kosheeva & Michail Liatos

program qregpd2
	version 11.2
	if replay() {
		if (`"`e(cmd)'"' != "qregpd2") error 301
		Replay `0'
	}
	else Estimate `0'
end
program Estimate, eclass

	syntax varlist [fw aw/] [if] [in] , 	///
		[OPTIMIZE(string)					///
			DRAWs(int 1000)					///
			BURN(int 0)						///
			ARATE(real .234)                ///
			SAMPLER(string)					///
			DAMPPARM(real 1)				///
			FROM(string)					///
			FROMVariance(string)			///
			SAVING(string)					///
			REPLACE							///
			APPEND							///
			THIN(int 1)						///
			JUMBLE							///
			NOISY							///
			grid1(numlist)					///
			grid2(numlist)					///
			NMSIMplex(string)				///
			ANALYTIC						///
			USEMAX							///
		IDentifier(varname)   ///
		INSTRuments(varlist)  ///
		FIX(varname)          ///
		Quantile(real 0.5)]

	marksample touse
	gettoken lhs rhs : varlist
	markout `touse' `lhs' `rhs'
	
	
	if `quantile' >= 1 {
			local tau = `quantile'/100
		}
		else	local tau "`quantile'"
		if `tau' <= 0 | `tau' >= 1 {
			di in red "quantiles(`quantile') out of range"
			exit 198
		}
	
	tempvar ones
	qui gen `ones'=1 if `touse'
	
	tempvar weightvar
	if `"`exp'"' != "" {
		qui gen double `weightvar' = `exp' if `touse'
	}
	else	{
		qui gen `weightvar' = `ones' if `touse'
	}
	
if ("`identifier'"=="") {
    di as error "Error: You must specify a panel identifier"
    exit
}

else if ("`fix'"=="") {
    di as error "Error: You must specify a fixed effect"
    exit
}

else {

    tempvar panel_nobs

    quietly bysort `identifier': ///
        egen long `panel_nobs' = total(`touse')

    quietly count if `touse' & `panel_nobs' == 1
    local N_singletons = r(N)
	
    if `N_singletons' > 0 {

        quietly replace `touse' = 0 ///
            if `touse' & `panel_nobs' == 1
    }

    quietly count if `touse'

    if r(N) == 0 {
        di as error ///
            "No observations remain after dropping singleton panel units."
        exit 2000
    }

    preserve

    sort `identifier' `fix'

    mata: id=st_data(.,"`identifier'","`touse'")
    mata: info=panelsetup(id,1)
    mata: fix=st_data(.,"`fix'","`touse'")
}

	if ("`optimize'" !="mcmc" & "`optimize'" !="grid") {
		local optimize "NM"
	}
	if ("`optimize'" !="mcmc") {
		if (`burn'!=0 | `draws'!=1000 | `arate'!=.234 | `dampparm'!=1 ///
		  | `thin'!=1 | "`sampler'"!="" ///
		  | "`saving'"!="" | "`replace'"=="replace" | "`append'"=="append" ///
		  | "`jumble'"=="jumble") {
			disp as err "An mcmc option has been set but estimation method is not mcmc. Option is superfluous." 
			disp as err "Please consult the help file for information on how to use mcmc options."
			disp as txt " "
		}
	}

	if ("`optimize'" !="grid") {
		if ("`grid1'"!="" | "`grid2'"!="") {
			disp as err "A grid search option has been set but estimation method is not grid search. Option is superfluous." 
			disp as err "Please consult the help file for information on how to use grid search options."
			disp as txt " "
		}
	}	
	
	if ("`instruments'" == "") {
		local instruments  "`rhs'"
	}
	local instruments : list retokenize instruments
	if wordcount("`instruments'") < wordcount("`rhs'") {
		di as error "Error: # of instruments must be >= # of right-hand side variables."
		exit
	}
	
	tempname excluded_inst exog_vars endog_vars
	local excluded_inst : list instruments - rhs
	local exog_vars : list instruments - excluded_inst
	local endog_vars : list rhs - instruments
	
	if ("`from'" == "" &  "`fromvariance'"=="") {
		qui qreg `lhs' `rhs' if `touse' [aw= `weightvar'], q(`tau')
		mata: beta_start=st_matrix("e(b)")
		mata: beta_start=beta_start[1..cols(beta_start)-1]
		mata: beta_start_var=st_matrix("e(V)")
		mata: beta_start_var=beta_start_var[1..(cols(beta_start)),1..(cols(beta_start))]
		mata: beta_start_se=sqrt(diagonal(beta_start_var))'
	}
	else if ("`from'" == "" & "`fromvariance'"!="" ) {
		mata: beta_start_var=st_matrix("`fromvariance'")
		mata: beta_start_se=sqrt(diagonal(beta_start_var))'
		qui qreg `lhs' `rhs' if `touse' [aw= `weightvar'], q(`tau')
		mata: beta_start=st_matrix("e(b)")
		mata: beta_start=beta_start[1..cols(beta_start)-1]
	}
	else if ("`from'" != "" & "`fromvariance'"=="" ) {
		mata: beta_start=st_matrix("`from'")
		qui qreg `lhs' `rhs' if `touse' [aw= `weightvar'], q(`tau')
		mata: beta_start_var=st_matrix("e(V)")
		mata: beta_start_var=beta_start_var[1..(cols(beta_start)),1..(cols(beta_start))]
		mata: beta_start_se=sqrt(diagonal(beta_start_var))'
	}
	else {
		mata: beta_start=st_matrix("`from'")
		mata: beta_start_var=st_matrix("`fromvariance'")
		mata: beta_start_se=sqrt(diagonal(beta_start_var))'
	}

	mata: tau     =strtoreal(st_local("tau"))
	mata: st_view(di=.,.,"`rhs'","`touse'")
	mata: st_view(zi=.,.,"`instruments'","`touse'")
	mata: st_view(yi=.,.,"`lhs'","`touse'")
	mata: st_view(wt=.,.,"`weightvar'","`touse'")

	mata: wt = wt:/colsum(wt)*rows(info)
	
	mata: mm_panels(id,ti=.)				/* Nx1 vec num of t per i,  */ 						
	mata: Ti=mm_expand(ti,ti,1,1) 			/* NTx1 vec of num of t per i */

	mata: Znew=demeanX(id,zi,wt)	
	
	mata: A=weightMatrix(beta_start,tau,yi,di,Znew,id,fix,zi,wt,"efficient")
	
	mata: M=moptimize_init()
	mata: moptimize_init_evaluator(M,&objectiveFunGmm())
	mata: moptimize_init_touse(M,"`touse'")	
	mata: moptimize_init_evaluatortype(M,"d0")
	mata: moptimize_init_technique(M,"nm")
	if "`nmsimplex'"=="" {
		mata: moptimize_init_nmsimplexdeltas(M,J(1,1,10))
	}
	else {
		mata: simplex=st_matrix("`nmsimplex'")
		mata: moptimize_init_nmsimplexdeltas(M,simplex)
	}
	mata: moptimize_init_depvar(M,1,yi)
	mata: moptimize_init_eq_indepvars(M,1,di)
	mata: moptimize_init_eq_cons(M,1, "off")
	mata: moptimize_init_eq_coefs(M,1, beta_start)

	mata: G=qrpdInfoInit(yi,Znew,di,A,info,fix,zi,wt,tau,id,Ti)
	mata: moptimize_init_userinfo(M,1,G)
	
/* MCMC optimization */
	if "`optimize'"=="mcmc" {
		di "Adaptive MCMC optimization"
		
		mata: junk=moptimize_evaluate(M)
		
		if (`burn'>=`draws') {
			di as error "Error: The # of burn in draws must be less than # draws"
			exit
		}	
				   
		if "`sampler'"!="mwg" {
			local sampler "global"
		}
		
		mata: noisy   =st_local("noisy")
		mata: sampler =st_local("sampler")	
		mata: arate   =strtoreal(st_local("arate"))
		mata: alginfo ="standalone",st_local("sampler")
		mata: draws   =strtoreal(st_local("draws"))
		mata: dampparm=strtoreal(st_local("dampparm"))
		mata: burn    =strtoreal(st_local("burn"))
		mata: thin    =strtoreal(st_local("thin"))
		mata: jumble  =st_local("jumble")
		
		mata: alginfo=sampler,"gf0","moptimize"
		mata: Draws = amcmc(alginfo,&objectiveFunGmm(),beta_start, beta_start_var,draws,burn,dampparm,arate,ar=.,vals=.,lam=.,.,M,noisy)
				
		tempname arates arates_mean arates_min arates_max
		mata: st_matrix("`arates'",ar)
		mata: st_numscalar("`arates_mean'",mean(ar'))
		mata: st_numscalar("`arates_min'",min(ar))
		mata: st_numscalar("`arates_max'",max(ar))
		
		mata: gqreg_process_draws(Draws,vals,jumble,thin,0,betaf=.,valsf=.)
		
		tempname dof
		mata: st_numscalar("`dof'",rows(betaf))
		local dof=`dof'
		
		tempname b V gamma valMean valMax valMin drawsKept
		
		if "`usemax'"!="" {
			mata: bmcmc=findRowMax(valsf,betaf)[1,.]
		}
		else {
			mata: bmcmc=mean(betaf)
		}
		
		if "`analytic'"!="" {
			mata: Vmcmc=qrpd_var(yi,di,bmcmc,tau,fix,info,zi,wt,beta_start_se,id,Znew)
		}
		else {
			mata: Vmcmc=variance(betaf)
		}
		
		mata: st_matrix("`b'",bmcmc)
		mata: st_matrix("`V'",Vmcmc)
		mata: st_numscalar("`valMean'",mean(valsf))
		mata: st_numscalar("`valMax'",max(valsf))
		mata: st_numscalar("`valMin'",min(valsf))
		mata: st_numscalar("`drawsKept'",rows(valsf))
	
		mata: gam=estimateGamma(bmcmc,yi,di,wt,tau,fix)
		mata: uniq_gam=uniqrows((fix,gam))						/* in "no constant" form */
		mata: constant=uniq_gam[rows(uniq_gam),2]				/* let constant be last value */

		mata: st_matrix("`gamma'", uniq_gam)
				
		tempname beta_start beta_start_var
		mata: st_matrix("`beta_start'",beta_start)
		mata: st_matrix("`beta_start_var'",beta_start_var)
	
		mat colnames `b' = `rhs' 
		mat colnames `V' = `rhs' 
		mat rownames `V' = `rhs' 
		mat colnames `gamma' = `fix' estimate
		
		tempname N_g N g_min g_max
		mata: st_numscalar("`N_g'",panelstats(info)[1])
		mata: st_numscalar("`N'",panelstats(info)[2])
		mata: st_numscalar("`g_min'",panelstats(info)[3])
		mata: st_numscalar("`g_max'",panelstats(info)[4])
		local N_g=`N_g'
		local N=`N'
		local g_min=`g_min'
		local g_max=`g_max'
	
		ereturn clear
		ereturn post `b' `V', esample(`touse') obs(`N') 
		
		ereturn local title "Quantile Regression for Panel Data (QRPD)"
		ereturn local cmd "qregpd"
		ereturn local indepvars "`rhs'"
		ereturn local instruments "`instruments'"
		ereturn local depvar "`lhs'"

		ereturn matrix gamma `gamma'
		ereturn matrix b_init=`beta_start'
		ereturn matrix V_init=`beta_start_var'
		
		ereturn local sampler `sampler'
		ereturn local saving `saving'
		ereturn local jumble `jumble'
		
		ereturn scalar N=`N'
		ereturn scalar N_singletons = `N_singletons'
		ereturn scalar N_g=`N_g'
		ereturn scalar g_min=`g_min'
		ereturn scalar g_max=`g_max'
		
		ereturn scalar draws=`draws'
		ereturn scalar burn=`burn'
		ereturn scalar thin=`thin'
		ereturn scalar damper=`dampparm'
	
		ereturn scalar opt_arate=`arate'
		ereturn scalar arates_mean=`arates_mean'
		ereturn scalar f_mean=`valMean'
		ereturn scalar f_max =`valMax'
		ereturn scalar f_min =`valMin'
	
		ereturn scalar draws_retained = `drawsKept'
		ereturn matrix arates=`arates'
	
if "`saving'"!="" {
    clear

    local exnames "`rhs'"
    getmata (`exnames')=betaf
    getmata (fun_val)=valsf

    gen t=_n
    if "`append'"=="append" {
        append using "`saving'"
    }
    else {
        qui save "`saving'", `replace'
    }

    restore

    ereturn repost, esample(`touse')
}		
		
		di _newline
		di as txt "`e(title)'" 
		di as txt "     Number of obs:        " as res %10.0f `e(N)'
		di as txt "     Number of groups:     " as res %10.0f `e(N_g)'
		di as txt "     Min obs per group:    " as res %10.0f `e(g_min)'
		di as txt "     Max obs per group:    " as res %10.0f `e(g_max)'
					
		ereturn display
		if ("`excluded_inst'"=="") {
			disp as text "No excluded instruments - standard QRPD estimation."
		}
		if ("`excluded_inst'"!="") {
			disp as text "Excluded instruments: " "`excluded_inst'"
		}
		
		di as txt "MCMC diagonstics:"
		di as txt "     Mean acceptance rate: " as res %10.3f `e(arates_mean)'
		di as txt "     Total draws:          " as res %10.0f `e(draws)'
		di as txt "     Burn-in draws:        " as res %10.0f `e(burn)'
		di as txt "     Draws retained:       " as res %10.0f `e(draws_retained)'
		di as txt "     Value of objective function:   "
		di as txt "             Mean:         " as res %12.4f `e(f_mean)'
		di as txt "             Min:          " as res %12.4f `e(f_min)'
		di as txt "             Max:          " as res %12.4f `e(f_max)'
		di as txt "MCMC notes:"
		if `e(thin)'!=1		{
		di as txt "     *One of every " %1.0f `e(thin)' " draws kept"
					} 
		if "`e(jumble)'"=="jumble" {
		di as txt "     *Draws Jumbled"
					} 
		if "`usemax'"!="" {
			di as txt  "     *Point estimates correspond to max value of obj function."
		}
		else {
			di as txt  "     *Point estimates correspond to mean of draws."
		}
		
        if "`analytic'"!="" {
            di as txt "     *Standard errors are derived analytically."
        }
        else {
            di as txt "     *Standard errors are derived from variance of draws."
        }
		capture mata: mata drop G M
    }

    /* grid search */    
    else if "`optimize'"=="grid" {
		di "Grid-search optimization"
		
		if wordcount("`rhs'") > 2 {
			di as error "Error: The # of RHS variables is > 2. Grid search not feasible."
			exit
		}

		if wordcount("`rhs'") == 2 {
			
			if "`grid1'"!="" {
				local grid1: subinstr local grid1 " " ", ", all
				mata: grid1=(`grid1')'
			}
			else {
				mata: min1=beta_start[1,1]-1.96*beta_start_se[1,1]
				mata: max1=beta_start[1,1]+1.96*beta_start_se[1,1]
				mata: intvl1=(min1-max1)/8
				mata: grid1=range(min1,max1,intvl1)
				mata: st_numscalar("min1",min1)
				mata: st_numscalar("max1",max1)
				mata: st_numscalar("intvl1",intvl1)
				local grid1 "`=min1'(`=intvl1')`=max1'"
				numlist "`grid1'"
			}
			
			if "`grid2'"!="" {
				local grid2: subinstr local grid2 " " ", ", all
				mata: grid2=(`grid2')'
			}
			else {
				mata: min2=beta_start[1,2]-1.96*beta_start_se[1,2]
				mata: max2=beta_start[1,2]+1.96*beta_start_se[1,2]
				mata: intvl2=(min2-max2)/8
				mata: grid2=range(min2,max2,intvl2)
				mata: st_numscalar("min2",min2)
				mata: st_numscalar("max2",max2)
				mata: st_numscalar("intvl2",intvl2)
				local grid2 "`=min2'(`=intvl2')`=max2'"
				numlist "`grid2'"
			}
		
			mata: expand_grid1=J(rows(grid2),1,grid1)
			mata: expand_grid2=vec(J(rows(grid1),1,grid2'))
			mata: grid12=(expand_grid1, expand_grid2)
			mata: st_numscalar("grid_points",rows(grid12))
			di as txt "Number of grid points: " `=grid_points'
			mata: grid12=grid12,J(rows(grid12),1,.)
			mata: grid_index=0
			
			local iter=0
			qui foreach num2 of numlist `grid2' { 
				qui foreach num of numlist `grid1' {
					if "`noisy'"!="" {
						local iter=`iter'+1
						noisily di as txt "." _c
						if mod(`iter',50)==0 {
							noisily di as txt "`iter'" 
						}
					}
					mata: grid_index=grid_index+1					
					mata: b_grid=grid12[grid_index,1..2]
					mata: gamma=estimateGamma(b_grid,yi,di,wt,tau,fix)
					//mata: below=wt:*(yi:<=(gamma:+di*b_grid')):/Ti
					mata: below=wt:*(yi:<=(gamma:+di*b_grid'))
					mata: g=cross(below,Znew)/sqrt(rows(info))
					mata: grid12[grid_index,3]=-0.5*g*A*g'
				}
			}
			mata: grid_max=findRowMax(grid12[.,3],grid12)
			mata: b1=(min(uniqrows(grid_max[.,1]))+max(uniqrows(grid_max[.,1])))/2
			mata: b2=(min(uniqrows(grid_max[.,2]))+max(uniqrows(grid_max[.,2])))/2
			mata: b=b1,b2
			mata: fval=grid_max[1,3]
		}	
		
		if wordcount("`rhs'") == 1 {
			if "`grid1'"!="" {
				local grid1: subinstr local grid1 " " ", ", all
				mata: grid1=(`grid1')'
			}
			else {
				mata: min1=beta_start[1,1]-1.96*beta_start_se[1,1]
				mata: max1=beta_start[1,1]+1.96*beta_start_se[1,1]
				mata: intvl1=(min1-max1)/80
				mata: grid1=range(min1,max1,intvl1)
				mata: st_numscalar("min1",min1)
				mata: st_numscalar("max1",max1)
				mata: st_numscalar("intvl1",intvl1)
				local grid1 "`=min1'(`=intvl1')`=max1'"
				numlist "`grid1'"
			}	
			
			mata: st_numscalar("grid_points",rows(grid1))
			di as txt "Number of grid points: " `=grid_points'
			mata: grid1=grid1,J(rows(grid1),1,.)
			mata: grid_index=0
			
			local iter=0
			qui foreach num of numlist `grid1' {
				if "`noisy'"!="" {
					local iter=`iter'+1
					noisily di as txt "." _c
					if mod(`iter',50)==0 {
						noisily di as txt "`iter'" 
					}
				}
				mata: grid_index=grid_index+1
				mata: b_grid=grid1[grid_index,1]
				mata: gamma=estimateGamma(b_grid,yi,di,wt,tau,fix)
				mata: below=wt:*(yi:<=(gamma:+di*b_grid'))
				mata: g=cross(below,Znew)/sqrt(rows(info))
				mata: grid1[grid_index,2]=-0.5*g*A*g'
			}
			mata: grid_max=findRowMax(grid1[.,2],grid1)
			mata: b=(min(uniqrows(grid_max[.,1]))+max(uniqrows(grid_max[.,1])))/2
			mata: fval=grid_max[1,2]
		}	
		
		tempname grid_betas
		mata: st_matrix("`grid_betas'",grid_max)
		mat colnames `grid_betas' = `rhs' obj_fun_value
		
		mata: deriv= beta_start_se		
		mata: V=qrpd_var(yi,di,b,tau,fix,info,zi,wt,deriv,id,Znew)
		
		mata: gam=estimateGamma(b,yi,di,wt,tau,fix)
		mata: uniq_gam=uniqrows((fix,gam))						/* in "no constant" form */
		mata: constant=uniq_gam[rows(uniq_gam),2]				/* let constant be last value */
		tempname b V gamma
		
		mata: st_matrix("`b'",b)
		mata: st_matrix("`V'",V)
		mata: st_matrix("`gamma'", uniq_gam)	
		
		mat colnames `b' = `rhs' 
		mat colnames `V' = `rhs' 
		mat rownames `V' = `rhs' 
		mat colnames `gamma' = `fix' estimate
		
		tempname N_g N g_min g_max fval
		mata: st_numscalar("`N_g'",panelstats(info)[1])
		mata: st_numscalar("`N'",panelstats(info)[2])
		mata: st_numscalar("`g_min'",panelstats(info)[3])
		mata: st_numscalar("`g_max'",panelstats(info)[4])
		local N_g=`N_g'
		local N=`N'
		local g_min=`g_min'
		local g_max=`g_max'
		mata: st_numscalar("`fval'",fval)
	
		ereturn clear
		ereturn post `b' `V', esample(`touse') obs(`N') 
		
		ereturn local title "Quantile Regression for Panel Data (QRPD)"
		ereturn local cmd "qregpd"
		ereturn local indepvars `rhs'
		ereturn local depvar `lhs'
		ereturn matrix solutions `grid_betas'
		ereturn matrix gamma `gamma'
		ereturn scalar N=`N'
		ereturn scalar N_g=`N_g'
		ereturn scalar g_min=`g_min'
		ereturn scalar g_max=`g_max'
		ereturn scalar fval=`fval'
		
		di _newline
		di as txt "`e(title)'" 
		di as txt "     Number of obs:        " as res %10.0f `e(N)'
		di as txt "     Number of groups:     " as res %10.0f `e(N_g)'
		di as txt "     Min obs per group:    " as res %10.0f `e(g_min)'
		di as txt "     Max obs per group:    " as res %10.0f `e(g_max)'
		ereturn display

		if ("`excluded_inst'"=="") {
			disp as text "No excluded instruments - standard QRPD estimation."
		}
		if ("`excluded_inst'"!="") {
			disp as text "Excluded instruments: " "`excluded_inst'"
		}		
		
		if rowsof(e(solutions)) > 1 {
			di as txt "Note: Alternative solutions exist. See e(solutions)."
		}
		di as txt "Value of objective function: " as res %12.8f `fval'
		
		
	} 
	
/* Nelder-Mead */
	else {
		di "Nelder-Mead optimization"
		mata: moptimize(M)	
		mata: b=moptimize_result_coefs(M)
		
		mata: deriv= beta_start_se	
		mata: V=qrpd_var(yi,di,b,tau,fix,info,zi,wt,deriv,id,Znew)
		
		mata: gam=estimateGamma(b,yi,di,wt,tau,fix)
		mata: uniq_gam=uniqrows((fix,gam))						/* in "no constant" form */
		mata: constant=uniq_gam[rows(uniq_gam),2]				/* let constant be last value */

		tempname b V gamma fval
		
		mata: st_matrix("`b'",b)
		mata: st_matrix("`V'",V)
		mata: st_matrix("`gamma'", uniq_gam)	
		mata: st_numscalar("`fval'",moptimize_result_value(M))
		
		mat colnames `b' = `rhs' 
		mat colnames `V' = `rhs' 
		mat rownames `V' = `rhs' 
		mat colnames `gamma' = `fix' estimate
		
		tempname N_g N g_min g_max
		mata: st_numscalar("`N_g'",panelstats(info)[1])
		mata: st_numscalar("`N'",panelstats(info)[2])
		mata: st_numscalar("`g_min'",panelstats(info)[3])
		mata: st_numscalar("`g_max'",panelstats(info)[4])
		local N_g=`N_g'
		local N=`N'
		local g_min=`g_min'
		local g_max=`g_max'
	
		ereturn clear
		ereturn post `b' `V', esample(`touse') obs(`N') 
		
		ereturn local title "Quantile Regression for Panel Data (QRPD)"
		ereturn local cmd "qregpd"
		ereturn local indepvars `rhs'
		ereturn local depvar `lhs'
		ereturn matrix gamma `gamma'
		ereturn scalar N=`N'
		ereturn scalar N_g=`N_g'
		ereturn scalar g_min=`g_min'
		ereturn scalar g_max=`g_max'
		ereturn scalar fval=`fval'

		di _newline
		di as txt "`e(title)'" 
		di as txt "     Number of obs:        " as res %10.0f `e(N)'
		di as txt "     Number of groups:     " as res %10.0f `e(N_g)'
		di as txt "     Min obs per group:    " as res %10.0f `e(g_min)'
		di as txt "     Max obs per group:    " as res %10.0f `e(g_max)'

		ereturn display
		if ("`excluded_inst'"=="") {
			disp as text "No excluded instruments - standard QRPD estimation."
		}
		if ("`excluded_inst'"!="") {
			disp as text "Excluded instruments: " "`excluded_inst'"
		}		
	
	}
end

	
program Replay
	syntax 
	ereturn display
end	

mata:
struct qrpdInfo {
	real matrix yi,Znew,di,A,info,fix,zi,wt,id,Ti
	real scalar tau
} 
struct qrpdInfo qrpdInfoInit(yi,Znew,di,A,info,fix,zi,wt,tau,id,Ti)
{
	struct qrpdInfo scalar G
	G.yi=yi
	G.Znew=Znew
	G.di=di
	G.A=A
	G.info=info
	G.fix=fix
	G.zi=zi
	G.wt=wt
	G.tau=tau
	G.id=id
	G.Ti=Ti
	return(G)
}

void objectiveFunGmm(M,todo,b,obj,h,H) 
{
	real scalar tau
	real matrix yi,Znew,di,info,fix,zi,wt,A,id,Ti
	struct qrpdInfo scalar G
	
	real matrix  gamma,below, g
	
	yi=moptimize_util_depvar(M,1)
	G=moptimize_util_userinfo(M,1)
	
	gamma=estimateGamma(b,yi,G.di,G.wt,G.tau,G.fix)
	below=G.wt:*(yi:<=(gamma:+G.di*b'))
	g=cross(below,G.Znew)/sqrt(rows(G.info))
	obj=-0.5*g*G.A*g'	
	
}


real matrix estimateGamma(	real matrix b,
							real matrix yi, 
							real matrix di, 
							real matrix wt, 
							real scalar tau, 
							real matrix fix)
{
	real matrix Gamma, u, fe_id, u_fe, wt_fe
	real scalar num_fe, gamma_fe, fe
	real matrix dum

	u = yi:-di*b'
	fe_id = uniqrows(fix)
	num_fe = rows(fe_id)
	Gamma = J(rows(di), 1, 0)
	dum = J(rows(di), num_fe, .)
	for (fe=1; fe<=num_fe; fe++) {
		u_fe=select(u, fix:==fe_id[fe])
		wt_fe=select(wt, fix:==fe_id[fe])
		gamma_fe=mm_quantile(u_fe, wt_fe, tau)
		dum[.,fe]=gamma_fe*(fix:==fe_id[fe])
		Gamma=Gamma+dum[.,fe]
	}
	return(Gamma)
}
	

real matrix weightMatrix(	real rowvector b,
							real scalar tau,
							real colvector yi,
							real matrix di,
							real matrix Znew,
							real colvector id,
							real matrix fix,
							real matrix zi,
							real colvector wt,
							string scalar type) 
{
	real scalar n,i
	real colvector info,gamma,below
	real matrix git,gi,A
	
	
	info=panelsetup(id,1)
	n=rows(info)
	
	if (type=="efficient") {
	
		gamma=estimateGamma(b,yi,di,wt,tau,fix)
		gi=J(rows(info),cols(zi),.)		
		for (i=1; i<=rows(info); i++) {
			dip=panelsubmatrix(di,i,info)
			zip=panelsubmatrix(zi,i,info)
			yip=panelsubmatrix(yi,i,info)
			wtp=panelsubmatrix(wt,i,info)
			gammap=panelsubmatrix(gamma,i,info)
			zim=panelsubmatrix(Znew,i,info)
			gi[i,.]=colsum(wtp:*zim:*(yip:<=(gammap:+dip*b')))
			//gi[i,.]=colsum(wtp:*zim:*(yip:<=(gammap:+dip*b')))/rows(dip)
		}
		A=invsym((1/rows(yi))*cross(gi,gi))
	}
	
	else {
		A=(1/(tau*(1-tau)))*invsym((1/rows(yi))*cross(Znew,Znew))	
	}
	
	return(A)
}
	

void gqreg_process_draws(real matrix Draws,
							 real matrix vals,
							 string scalar jumble,
							 real scalar thin,
							 real scalar burn,
							 transmorphic betaf,
							 transmorphic valsf)
{
	real scalar i
	real matrix ordVec

	Draws = Draws[burn+1::rows(Draws),]

	if (thin!=1) {
		keep=J(0,1,.)
		for (i=thin;i<=rows(Draws);i=i+thin) keep=keep \ i
	}
	else keep=(1::rows(Draws))

	if (jumble=="jumble") _jumble(keep)
	betaf=Draws[keep,]
	valsf=vals[keep,]
}	

// QRPD analytic standard errors.
real matrix qrpd_var(yi,di,beta,tau,fix,info,zi,wt,deriv,id,Znew)
{
	real scalar bdw,i
	real matrix dib,gamma,uu,iqr,meanuu,sqruu,pick,stderr,meansqr 
	real matrix beta_plus,gamma_psi_plus,beta_minus, gamma_psi_minus
	real matrix yearderiv,x_plus_dgdd,abs_u,abs_u_below_bdw
	real matrix G1,G2,G3,G,below,S1,S2,g0,g2,Sigma,variance,A

	dib=di*beta'
	
	info=panelsetup(id,1)

	gamma=estimateGamma(beta,yi,di,wt,tau,fix)		 					

	uu = yi:-dib:-gamma
	iqr=mm_iqrange(uu, wt)
	meanuu=colsum(uu:*wt)/colsum(wt)
	sqruu=((uu:*wt):-meanuu):^2
	meansqr=(colsum(sqruu)/(colsum(wt)-1))^.5
	pick=(meansqr,iqr/1.349)
	stderr=rowmin(pick)
	bdw=1.0589*stderr*(rows(yi))^(-.2)
	
	yearderiv=J(rows(di),cols(di),0)
	for (l=1; l<=cols(di); l++) {
		beta_plus=beta
		beta_plus[l]=beta[l]:+deriv[l]
		gamma_psi_plus=estimateGamma(beta_plus, yi,di,wt,tau,fix)
		beta_minus=beta
		beta_minus[l]=beta[l]:-deriv[l]
		gamma_psi_minus=estimateGamma(beta_minus,yi,di,wt,tau,fix)
		yearderiv[.,l]=(gamma_psi_plus:-gamma_psi_minus):/(2:*deriv[l])
	}

	x_plus_dgdd=di:+yearderiv	
	
	abs_u = abs(yi:-dib:-gamma)
	abs_u_below_bdw=(abs_u:<= bdw)
	
	A=weightMatrix(beta,tau,yi,di,Znew,id,fix,zi,wt,"efficient")
	
	wtZnew=wt:*Znew
	
	G1=panelsubmatrix(x_plus_dgdd,1,info)			
	G2=panelsubmatrix(wtZnew,1,info)
	G3=panelsubmatrix(abs_u_below_bdw,1,info)	
	G=cross(G2,G3,G1)
	for (i=2; i<=rows(info); i++) {
		G1=panelsubmatrix(x_plus_dgdd, i, info)			
		G2=panelsubmatrix(wtZnew, i, info)
		G3=panelsubmatrix(abs_u_below_bdw,i,info)
		G=cross(G2,G3,G1)+G
		//G=cross(G2,G3,G1)/rows(G1)+G 
	}	
	G=((1/2)*(1/bdw)):*G:*(1/rows(info))

	below=(yi:<=(gamma:+di*beta'))
	S1=panelsubmatrix(below, 1, info)
	S2=panelsubmatrix(wtZnew,1,info)
	g0=(cross(S1,S2))
	g2=cross(g0,g0)
	for (i=2; i<=rows(info); i++) {
		S1=panelsubmatrix(below, i, info)			
		S2=panelsubmatrix(wtZnew, i, info)					
		g0=(cross(S1,S2))
		g2=cross(g0,g0)+g2
	}
	
	Sigma=g2/rows(info)
	
	variance=invsym(G'*A*G)*(G'*A*Sigma*A*G)*invsym(G'*A*G)/rows(info)
	
	return(variance)
}

real matrix demeanX(real colvector id,
					real matrix X,
					real colvector wt)
{
	real matrix Xmean,Xdemean
	real colvector Ti			

	mm_panels(id,Ti=.)				/* num of t per i */
	Xmean=mm_collapse(X,wt,id)		/* collapses into a Nx1 vec */
	Xmean=mm_expand(Xmean,Ti,1,1) 	/* expands into NTx1 */
	Xmean=Xmean[.,2..cols(Xmean)]	/* drop ids in first col */
	Xdemean=X:-Xmean
	return(Xdemean)
}

real matrix findRowMax(	real colvector X1,
						real matrix X2)
{
	real colvector row_index
	real matrix trash
	
	row_index=.; trash=.; maxindex(X1,1,row_index,trash)
	RowMax=X2[row_index,.]
	
	return(RowMax)
}

end