rolling window multiplication matrix

[Sofia Frascalo]

Hi everyone,



I am trying to use a cumulative multiplication formula over a 12 month window for 47 industries.
I am trying to create a loop to do this all at once in a rolling window.
I tried to look up a stata formula for cumulative multiplication but I could not find anything.
so
N12=Nk,t:t-11= [(1+Nk,t)*(1+Nk,t-1).........(1+Nk,t-11)]-1
where k indicates industry k and t indicates the specific month.

This is what I have:

foreach k of varlist N_Agric N_Food N_Beer N_Smoke N_Toys N_Fun N_Books N_Hshld N_Clths N_MedEq N_Drugs N_Chems N_Txtls N_BldMt N_Cnstr N_Steel N_Mach N_ElcEq N_Autos N_Aero N_Ships N_Mines N_Coal N_Oil N_Util N_Telcm N_PerSv N_BusSv N_Hardw N_Chips N_LabEq N_Boxes N_Trans N_Whlsl N_Rtail N_Meals N_Banks N_Insur N_RlEst N_Fin N_Other{
gen N12 = (`k'[_n] * `k'[_n-1]* `k'[_n-2]* `k'[_n-3]* `k'[_n-4]* `k'[_n-5]* `k'[_n-6]* `k'[_n-7]* `k'[_n-8]* `k'[_n-9]* `k'[_n-10]* `k'[_n-11])-1
}

does anyone know how to do this?

Thank you in advance!

[Clyde Schechter]

The problem with your code is that it tries to generate the same variable N12 on each iteration of the loop. That will work the first time, but then the second time Stata will complain that N12 already exists. So you have to create a different variable, whose name refers to `k', sometime in the loop. So changing -gen N12 = - to -gen N12`k' =- is the first step. The second issue is that if you are looking at, say the second observation for the second industry, `k'[_n-3] is going to refer to an observation in the preceding industry. Similarly, if there are some gaps in the dates, you will be extending your rolling product too far back in time. So I would use a safer approach:

xtset industry month
foreach k of varlist N_Agric N_Food /*etc.*/ {
    by industry (month), sort: gen double running_product_`k' = `k' if _n == 1
    by industry (month): replace running_product_`k' = running_product_`k'[_n-1] * `k' if _n > 1
    by industry (month): gen double rolling_product_`k' = running_product_`k'/L11.running_product_`k'
    gen N12`k' = rolling_product_`k' - 1
}

[Attaullah Shah]

The product option of asrol (can be downloaded from SSC) implements a log transformation method of finding products of values. Here is one example how the product and log transformation work.


Product method

 
Values
cumulative product

0.5
0.5

0.6
0.3

0.3
0.09

The log transformation method

       values          log(values)
---------------------------------- 
0.5
-0.693147181

0.6
-0.510825624

0.3
-1.203972804

sum
-2.407945609

=exp(-2.407)
0.09

asrol example

clear
input float values
.5
.6
.3
end
ssc install asrol
asrol values, stat(product)


     +--------------------+
     | values   product~s |
     |--------------------|
  1. |     .5   .09000001 |
  2. |     .6   .09000001 |
  3. |     .3   .09000001 |
     +--------------------+

You can also consider option add() of asrol.

Therefore, in your example, the asrol code could be:

bys industries: asrol N_Agric N_Food N_Beer N_Smoke N_Toys N_Fun N_Books N_Hshld N_Clths N_MedEq N_Drugs N_Chems N_Txtls N_BldMt N_Cnstr N_Steel N_Mach N_ElcEq N_Autos N_Aero N_Ships N_Mines N_Coal N_Oil N_Util N_Telcm N_PerSv N_BusSv N_Hardw N_Chips N_LabEq N_Boxes N_Trans N_Whlsl N_Rtail N_Meals N_Banks N_Insur N_RlEst N_Fin N_Other, window(month 16) stat(product) 

Note: You are using the _n operator, which disregards any gap in the time series. However, asrol is aware of any time gap and accordingly adjusts the rolling window calculations.

[Nick Cox]

