function y = bbs3(myy, myx, num, k)
% non-overlapped block bootstrap
% num: times 
% k: block length
% assume a random walk for the DGP of myx
% 198301 - 202203 

T = size(myy,1);
j = mod(T,k);   % first j observations are fixed
n = (T-j)/k;    % the number of blocks
tboots = nan(num,1);

ymean = mean(myy,1);
res = myy - ymean;

res1 = reshape(res(j+1:end,1),k,n);
res2 = reshape(res(j+1:end,2),k,n);
res3 = reshape(res(j+1:end,3),k,n);
ypred = myy - res;

xdiff = diff(myx(j:end,:),1);
xdiff1 = reshape(xdiff(:,1),k,n);
xdiff2 = reshape(xdiff(:,2),k,n);
xdiff3 = reshape(xdiff(:,3),k,n);
xboots = nan(T,3);
xboots(1:j,:) = myx(1:j,:);
xcell = cell(T,1);

for i = 1:num
    ind = randi(n,n,1);
    x = [xboots(j,:); reshape(xdiff1(:,ind),[],1) reshape(xdiff2(:,ind),[],1) reshape(xdiff3(:,ind),[],1)];
    xboots(j:end,:) = cumsum(x);
    for ii = 1:T
        xcell{ii} = [eye(3), xboots(ii,:)'];
    end
    myresid = [reshape(res1(:,ind),[],1),reshape(res2(:,ind),[],1),reshape(res3(:,ind),[],1)];
    yb = ypred + [res(1:j,:) ;myresid];
    
    [b,~,~,CovB] = mvregress(xcell,yb, 'algorithm','mvn');
    tboots(i) = b(4)/sqrt(CovB(4,4));
    if mod(i,1000)==0
        disp(i)
    end
end

y = tboots;

end 