Re: st: looping in MATA

[wgould@stata.com (William Gould, StataCorp LP)]

> I need to create a rolling loop in stata mata language, I am using
> some simple matrix algebra and it is easier to do in Mata than in
> matrix functions.
> 
> I need to do something like this (within a do file):
> *****************************************
> mata
> 
> foreach y in 1970/2000 {
> 
> x`y' = st_data(., (" g`y'", " g`y+1'", " g`y+2'"," g`y+3'", " g`y+4'"), 0)
> 
> r=rows(x`y')
> vc=x`y'*x`y''
> vc=vc/25
> var`y'=sqrt(trace(vc)/rows(vc))
> tr=vc-diag(vc)
> cov`y'=sum(tr)
> cov`y'=sqrt(cov`y'/(r^2-r)/2)
>
> .end
> *************************
>
> Ideally I would also like to save these cov`y' and var`y' data into a matrix.

Here is a solution:

        mata:

        result = J(0, 3, .)
        for (y=1970; y<=2000; y++) {
                v0 = sprintf("g%f", y)
                v1 = sprintf("g%f", y+1)
                v2 = sprintf("g%f", y+2)
                v3 = sprintf("g%f", y+3)
                v4 = sprintf("g%f", y+4)

                x  = st_data(., (v0, v1, v2, v3, v4))
                r  = rows(x)
                vc = (x*x')/25        // I think 25 should be r
                var= sqrt(trace(vc)/25)
                tr = vc - diag(vc)
                cov= sum(tr)
                cov= sqrt(cov/(r^2-r)/2)

                result = result \ (y,var,cov)
        }
        result

        end

N.B., I did not use macros!  As Nick Cox <n.j.cox@durham.ac.uk> Cox said in
his reply, "The key principle here is that local macro references are not part
of Mata. So, they must be interpreted before Mata sees them."  Nick went from
there to suggesting a a "clumsy" (his word) solution that called Mata line by
line.

I eliminated the macros altogether.  Whenever you're thinking macros and 
Mata simultaneously, you are thinking wrong, and the solution is not to 
try to figure out how to make the macros work, but instead to think about 
how to eliminate the macros altogether.

That was easy.

Olga thought, "I need to use use variables named g`y', g`=y+1', ..., where 
y = 1970, 1971, ...  That's the right way to think in Stata, but not 
in Mata.  Whereever you use macros in Stata, you need yet another variable
in Mata.  Variables in Mata are used for everything.  I created a sequence
of Mata string variables to hold the names I wanted to access:

                v0 = sprintf("g%f", y)
                v1 = sprintf("g%f", y+1)
                v2 = sprintf("g%f", y+2)
                v3 = sprintf("g%f", y+3)
                v4 = sprintf("g%f", y+4)

I used sprintf() to assemble the name, but I could just as well have coded

                v0 = "g" + strofreal(y) 
                v1 = "g" + strofreal(y+1) 
                v2 = "g" + strofreal(y+2) 
                v3 = "g" + strofreal(y+3) 
                v4 = "g" + strofreal(y+4) 

Code it however you wish, but understand, variable y is a real scalar, 
and I'm forming variable names and stroging them in variables v0, v1, ...,
v4, which are string scalars.  Once I had the variables containing the 
names of the Stata variables, I was ready to form matrix x:

                x  = st_data(., (v0, v1, v2, v3, v4))

By the way, I collected the results together in a single matrix, just 
as Olga requested.

By the way, I think Olga needs to check that she has the correct formulas 
for the calculation.  There is a 25 in the code that I believe should be 
r, and I note that the calculation as written is scale dependent.  To 
wit, I tested the above code on some literally (psuedo) random data.
Then I took the data and added 10 to to every variable.  Answers changed.

The scale dependencies can be traced to the line, 

                vc = (x*x')/25        // I think 25 should be r

I suspect Olga may be thinking of matrix x as a matrix of values recorded
in deviation from the mean form.  If so, changing absolute x into deviation 
form is easy enough: 

               x = x :- (colsum(x):/rows(x)) 

or, if you prefer:

               mean = colsum(x) :/ rows(x)
               x    = x :- mean

-- Bill
wgould@stata.com
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/