有谁知道 如何用宏语言实现bootstrap resampling?

jingju11 发表于 2010-4-7 08:59
另外，你可以把你bootstrap 分析 生成的data的过程简略地说一下吗？或者是sas程序。我对bootstrap的原理一直不很理解，许多paper对bootstrap的应用也很牵强，似是而非的。

其实我想做的boostrap属于残值再取样的，具体步骤是这样的：

1, 对每家公司，按公式：Yi=ai+bi*Bi,t+ci*Ci,t 做时间序列回归。得出每家公司的ai,bi,ci,以及一组残值ei,t;

2, 然后对每家公司的残值进行有放回的再抽（抽取的数目与原来的数目一样），得出bootstrapping 残值（ebi,t）;

3, 假设a=0, 将ebi,t带回公式Ybi= bi*Bi,t+ci*Ci,t+ebi,t,得出一个新的Ybi （bootstrapping Y）;

4; 再用这个新生成的Ybi （bootstrapping Y）去做 Ybi =ai+bi*Bi,t+ci*Ci,t 时间序列回归，得出一个新的ai （我称之为ab）,boostrapping alpha.

5. 步骤2－4，重复1000次。也就是说每家公司最后有1000个boostrapping alpha.

这个就是我想要做的，不知您觉得什么语句对我来说更合适些呢？其实对SAS，我了解很少，主要是个工具来运用的。

如果哪位有非常好的语句，欢迎指教！

1. data mydata;
   
2.         do co=10001 to 10020 by 1;
   
3.                 nobs=ceil(40+ranuni(123)*20);
   
4.                 do i=0 to nobs-1 by 1;
   
5.                         date=intnx("MONTH","01Jan2001"d,i,"END");
   
6.                         resid=ranuni(123);
   
7.                         output;
   
8.                 end;
   
9. 
   
10.         end;
    
11.         drop i nobs;
    
12.         format date yymmn6.;
    
13. run;
    
14. 
    
15. proc sql;
    
16.         create table t1 as
    
17.                 select
    
18.                         co
    
19.                         ,date
    
20.                         ,resid
    
21.                         ,count(*) as Cnt
    
22.                 from mydata
    
23.                 group by co
    
24.                 order by co,date
    
25.         ;
    
26. quit;
    
27. 
    
28. %MACRO T;
    
29. data t2;
    
30.         set t1;
    
31.         by co date;
    
32.         if first.co then nobs=0;
    
33.         nobs+1;
    
34.         %do i=1 %to 1000 %by 1;
    
35.                 resid_boot&i=input(compress(put(co,z5.)||put(ceil(ranuni(123)*Cnt),z3.)),8.);       
    
36.         %end;
    
37. 
    
38. run;
    
39. 
    
40. data t3;
    
41.         set t2(keep=co nobs resid);
    
42.         fmtname="REPLACE";
    
43.         start=input(compress(put(co,z5.)||put(nobs,z3.)),8.);
    
44.         end=start;
    
45.         label=resid;
    
46.         keep fmtname start label;
    
47. 
    
48. run;
    
49. 
    
50. proc format cntlin=t3;run;quit;
    
51. 
    
52. data t4;
    
53.         set t2;
    
54.         %do i=1 %to 1000 %by 1;
    
55.         resid_boot&i=input(put(resid_boot&i,replace.),best12.);
    
56.         %end;
    
57. run;
    
58. 
    
59. %MEND T;
    
60. %T;

复制代码

提供一种思路
写得很仓促，FORMAT替换很粗糙，希望能改进
思路是
放回随机抽取-->生成1-每个co下记录个数的随机整数-->替换整数为原值

希望是对的，对你有帮助

15# jingju11

用上面语句，是能实现resid的再抽样，但是它的排序是按照时间来排的，而且有的值重复抽取，会把这个重复抽取的值放在一起。

但是我想要得到的是：在resid中任抽一个放在第一个，放回，再抽一个放在第二个，放回。。。。。所以最后抽取出的结果，时间上一定是打乱的；而且重复抽到的值，不一定会在一起（有可能第5次抽到这个值，第10次再抽到这个值，但它们应该分布放在第5和第10的位置上）。

不知如何能够改进能够实现这个目的？

其实我原来的bootstrap语句(在13楼)能够做到这一点，但就是有错误，哪位帮忙看看，错在哪儿了？怎么改？


万分感激！

