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

shiyiming

如果y=x仅一个变量我可以做到
real time 1.71s
cpu time 0.76s
我只运算coefficient, that is, beta1
如果是老老实实做regression（将近500k次），可以报告se， p等等，我差不多可以控制在5分钟之内
电脑cpu2.33G， memory4.0G, network processing

shiyiming

to jingju11
好，比我强
代码是啥？

shiyiming

to jeozu
要的就是代码啊

shiyiming

[code:1y9thflp]data had;
        keep x1 y;
        array x[10];
        do i =1 to 5e5;
                do j =1 to 10;
                        x[j] =round(ranuni(111)*10);
                end;
                y =round(rannor(11),0.01);
                output;
        end;
run;
**创建数据的时间应该不算吧;
sasfile had load;
data have;
        set had;
        array t1[0:19] _temporary_; array t2[0:19] _temporary_; array t3[0:19] _temporary_; array t4[0:19] _temporary_;
        _dim =mod(_n_, 20);
        t1[_dim] =x1; t2[_dim] =y; t3[_dim] =x1*y; t4[_dim] =x1*x1;
        if _n_ >=20 then b1 =(sum(of t3[*])-sum(of t1[*])*sum(of t2[*]) /20)        /(sum(of t4[*])-20*mean(of t1[*])**2);
run;
sasfile had close;[/code:1y9thflp]

shiyiming

本想用多维数组，最后没有找到reference的方法。所以程序很难看。
如果做regression的话，白天写了个程序用by processing,时间大概5.5分钟。程序懒的再编写了。白天network比较拥挤。如果是现在。我想应该在5分钟之内。

shiyiming

to jingju11
嗯，直接用公式。扩展到多变量就把第二个解决了
多维数组：
def: array _x{4,20} _temporary_;
ref: _x[3,2]=....
fun:  dim(_x, 2) ...

shiyiming

说实话。我只能依靠_temporary_数组的记忆功能。如果要loop的话，我也没有哪个本事写出来。
我想不明白如何推广到多变量的形式。不过既然每个x形成的面是互相垂直的，说不定简单的相加即可。到此为止。
我现在的兴趣是探讨nlmixed过程的converge问题。希望有兴趣的同志多多探讨一下。谢谢。

京剧

shiyiming

好吧，既然大家没啥兴趣，就算了吧
[code:2krjs26y]
/* Method 1.*/
data rest0;
        set test;
                array _x{4} _temporary_;
                array _a{2,20}  _temporary_;
                m=mod(_n_-1, 20)+1;
                _a[1, m]=x1; _a[2,m]=y;
                link filler;

                m2=mod(_n_-20, 20)+1;
                if _n_>=20 then do;
                  if _n_>20 then do;
              link deduct;
                  end;
                  beta=(_x[2]-_x[1]*_x[4]/20)/(_x[3]-_x[1]**2/20);
                  intercept=_x[4]/20 - beta*_x[1]/20;
                  keep  seq   intercept  beta  ;
                  output;
                end;
                return;      
filler:
     _x[1]+x1;
         _x[2]+x1*y;
         _x[3]+x1**2;
         _x[4]+y;
return;
deduct:
     _x[1]=_x[1]-_a[1,m2];
         _x[2]=_x[2]-_a[1,m2]*_a[2,m2];
         _x[3]=_x[3]-_a[1,m2]**2;
         _x[4]=_x[4]-_a[2,m2];
return;
run;




/* Method 2.*/

%macro wrap;
%let window=20;
data testv/view=testv;
     set test;
       xy=x1*y;  
run;

proc expand data=testv  method=none  out=summary(keep=seq sumxy  sumx1  sumy  ussx1  may  max);
       convert  x1=sumx1/transformout=(movsum &window);
       convert  xy=sumxy/transformout=(movsum &window);
       convert  x1=ussx1/transformout=(movuss &window);
       convert  y =sumy /transformout=(movsum &window);
       convert  y =may / transformout=(movave &window);
       convert  x1 =max / transformout=(movave &window);  
run;

data result1;
     set summary(firstobs=&window);
       beta = (sumxy - sumx1*sumy/&window)/(ussx1 - sumx1/&window.*sumx1);  
       alpha= may - beta*max;
       keep seq  beta  alpha;       
run;
%mend;

%let t0=%sysfunc(datetime(), datetime24.);
*options nosource nonotes;
%wrap;
options source notes;
%let t1=%sysfunc(datetime(), datetime24.);
%put Start @ &t0;
%put End   @ &t1;



/* Method 3.*/
%let t0=%sysfunc(datetime(), datetime.);

data test2v/view=test2v;
       set test;
       array _x{2, 20} _temporary_ (20*0 20*0);
       k=mod(_n_-1, 20)+1;
       _x[1, k]=x1; _x[2, k]=y;
       if _n_>=20 then do;
          do j=1 to dim(_x, 2);
               x=_x[1, j]; y=_x[2, j];
               output;
               keep time x y;
            end;
       end;
run;

ods select none;
proc  reg data=test2v  outest=res2(keep=time x intercept);
         by time;
         model y = x;
run;quit;
ods select all;

%let t1=%sysfunc(datetime(), datetime.);
%put Start @ &t0;
%put End   @ &t1;

[/code:2krjs26y]