st: solving system of equations with Mata

[Nicola Orsini <nicola.orsini@ki.se>]

Dear Mata users,

I programmed the below iterative method in Mata
to solve a system of equations with one constraint
but I'm wondering if there is a more efficient way to do that.

Any thought is welcome!

Best,
Nicola

<--- Do-file start here --->

// Worked example

clear
version 9.1
mata

// suppose we have 3 column vectors (4 by 1)

A = (165\ 74\ 90\ 122)
N = (337\167\186\212)
Y = (1\0.8\1.16\1.57)

(A, N, Y)

// Goal: Solve a system of 4 equations for A

// log(Y_i)+log(M1- AP)+log(N_i - A_i)-log(A_i)-log(N0-M1+AP)=0

// where the index i denotes each row with i = 1,2,3,4
// and
// M1 = colsum(A)
// AP = colsum(A[|2,1 \ ., .|])
// N0 = N[1]
//
// with the constraint that M1 has to be kept constant

// what I did is the following

    M1 = colsum(A)
    NS = colsum(N)
    N0 = N[1]

    // Loop until converence (no more changes of A)

    SINCR = 1

while (abs(SINCR)>1e-5) {

    AP = colsum(A[|2,1 \ ., .|])
    A0 = M1 - AP
     ex = J(rows(A), 1,.)
    cx = (1:/A) + 1:/(N:-A)

for (i=2; i<=rows(A); i++) {
     ex[i] = log(Y[i]) + log(A0)+log(N[i]-A[i])-log(A[i])-log(N0-A0)
}

      H1 = J(rows(A)-1,rows(A)-1,cx[1])
      H2 = diag(cx)
      H = H1 + H2[|2,2\. , .|]

      AP = A[|2,1 \ ., .|] + invsym(H)*ex[|2,1 \ ., .|]

    INCR = invsym(H)*ex[|2,1 \ ., .|]
    SINCR = colsum(INCR)
     DCA = M1-sum(AP)

    // New colvector of fitted A
    A  = (DCA \ AP)
}

// The solution is the vector A

A

// To check if the vector A is the correct solution

 AP = colsum(A[|2,1 \ ., .|])
 X = log(Y) :+ log(M1 - AP) :+ log(N-A) :- log(A) :- log(N0 - M1 + AP)
 X
end

<--- Do-file end here --->

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/