function y = panel_bootstrap(myy, myx, num, myseed)
% bootstrap for the panel regression of r* and macro-variables
% non-overlapped block bootstrap
% assume a random walk for r*(myy) 
% assume an AR(1) for macro-variables(myx)
%
% num: times 
% myseed: seed

[T,N] = size(myy);
% j = mod(T,k); 
% n = (T-j)/k;    % the number of blocks
tboots = nan(num,1);

rng(myseed)

if N == 3
    % Three countries
    ydiff = diff(myy,1);
    std1 = std(ydiff(:,1));
    std2 = std(ydiff(:,2));
    std3 = std(ydiff(:,3));
    y_innov1 = std1*randn(T-1,num);
    y_innov2 = std2*randn(T-1,num);
    y_innov3 = std3*randn(T-1,num);

    % ydiff1 = reshape(ydiff(:,1),k,n);
    % ydiff2 = reshape(ydiff(:,2),k,n);
    % ydiff3 = reshape(ydiff(:,3),k,n);
    % yboots = nan(T,3);
    % yboots(1:j,:) = myy(1:j,:);  % first k observations are fixed

    Mdl = arima(1,0,0);
    Est1 = estimate(Mdl,myx(:,1),'Display','off');
    Est2 = estimate(Mdl,myx(:,2),'Display','off');
    Est3 = estimate(Mdl,myx(:,3),'Display','off');
    % const1 = Est1.Constant;
    % const2 = Est2.Constant;
    % const3 = Est3.Constant;
    % ar_coeff1 = Est1.AR{:};
    % ar_coeff2 = Est2.AR{:};
    % ar_coeff3 = Est3.AR{:};
    % std1 = sqrt(Est1.Variance);
    % std2 = sqrt(Est2.Variance);
    % std3 = sqrt(Est3.Variance);
    xboots = nan(T,3);
    % xboots(1:j,:) = myx(1:j,:);
    xcell = cell(T,1);

    for i = 1:num
        % ind = randi(n,n,1);
        % y = [myy(j,:); reshape(ydiff1(:,ind),[],1) reshape(ydiff2(:,ind),[],1) reshape(ydiff3(:,ind),[],1)];
        % yboots(j:end,:) = cumsum(y);
        y = [myy(1,:); y_innov1(:,i) y_innov2(:,1) y_innov3(:,i)];
        yboots = cumsum(y);
        xboots(:,1) = simulate(Est1,T);
        xboots(:,2) = simulate(Est2,T);
        xboots(:,3) = simulate(Est3,T);
        for ii = 1:T
            xcell{ii} = [eye(3), xboots(ii,:)'];
        end
        [b,~,~,CovB] = mvregress(xcell,yboots, 'algorithm','cwls');
        tboots(i) = b(4)/sqrt(CovB(4,4));
        if mod(i,1000)==0
            disp(i)
        end
    end
else
    % Two countries
    ydiff = diff(myy,1);
    std1 = std(ydiff(:,1));
    std2 = std(ydiff(:,2));
    y_innov1 = std1*randn(T-1,num);
    y_innov2 = std2*randn(T-1,num);

    Mdl = arima(1,0,0);
    Est1 = estimate(Mdl,myx(:,1),'Display','off');
    Est2 = estimate(Mdl,myx(:,2),'Display','off');
    xboots = nan(T,2);
    xcell = cell(T,1);

    for i = 1:num
        % ind = randi(n,n,1);
        % y = [myy(j,:); reshape(ydiff1(:,ind),[],1) reshape(ydiff2(:,ind),[],1) reshape(ydiff3(:,ind),[],1)];
        % yboots(j:end,:) = cumsum(y);
        y = [myy(1,:); y_innov1(:,i) y_innov2(:,1)];
        yboots = cumsum(y);
        xboots(:,1) = simulate(Est1,T);
        xboots(:,2) = simulate(Est2,T);
        for ii = 1:T
            xcell{ii} = [eye(2), xboots(ii,:)'];
        end
        [b,~,~,CovB] = mvregress(xcell,yboots, 'algorithm','cwls');
        tboots(i) = b(3)/sqrt(CovB(3,3));
        if mod(i,1000)==0
            disp(i)
        end
    end
end

    y = tboots;

end


