Friday, February 03, 2006

st: solving system of equations with Mata

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/


Tag:




<< Home

This page is powered by Blogger. Isn't yours?