*******************************
* Mata funcitons: elasticities
*******************************

mata:

mata clear

// simple logit elasticities and diversion ratios
void elasticities_logit(
	string scalar market0,
	string scalar brand0,
	string scalar p0,
	string scalar q0,
	string scalar share0,
	string scalar alpha0,
	string scalar touse)
{

	// declarations
	real colvector market,brand,p,q,s,marketm,pm,qm,sm,brandm,alpha,alpham,sbm,pbm,missing_index
	real matrix amsumb,el,dr,weight,msumb,elm,qbm,weightm,elm0
	real scalar b,bb,m,elmij,j,i,_alpha

	// load data
	st_view(market, .,tokens(market0),touse)
	st_view(brand, .,tokens(brand0),touse)
	st_view(p, .,tokens(p0),touse)
	st_view(q, .,tokens(q0),touse)
	st_view(s, .,tokens(share0),touse)
	st_view(alpha, .,tokens(alpha0),touse)
	
	// reindexing group variable (categorical variables should start from 1 and increment by 1)
	market=reindex(market)
	brand=reindex(brand)

	// matrix for summations
	amsumb=J(rows(market),colmax(brand)-colmin(brand)+1,0)
	b=0
	for (bb=colmin(brand); bb<=colmax(brand); bb++) {
		b=b+1
		amsumb[.,b]=(brand:==bb)
	}

	el=J(colmax(brand),colmax(brand),0)												// brand level elasticity matrix
	dr=J(colmax(brand),colmax(brand),0)												// brand level diversion ratio matrix
	weight=J(colmax(brand),colmax(brand),0)											// weight matrix

	for (m=colmin(market); m<=colmax(market); m++) {								// calcualtions by markets

		// variables
		marketm=select(market,market:==m)
		pm=select(p,market:==m)
		qm=select(q,market:==m)
		sm=select(s,market:==m)
		brandm=select(brand,market:==m)
		alpham=select(alpha,market:==m)
		_alpha=mean(alpham)
	
		// brand matrix for summations
		msumb=select(amsumb,market:==m)

		// calculation of brand level elasticity matrix
		sbm=msumb'*sm																// brand level market shares
		pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')								// brand level weighted average prices
		qbm=msumb'*qm																// brand level quantities
		elm=J(colmax(brandm),colmax(brandm),1):*sbm'								// brand level elasticity matrix of market m
		elm=elm :- diag( J(rows(elm),1,1) )
		elm=_alpha*elm
		elm=elm:*pbm'															// transforming into elasticities: multiplying by derivating prices
		// changing missing values to zero (in unbalanced panels elasticities involving product(s) not present at the given market are missing)
		if (rows(elm)>1) {
			missing_index=(elm[.,1]:==.):*runningsum(J(rows(elm),1,1))
			missing_index=select(missing_index,missing_index:>0)
			elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
			elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
			missing_index=(elm[1,.]':==.):*runningsum(J(rows(elm),1,1))
			missing_index=select(missing_index,missing_index:>0)
			elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
			elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
		}

		// weigth matrix for market m
		weightm=qbm#J(1,rows(qbm),1)

		// adding up weights and weighted elasticities
		weight=weight:+weightm
		el=el:+(elm:*weightm)

	}	// end of m (markets) loop 

	el=el:/weight																	// normalizing: final, weighted average (across markets) elasticity matrix

	// calculation of loss diversations for average elasticites
	for (j=1; j<=colmax(brand); j++) {
		if (sum(brand:==j)!=0) {		// only if brand j is present
			for (i=1; i<=colmax(brand); i++) {
				if (sum(brand:==i)!=0) {		// only if brand i is present
					dr[i,j]=100*sum(q:*(brand:==i))*el[i,j]/abs(sum(q:*(brand:==j))*el[j,j])
				}
			}
		}
	}
	
	// storing in Stata
	st_matrix("el", el)
	st_matrix("dr", dr)
	
}	// end of elasticities_logit function
mata mlib add lrcl elasticities_logit()


// one-level nested logit elasticities and diversion ratios
void elasticities_nlogit(
	string scalar market0,
	string scalar brand0,
	string scalar g0,
	string scalar p0,
	string scalar q0,
	string scalar share0,
	string scalar sjg0,
	string scalar alpha0,
	string scalar sigmag0,
	string scalar touse)
{

	// declarations
	real colvector market,brand,g,p,q,s,sjg,marketm,pm,qm,sm,brandm,alpha,sigmag,alpham,sigmagm,gm,sjgm,sbm,pbm,scg,sjscg,missing_index,im
	real matrix amsumb,el,dr,weight,msumb,msumg,elm,qbm,weightm,elm0
	real scalar b,bb,m,elmij,j,i,gg,_alpha,_sigmag

	// load data
	st_view(market, .,tokens(market0),touse)
	st_view(brand, .,tokens(brand0),touse)
	st_view(g, .,tokens(g0),touse)
	st_view(p, .,tokens(p0),touse)
	st_view(q, .,tokens(q0),touse)
	st_view(s, .,tokens(share0),touse)
	st_view(sjg, .,tokens(sjg0),touse)
	st_view(alpha, .,tokens(alpha0),touse)
	st_view(sigmag, .,tokens(sigmag0),touse)
	
	// reindexing group variable (categorical variables should start from 1 and increment by 1)
	market=reindex(market)
	brand=reindex(brand)
	g=reindex(g)

	// matrix for summations
	amsumb=J(rows(market),colmax(brand)-colmin(brand)+1,0)
	b=0
	for (bb=colmin(brand); bb<=colmax(brand); bb++) {
		b=b+1
		amsumb[.,b]=(brand:==bb)
	}

	el=J(colmax(brand),colmax(brand),0)												// brand level elasticity matrix
	dr=J(colmax(brand),colmax(brand),0)												// brand level diversion ratio matrix
	weight=J(colmax(brand),colmax(brand),0)											// weight matrix

	for (m=colmin(market); m<=colmax(market); m++) {								// calcualtions by markets

		// variables
		marketm=select(market,market:==m)
		pm=select(p,market:==m)
		qm=select(q,market:==m)
		sm=select(s,market:==m)
		brandm=select(brand,market:==m)
		alpham=select(alpha,market:==m)
		sigmagm=select(sigmag,market:==m)
		gm=select(g,market:==m)
		sjgm=select(sjg,market:==m)
		_alpha=mean(alpham)
		_sigmag=mean(sigmagm)
	
		// matrices for summations
		msumb=select(amsumb,market:==m)
		im=runningsum(J(rows(gm),1,1))
		msumg=amsumf(im,gm)

		// calculation of brand level elasticity matrix
		qbm=msumb'*qm																// brand level quantities
		if (rows(uniqrows(brandm))==rows(brandm)) {									// if no grouping (no aggregation needed)
			elm0=J(rows(brandm),rows(brandm),1):*sm'
			elm0=elm0 :+ (_sigmag/(1-_sigmag))*(msumg):*(sjgm')
			elm0=elm0 :- diag( J(rows(elm0),1,(1/(1-_sigmag))) )
			elm0=_alpha*elm0
			elm0=elm0:*pm'															// transforming into elasticities: multiplying by derivating prices
			elm=J(rows(el),cols(el),0)
			elm[brandm,brandm]=elm0													// elasticity matrix of market m
		}
		if (rows(uniqrows(brandm))!=rows(brandm) & rows(uniqrows(brandm))!=1) {		// if grouping (aggregation needed)
			pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')							// brand level weighted average prices
			sbm=msumb'*sm															// brand level market shares
			elm=J(rows(brandm),1,1)*(sbm')											// brand level elasticity matrix of market m
			elm=elm :+ (_sigmag/(1-_sigmag))*(msumg:*(sjgm'))*msumb
			elm=elm :- msumb*(1/(1-_sigmag))
			elm=_alpha*elm:*sm
			elm=msumb'*elm
			elm=elm:/sbm															// transforming into semi elasticities: dividing by shares to be derivated
			elm=elm:*pbm'															// transforming into elasticities: multiplying by derivating prices 
			// changing missing values to zero (in unbalanced panels elasticities involving product(s) not present at the given market are missing)
			if (rows(elm)>1) {
				missing_index=(elm[.,1]:==.):*runningsum(J(rows(elm),1,1))
				missing_index=select(missing_index,missing_index:>0)
				elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
				elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
				missing_index=(elm[1,.]':==.):*runningsum(J(rows(elm),1,1))
				missing_index=select(missing_index,missing_index:>0)
				elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
				elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
			}
		}
		if (rows(uniqrows(brandm))!=rows(brandm) & rows(uniqrows(brandm))==1) {		// if market level elasticity
			pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')							// market level weighted average price
			elm=-_alpha*(1-sum(sm))*pbm
		}

		// weigth matrix for market m
		weightm=qbm#J(1,rows(qbm),1)

		// adding up weights and weighted elasticities
		weight=weight:+weightm
		el=el:+(elm:*weightm)

	}	// end of m (markets) loop 

	el=el:/weight																	// normalizing: final, weighted average (across markets) elasticity matrix

	// calculation of loss diversations for average elasticites
	for (j=1; j<=colmax(brand); j++) {
		if (sum(brand:==j)!=0) {		// only if brand j is present
			for (i=1; i<=colmax(brand); i++) {
				if (sum(brand:==i)!=0) {		// only if brand i is present
					dr[i,j]=100*sum(q:*(brand:==i))*el[i,j]/abs(sum(q:*(brand:==j))*el[j,j])
				}
			}
		}
	}

	// storing in Stata
	st_matrix("el", el)
	st_matrix("dr", dr)
	
}	// end of elasticities_nlogit function
mata mlib add lrcl elasticities_nlogit()


// two-level nested logit elasticities and diversion ratios
void elasticities_nlogit2(
	string scalar market0,
	string scalar brand0,
	string scalar g0,
	string scalar h0,
	string scalar p0,
	string scalar q0,
	string scalar share0,
	string scalar sjg0,
	string scalar sjh0,
	string scalar alpha0,
	string scalar sigmag0,
	string scalar sigmah0,
	string scalar touse)
{

	// declarations
	real colvector market,brand,g,h,p,q,s,sjg,sjh,alpha,sigmag,sigmah,marketm,pm,qm,sm,brandm,alpham,sigmagm,sigmahm,gm,hm,sjgm,sjhm,sbm,pbm,scg,sjscg,schg,sjschg,missing_index,im
	real matrix amsumb,el,dr,weight,msumb,msumg,msumhg,elm,qbm,weightm,elm0
	real scalar b,bb,m,elmij,j,i,gg,hh,_alpha,_sigmag,_sigmah

	// load data
	st_view(market, .,tokens(market0),touse)
	st_view(brand, .,tokens(brand0),touse)
	st_view(g, .,tokens(g0),touse)
	st_view(h, .,tokens(h0),touse)
	st_view(p, .,tokens(p0),touse)
	st_view(q, .,tokens(q0),touse)
	st_view(s, .,tokens(share0),touse)
	st_view(sjg, .,tokens(sjg0),touse)
	st_view(sjh, .,tokens(sjh0),touse)
	st_view(alpha, .,tokens(alpha0),touse)
	st_view(sigmag, .,tokens(sigmag0),touse)
	st_view(sigmah, .,tokens(sigmah0),touse)
	
	// reindexing group variable (categorical variables should start from 1 and increment by 1)
	market=reindex(market)
	brand=reindex(brand)
	g=reindex(g)
	h=reindex(h)

	// matrix for summations
	amsumb=J(rows(market),colmax(brand)-colmin(brand)+1,0)
	b=0
	for (bb=colmin(brand); bb<=colmax(brand); bb++) {
		b=b+1
		amsumb[.,b]=(brand:==bb)
	}

	el=J(colmax(brand),colmax(brand),0)												// brand level elasticity matrix
	dr=J(colmax(brand),colmax(brand),0)												// brand level diversion ratio matrix
	weight=J(colmax(brand),colmax(brand),0)											// weight matrix
	for (m=colmin(market); m<=colmax(market); m++) {								// calcualtions by markets

		// variables
		marketm=select(market,market:==m)
		pm=select(p,market:==m)
		qm=select(q,market:==m)
		sm=select(s,market:==m)
		brandm=select(brand,market:==m)
		alpham=select(alpha,market:==m)
		sigmagm=select(sigmag,market:==m)
		sigmahm=select(sigmah,market:==m)
		gm=select(g,market:==m)
		hm=select(h,market:==m)
		sjgm=select(sjg,market:==m)
		sjhm=select(sjh,market:==m)
		_alpha=mean(alpham)
		_sigmag=mean(sigmagm)
		_sigmah=mean(sigmahm)
	
		// matrices for summations
		msumb=select(amsumb,market:==m)
		im=runningsum(J(rows(gm),1,1))
		msumg=amsumf(im,gm)
		msumhg=amsumf(im,hm):*msumg

		// calculation of brand level elasticity matrix
		qbm=msumb'*qm																// brand level quantities
		elm=J(colmax(brand),colmax(brand),0)										// brand level elasticity matrix of market m
		if (rows(uniqrows(brandm))==rows(brandm)) {									// if no grouping (no aggregation needed)
			elm0=J(rows(brandm),rows(brandm),1):*sm'
			elm0=elm0 :+ (_sigmag/(1-_sigmag))*(msumg):*(sjgm')
			elm0=elm0 :+ ( (1/(1-_sigmah)) - (1/(1-_sigmag)) )*(msumhg):*(sjhm')
			elm0=elm0 :- diag( J(rows(elm0),1,(1/(1-_sigmah))) )
			elm0=_alpha*elm0
			elm0=elm0:*pm'															// transforming into elasticities: multiplying by derivating prices
			elm=J(rows(el),cols(el),0)
			elm[brandm,brandm]=elm0													// elasticity matrix of market m 
		}
		if (rows(uniqrows(brandm))!=rows(brandm) & rows(uniqrows(brandm))!=1) {		// if grouping (aggregation needed)
			pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')							// brand level weighted average prices
			sbm=msumb'*sm															// brand level market shares
			elm=J(rows(brandm),1,1)*(sbm')											// brand level elasticity matrix of market m
			elm=elm :+ (_sigmag/(1-_sigmag))*(msumg:*(sjgm'))*msumb
			elm=elm :+ ( (1/(1-_sigmah)) - (1/(1-_sigmag)) )*(msumhg:*(sjhm'))*msumb
			elm=elm :- msumb*(1/(1-_sigmah))
			elm=_alpha*elm:*sm
			elm=msumb'*elm
			elm=elm:/sbm															// transforming into semi elasticities: dividing by shares to be derivated
			elm=elm:*pbm'															// transforming into elasticities: multiplying by derivating prices
			// changing missing values to zero (in unbalanced panels elasticities involving product(s) not present at the given market are missing)
			if (rows(elm)>1) {
				missing_index=(elm[.,1]:==.):*runningsum(J(rows(elm),1,1))
				missing_index=select(missing_index,missing_index:>0)
				elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
				elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
				missing_index=(elm[1,.]':==.):*runningsum(J(rows(elm),1,1))
				missing_index=select(missing_index,missing_index:>0)
				elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
				elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
			}
		}
		if (rows(uniqrows(brandm))!=rows(brandm) & rows(uniqrows(brandm))==1) {		// if market level elasticity
			pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')							// market level weighted average price
			elm=-_alpha*(1-sum(sm))*pbm
		}

		// weigth matrix for market m
		weightm=qbm#J(1,rows(qbm),1)

		// adding up weights and weighted elasticities
		weight=weight:+weightm
		el=el:+(elm:*weightm)
	
	}	// end of m (markets) loop 

	el=el:/weight																	// normalizing: final, weighted average (across markets) elasticity matrix

	// calculation of loss diversations for average elasticites
	for (j=1; j<=colmax(brand); j++) {
		if (sum(brand:==j)!=0) {		// only if brand j is present
			for (i=1; i<=colmax(brand); i++) {
				if (sum(brand:==i)!=0) {		// only if brand i is present
					dr[i,j]=100*sum(q:*(brand:==i))*el[i,j]/abs(sum(q:*(brand:==j))*el[j,j])
				}
			}
		}
	}

	// storing in Stata
	st_matrix("el", el)
	st_matrix("dr", dr)
	
}	// end of elasticities_nlogit2 function
mata mlib add lrcl elasticities_nlogit2()


// three-level nested logit elasticities and diversion ratios
void elasticities_nlogit3(
	string scalar market0,
	string scalar brand0,
	string scalar g0,
	string scalar h0,
	string scalar k0,
	string scalar p0,
	string scalar q0,
	string scalar share0,
	string scalar sjg0,
	string scalar sjh0,
	string scalar sjk0,
	string scalar alpha0,
	string scalar sigmag0,
	string scalar sigmah0,
	string scalar sigmak0,
	string scalar touse)
{

	// declarations
	real colvector market,brand,g,h,k,p,q,s,sjg,sjh,sjk,alpha,sigmag,sigmah,sigmak,marketm,pm,qm,sm,brandm,alpham,sigmagm,sigmahm,sigmakm,gm,hm,km,sjgm,sjhm,sjkm,sbm,pbm,scg,sjscg,schg,sjschg,sckhg,sjsckhg,missing_index,im
	real matrix amsumb,el,dr,weight,msumb,msumg,msumhg,msumkhg,elm,qbm,weightm,elm0
	real scalar b,bb,m,elmij,j,i,gg,hh,kk,_alpha,_sigmag,_sigmah,_sigmak

	// load data
	st_view(market, .,tokens(market0),touse)
	st_view(brand, .,tokens(brand0),touse)
	st_view(g, .,tokens(g0),touse)
	st_view(h, .,tokens(h0),touse)
	st_view(k, .,tokens(k0),touse)
	st_view(p, .,tokens(p0),touse)
	st_view(q, .,tokens(q0),touse)
	st_view(s, .,tokens(share0),touse)
	st_view(sjg, .,tokens(sjg0),touse)
	st_view(sjh, .,tokens(sjh0),touse)
	st_view(sjk, .,tokens(sjk0),touse)
	st_view(alpha, .,tokens(alpha0),touse)
	st_view(sigmag, .,tokens(sigmag0),touse)
	st_view(sigmah, .,tokens(sigmah0),touse)
	st_view(sigmak, .,tokens(sigmak0),touse)

	// reindexing group variable (categorical variables should start from 1 and increment by 1)
	market=reindex(market)
	brand=reindex(brand)
	g=reindex(g)
	h=reindex(h)
	k=reindex(k)

	// matrix for summations
	amsumb=J(rows(market),colmax(brand)-colmin(brand)+1,0)
	b=0
	for (bb=colmin(brand); bb<=colmax(brand); bb++) {
		b=b+1
		amsumb[.,b]=(brand:==bb)
	}

	el=J(colmax(brand),colmax(brand),0)												// brand level elasticity matrix
	dr=J(colmax(brand),colmax(brand),0)												// brand level diversion ratio matrix
	weight=J(colmax(brand),colmax(brand),0)											// weight matrix

	for (m=colmin(market); m<=colmax(market); m++) {								// calcualtions by markets

		// variables
		marketm=select(market,market:==m)
		pm=select(p,market:==m)
		qm=select(q,market:==m)
		sm=select(s,market:==m)
		brandm=select(brand,market:==m)
		alpham=select(alpha,market:==m)
		sigmagm=select(sigmag,market:==m)
		sigmahm=select(sigmah,market:==m)
		sigmakm=select(sigmak,market:==m)
		gm=select(g,market:==m)
		hm=select(h,market:==m)
		km=select(k,market:==m)
		sjgm=select(sjg,market:==m)
		sjhm=select(sjh,market:==m)
		sjkm=select(sjk,market:==m)
		_alpha=mean(alpham)
		_sigmag=mean(sigmagm)
		_sigmah=mean(sigmahm)
		_sigmak=mean(sigmakm)

		// matrices for summations
		msumb=select(amsumb,market:==m)
		im=runningsum(J(rows(gm),1,1))
		msumg=amsumf(im,gm)
		msumhg=amsumf(im,hm):*msumg
		msumkhg=amsumf(im,km):*msumhg

		// calculation of brand level elasticity matrix
		qbm=msumb'*qm																// brand level quantities
		if (rows(uniqrows(brandm))==rows(brandm)) {									// if no grouping (no aggregation needed)
			elm0=J(rows(brandm),rows(brandm),1):*sm'
			elm0=elm0 :+ (_sigmag/(1-_sigmag))*(msumg):*(sjgm')
			elm0=elm0 :+ ( (1/(1-_sigmah)) - (1/(1-_sigmag)) )*(msumhg):*(sjhm')
			elm0=elm0 :+ ( (1/(1-_sigmak)) - (1/(1-_sigmah)) )*(msumkhg):*(sjkm')
			elm0=elm0 :- diag( J(rows(elm0),1,(1/(1-_sigmak))) )
			elm0=_alpha*elm0
			elm0=elm0:*pm'															// transforming into elasticities: multiplying by derivating prices
			elm=J(rows(el),cols(el),0)
			elm[brandm,brandm]=elm0													// elasticity matrix of market m 
		}
		if (rows(uniqrows(brandm))!=rows(brandm) & rows(uniqrows(brandm))!=1) {		// if grouping (aggregation needed)
			pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')							// brand level weighted average prices
			sbm=msumb'*sm															// brand level market shares
			elm=J(rows(brandm),1,1)*(sbm')											// brand level elasticity matrix of market m
			elm=elm :+ (_sigmag/(1-_sigmag))*(msumg:*(sjgm'))*msumb
			elm=elm :+ ( (1/(1-_sigmah)) - (1/(1-_sigmag)) )*(msumhg:*(sjhm'))*msumb
			elm=elm :+ ( (1/(1-_sigmak)) - (1/(1-_sigmah)) )*(msumkhg:*(sjkm'))*msumb
			elm=elm :- msumb*(1/(1-_sigmak))
			elm=_alpha*elm:*sm
			elm=msumb'*elm
			elm=elm:/sbm															// transforming into semi elasticities: dividing by shares to be derivated
			elm=elm:*pbm'															// transforming into elasticities: multiplying by derivating prices
			// changing missing values to zero (in unbalanced panels elasticities involving product(s) not present at the given market are missing)
			if (rows(elm)>1) {
				missing_index=(elm[.,1]:==.):*runningsum(J(rows(elm),1,1))
				missing_index=select(missing_index,missing_index:>0)
				elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
				elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
				missing_index=(elm[1,.]':==.):*runningsum(J(rows(elm),1,1))
				missing_index=select(missing_index,missing_index:>0)
				elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
				elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
			}
		}
		if (rows(uniqrows(brandm))!=rows(brandm) & rows(uniqrows(brandm))==1) {		// if market level elasticity
			pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')							// market level weighted average price
			elm=-_alpha*(1-sum(sm))*pbm
		}

		// weigth matrix for market m
		weightm=qbm#J(1,rows(qbm),1)

		// adding up weights and weighted elasticities
		weight=weight:+weightm
		el=el:+(elm:*weightm)
	
	}	// end of m (markets) loop 

	el=el:/weight																	// normalizing: final, weighted average (across markets) elasticity matrix

	// calculation of loss diversations for average elasticites
	for (j=1; j<=colmax(brand); j++) {
		if (sum(brand:==j)!=0) {		// only if brand j is present
			for (i=1; i<=colmax(brand); i++) {
				if (sum(brand:==i)!=0) {		// only if brand i is present
					dr[i,j]=100*sum(q:*(brand:==i))*el[i,j]/abs(sum(q:*(brand:==j))*el[j,j])
				}
			}
		}
	}
	
	// storing in Stata
	st_matrix("el", el)
	st_matrix("dr", dr)
	
}	// end of elasticities_nlogit3 function
mata mlib add lrcl elasticities_nlogit3()


// BLP elasticity and diversion ratio matrices (averages over markets)
// elvar: variable specifying the groups to which the aggregation is to be done (product, brand, firm, etc.)
void elasticities_rcl(
	real colvector theta1,
	real rowvector theta2,
	real colvector delta,
	real colvector elvar)
{

	// declarations
	external real colvector market,p,s,msize
	external real matrix rc,simdraws,msumm,prices,rc_prices
	external real scalar _is_rc_on_p
	real colvector brand,q,sigmas,pm,qm,sm,sbm,qbm,pbm,missing_index
	real matrix amsumb,alphai,mu,shati,el,dr,weight,shatim,msumb,shatimb,elm,weightm
	real scalar b,bb,alpha,m,j,i
	
	// reindexed group variable (henceforth, it is called "brand")
	brand=reindex(elvar)

	// quantities
	q=s:*msize

	// matrix for summations
	amsumb=J(rows(market),colmax(brand)-colmin(brand)+1,0)
	b=0
	for (bb=colmin(brand); bb<=colmax(brand); bb++) {
		b=b+1
		amsumb[.,b]=(brand:==bb)
	}

	// parameters
	sigmas=theta2'																	// vector of random coefficient parameters
	alpha=-theta1[1,1]																// mean price coefficient
	if (_is_rc_on_p==1) {
		alphai=alpha:-(sigmas[1,1]*simdraws[1,.])									// individual price coefficients (if there is random coefficient on price)
	}
	if (_is_rc_on_p==0) {
		alphai=alpha:*J(1,cols(simdraws),1)											// individual price coefficients (if there is no random coefficient on price)
	}
	mu=rc*(sigmas:*simdraws[1..rows(simdraws)-1,.])									// matrix of observed consumer heterogeneity (separate column for each consumer)

	// predicted market shares
	shati=exp(delta:+mu)															// predicted individual choice probabilities
	shati=shati:/(1:+msumm*((shati'*msumm)'))

	el=J(colmax(brand),colmax(brand),0)												// brand level elasticity matrix
	dr=J(colmax(brand),colmax(brand),0)												// brand level diversion ratio matrix
	weight=J(colmax(brand),colmax(brand),0)											// weight matrix

	for (m=colmin(market); m<=colmax(market); m++) {								// calcualtions by markets

		// variables
		pm=select(p,market:==m)
		qm=select(q,market:==m)
		sm=select(s,market:==m)
		brandm=select(brand,market:==m)
		shatim=select(shati,market:==m)

		// brand matrix for summations
		msumb=select(amsumb,market:==m)

		/* calculation of brand level elasticity matrix */

		// auxiliary objects
		shatimb=msumb'*shatim														// individual choice probabilities of brands (sums of product level probabilities for a given consumer)
		sbm=msumb'*sm																// brand level market shares
		qbm=msumb'*qm																// brand level quantities
		pbm=((msumb:*qm)'*pm):/(colsum((msumb:*qm))')								// brand level weighted average prices

		// brand level elasticity matrix
		// rows: brands whose shares respond (derivated shares)
		// columns: brands whose prices increase (derivating prices)
		elm=((alphai:*shatimb):*(simdraws[rows(simdraws),.]))*shatimb'				// general (to all elasticities) cross terms
		elm=elm:-diag((alphai:*shatimb)*simdraws[rows(simdraws),.]')				// terms specific to own-price elasticites
		elm=elm:/sbm																// transforming into semi elasticities: dividing by shares to be derivated
		elm=elm:*pbm'																// transforming into elasticities: multiplying by derivating prices

		// changing missing values to zero (in unbalanced panels elasticities involving product(s) not present at the given market are missing)
		if (rows(elm)>1) {
			missing_index=(elm[.,1]:==.):*runningsum(J(rows(elm),1,1))
			missing_index=select(missing_index,missing_index:>0)
			elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
			elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
			missing_index=(elm[1,.]':==.):*runningsum(J(rows(elm),1,1))
			missing_index=select(missing_index,missing_index:>0)
			elm[missing_index,.]=J(rows(elm[missing_index,.]),cols(elm[missing_index,.]),0)
			elm[.,missing_index]=J(rows(elm[.,missing_index]),cols(elm[.,missing_index]),0)
		}

		// weigth matrix for market m
		weightm=qbm#J(1,rows(sbm),1)

		// adding up weights and weighted elasticities
		weight=weight:+weightm
		el=el:+(elm:*weightm)

	}	// end of m (markets) loop 

	el=el:/weight																	// normalizing: final, weighted average (across markets) elasticity matrix

	// calculation of loss diversations for average elasticites
	for (j=1; j<=colmax(brand); j++) {
		if (sum(brand:==j)!=0) {		// only if brand j is present
			for (i=1; i<=colmax(brand); i++) {
				if (sum(brand:==i)!=0) {		// only if brand i is present
					dr[i,j]=100*sum(q:*(brand:==i))*el[i,j]/abs(sum(q:*(brand:==j))*el[j,j])
				}
			}
		}
	}

	
	// storing in Stata
	st_matrix("el", el)
	st_matrix("dr", dr)
	
}	// end of elasticities_rcl function
mata mlib add lrcl elasticities_rcl()

end