%% Wild bootstrap with volatility scaling using cumulative returns (Table 4, by asset class) %%

%% Commodity %%
clear;clc
addpath 'E:\RESEARCH\TSMOM\Codes_new';
addpath 'E:\Dropbox\My toolbox';
load Return.mat;

Ns = 1; Ne = 24; 
Return = Return(Return(:,1)>=198501,:);
Return = Return(Return(:,6)<=Ne,:);

Rs = Return(:,2);
rf = Return(:,4)+1;
exVol = Return(:,3);
Rcum = Func_GenCum(rf, Rs);
Rcumvol = Rcum(2:end,:)./exVol(1:end-1);
Rcumvol = [nan(1,12); Rcumvol];
Rcumvol = [Rcumvol, Return(:,6)];
% Drop first 12 observations in case of NaN %
for ss = Ns:Ne
      aa = find(Rcumvol(:,13)==ss);
      Rcumvol(aa(1):aa(12),:) = [];
      bb = find(Return(:,6)==ss);
      Return(bb(1):bb(12),:) = [];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                       %            1       2       3   4    5         6-17
z0 = [z0, Rcumvol(:,1:12)];                           %  z0 =  [ym, xret, rf, id, vol, xcum/vol]

% Run predictive regression %
res = []; res_sign = [];
for h = 1:12
      yy = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % [ym, xret, rf, id, vol, xcum/vol]                                                                           
      yy = [yy; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % yy = [y, x,  ym, ID]; 
    end
      b = clusterreg(yy(:,1), [yy(:,2), ones(size(yy(:,2)))], yy(:,3));
      res = [res; b(1, [1, 3]), b(2, [1, 3])];      % res = [b, b_tstat, cons, cons_tstat]
      bsign = clusterreg(yy(:,1), [sign(yy(:,2)), ones(size(yy(:,2)))], yy(:,3));
      res_sign = [res_sign; bsign(1, [1, 3]), bsign(2, [1, 3])];
end

% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),2); Exhat_sign = nan(length(z0),2);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = [z0(i,5+h), 1]*[res(h,1), res(h,3)]';
        Exhat_sign(i+1,1) = [sign(z0(i,5+h)), 1]*[res_sign(h,1), res_sign(h,3)]';
    end
    for i = 2:length(z0)
        Exhat(i,2) = z0(i,6) - Exhat(i,1);                            
        Exhat_sign(i,2) = z0(i,6) - Exhat_sign(i,1);
    end
Pedhat = [Pedhat, Exhat(:,1)];
Pedhat_sign = [Pedhat_sign, Exhat(:,1)];
Reshat = [Reshat, Exhat(:,2)];
Reshat_sign = [Reshat_sign, Exhat_sign(:,2)];
end
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_sign, z0(:,4)];
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
for s = Ns:Ne
    line1 = find(Reshat(:,13)==s);
    Reshat(line1(1),:) = [];
    line2 = find(Reshat_sign(:,13)==s);
    Reshat_sign(line2(1),:) = [];
    line3 = find(Pedhat(:,13)==s);
    Pedhat(line3(1),:) = [];
    line4 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line4(1),:) = [];
    line5 = find(z0(:,4)==s);
    z0(line5(1),:) = [];
end

