function y = bbs4(myy, myx1, myx2, num, k, ix)
% non-overlapped block bootstrap： constant + 2 explainary variables
% num: times 
% k: block length
% assume a random walk for the DGP of myx
% ix: the explainary variable we want to test 

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);

xinit11 = myx1(1:j,1);
xinit12 = myx1(1:j,2);
xinit13 = myx1(1:j,3);
xinit21 = myx2(1:j,1);
xinit22 = myx2(1:j,2);
xinit23 = myx2(1:j,3);

xdiff1 = diff(myx1(j:end,:),1);
xdiff2 = diff(myx2(j:end,:),1);

xdiff11 = reshape(xdiff1(:,1),k,n);
xdiff12 = reshape(xdiff1(:,2),k,n);
xdiff13 = reshape(xdiff1(:,3),k,n);

xdiff21 = reshape(xdiff2(:,1),k,n);
xdiff22 = reshape(xdiff2(:,2),k,n);
xdiff23 = reshape(xdiff2(:,3),k,n);

xboots1 = nan(T,3);
xboots1(1:j,:) = [xinit11 xinit12 xinit13];
xboots2 = nan(T,3);
xboots2(1:j,:) = [xinit21 xinit22 xinit23];
xcell = cell(T,1);

if ix == 1
    bootsx = cell(T,1);
    for i = 1:T
        bootsx{i} = [eye(3), myx2(i,:)'];
    end
    [bhat,~,E,~] = mvregress(bootsx,myy, 'algorithm','mvn');
    
    res1 = reshape(E(j+1:end,1),k,n);
    res2 = reshape(E(j+1:end,2),k,n);
    res3 = reshape(E(j+1:end,3),k,n);
    
    for i = 1:num
        ind = randi(n,n,1);
        
        xboots1(j+1:end,1) = reshape(xdiff11(:,ind),[],1);
        xboots1(j+1:end,2) = reshape(xdiff12(:,ind),[],1);
        xboots1(j+1:end,3) = reshape(xdiff13(:,ind),[],1);
        xboots1(j:end,:) = cumsum(xboots1(j:end,:),1);
        
        xboots2(j+1:end,1) = reshape(xdiff21(:,ind),[],1);
        xboots2(j+1:end,2) = reshape(xdiff22(:,ind),[],1);
        xboots2(j+1:end,3) = reshape(xdiff23(:,ind),[],1);
        xboots2(j:end,:) = cumsum(xboots2(j:end,:),1);
        
        for ii = 1:T
            xcell{ii} = [eye(3), xboots1(ii,:)', xboots2(ii,:)'];
        end
        
        ypred = nan(T,3);
        myresid = nan(T,3);
        myresid(1:j,:) = E(1:j,:);
        ypred(:,1) = [ones(T,1),xboots2(:,1)]*[bhat(1); bhat(4)];
        ypred(:,2) = [ones(T,1),xboots2(:,2)]*[bhat(2); bhat(4)];
        ypred(:,3) = [ones(T,1),xboots2(:,3)]*[bhat(3); bhat(4)];
        myresid(j+1:end,:) = [reshape(res1(:,ind),[],1),reshape(res2(:,ind),[],1),reshape(res3(:,ind),[],1)];
        yb = ypred + 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 
elseif ix==2
    bootsx = cell(T,1);
    for i = 1:T
        bootsx{i} = [eye(3), myx1(i,:)'];
    end
    [bhat,~,E,~] = mvregress(bootsx,myy, 'algorithm','mvn');
    
    res1 = reshape(E(j+1:end,1),k,n);
    res2 = reshape(E(j+1:end,2),k,n);
    res3 = reshape(E(j+1:end,3),k,n);
    
    for i = 1:num
        ind = randi(n,n,1);
        
        xboots1(j+1:end,1) = reshape(xdiff11(:,ind),[],1);
        xboots1(j+1:end,2) = reshape(xdiff12(:,ind),[],1);
        xboots1(j+1:end,3) = reshape(xdiff13(:,ind),[],1);
        xboots1(j:end,:) = cumsum(xboots1(j:end,:),1);
        
        xboots2(j+1:end,1) = reshape(xdiff21(:,ind),[],1);
        xboots2(j+1:end,2) = reshape(xdiff22(:,ind),[],1);
        xboots2(j+1:end,3) = reshape(xdiff23(:,ind),[],1);
        xboots2(j:end,:) = cumsum(xboots2(j:end,:),1);
        
        for ii = 1:T
            xcell{ii} = [eye(3), xboots1(ii,:)', xboots2(ii,:)'];
        end
        
        ypred = nan(T,3);
        myresid = nan(T,3);
        myresid(1:j,:) = E(1:j,:);
        ypred(:,1) = [ones(T,1),xboots1(:,1)]*[bhat(1); bhat(4)];
        ypred(:,2) = [ones(T,1),xboots1(:,2)]*[bhat(2); bhat(4)];
        ypred(:,3) = [ones(T,1),xboots1(:,3)]*[bhat(3); bhat(4)];
        myresid(j+1:end,:) = [reshape(res1(:,ind),[],1),reshape(res2(:,ind),[],1),reshape(res3(:,ind),[],1)];
        yb = ypred + myresid;
        
        [b,~,~,CovB] = mvregress(xcell,yb, 'algorithm','mvn');
        tboots(i) = b(5)/sqrt(CovB(5,5));
        if mod(i,1000)==0
            disp(i)
        end
    end
end
    
y = tboots;

end 