23# funwin
时间也要一起抽啊，多做个format了要......

funwin 发表于 2010-4-7 21:46
15# jingju11

用上面语句，是能实现resid的再抽样，但是它的排序是按照时间来排的，而且有的值重复抽取，会把这个重复抽取的值放在一起。

但是我想要得到的是：在resid中任抽一个放在第一个，放回，再抽一个放在第二个，放回。。。。。所以最后抽取出的结果，时间上一定是打乱的；而且重复抽到的值，不一定会在一起（有可能第5次抽到这个值，第10次再抽到这个值，但它们应该分布放在第5和第10的位置上）。

不知如何能够改进能够实现这个目的？

其实我原来的bootstrap语句(在13楼)能够做到这一点，但就是有错误，哪位帮忙看看，错在哪儿了？怎么改？


万分感激！

1. 
   
2. %macro BootMcr(b = 10);
   
3. PROC SQL NOPRINT;
   
4.   SELECT COUNT(DISTINCT CO) INTO :NCOs FROM Original ORDER BY CO;
   
5.   SELECT DISTINCT CO, COUNT(*) AS COUNT INTO :COs SEPARATED BY ' ' , :Nobs SEPARATED BY ' ' FROM Original GROUP BY CO;
   
6. QUIT;
   
7. Data %DO t = 1 %TO &NCOs; %LET CO = %SCAN(&COs, &t);  boots&t(WHERE = (CO = &CO)) %END;;
   
8.   Set Original;
   
9. Run;
   
10. Data bootstrap;
    
11.   Do sample = 1 to &b;
    
12.    %do t=1 %to &NCOs;
    
13.     %LET Nobs_= %SCAN(&Nobs, &t);
    
14.       Do i=1 to &Nobs_;
    
15.          x=CEIL(ranuni(&t)*&nobs_);
    
16.          Set boots&t point=x;
    
17.          Output;
    
18.        End;
    
19.      %End;
    
20.   end;
    
21.   Stop;
    
22. Run;
    
23. %MEND BootMcr;
    
24. %BootMcr
    
25. ;

复制代码

25# jingju11

好像是我想要的！容我再好好研究一下！
非常感谢您的帮助！！

楼主是不是想求回归系数的置信区间？如果这样，对重复抽取的数据你如何处理？因为bootstrapping的本质就是有放回随机抽样，如果设定抽取的样本数等于原来的样本数，那么很大可能有重复的。
比如，你第一次回归后从残差如下：
1月 0.1
2月 0.2
3月 0.3
4月 0.4
，bootstrapping抽样后，会出现这种情况：
1月 0.1
1月 0.1
3月 0.3
4月 0.4
，其中1月份重复抽取了，你想如何处理上面的数据？

重复没有关系，因为抽取后的残值，我代回原公式，但是假设alpha=0，重新得到一组新的bootstrapping y。然后再用时间序列回归，算出bootstrapping的alpha. 这个是bootstrap的残值再取样的方法。

可能版主想问的是我怎么把残值代回公式。
任意一个抽取的残值放在第一月，就作为第一个月的新残值。然后再抽放在第二个月，依此类推。比如抽取出来的前三个是5月，12月，5月的，那么就把这三个作为1，2，3月的新残值来用。。。。。
只是我对bootstrap残值再取样的理解.

25# jingju11

Jingju11,您的语言非常好，实现了我想要的bootstrap residual resampling。非常感谢！现在还有个问题，不知用sas宏语句（或者其他语句）如何能完成下面这些目的？ 之前的求了100次的残值为
_1, _2,…, _100

co	date	X1	X2	t	pre	_1	_2	…..	_100
10001	200501			0.89		0.023	0.0025	0.023	…
…..	…..			0.89		0.0025	-0.0154	0.023	…
10001	200812			0.89		-0.0154	-0.0154	0.0025	…
10003
……

对每家公司：得到新的y：
y1=pre+_1, y2=pre+_2……y100=pre+_100; (同时为了节省空间把残值_1,_2…._100去掉). 用新的y,
对每家公司每个y值分别做时间序列回归，proc reg; model y1=x1 x2;
…….  
model y100=x1 x2; 每家公司得出100个intercept, intercept的t-statistics; 再得出这100个t-statistics中有多少个是大于等于t (原来的t-value)的，得出一个p_boot=大于等于t 的个数/抽样数100 . 希望再帮忙解决一下我的问题。万分感激！