练手系列6: OLS regression in a rolling window

shiyiming

Can't input Chinese at workstation, but here is an interesting question. I have a solution I think that is the fastest so far, but I would like to see how you guys come up with noval solutions:

1. Easy One:
    I have a dataset with 500K observations, one independent variable X and one dependent variable Y. I want to conduct an OLS regression sequentially in a 20-observation rolling window and obtain the slope coefficient for X, say \beta. Note that this Linear Normal regression model should have an intercept, but I am not interested in the intercept term.

2. Difficult One:
    I have a dataset with 500K observations, M independent variables X1-Xm where m ranges from 2 to 10 and one dependent variable Y. I want to conduct a OLS regression in a 20 observation rolling window and obtain the slope coefficients for X1-Xm, a.k.a \beta_1 to \beta_m. Note that your linear normal regression model should have an intercept term but it is of no interets to me.

I will post my solutions to both questions after 5 followup posts or after this weekend, whichever comes first.

All solutions will be gauged on my PC or Shiyiming's machine based on total real time and memory consumption. Folks, fire up you knowledge on Stat and SAS! <!-- s:D --><img src="{SMILIES_PATH}/icon_biggrin.gif" alt=":D" title="Very Happy" /><!-- s:D -->

shiyiming

好难啊，怎么又扯上我的机器了，哈哈哈。 <!-- s:lol: --><img src="{SMILIES_PATH}/icon_lol.gif" alt=":lol:" title="Laughing" /><!-- s:lol: -->

shiyiming

第一个还好吧，我给了两个方法，第一种方法在我公司的机器上总计用时2秒，在我家的烂PC上总计用时7秒；第二种方法应该是一个比较有经验的SAS programmer能想到的，在我公司的机器上总计用时71秒，在我家总计用时104秒。当然还有用macro来慢慢loop的，那就慢慢等上几个小时吧。。。两种方法都不需要很难的coding。

第二个题需要有点线性代数的底子，coding稍微难点

shiyiming

50万的记录即使读入数据在我的机子上都要超过7秒。差不多50万次运算（regression），7秒真是难以想象。我是2g的内存，cpu速度忘了，用时5.5×60秒。

shiyiming

to jingju11
1.86G的C2D啊，3GB内存，不过2GB也够了，这个方法用了53MB的内存，比较耗。
还有更快的，就是要多写点代码，数据只要扫一遍。今晚调一下这个代码。
你的硬盘估计不是很给力

shiyiming

你指的是moving window，而不是每20个记录，right？也就是说如果一个一个的fit model要差不多500k个是吧？

shiyiming

to jingju11
是moving window，也就是比如说记录1-20算一组，做OLS回归； 记录2-21算一组，做OLS回归，记录3-23算一组，做回归，以此类推。
所以如果要每一组都fit model的话，确实要499981个回归

我弄这个练手系列6就是要说明搞统计的，一定要对算法有了解

你那个5.5X60秒是怎么来的？

下面是我的log。 我故意把notes和源代码都隐藏了。

2053  %let t0=%sysfunc(datetime(), datetime24.);
2054  options nosource nonotes;
2057  %let t1=%sysfunc(datetime(), datetime24.);
2058  %put Start @ &amp;t0;
Start @ 21DEC2010:22:09:14
2059  %put End   @ &amp;t1;
End   @ 21DEC2010:22:09:21

shiyiming

to oloolo
能用OLS的计算公式计算吗? 这个是最快的，一遍搞定～～一个array{20,2}. FIFO处理方式。
原来学计量的时候都是用gauss矩阵算所有的。

shiyiming

to jeozu
当然可以，不用IML，直接在data step里面写
等你的代码了

shiyiming

to oloolo
逻辑都全部说出来了也就不用再写了吧～～
对了， IML实现那种算法反而很蹩脚～～
我现在挑灯夜战写开题报告.一大把年纪回到学校上课不容易，更不容易的是还要跟一群小青年比体力～～
献个丑，1.8G, 900M, 1CPU 虚拟机 0.49m
[code:22kfj5cg]
%let N=50;
data a;
do i=1 to 500000;
   X=ranuni(0);
   Y=0&#46;25*X+rannor(0)+1;
   output;
end;
run;
%macro A(t);
ifn(_n_&lt;=&amp;N, 0, &amp;t{i})
%mend;
options mprint ;
data b;
   array R{0&#58;%eval(&amp;N-1)} _temporary_;
   array L{0&#58;%eval(&amp;N-1)} _temporary_;
   set a;
   i=mod(_n_,&amp;N);
   sum_X+X-%A(R);
   sum_Y+Y-%A(L);
   sum_XY+X*Y-%A(R) * %A(L);
   sum_X2+X**2-%A(R)**2;
   if _n_ &gt;=&amp;N then do;
      a=(&amp;n*sum_XY-sum_X*sum_Y)/(&amp;n*sum_X2-sum_X**2);
      b=sum_Y/&amp;N -a*sum_X/&amp;N;
      output;
   end;
   R{i}=X;
   L{i}=Y;
keep a b i;
run;
proc sql noprint;
create table xx as
select mean(a) as a, mean(b) as b
from b;
quit;
[/code:22kfj5cg]