*! Version 1.0.2 08august2022
/*
	Syntax:
	hmindex, <Price> <Quantity> [<op1> <op2> <op3> <op4> <op5> <op6>]  
		
		<op1> := AXiom[(string)]
			eWGARP	   (default)
			eWARP
		<op2> := power
		<op3> := SEED[(real 12345)]
		<op4> := SIMulations[(real 1000)]
		<op5> := IDvar[(string)]
		<op6> := GENerate[(string)]
*/
/* ================================================================== */
* Main program
capture program drop hmindex
program hmindex, rclass sortpreserve
	version 15.1
	
	syntax, Price(string) Quantity(string) ///
			[AXiom(string)			///
			 DISTribution   		///
			 SIMulations(real 1000) ///
			 SEED(real 12345)		///
			 GENerate(string)		///
			 ID(string)]
	
	
				*******************************
				*** Checking data structure ***
				***   and syntax validity   ***
				*******************************

	* Is there at least one product with a non-zero quantity in every given period?
	mata: st_numscalar("r(nrz)", (min(rowsum(st_matrix("`quantity'")))==0 ? 1 : 0))
	if `r(nrz)' == 1 {
		display as error 	" The data contains observations" ///
		as error 			" with zero (or missing) quantities."
		exit 459 /* "Something that should be true of your data is not" error */
	}
	* Are prices and quantity data of the same dimension?
	mata: st_numscalar("r(OK)", (rows(st_matrix("`price'")) == rows(st_matrix("`quantity'")) ///
		& cols(st_matrix("`price'")) == cols(st_matrix("`quantity'")) ? 1 : 0))
	if `r(OK)' == 0 {
		display as error 	"Invalid matrix dimensions."
		exit 459 /* "Something that should be true of your data is not" error */
	}
	* Are prices strictly positive?
	mata: st_numscalar("r(PNSP)", (min(st_matrix("`price'")) < 0 ? 1 : 0))
	if `r(PNSP)' > 0 {
		display as error 	" The price matrix contains non-positive values."
		exit 411 /* "Nonpositive values encountered" error */
	}
	
	* Which axiom(s) does the user want to check?
	* And creating necessary scalars, vectors and tempnames accordingly.
	local tempname_prefix hm min_hm mean_hm median_hm max_hm std_hm q1_hm q3_hm sim rawResults
	local axiom = lower("`axiom'")
	if ("`axiom'"=="")		local axiom wgarp	/* WGARP set to default */
	
	if ("`axiom'"=="all")	local axiom wgarp warp
	tokenize `axiom'
	local axioms "`1' `2'"
	foreach ax of local axioms {
		if !inlist("`ax'", "wgarp", "warp") {
			display as error 	" Axiom() must be either WGARP or WARP; case-insensitive."
			display as error 	" If not specified, the default setting for" ///
								" Axiom() is WGARP."
			exit 198 /* "Invalid syntax --> range invalid" error */
		}
		
		else if inlist("`ax'", "wgarp", "warp") {
			
			foreach temp of local tempname_prefix {
				tempname `temp'_`ax'
			}
 		}

 		if "`distribution'" != "" {
 			matrix `sim_`ax'' = J(`simulations', 2,.)
 			matrix colname `sim_`ax'' = HM_NUM HM_FRAC

 		}

	}
						**************
						*** AXIOMS ***
						**************
	local goods `=colsof(`price')'
	local obs 	`=rowsof(`price')'
	
	local first_ax = 1

 	local allAxiomsDisplay ""

 	tempname rawResults generalInfoTable sumStatsTable ind	/* Temp names for tables */

	foreach ax of local axioms {
		capture matrix drop pnondrop
		capture matrix drop qnondrop
		
		* Creating Quantity Matrix with a sequence vector
		mata: st_matrix("seq", range(1,`obs',1))
		matrix qnondrop = (seq, `quantity')
		
		* Creating Price Matrix with a zero vector
		matrix pnondrop = (J(`obs',1,0), `price')

		mata: test_`ax'("pnondrop", "qnondrop")
		
		if `r(PASS)' == 1 {
		
			return scalar HM_`ax' = `=rowsof(pnondrop)'
			
		}
		
		else if `r(PASS)' != 1{
		
			matrix pdrop = J(1,`goods' + 1,.)
			matrix qdrop = J(1,`goods' + 1,.)
	
			local ind = 0
			while `ind' == 0 {
	
				mata: adjmat_`ax'("pnondrop", "qnondrop")

				mata: remove("AdjMat", "pnondrop", "qnondrop")
				
				
				* TO DO: capture only if pd doesn't exist
				capture matrix pdrop = (pdrop\pd)
				capture matrix qdrop = (qdrop\qd)

				mata: test_`ax'("pnondrop", "qnondrop")

				if `r(PASS)' == 1 {
					local ind = 1
				}
				
			} /* Ends while loop */
		
		} /* Ends else if */
		
		capture matrix drop pd
		capture matrix drop qd
		capture matrix drop AdjMat
		
		* Sorting matrices
		matrix dnonsort = [qnondrop, pnondrop]
		mata : st_matrix("dnonsort", sort(st_matrix("dnonsort"), 1))
		matrix qnondrop = dnonsort[1...,2..`goods' + 1]
		matrix pnondrop = dnonsort[1...,`goods' + 3..`=colsof(dnonsort)']
			
		* If matrix pdrop exists (because at least 1 observation was dropped),
		* then clean and sort pdrop and qdrop matrices. Otherwise, skip.
		capture confirm matrix pdrop
		
		if _rc == 0 {
			matrix pdrop = pdrop[2...,1...]
			matrix qdrop = qdrop[2...,1...]
			matrix dsort = [qdrop, pdrop]
			mata : st_matrix("dsort", sort(st_matrix("dsort"), 1))
			matrix qdrop = dsort[1...,2..`goods' + 1]
			matrix pdrop = dsort[1...,`goods' + 3..`=colsof(dsort)']
			
			* Creating column vector with identifier for dropped observations
			matrix obsdrop = dsort[1..., 1]

			* Creating indicator column vector
			matrix I = J(`obs',1, 1)
			forvalues i = 1/`=rowsof(obsdrop)' {
				matrix I[obsdrop[`i',1], 1] = 0
			}
		}
	

		local FC_`ax': di %3.2f scalar((`=rowsof(pnondrop)'/`obs'))
		local FRAC_HM_`ax' = `FC_`ax''
		
		** Creating output & return list tables
		local axiomDisplay = upper("`ax'")
		local allAxiomsDisplay "`allAxiomsDisplay' `axiomDisplay'"
			
		* Creating output table
		matrix `rawResults_`ax'' = `=rowsof(pnondrop)', `FRAC_HM_`ax''
		matrix rowname `rawResults_`ax'' = "`axiomDisplay'"
		if 		`first_ax' == 1		matrix `rawResults' = `rawResults_`ax''
		else if `first_ax'  > 1		matrix `rawResults' = `rawResults' \ `rawResults_`ax''
			
		local first_ax = `first_ax' + 1
		
		
		* Return list for several axioms
		return scalar HM_NUM_`axiomDisplay' 		= `=rowsof(pnondrop)'
		return scalar HM_FRAC_`axiomDisplay'		= `FRAC_HM_`ax''
		
		* If at least one observation was drop, return violater sets
		if `=rowsof(pnondrop)' < `obs' {
			return matrix VS_price_`axiomDisplay'		= pdrop
			return matrix VS_quantity_`axiomDisplay'	= qdrop
			
			return matrix OBSDROP_`axiomDisplay'		= obsdrop
			return matrix INDICATOR_`axiomDisplay'		= I
		} 
		
		return matrix CS_price_`axiomDisplay'		= pnondrop
		return matrix CS_quantity_`axiomDisplay'	= qnondrop				

 	} /* Ends foreach ax */

 	********************
 	*** DISTRIBUTION ***
 	********************
 	if "`distribution'" != "" {
		
 		tempname gmat te
 		mata: newGMAT("`price'","`quantity'", `seed', `simulations')
	
 		matrix `gmat' = gamma_matrix
 		matrix `te' = total_expenditure
 		local rK = r(K)
 		local rT = r(T)

 		forvalues i = 1(1)`simulations' {

 			mata: genXS("`gmat'", `rT', `rK', "`te'", "`price'", `i')

 			foreach ax of local axioms {

 				local axiomDisplay = upper("`ax'")

 				** HMINDEX results
 				quietly hmindex, price("`price'") ///
 								 quantity("simulated_quantities") ///
 								 axiom("`ax'") 
 				quietly return list
				
 				* Number & fraction of violations
 				local HM_NUM = r(HM_NUM_`axiomDisplay')
 				local HM_FRAC = r(HM_FRAC_`axiomDisplay')

 				matrix `sim_`ax''[`i',1] = `HM_NUM'
 				matrix `sim_`ax''[`i',2] = `HM_FRAC'

 			} /* Ends foreach ax */

 		} /* Ends forvalues i */

 	} /* Ends if statement */

 	* Return list common for all axioms
 	return scalar OBS						= `obs'
 	return scalar GOODS						= `goods'
 	return local  AXIOM						"`allAxiomsDisplay'"

 	* Displaying output table
 	if "`distribution'" == "" {
 		matrix `generalInfoTable' = `obs', `goods'
 		matrix `generalInfoTable' = `generalInfoTable''
 		matrix rowname `generalInfoTable' = "	Number of obs		= " ///
 											"	Number of goods		= " 
 	}
 	else if "`distribution'" != "" {
 		matrix `generalInfoTable' = `obs', `goods', `simulations'
 		matrix `generalInfoTable' = `generalInfoTable''
 		matrix rowname `generalInfoTable' = "	Number of obs		= " ///
 											"	Number of goods		= " ///
 											"	Simulations 		= " 
 	}

 	matrix colname `generalInfoTable' = " "
 	matlist `generalInfoTable', border(none) lines(none) ///
 			format(%7.2g) names(rows) left(0) twidth(30)

 	matrix colnames `rawResults' = #HM %HM
 	matlist `rawResults', border(top bottom) rowtitle("Axiom")

 	if "`distribution'" != "" {

 		di " "
 		di as text "Summary statistics for simulations:"	

 		local allAxiomsDisplay ""
 		foreach ax of local axioms {

 			local axiomDisplay = upper("`ax'")

 			local allAxiomsDisplay "`allAxiomsDisplay' `axiomDisplay'"

 			* Summary stats table
 			tempvar HM_NUM HM_FRAC

 			mata: A = st_matrix("`sim_`ax''")

 			getmata (`HM_NUM' `HM_FRAC') = A, force

 			quietly tabstat `HM_NUM' `HM_FRAC', stat(mean sd min p25 ///
 					median p75 max) save

 			quietly return list 

 			matrix `sumStatsTable' = r(StatTotal)
 			matrix colnames `sumStatsTable' =	"#HM" "%HM"
 			matrix rownames `sumStatsTable' =	Mean "Std. Dev." Min ///
 													Q1 Median Q3 Max


 			matlist `sumStatsTable',	border(top bottom) ///
 										rowtitle("`axiomDisplay'")

 			* Return list for several axioms
 			return scalar SIM						= `simulations'
 			return local  AXIOM						"`allAxiomsDisplay'"
 			return matrix SIMRESULTS_`axiomDisplay'	= `sim_`ax''
 			return matrix SUMSTATS_`axiomDisplay'	= `sumStatsTable'

 		} /* Ends foreach ax */

 	}	/* Ends if statement */

 end

  
mata
/* ================================================================== */
					/* AXIOM 1: WGARP */
function adjmat_wgarp(string P_temp, string X_temp)
{
	P_mat = st_matrix(P_temp)
	X_mat = st_matrix(X_temp)
	
	T = rows(P_mat)
	
	// Create empty matrices DRP and SDRP of size T x T
	DRP = J(T, T, 0)
	SDRP = J(T, T, 0)
		// Looping over i and j
	for (i=1; i<= T; i++) {
		for (j=1; j<= T; j++) {
			if (P_mat[i,.] * (X_mat[i,.])' >= P_mat[i,.] * (X_mat[j,.])') {
			
				DRP[i,j] = 1
				
			  if (P_mat[i,.] * (X_mat[i,.])' > P_mat[i,.] * (X_mat[j,.])') {
					
					SDRP[i,j] = 1
				
				}		/* Ends if > */
			
			}		/* Ends if >= */
			
		}		/* Ends for j */
		
	}		/* Ends for i */
	
	
	// Computing edges
	A = J(T, T, 0)
	
	for (i=1; i<= T; i++) {
		
		for (j=1; j<= T; j++) {
			
			if (DRP[i,j]==1 & SDRP[j,i]==1) {
				
				A[i,j] = 1
				
			} /* Ends if */
		
		} /* Ends for j */
		
	}	/* Ends for i */

		
	st_matrix("AdjMat", A)
	
}
					/* AXIOM 1: WGARP */
/* ================================================================== */
/* ================================================================== */
					/* AXIOM 2: WARP */
function adjmat_warp(string P_temp, string X_temp)
{
	P_mat = st_matrix(P_temp)
	X_mat = st_matrix(X_temp)
	T = rows(P_mat)
	// Create empty matrices DRP of size T x T
	DRP = J(T, T, 0)		
	
	// Looping over i and j
	for (i=1; i<= T; i++) {
		for (j=1; j<= T; j++) {
			if (P_mat[i,.] * (X_mat[i,.])' >= P_mat[i,.] * (X_mat[j,.])') {
			
				DRP[i,j] = 1
				
			}		/* Ends if >= */
		}		/* Ends for j */
		
	}		/* Ends for i */
	// Computing edges
	A = J(T, T, 0)
	for (i=1; i<= T; i++) {
		
		for (j=1; j<= T; j++) {
		
			if (X_mat[i,.] != X_mat[j,.]) {
				if (DRP[i,j]==1 & DRP[j,i]==1) {
					
					A[i,j] = 1
					
				} /* Ends if */
			
			} /* Ends if */
			
		} /* Ends for j */
		
	}	/* Ends for i */
	st_matrix("AdjMat", A)
	
}
					/* AXIOM 2: WARP */
/* ================================================================== */
/* ================================================================== */
					/* TEST WGARP*/
void test_wgarp(string P_temp, string X_temp)
{
	real scalar PASS
	P_mat = st_matrix(P_temp)
	X_mat = st_matrix(X_temp)
	T = rows(P_mat)

	
	// Create empty matrices DRP and SDRP of size T x T
	DRP = J(T, T, 0)		/* J(rows, columns, values) */
	SDRP = J(T, T, 0)
	
	// Looping over i and j
	for (i=1; i<= T; i++) {
		for (j=1; j<= T; j++) {
			if (P_mat[i,.] * (X_mat[i,.])' >= P_mat[i,.] * (X_mat[j,.])') {
			
				DRP[i,j] = 1
				
			  if (P_mat[i,.] * (X_mat[i,.])' > P_mat[i,.] * (X_mat[j,.])') {
					
					SDRP[i,j] = 1
				
				}		/* Ends if > */
			
			}		/* Ends if >= */
			
		}		/* Ends for j */
		
	}		/* Ends for i */
	
	
	// Search for WARP violations
	PASS = 1
	for (i=1; i<= T; i++) {
		for (j=1; j<= T; j++) {
		
			if (j > i) {
				if (DRP[i,j] == 1 && SDRP[j,i] == 1) {
				
					PASS = 0
					
					break
					
					}		/* Ends if */
					
				else if (DRP[j,i] == 1 && SDRP[i,j] == 1) {
					
					PASS = 0
					
					break
					
					}		/* Ends if */
					
				}		/* Ends if */
			
			}		/* Ends for j */
	
		}		/* Ends for i */
		
	// Returning results
	st_numscalar("r(PASS)", PASS)
}
		
			
/* ================================================================== */
/* ================================================================== */
					/* TEST WARP*/
void test_warp(matrix P_temp, matrix X_temp)
{
	real scalar PASS
	P_mat = st_matrix(P_temp)
	X_mat = st_matrix(X_temp)
	T = rows(P_mat)
	
	// Create empty matrices DRP and SDRP of size T x T
	DRP = J(T, T, 0)		/* J(rows, columns, values) */
	SDRP = J(T, T, 0)
	
	// Looping over i and j
	for (i=1; i<= T; i++) {
		for (j=1; j<= T; j++) {
			if (P_mat[i,.] * (X_mat[i,.])' >= P_mat[i,.] * (X_mat[j,.])') {
			
				DRP[i,j] = 1
			}		/* Ends if >= */
			
		}		/* Ends for j */
		
	}		/* Ends for i */
		
	// Search for WARP violations	
	PASS = 1
	for (i=1; i<= T; i++) {
		for (j=1; j<= T; j++) {
		
			if (j > i) {
				
				if (X_mat[i,.] != X_mat[j,.]) {
					if (DRP[i,j] == 1 && DRP[j,i] == 1) {
					
						PASS = 0
						
						break
						
					}	/* Ends if */
				
				}	/* Ends if */
					
			}	/* Ends if */
			
		}	/* Ends for j */
	
	}	/* Ends for i */
	
	// Returning results
	st_numscalar("r(PASS)", PASS)
	
}				
		
/* ================================================================== */
/* ================================================================== */
					/* Submodule 1: Remove */
function remove(string AdjMat, string P_temp, string X_temp)
{
	
	P_mat = st_matrix(P_temp)
	X_mat = st_matrix(X_temp)
	
	T = rows(P_mat)
	Columns = cols(P_mat)
	
	AM = st_matrix(AdjMat)
	
	degs = rowsum(AM)
	
	
	maxdegs = max(degs)
	
	d = selectindex((degs:==maxdegs))
		
	dropV = J(T, 1, 0)
	for (i=1; i<=rows(d); i++) {
		j = d[i]
		
		Ai = selectindex(AM[j,.]:>0)
		
		Aidegs = degs[Ai]
		Ai1 = Aidegs:==1

		/* 
		Procedure 1: This procedure looks at the case when the set of maximal
		nodes is not a singleton. Consider node i. If the degree of node i is 
		greater than the degrees of all nodes adjacent to i, 
		[i.e., degs(i)>max(Aidegs)], then we drop node i.
		*/
		
		if (degs[j]>max(Aidegs)) {	 // if (a)
			
			dropV[j] = 1
			
		}
		/*
		Procedure 2: This procedure looks at the case when the set of maximal
		nodes is not a singleton. Consider node i. If the degree of node i is
		equal to the degree of a node h adjacent to node i, 
		[i.e., degs(i)==degs(h)], then we have three cases; see below.
		*/
		
		else {
		
			AiRows = rows(Ai')
			for (h=1; h<=AiRows; h++) {
	
				k = Ai[h]'
				if (degs[j]==degs[k] && j<k) {	// if (b)
					
					Ah = selectindex(AM[k,.]:>0)
					
					Ahdegs = degs[Ah]
					Ah1 = Ahdegs:==1
					
					if (any(Ah:==j)) {			// if (c)

						if (rows(Ai1) != 0) {		// if (d)
							
							dropV[j] = 1
							
						}
						else if (rows(Ah1)!=0) {
							dropV[k] = 1
						
						}
						
						else if (rows(Ai1)==0 && rows(Ah1)==0) {
							
							dropV[j] = 1
							
						}	// Ends if (d)
					
					}	// Ends if (c)
											
				}	// Ends if (b)
			}	// Ends for h
			
		}	// Ends if (a)
	
	}	// Ends for i
	
	
	// Matrices
	pd = P_mat[selectindex(dropV:==1),.]
	qd = X_mat[selectindex(dropV:==1),.]
		
	pnondrop = P_mat[selectindex(dropV:==0),.]
	qnondrop = X_mat[selectindex(dropV:==0),.]
	
	st_matrix("pd", pd)
	st_matrix("qd", qd)
	st_matrix("pnondrop", pnondrop)	
	st_matrix("qnondrop", qnondrop)

}

					/* Submodule 1: Remove */				
/* ================================================================== */
/* ================================================================== */
					/* powerCalc */
					
function newGMAT(string P_temp, string X_temp, scalar seed, scalar S)					
{
	p = st_matrix(P_temp)
	x = st_matrix(X_temp)
	rseed(seed)							// setting the random seed 
	T = rows(p)							// # observations
	K = cols(p)							// # goods
	TE = (p:*x)*J(K, 1, 1);				// total expenditure
	GMAT = rgamma(T*K,S,1,1)			// (T*K)xS matrix of Gamma(1,1) random numbers
	
	st_matrix("gamma_matrix", GMAT)
	st_numscalar("r(T)", T)
	st_numscalar("r(K)", K)
	st_matrix("total_expenditure", TE)
}
	
function genXS(matrix GMAT, scalar T, scalar K, matrix TE, matrix p, scalar s)
{
	GMAT = st_matrix(GMAT)
	TE = st_matrix(TE)
	p = st_matrix(p)
	
	G = rowshape(GMAT[.,s],T);          // making a TxK matrix 
	D = G:/J(1, K, G*J(K,1,1))
	
	x_S = D:*(J(1, K, TE):/p)           // simulated quantities
	
	st_matrix("simulated_quantities", x_S)
}
					/* powerCalc */
/* ================================================================== */
end