function [bsdata, indices]=stationary_bootstrap(data,w)
% Implements the stationary bootstrap for bootstrapping stationary, dependent series
% 
% USAGE:
% [BSDATA, INDICES] = stationary_bootstrap(DATA,B,W)
% 
% INPUTS:
% DATA - T by 1 vector of data to be bootstrapped
% B - Number of bootstraps
% W - Average block length. P, the probability of starting a new block is defined P=1/W
% 
% OUTPUTS:
% BSDATA - T by B matrix of bootstrapped data
% INDICES - T by B matrix of locations of the original BSDATA=DATA(indexes);
% 
% COMMENTS:
% To generate bootstrap sequences for other uses, such as bootstrapping vector processes,
% set DATA to (1:N)'. 
%
% See also block_bootstrap
% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 2 Date: 12/31/2001
t=size(data,1);
p=1/w; % Define the probability of a new block
% Large p ==> large bias small variance
indices=zeros(t,1); % Set up the bsdata and indices
indices(1,:)=ceil(t*rand(1,1)); % Initial positions
% Set up the random numbers
select=rand(t,1)<p;
indices(select)=ceil(rand(1,sum(sum(select)))*t);
for i=2:t % Determine whether we stay (rand>p) or move to a new starting value (rand<p)
indices(i,~select(i,:))=indices(i-1,~select(i,:))+1;
end
indices(indices>t) = indices(indices>t)-t;
% The indices make finding the bsdata simple
bsdata=data(indices,:);
end