% Wild bootstrap for 1000 times %
B = 1000; 
beta_all = nan(B,12); tstat_all = nan(B,12);
beta_sign_all = nan(B,12); tstat_sign_all = nan(B,12);
for i = 1:B
    tic
    Boot_ret = [];y_boot1 = []; y_boot2 = [];
    for s = Ns:Ne
          rng(s*i);                                                  % set seed
          tY = z0(z0(:,4)==s,1);                           % z0 = [ym, xret, rf, id, vol, xcum/vol]    
          T = length(tY);
          
          Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
          Pedhat_sign_s = Pedhat_sign(Pedhat_sign(:,13)==s,:);
          Reshat_s = Reshat(Reshat(:,13)==s,:);
          Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);

          rand_b = [rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1)];
          
          y_boot_all = nan(T,12);
          y_boot_all_sign = nan(T,12);
          % Construct pseudo sample path %
          for t = 1:T
            vt = ((rand_b(t+100,:)>=0.5)*(1) + (rand_b(t+100,:)<0.5)*(-1));
            y_boot_all(t,:) = Pedhat_s(t,1:12) + Reshat_s(t,1:12).*vt;                        
            y_boot_all_sign(t,:) = Pedhat_sign_s(t,1:12) + Reshat_sign_s(t,1:12).*vt;     
          end
          y_boot1 = [y_boot1; y_boot_all];
          y_boot2 = [y_boot2; y_boot_all_sign];
   end
   z1 = [z0(:,1), z0(:,4), y_boot1];
   z2 = [z0(:,1), z0(:,4), y_boot2];
     
   res_i = []; res_sign_i = [];
  for h = 1:12
       yy1 = []; yy2 = [];
       for n = Ns:Ne 
       tt = find(z1(:,2) == n);                             %    1         2         3-14       
       a = z1(tt,:);                                                %   [ym      ID       cum/vol]  
            
       ta = find(z0(:,4)==n);                              %     1          2          3         4         5         6                                         
       d = z0(ta,:);                                               %   [ym       xret      rf        ID     vol   xcum/vol]
       yy1 = [yy1; a(2:end,2+h), d(1:end-1,5+h), a(2:end,[1, 2])];  %  use original independent variables and bootstrapped  dependent variables 
            
       td = find(z2(:,2)==n);
       c = z2(td,:);
       yy2 = [yy2; c(2:end,2+h), d(1:end-1,5+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
       end
       b = clusterreg(yy1(:,1), [yy1(:,2), ones(size(yy1(:,2)))], yy1(:,3));
       res_i = [res_i; b(1, [1, 3]), b(2, [1, 3])];
  
       bsign = clusterreg(yy2(:,1), [sign(yy2(:,2)), ones(size(yy2(:,2)))], yy2(:,3));
       res_sign_i = [res_sign_i; bsign(1, [1, 3]), bsign(2, [1, 3])];
   end
      beta_all(i,:) = res_i(:,1)';
      tstat_all(i,:) = res_i(:,2)';
      beta_sign_all(i,:) = res_sign_i(:,1)';
      tstat_sign_all(i,:) = res_sign_i(:,2)';
      toc
      i
end

save('E:\RESEARCH\TSMOM\Codes_new\Tab4_WildCom.mat','beta_all', 'beta_sign_all', 'tstat_all', 'tstat_sign_all')

%% Equity %%
clear;clc
addpath 'E:\RESEARCH\TSMOM\Codes_new';
addpath 'E:\Dropbox\My toolbox';
load Return.mat;

Ns = 25; Ne = 33; 
Return = Return(Return(:,1)>=198501,:);
Return = Return(Return(:,6)>=Ns & Return(:,6)<=Ne,:);

Rs = Return(:,2);
rf = Return(:,4)+1;
exVol = Return(:,3);
Rcum = Func_GenCum(rf, Rs);
Rcumvol = Rcum(2:end,:)./exVol(1:end-1);
Rcumvol = [nan(1,12); Rcumvol];
Rcumvol = [Rcumvol, Return(:,6)];
% Drop first 12 observations in case of NaN %
for ss = Ns:Ne
      aa = find(Rcumvol(:,13)==ss);
      Rcumvol(aa(1):aa(12),:) = [];
      bb = find(Return(:,6)==ss);
      Return(bb(1):bb(12),:) = [];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                       %            1       2       3   4    5         6-17
z0 = [z0, Rcumvol(:,1:12)];                           %  z0 =  [ym, xret, rf, id, vol, xcum/vol]

% Run predictive regression %
res = []; res_sign = [];
for h = 1:12
      yy = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % [ym, xret, rf, id, vol, xcum/vol]                                                                           
      yy = [yy; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % yy = [y, x,  ym, ID]; 
    end
      b = clusterreg(yy(:,1), [yy(:,2), ones(size(yy(:,2)))], yy(:,3));
      res = [res; b(1, [1, 3]), b(2, [1, 3])];      % res = [b, b_tstat, cons, cons_tstat]
      bsign = clusterreg(yy(:,1), [sign(yy(:,2)), ones(size(yy(:,2)))], yy(:,3));
      res_sign = [res_sign; bsign(1, [1, 3]), bsign(2, [1, 3])];
end

% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),2); Exhat_sign = nan(length(z0),2);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = [z0(i,5+h), 1]*[res(h,1), res(h,3)]';
        Exhat_sign(i+1,1) = [sign(z0(i,5+h)), 1]*[res_sign(h,1), res_sign(h,3)]';
    end
    for i = 2:length(z0)
        Exhat(i,2) = z0(i,6) - Exhat(i,1);                            
        Exhat_sign(i,2) = z0(i,6) - Exhat_sign(i,1);
    end
Pedhat = [Pedhat, Exhat(:,1)];
Pedhat_sign = [Pedhat_sign, Exhat(:,1)];
Reshat = [Reshat, Exhat(:,2)];
Reshat_sign = [Reshat_sign, Exhat_sign(:,2)];
end
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_sign, z0(:,4)];
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
for s = Ns:Ne
    line1 = find(Reshat(:,13)==s);
    Reshat(line1(1),:) = [];
    line2 = find(Reshat_sign(:,13)==s);
    Reshat_sign(line2(1),:) = [];
    line3 = find(Pedhat(:,13)==s);
    Pedhat(line3(1),:) = [];
    line4 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line4(1),:) = [];
    line5 = find(z0(:,4)==s);
    z0(line5(1),:) = [];
