编写stationary bootstrap 程序

以下程序是我从网上找到的用MATLAB编写的实现该要求的程序：
function [bsdata, indexes]=stationary_bootstrap(data,w,B);
% PURPOSE:
%mplements Politis' continuous bootstrap for bootstrapping unit root series
%
% USAGE:
%     [bsdata, indexes]=stationary_bootstrap(data,w,B);
%
% INPUTS:
%     data: T by 1 matrix of data to be bootstrapped(shoudl be unit root)
%     w:    Desired average windows length
%     B:    Number fo bootstraps
%
% OUTPUTS:
%     bsdata: T x B matrix of bootstrapped data
%     indexes: T by B matrix of location sof the original data(data(indexes)=bsdata;
%
% COMMENTS:
%
% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 2    Date: 12/31/2001


p=1/w;
[t,k]=size(data);
bsdata=zeros(t,B);
indexes=zeros(t,B);
indexes(1,:)=ceil(t*rand(1,B));
bsdata(1,:)=data(indexes(1,:))';

for i=2:t
    select=rand(1,B)<p;
    indexes(i,:)=ceil((t-1)*rand(1,B)).*select+(mod(~select.*(indexes(i-1,:)),t)+1);
    bsdata(i,:)=data(indexes(i,:))';
end
希望给楼上几位以参考和启发，看看stationary bootstrap 到底是怎么回事。希望能够用IML给出相应的程序，谢谢！

harlon1976 发表于 2011-8-18 21:14
以下程序是我从网上找到的用MATLAB编写的实现该要求的程序：
function =stationary_bootstrap(data,w,B);
...

I cannot understand some part of the matlab code

(indexes(i-1,:)

I though index(i,.) needs a number

I search the methods and find they are sometimes different.
one was like that:

Algorithm (Stationary Bootstrap)
http://www.kevinsheppard.com/images/5/55/Week_1_slides.pdf

that is way much easier than what we did.
JingJu

indexes(i-1,:)就像IML中的矩阵表示第i-1行的所有元素，谢谢你提供的这些程序，我也不懂MATLAB,只要再学习一下，然后利用它的思路用SAS/IML来完成了。从这个MATLAB程序看，似乎以概率p取所有元素（而不是排除我原来理解的下一个元素），而以1-p取下一个元素，我要再查查这个概念。
这个MATLAB程序的过程是这样的：首先从1-T任意选取第一个观测号，不妨记为u(i-1)，然后产生一个【0,1】上的均匀分布随机数v，对于事先给定的概率p=1/m，如果v大于等于p，显然这个概率为1-p,则第二个观测号为u(i)=u(i-1)+1,为了防止u(i)大于T，一旦出现这种情况就用它u(i)-T来控制；如果v小于等于p（这个概率为p），则从1-T中任意去一个即可。
因此从这个程序来看，应该以概率p取所有元素，包括u(i)元素，不知我的理解是否正确。以上是其算法，但这个程序本身并没有给出u(i)-T语句来控制出现超过最大观测号的情况，这是为什么呢？

而我所找到的这个程序中似乎有个问题就是这句：
indexes(i,:)=ceil((t-1)*rand(1,B)).*select+(mod(~select.*(indexes(i-1,:)),t)+1);
这里的ceil((t-1)而不是ceil((t)，这里的t就是T,似乎表明是从剩下的T-1中来选取，而且后面的MOD语句中的最外层括号里面还加1，这个MOD函数就是为了防止观测好超过T而设置的，一旦某次的观测好为T，则下一个以概率1-p为第一个观测，但这个外层括号似乎表明这个1是和MOD在一起的，但这个括号似乎又能去掉，此时加1似乎又能和前面的连在一起，弥补前面的T-1而不是T的影响，所以问题就在这里。

harlon1976 发表于 2011-8-19 06:51
而我所找到的这个程序中似乎有个问题就是这句：
indexes(i,:)=ceil((t-1)*rand(1,B)).*select+(mod(~selec ...

我想这一部分是可以理解的：
如果select =1， 那么这一句其实就是
indexes(i,:)=ceil((t-1)*rand(1,B))+1;选取从2 到T的任一值。但是select =0
的时候我有些看不懂
mod((indexes(i-1,:)),t)+1=i？
京剧

jingju11 发表于 2011-8-19 08:24
我想这一部分是可以理解的：
如果select =1， 那么这一句其实就是
indexes(i,:)=ceil((t-1)*rand(1,B)) ...

mod((),t)+1这个实际上就是表示取上次位置上的相邻的下一个位置，就是以概率1-p取的，mod是取余数，这个SAS也有，就是indexes(i-1,:)[上次的位置]除以T后的余数，例如上次为8，如果T=10，则余数为8，再加1，则就是9 了。我的不是当select =1时，为什么ceil((t-1)而不是ceil((t)，为什么是选取从2 到T的任一值，而不是从1-T中取值呢？

harlon1976 发表于 2011-8-19 10:02
mod((),t)+1这个实际上就是表示取上次位置上的相邻的下一个位置，就是以概率1-p取的，mod是取余数，这个S ...

因为此时第一个位置已经有值了。

凑个热闹,SAS不懂,Matlab略懂

假设:

data=[8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.0]'
B=2
w=5
p=1/w
% Initial positions
indices =
     2     2
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
     0     0
% Set up the random numbers

    0.9991    0.1564
    0.1711    0.8555
    0.0326    0.6448
    0.5612    0.3763
    0.8819    0.1909
    0.6692    0.4283
    0.1904    0.4820
    0.3689    0.1206
    0.4607    0.5895
    0.9816    0.2262
%select=rand(t,B)<p
select =

     0     1
     1     0
     1     0
     0     0
     0     1
     0     0
     1     0
     0     1
     0     0
     0     0
%indices(select)=ceil(rand(1,sum(sum(select)))*t)
indices =

     2     3
     4     0
     6     0
     0     0
     0     7
     0     0
     3     0
     0     3
     0     0
     0     0
%when i=2
indices =

     2     3
     4     4
     6     0
     0     0
     0     7
     0     0
     3     0
     0     3
     0     0
     0     0
%when i=3
indices =

     2     3
     4     4
     6     5
     0     0
     0     7
     0     0
     3     0
     0     3
     0     0
     0     0
%%%%
.....
.....
%%%%when i=10
indices =

     2     3
     4     4
     6     5
     7     6
     8     7
     9     8
     3     9
     4     3
     5     4
     6     5
bsdata =

    8.2000    8.3000
    8.4000    8.4000
    8.6000    8.5000
    8.7000    8.6000
    8.8000    8.7000
    8.9000    8.8000
    8.3000    8.9000
    8.4000    8.3000
    8.5000    8.4000
    8.6000    8.5000

不知哪位能够懂这两门语言呢？楼上的兄弟的程序与我给的程序好像不一样啊！

jingju11 发表于 2011-8-19 03:04
I cannot understand some part of the matlab code
I though index(i,.) needs a number

你好，你上面的链接我打不开，能给我那些代码吗