Just in case you need it, here is a reminder that converting to logarithms requires that all terms in the product are positive. A product (a) is necessarily 0 if any term is 0 and (b) otherwise could be positive or negative if there is a mix of positive and negative terms, but using logarithms would lead to a missing result in either case. That said, almost all posts asking for this ask when terms are all positive. But that's not explicit in #1.

I don't know if asrol checks for (a) or (b) and gives a warning in either case. That would be a good feature if so, and otherwise is a suggestion for a future revision.

Note also this solution from 1999:

STB-51 dm71 . . . . . . . . . . . . Calculating the product of observations
(help prod if installed) . . . . . . . . . . . . . . . . . . P. Ryan
9/99 pp.3--4; STB Reprints Vol 9, pp.45--48
extension to egen for producing the product of observations

which also uses the logarithm method

*! version 1.0.0                    [STB51: dm71]
*! Calculate the product of observations
*! extension to egen
*! Philip Ryan  August 6 1999  v1.3.4
*! with coding suggestions by Michael Blasnik 
program define _gprod
version 6
#delimit ;
syntax newvarname = /exp [if] [in] [, BY(varlist) PMiss(string)];
if "`pmiss'" == "" {;
local pmiss "ignore";
};
if "`pmiss'" != "missing" & "`pmiss'" != "1" & "`pmiss'" != "ignore" {;
displ in re "Syntax error: specify " in wh "pmiss(missing)" in re " or "
in wh "pmiss(1)"  in re " or " in wh "pmiss(ignore)";
exit 198;
};
marksample touse, novarlist;
tempvar absvar oldsrt sumlog negsign;
quietly{;
local sortby: sortedby;
gen long `oldsrt' = _n;
gen double `absvar' = abs(`exp');
sort `touse' `by' `absvar';
gen double `sumlog' = log(`absvar') if `touse';
by `touse' `by': replace `sumlog'=sum(`sumlog') if `touse';
replace `sumlog'=exp(`sumlog');
by `touse' `by': gen `negsign' = sum(`exp'<0) if `touse';
by `touse' `by': replace `sumlog'=`sumlog'[_N]*
                 cond(mod(`negsign'[_N],2)==0,1,-1)*
                 (`absvar'[1]>0 );
if "`pmiss'" == "ignore" {;
by `touse' `by': replace `sumlog'= . if `absvar'[1]==.;
};
if "`pmiss'" == "missing" {;
by `touse' `by': replace `sumlog'= . if `absvar'[_N]==.;
};

rename `sumlog' `varlist' ;
sort `sortby' `oldsrt';
};
#delimit cr
end
exit

[Attaullah Shah]

Thanks, Nick for the suggestion. In fact, to get around the problem of zero and negative values, asrol offers the option add(). The help file of asrol has details on option add(#). These details are lso copied below.

add(#) adds the value # to each value in a given window before computing the geometric mean or products of values. This option is useful when analyzing variables with non-positive values. Since log of negative numbers is undefined, one workaround is to add a positive number that is larger enough to convert non-positive values into positive values before log transformation. Option add(#) can be used for this purpose. asrol adds the positive number supplied through add(#) option to each of the value in a given window, finds the required statistic, and then subtracts back the positive number.

Please note that the results (specifically the geometric mean) are usually sensitive to the selection of the value that is added before calculation.

[Nick Cox]

I wasn't suggesting a way "to get around" zero and negative values. The mathematical facts from early education are that

1. even one zero in a set of values means that their product is zero, or in colourful terminology sometimes used, it annihilates the product.

2. the sign of a product is negative if and only if the number of negative signs in the set is odd.

These cases (any zeros or negative values) are fatal if a product is calculated through a sum of logarithms as the logarithm of any zero or negative value is returned as missing by Stata and the product is in turn missing.

I was suggesting that the problem -- for your program asrol and any other that uses that method -- merits a warning to users, if not bailing out with an error.

I know that there are other good reasons e.g. to add 1 as expressions like log(1 + x) often arise naturally, but I don't see that that can be a general solution. (Adding another constant sounds a lot weirder.)

Suppose I have 0 1 2 and add 1 to get 1 2 3 and a product of 8. 1 1 1 would give the same product of 8 if bumped up. There is no way to recover the fact for the first case that one term was 0 from the product, nor a way of knowing that the correct product should be zero. Similarly suppose that I have -3 -2 -1 and add 4 to get 1 2 3 and a product of 6 where is the information that enables me to reconstruct that the product should be negative?

This is more than the fact that you can often get the same product in different ways. It's that a fudge can't give you the correct answer at all.

I don't ever use a geometric mean unless all values are positive. The geometric mean of several values + constant is just some other beast. It could be a crude approximation to the geometric mean in some cases, but otherwise what is its use or interest?

I'll flag that rangestat (SSC) doesn't support products. If users of rangestat want to do this, they can take logs first and then exponentiate afterwards and they won't ever get results that are not justifiable. Or they can turn to rangerun (SSC) and code up a solution using Clyde's method, which always gives the right answer. Or, naturally, use code like Clyde's.

[Attaullah Shah]

Nick, As you mentioned, calculation of products is always tricky due to the reasons you mentioned in post #4. User's should be alert to these pitfalls. I did flag

Please note that the results (specifically the geometric mean) are usually sensitive to the selection of the value that is added before calculation.

.
On a side note, I was expecting that you would comment on the use of _n in rolling window estimation.

[Nick Cox]

I am not following. Either

(A) asrol has a work-around for negative or zero terms that makes sense (I don't understand what it could be)

or

(B) it can't solve the problem -- and users of your command deserve a signal to that effect

or

(C) it is useful to add a constant first and calculate something else (but why?) .

Let's make this concrete: What is the asrol way to calculate the products 0 1 2 or -1 -2 -3 (correct answers 0 and -6 respectively)?

Which is it? Or what I am misunderstanding?

Are you suggesting that given terms like a b c d then the product

a b c d ... = exp[log(a + k) + log(b + k) + log(c + k) + log(d + k) + ...] - k

because that is generally valid only if k = 0 and a b c d ... are all positive, in which case it is pointless.

P.S.

I was expecting that you would comment on the use of _n in rolling window estimation.

I don't understand that. Observation numbers may or may not be a metric to count in.

[Attaullah Shah]

Well, adding a warning message is already on my agenda in my next version of asrol. Thanks for this suggestion. If one searches this forum, s/he can easily find that log transformation for finding products is a matter of discussion for long. While posting my asrol solution, I showed an example, clearly highlighted that asrol uses this method, and then it is up to the user to use this log method or not. I did not develop the log method, it is there since long and has its merits and demerits. I cannot do much about it. In asrol, it is just a matter of option.

I have great respect for Clyde but omissions do happen. Why is it so that the Clyde code does not adjust for zeros, and you ignored that. you said.

Or they can turn to rangerun (SSC) and code up a solution using Clyde's method, which always gives the right answer. Or, naturally, use code like Clyde's.

There is hardly any post related to asrol that you do not pick on me or my programs. This was the attitude I left this platform a while ago but then thought why not try again.

[Nick Cox]

Thanks for your reply. As far as I can see most of my questions from #8 didn't get answered.

Everyone is allowed to bring in other issues, but here it just looks like a diversionary tactic. Let's divide the issues:

1. What is exactly right for Sofia's data can only be established with access to Sofia's data.

2. The essence of Clyde's method is

a. initialise product = first value
b. for each later value : updated product = previous product * later value.

I am emphasising that this method gives the correct answer for zero or negative terms. The idea that Clyde needs to adjust for zeros (???) appears to mean that Clyde is at fault for not changing the definition of a product. I don't recommend that he does that.

You're right that if Sofia wants a window that is in terms of a time variable but has gaps in her data then observation numbers are not the right metric, but you have already pointed this out clearly in #3. I didn't see reason to comment on it yet again. Sorry, but I can't discuss all issues at once. But what you quote from me in #9 does carry the rider -- provided that observation numbers are not wrong for the data being used.

On the main points in recent posts: I regret that you are so sensitive to direct questions about what your code does. I am interested only in getting answers to technical questions. As you've cited yourself in various places I've been writing commands in this territory for some long time and it interests me.

If asrol is wrong you should want to know that and to fix it. If asrol is right and I am wrong you can explain why and many people would find that entertaining or instructive. You owe that to people in the community. Whoever is asking about this is just whoever it is. It's always true that people can choose whether or not to use your program, or any other, but whether the program gives correct and defensible results is germane to that choice.

Of course using logarithms as a way to get products of positive terms is a long-established method. I already pointed out that Philip Ryan used this method in Stata code in 1999 and I know him well enough to imagine that he would laugh out loud at the idea that he invented it.

But that doesn't mean that it extends to getting the right answer when zeros or negative values are present.

[Attaullah Shah]

Here is the problem I was referring to

a. initialise product = first value
b. for each later value : updated product = previous product * later value.

I am emphasising that this method gives the correct answer for zero or negative terms. The idea that Clyde needs to adjust for zeros (???) appears to mean that Clyde is at fault for not changing the definition of a product. I don't recommend that he does that.

clear
set obs 20
gen industry = "A"
gen month = _n+570
format month %tm
gen values = uniform()
replace value = 0 in 1
tsset
by industry (month), sort: gen double running_product = values if _n == 1
by industry (month): replace running_product = running_product[_n-1] * values if _n > 1
tsset
by industry (month): gen double rolling_product = running_product/L11.running_product

* fails
not sorted
r(5);

* since I have one industry, therefore removing by
tsset

gen double rolling_product = running_product/L11.running_product


list

     +-----------------------------------------------------+
     | industry     month     values   runnin~t   rollin~t |
     |-----------------------------------------------------|
  1. |        A    2007m8          0          0          . |
  2. |        A    2007m9   .8077623          0          . |
  3. |        A   2007m10   .6967398          0          . |
  4. |        A   2007m11    .124575          0          . |
  5. |        A   2007m12   .0093527          0          . |
     |-----------------------------------------------------|
  6. |        A    2008m1   .7423609          0          . |
  7. |        A    2008m2   .6246868          0          . |
  8. |        A    2008m3   .3707404          0          . |
  9. |        A    2008m4   .9758291          0          . |
 10. |        A    2008m5   .6274327          0          . |
     |-----------------------------------------------------|
 11. |        A    2008m6     .68028          0          . |
 12. |        A    2008m7   .9618542          0          . |
 13. |        A    2008m8   .3463633          0          . |
 14. |        A    2008m9   .2641014          0          . |
 15. |        A   2008m10   .7256252          0          . |
     |-----------------------------------------------------|
 16. |        A   2008m11   .3859604          0          . |
 17. |        A   2008m12   .8343201          0          . |
 18. |        A    2009m1   .9593934          0          . |
 19. |        A    2009m2   .7894724          0          . |
 20. |        A    2009m3   .9936081          0          . |
     +-----------------------------------------------------+

.

Anyway, thanks again. I shall add another method of finding products as well in addition to log transformation. Hopefully, I shall be able to do that in the summers.

[Nick Cox]

That's concrete but the problem is not in any code. The code is right: the product is zero and it's not useful to divide by it. What that means substantively and whether people would be better off calculating something else are different questions.

[Attaullah Shah]

As said before, I have great respect for Clyde and he happens to be in this post, I do not want to discuss his code, but I am referring to it just because I have to due to the circumstances created. Even if you think that adjusting for zeros is not warranted, what about missing values. My point is that you opt to ignore all these subtilities, and you love to just pick on me.

see this example

clear
set obs 20
gen industry = "A"
gen months = _n+570
format months %tm
gen values = uniform()
replace value = . in 3
encode industry, gen(id)
tsset id month
by industry (month), sort: gen double running_product = values if _n == 1
by industry (month): replace running_product = running_product[_n-1] * values if _n > 1
tsset

gen double rolling_product = running_product/L11.running_product



 +-----------------------------------------------------------+
     | industry    months     values   id   running~t   rollin~t |
     |-----------------------------------------------------------|
  1. |        A    2007m8   .2014603    A   .20146026          . |
  2. |        A    2007m9   .9207765    A   .18549987          . |
  3. |        A   2007m10          .    A           .          . |
  4. |        A   2007m11   .2893776    A           .          . |
  5. |        A   2007m12   .9847161    A           .          . |
     |-----------------------------------------------------------|
  6. |        A    2008m1   .9066443    A           .          . |
  7. |        A    2008m2   .1602792    A           .          . |
  8. |        A    2008m3   .2316082    A           .          . |
  9. |        A    2008m4   .7238458    A           .          . |
 10. |        A    2008m5   .6402453    A           .          . |
     |-----------------------------------------------------------|
 11. |        A    2008m6   .7056089    A           .          . |
 12. |        A    2008m7   .2857377    A           .          . |
 13. |        A    2008m8   .9954505    A           .          . |
 14. |        A    2008m9   .0044654    A           .          . |
 15. |        A   2008m10    .776125    A           .          . |
     |-----------------------------------------------------------|
 16. |        A   2008m11   .0321118    A           .          . |
 17. |        A   2008m12    .371257    A           .          . |
 18. |        A    2009m1   .0281739    A           .          . |
 19. |        A    2009m2   .1809183    A           .          . |
 20. |        A    2009m3   .9615714    A           .          . |
     +-----------------------------------------------------------+

[Sofia Frascalo]

Hi everyone,


Thank you for the suggestions and examples.
Clyde, when I run your code I get the following error message:

time variable not set

maybe this has to do with the fact that my industry variables are separate vertical columns. it looks more like a matrix.
Please find a screenshot attached where you can see how my data looks like.
Furthermore, N_i is the sum of the intercept and the residual estimated from a rolling window so the values are non-zero.

Do you know how to solve this error?

[Nick Cox]

Sofia: Please provide a readable data example using dataex as requested in the FAQ Advice. What's evident in your image is that your raw data are often less than -1 so even adding 1 won't be sufficient to make all values positive. Therefore any approach to such products based on taking logarithms will fail.

I address the question of cumulative products using what appears to be a similar dataset. The numbers are probably silly but I wish to show that what follows always produces the correct results if there are any zeros or negative values.

Without trying to discuss all the complications that might arise for your data, most of which we can't see, I note the following.

1. The role here of tsset is to lay down a framework in which tsfill will make gaps explicit as missing values. Otherwise rangestat does not assume tsset or xtset.

2. If there are any missing values, because of gaps or already in the data, they won't be ignored. For financial calculations I believe this to be the better approach, but code could be written to ignore missings.

3. There are two sets of results, one for the cumulative product and one for the cumulative product of 1 + the data. As you wish to subtract 1 afterwards you can either do that yourself or edit the Mata code.

This is meant to be indicative, not definitive.

clear
set obs 24
gen dm = ym(1951, 12) + _n
format dm %tm
tsset dm
tsfill
local j = 1
foreach v in N_agric N_food N_beer {
    gen `v' = `j' * (-13 + _n )
    local ++j
}

* define a Mata function that produces a running product
mata:  
mata clear

real rowvector product(real matrix X) {
    real rowvector product
    real scalar i
    product = X[1,]
    if (rows(X) == 1) return(product)
    
    for (i = 2; i <= rows(X); i++) {
        product = product :* X[i,]
    }
    return(product)
    }
end

rangestat (product) N_* , interval(dm -11 0)

. list dm N_* product?

     +----------------------------------------------------------------------------+
     |      dm   N_agric   N_food   N_beer     product1     product2     product3 |
     |----------------------------------------------------------------------------|
  1. |  1952m1       -12      -24      -36          -12          -24          -36 |
  2. |  1952m2       -11      -22      -33          132          528         1188 |
  3. |  1952m3       -10      -20      -30        -1320       -10560       -35640 |
  4. |  1952m4        -9      -18      -27        11880       190080       962280 |
  5. |  1952m5        -8      -16      -24       -95040     -3041280    -23094720 |
     |----------------------------------------------------------------------------|
  6. |  1952m6        -7      -14      -21       665280     42577920    4.850e+08 |
  7. |  1952m7        -6      -12      -18     -3991680   -5.109e+08   -8.730e+09 |
  8. |  1952m8        -5      -10      -15     19958400    5.109e+09    1.309e+11 |
  9. |  1952m9        -4       -8      -12    -79833600   -4.087e+10   -1.571e+12 |
 10. | 1952m10        -3       -6       -9    2.395e+08    2.452e+11    1.414e+13 |
     |----------------------------------------------------------------------------|
 11. | 1952m11        -2       -4       -6   -4.790e+08   -9.810e+11   -8.485e+13 |
 12. | 1952m12        -1       -2       -3    4.790e+08    1.962e+12    2.546e+14 |
 13. |  1953m1         0        0        0            0            0            0 |
 14. |  1953m2         1        2        3            0            0            0 |
 15. |  1953m3         2        4        6            0            0            0 |
     |----------------------------------------------------------------------------|
 16. |  1953m4         3        6        9            0            0            0 |
 17. |  1953m5         4        8       12            0            0            0 |
 18. |  1953m6         5       10       15            0            0            0 |
 19. |  1953m7         6       12       18            0            0            0 |
 20. |  1953m8         7       14       21            0            0            0 |
     |----------------------------------------------------------------------------|
 21. |  1953m9         8       16       24            0            0            0 |
 22. | 1953m10         9       18       27            0            0            0 |
 23. | 1953m11        10       20       30            0            0            0 |
 24. | 1953m12        11       22       33            0            0            0 |
     +----------------------------------------------------------------------------+

.


mata :

real rowvector add1product(real matrix X) {
    real rowvector product
    real scalar i
    product = 1 :+ X[1,]
    if (rows(X) == 1) return(product)
    
    for (i = 2; i <= rows(X); i++) {
        product = product :* (1 :+ X[i,])
    }
    return(product)
    }
end

rangestat (add1product) N_* , interval(dm -11 0)

. list dm N_* add1product?

     +---------------------------------------------------------------------------+
     |      dm   N_agric   N_food   N_beer   add1pro~1   add1prod~2   add1prod~3 |
     |---------------------------------------------------------------------------|
  1. |  1952m1       -12      -24      -36         -11          -23          -35 |
  2. |  1952m2       -11      -22      -33         110          483         1120 |
  3. |  1952m3       -10      -20      -30        -990        -9177       -32480 |
  4. |  1952m4        -9      -18      -27        7920       156009       844480 |
  5. |  1952m5        -8      -16      -24      -55440     -2340135    -19423040 |
     |---------------------------------------------------------------------------|
  6. |  1952m6        -7      -14      -21      332640     30421755    3.885e+08 |
  7. |  1952m7        -6      -12      -18    -1663200   -3.346e+08   -6.604e+09 |
  8. |  1952m8        -5      -10      -15     6652800    3.012e+09    9.245e+10 |
  9. |  1952m9        -4       -8      -12   -19958400   -2.108e+10   -1.017e+12 |
 10. | 1952m10        -3       -6       -9    39916800    1.054e+11    8.136e+12 |
     |---------------------------------------------------------------------------|
 11. | 1952m11        -2       -4       -6   -39916800   -3.162e+11   -4.068e+13 |
 12. | 1952m12        -1       -2       -3           0    3.162e+11    8.136e+13 |
 13. |  1953m1         0        0        0           0   -1.375e+10   -2.325e+12 |
 14. |  1953m2         1        2        3           0    1.964e+09    2.906e+11 |
 15. |  1953m3         2        4        6           0   -5.169e+08   -7.014e+10 |
     |---------------------------------------------------------------------------|
 16. |  1953m4         3        6        9           0    2.128e+08    2.698e+10 |
 17. |  1953m5         4        8       12           0   -1.277e+08   -1.525e+10 |
 18. |  1953m6         5       10       15           0    1.081e+08    1.220e+10 |
 19. |  1953m7         6       12       18           0   -1.277e+08   -1.363e+10 |
 20. |  1953m8         7       14       21           0    2.128e+08    2.142e+10 |
     |---------------------------------------------------------------------------|
 21. |  1953m9         8       16       24           0   -5.169e+08   -4.869e+10 |
 22. | 1953m10         9       18       27           0    1.964e+09    1.704e+11 |
 23. | 1953m11        10       20       30           0   -1.375e+10   -1.057e+12 |
 24. | 1953m12        11       22       33   4.790e+08    3.162e+11    1.796e+13 |
     +---------------------------------------------------------------------------+