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/