end

% Wild bootstrap for 1000 times %
B = 1000; 
beta_all = nan(B,12); tstat_all = nan(B,12);
beta_sign_all = nan(B,12); tstat_sign_all = nan(B,12);
for i = 1:B
    tic
    Boot_ret = [];y_boot1 = []; y_boot2 = [];
    for s = Ns:Ne
          rng(s*i);                                                  % set seed
          tY = z0(z0(:,4)==s,1);                           % z0 = [ym, xret, rf, id, vol, xcum/vol]    
          T = length(tY);
          
          Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
          Pedhat_sign_s = Pedhat_sign(Pedhat_sign(:,13)==s,:);
          Reshat_s = Reshat(Reshat(:,13)==s,:);
          Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);

          rand_b = [rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1)];
          
          y_boot_all = nan(T,12);
          y_boot_all_sign = nan(T,12);
          % Construct pseudo sample path %
          for t = 1:T
            vt = ((rand_b(t+100,:)>=0.5)*(1) + (rand_b(t+100,:)<0.5)*(-1));
            y_boot_all(t,:) = Pedhat_s(t,1:12) + Reshat_s(t,1:12).*vt;                        
            y_boot_all_sign(t,:) = Pedhat_sign_s(t,1:12) + Reshat_sign_s(t,1:12).*vt;     
          end
          y_boot1 = [y_boot1; y_boot_all];
          y_boot2 = [y_boot2; y_boot_all_sign];
   end
   z1 = [z0(:,1), z0(:,4), y_boot1];
   z2 = [z0(:,1), z0(:,4), y_boot2];
     
   res_i = []; res_sign_i = [];
  for h = 1:12
       yy1 = []; yy2 = [];
       for n = Ns:Ne 
       tt = find(z1(:,2) == n);                             %    1         2         3-14       
       a = z1(tt,:);                                                %   [ym      ID       cum/vol]  
            
       ta = find(z0(:,4)==n);                              %     1          2          3         4         5         6                                         
       d = z0(ta,:);                                               %   [ym       xret      rf        ID     vol   xcum/vol]
       yy1 = [yy1; a(2:end,2+h), d(1:end-1,5+h), a(2:end,[1, 2])];  %  use original independent variables and bootstrapped  dependent variables 
            
       td = find(z2(:,2)==n);
       c = z2(td,:);
       yy2 = [yy2; c(2:end,2+h), d(1:end-1,5+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
       end
       b = clusterreg(yy1(:,1), [yy1(:,2), ones(size(yy1(:,2)))], yy1(:,3));
       res_i = [res_i; b(1, [1, 3]), b(2, [1, 3])];
  
       bsign = clusterreg(yy2(:,1), [sign(yy2(:,2)), ones(size(yy2(:,2)))], yy2(:,3));
       res_sign_i = [res_sign_i; bsign(1, [1, 3]), bsign(2, [1, 3])];
   end
      beta_all(i,:) = res_i(:,1)';
      tstat_all(i,:) = res_i(:,2)';
      beta_sign_all(i,:) = res_sign_i(:,1)';
      tstat_sign_all(i,:) = res_sign_i(:,2)';
      toc
      i
end

save('E:\RESEARCH\TSMOM\Codes_new\Tab4_WildEqt.mat','beta_all', 'beta_sign_all', 'tstat_all', 'tstat_sign_all')

%% Bond %%
clear;clc
addpath 'E:\RESEARCH\TSMOM\Codes_new';
addpath 'E:\Dropbox\My toolbox';
load Return.mat;

Ns = 34; Ne = 46; 
Return = Return(Return(:,1)>=198501,:);
Return = Return(Return(:,6)>=Ns & Return(:,6)<=Ne,:);

Rs = Return(:,2);
rf = Return(:,4)+1;
exVol = Return(:,3);
Rcum = Func_GenCum(rf, Rs);
Rcumvol = Rcum(2:end,:)./exVol(1:end-1);
Rcumvol = [nan(1,12); Rcumvol];
Rcumvol = [Rcumvol, Return(:,6)];
% Drop first 12 observations in case of NaN %
for ss = Ns:Ne
      aa = find(Rcumvol(:,13)==ss);
      Rcumvol(aa(1):aa(12),:) = [];
      bb = find(Return(:,6)==ss);
      Return(bb(1):bb(12),:) = [];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                       %            1       2       3   4    5         6-17
z0 = [z0, Rcumvol(:,1:12)];                           %  z0 =  [ym, xret, rf, id, vol, xcum/vol]

% Run predictive regression %
res = []; res_sign = [];
for h = 1:12
      yy = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % [ym, xret, rf, id, vol, xcum/vol]                                                                           
      yy = [yy; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % yy = [y, x,  ym, ID]; 
    end
      b = clusterreg(yy(:,1), [yy(:,2), ones(size(yy(:,2)))], yy(:,3));
      res = [res; b(1, [1, 3]), b(2, [1, 3])];      % res = [b, b_tstat, cons, cons_tstat]
      bsign = clusterreg(yy(:,1), [sign(yy(:,2)), ones(size(yy(:,2)))], yy(:,3));
      res_sign = [res_sign; bsign(1, [1, 3]), bsign(2, [1, 3])];
end

% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),2); Exhat_sign = nan(length(z0),2);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = [z0(i,5+h), 1]*[res(h,1), res(h,3)]';
        Exhat_sign(i+1,1) = [sign(z0(i,5+h)), 1]*[res_sign(h,1), res_sign(h,3)]';
    end
    for i = 2:length(z0)
        Exhat(i,2) = z0(i,6) - Exhat(i,1);                            
        Exhat_sign(i,2) = z0(i,6) - Exhat_sign(i,1);
    end
Pedhat = [Pedhat, Exhat(:,1)];
Pedhat_sign = [Pedhat_sign, Exhat(:,1)];
Reshat = [Reshat, Exhat(:,2)];
Reshat_sign = [Reshat_sign, Exhat_sign(:,2)];
end
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_sign, z0(:,4)];
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
for s = Ns:Ne
    line1 = find(Reshat(:,13)==s);
    Reshat(line1(1),:) = [];
    line2 = find(Reshat_sign(:,13)==s);
    Reshat_sign(line2(1),:) = [];
    line3 = find(Pedhat(:,13)==s);
    Pedhat(line3(1),:) = [];
    line4 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line4(1),:) = [];
    line5 = find(z0(:,4)==s);
    z0(line5(1),:) = [];
