*! version 1.0.1  15dec2005

/* Forms the derivative matrix needed to apply the delta method to obtain
   the variance of the full vector of parameters from the estimated vector
   of parameters.
*/

program quaids_delta2, rclass

	syntax , EQuations(integer)
	
	if `equations' < 2 {
		di as error "equations() must be at least two"
		exit 198
	}
	
	tempname Delta
	mata: _quaids_delta2(`equations', "`Delta'")
	return matrix Delta = `Delta'
		
end

mata:
mata set matastrict on
void _quaids_delta2(real scalar ng, string scalar Dmat)
{

	real scalar i, j, ic, jc, m, n		/* counters */
	real scalar nparest, npartot, ngm1
	real vector mi
	real matrix Delta

	nparest = 3*(ng-1) + 0.5*ng*(ng-1)
	npartot = 3*ng + 0.5*(ng+1)*ng
	
	Delta = J(npartot, nparest, 0)
	
	ngm1 = ng - 1
	mi = J(1, ngm1, -1)
	mi = I(ngm1) \ mi
	Delta[|1,1 \ ng, ngm1|] = mi				/* Alpha */
	Delta[|ng+1, ng \ 2*ng, 2*(ng-1)|] = mi			/* Beta  */
	Delta[|npartot-ng+1, nparest-ngm1+1 \ ., .|] = mi	/* Lambda */
	
	m = 2*ng+1
	for(j=1; j<=ng; ++j) {
		for(i=j; i<=ng; ++i) {
			n = 2*ng-1
			for(jc=1; jc<ng; ++jc) {
				for(ic=jc; ic<ng; ++ic) {
					if (j < ng & i < ng) {     /* Case 1 */
						if (jc==j && ic==i) 
							Delta[m,n] = 1
					}
					else if (j < ng & i==ng) { /* Case 2 */
						if (jc==j || ic==j)
							Delta[m,n]=-1
					}
					else if (j==ng & i==ng) {  /* Case 3 */
						Delta[m,n] = Delta[m,n]+1
						if (ic != jc) 
							Delta[m,n] =
								Delta[m,n] + 1
					}
					++n
				}
			}
			++m
		}
	}

	st_matrix(Dmat, Delta)
	
}
	
end


/*
	Derivatives for the Gamma parameters.

	Let N = number of goods

	Case 1: We have Gamma_{i,j} for i,j < N.  The derivative is
		simply unity
		
	Case 2: We have Gamma_{i,j} for i==N and j < N.  In this case
		
			Gamma_{N,j} = 1 - Gamma_{1,j} - Gamma_{2,j} - ...
		
		So the required derivative is minus one.
		
	Case 3: We are at Gamma_{N,N}.
	
			Gamma_{N,N} = 1 - Gamma_{1,N} - Gamma_{2,N} - ...
		(Slutsky symm.)	    = 1 - Gamma_{N,1} - Gamma_{N,2} - ...
				    = 1 - (1 - Sum_{j=1}^{j=N-1} Gamma_{1,j})
				    	- (1 - Sum_{j=1}^{j=N-1} Gamma_{2,j})
				    	- ...                           ^
				    	                    Call that i |
				    	
		So the required derivative is one if i==j and two if i!=j.
		
*/