*! 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