end

% Wild bootstrap for 1000 times %
B = 1000; 
beta_all = nan(B,12); tstat_all = nan(B,12);
beta_sign_all = nan(B,12); tstat_sign_all = nan(B,12);
for i = 1:B
    tic
    Boot_ret = [];y_boot1 = []; y_boot2 = [];
    for s = Ns:Ne
          rng(s*i);                                                  % set seed
          tY = z0(z0(:,4)==s,1);                           % z0 = [ym, xret, rf, id, vol, xcum/vol]    
          T = length(tY);
          
          Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
          Pedhat_sign_s = Pedhat_sign(Pedhat_sign(:,13)==s,:);
          Reshat_s = Reshat(Reshat(:,13)==s,:);
          Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);

          rand_b = [rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1)];
          
          y_boot_all = nan(T,12);
          y_boot_all_sign = nan(T,12);
          % Construct pseudo sample path %
          for t = 1:T
            vt = ((rand_b(t+100,:)>=0.5)*(1) + (rand_b(t+100,:)<0.5)*(-1));
            y_boot_all(t,:) = Pedhat_s(t,1:12) + Reshat_s(t,1:12).*vt;                        
            y_boot_all_sign(t,:) = Pedhat_sign_s(t,1:12) + Reshat_sign_s(t,1:12).*vt;     
          end
          y_boot1 = [y_boot1; y_boot_all];
          y_boot2 = [y_boot2; y_boot_all_sign];
   end
   z1 = [z0(:,1), z0(:,4), y_boot1];
   z2 = [z0(:,1), z0(:,4), y_boot2];
     
   res_i = []; res_sign_i = [];
  for h = 1:12
       yy1 = []; yy2 = [];
       for n = Ns:Ne 
       tt = find(z1(:,2) == n);                             %    1         2         3-14       
       a = z1(tt,:);                                                %   [ym      ID       cum/vol]  
            
       ta = find(z0(:,4)==n);                              %     1          2          3         4         5         6                                         
       d = z0(ta,:);                                               %   [ym       xret      rf        ID     vol   xcum/vol]
       yy1 = [yy1; a(2:end,2+h), d(1:end-1,5+h), a(2:end,[1, 2])];  %  use original independent variables and bootstrapped  dependent variables 
            
       td = find(z2(:,2)==n);
       c = z2(td,:);
       yy2 = [yy2; c(2:end,2+h), d(1:end-1,5+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
       end
       b = clusterreg(yy1(:,1), [yy1(:,2), ones(size(yy1(:,2)))], yy1(:,3));
       res_i = [res_i; b(1, [1, 3]), b(2, [1, 3])];
  
       bsign = clusterreg(yy2(:,1), [sign(yy2(:,2)), ones(size(yy2(:,2)))], yy2(:,3));
       res_sign_i = [res_sign_i; bsign(1, [1, 3]), bsign(2, [1, 3])];
   end
      beta_all(i,:) = res_i(:,1)';
      tstat_all(i,:) = res_i(:,2)';
      beta_sign_all(i,:) = res_sign_i(:,1)';
      tstat_sign_all(i,:) = res_sign_i(:,2)';
      toc
      i
end

save('E:\RESEARCH\TSMOM\Codes_new\Tab4_WildBnd.mat','beta_all', 'beta_sign_all', 'tstat_all', 'tstat_sign_all')

%% Currency %%
clear;clc
addpath 'E:\RESEARCH\TSMOM\Codes_new';
addpath 'E:\Dropbox\My toolbox';
load Return.mat;

Ns = 47; Ne = 55; 
Return = Return(Return(:,1)>=198501,:);
Return = Return(Return(:,6)>=Ns & Return(:,6)<=Ne,:);

Rs = Return(:,2);
rf = Return(:,4)+1;
exVol = Return(:,3);
Rcum = Func_GenCum(rf, Rs);
Rcumvol = Rcum(2:end,:)./exVol(1:end-1);
Rcumvol = [nan(1,12); Rcumvol];
Rcumvol = [Rcumvol, Return(:,6)];
% Drop first 12 observations in case of NaN %
for ss = Ns:Ne
      aa = find(Rcumvol(:,13)==ss);
      Rcumvol(aa(1):aa(12),:) = [];
      bb = find(Return(:,6)==ss);
      Return(bb(1):bb(12),:) = [];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                       %            1       2       3   4    5         6-17
z0 = [z0, Rcumvol(:,1:12)];                           %  z0 =  [ym, xret, rf, id, vol, xcum/vol]

% Run predictive regression %
res = []; res_sign = [];
for h = 1:12
      yy = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % [ym, xret, rf, id, vol, xcum/vol]                                                                           
      yy = [yy; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % yy = [y, x,  ym, ID]; 
    end
      b = clusterreg(yy(:,1), [yy(:,2), ones(size(yy(:,2)))], yy(:,3));
      res = [res; b(1, [1, 3]), b(2, [1, 3])];      % res = [b, b_tstat, cons, cons_tstat]
      bsign = clusterreg(yy(:,1), [sign(yy(:,2)), ones(size(yy(:,2)))], yy(:,3));
      res_sign = [res_sign; bsign(1, [1, 3]), bsign(2, [1, 3])];
end

% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),2); Exhat_sign = nan(length(z0),2);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = [z0(i,5+h), 1]*[res(h,1), res(h,3)]';
        Exhat_sign(i+1,1) = [sign(z0(i,5+h)), 1]*[res_sign(h,1), res_sign(h,3)]';
    end
    for i = 2:length(z0)
        Exhat(i,2) = z0(i,6) - Exhat(i,1);                            
        Exhat_sign(i,2) = z0(i,6) - Exhat_sign(i,1);
    end
Pedhat = [Pedhat, Exhat(:,1)];
Pedhat_sign = [Pedhat_sign, Exhat(:,1)];
Reshat = [Reshat, Exhat(:,2)];
Reshat_sign = [Reshat_sign, Exhat_sign(:,2)];
end
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_sign, z0(:,4)];
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
for s = Ns:Ne
    line1 = find(Reshat(:,13)==s);
    Reshat(line1(1),:) = [];
    line2 = find(Reshat_sign(:,13)==s);
    Reshat_sign(line2(1),:) = [];
    line3 = find(Pedhat(:,13)==s);
    Pedhat(line3(1),:) = [];
    line4 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line4(1),:) = [];
    line5 = find(z0(:,4)==s);
    z0(line5(1),:) = [];
end

% Wild bootstrap for 1000 times %
B = 1000; 
beta_all = nan(B,12); tstat_all = nan(B,12);
beta_sign_all = nan(B,12); tstat_sign_all = nan(B,12);
for i = 1:B
    tic
    Boot_ret = [];y_boot1 = []; y_boot2 = [];
    for s = Ns:Ne
          rng(s*i);                                                  % set seed
          tY = z0(z0(:,4)==s,1);                           % z0 = [ym, xret, rf, id, vol, xcum/vol]    
          T = length(tY);
          
          Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
          Pedhat_sign_s = Pedhat_sign(Pedhat_sign(:,13)==s,:);
          Reshat_s = Reshat(Reshat(:,13)==s,:);
          Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);

          rand_b = [rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1), ...
              rand(T+100,1), rand(T+100,1), rand(T+100,1), rand(T+100,1)];
          
          y_boot_all = nan(T,12);
          y_boot_all_sign = nan(T,12);
          % Construct pseudo sample path %
          for t = 1:T
            vt = ((rand_b(t+100,:)>=0.5)*(1) + (rand_b(t+100,:)<0.5)*(-1));
            y_boot_all(t,:) = Pedhat_s(t,1:12) + Reshat_s(t,1:12).*vt;                        
            y_boot_all_sign(t,:) = Pedhat_sign_s(t,1:12) + Reshat_sign_s(t,1:12).*vt;     
          end
          y_boot1 = [y_boot1; y_boot_all];
          y_boot2 = [y_boot2; y_boot_all_sign];
   end
   z1 = [z0(:,1), z0(:,4), y_boot1];
   z2 = [z0(:,1), z0(:,4), y_boot2];
     
   res_i = []; res_sign_i = [];
  for h = 1:12
       yy1 = []; yy2 = [];
       for n = Ns:Ne 
       tt = find(z1(:,2) == n);                             %    1         2         3-14       
       a = z1(tt,:);                                                %   [ym      ID       cum/vol]  
            
       ta = find(z0(:,4)==n);                              %     1          2          3         4         5         6                                         
       d = z0(ta,:);                                               %   [ym       xret      rf        ID     vol   xcum/vol]
       yy1 = [yy1; a(2:end,2+h), d(1:end-1,5+h), a(2:end,[1, 2])];  %  use original independent variables and bootstrapped  dependent variables 
            
       td = find(z2(:,2)==n);
       c = z2(td,:);
       yy2 = [yy2; c(2:end,2+h), d(1:end-1,5+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
       end
       b = clusterreg(yy1(:,1), [yy1(:,2), ones(size(yy1(:,2)))], yy1(:,3));
       res_i = [res_i; b(1, [1, 3]), b(2, [1, 3])];
  
       bsign = clusterreg(yy2(:,1), [sign(yy2(:,2)), ones(size(yy2(:,2)))], yy2(:,3));
       res_sign_i = [res_sign_i; bsign(1, [1, 3]), bsign(2, [1, 3])];
   end
      beta_all(i,:) = res_i(:,1)';
      tstat_all(i,:) = res_i(:,2)';
      beta_sign_all(i,:) = res_sign_i(:,1)';
      tstat_sign_all(i,:) = res_sign_i(:,2)';
      toc
      i
end

save('E:\RESEARCH\TSMOM\Codes_new\Tab4_WildCur.mat','beta_all', 'beta_sign_all', 'tstat_all', 'tstat_sign_all')