%% Fixed-design 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

% Generate demeaned vol-scaled cumulative returns %
Rcumvol_dm = []; Rcumvol_sign_dm = [];
for m = Ns:Ne
    Rcumvol_m = Rcumvol(Rcumvol(:,13)==m,1:12);
    Rcumvol_n = Rcumvol_m-mean(Rcumvol_m);
    Rcumvol_dm = [Rcumvol_dm; Rcumvol_n];
    
    Rcumvol_sign = sign(Rcumvol_m) - mean(sign(Rcumvol_m));
    Rcumvol_sign_dm = [Rcumvol_sign_dm; Rcumvol_sign];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                                                                                                          %            1     2       3   4    5         6-17                  18-29           30-41                    42-53
z0 = [z0, Rcumvol_dm, Rcumvol(:,1:12), Rcumvol_sign_dm, sign(Rcumvol(:,1:12))];      % z0 =  [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol, sign-dm-xcum/vol, sign-xcum/vol]

% Run predictive pool regressions with demeaned values %
res = []; res_sign = [];
for h = 1:12
      yy1 = []; yy2 = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % a = [ym, xret, rf, id, vol, dm-xcum/vol]                                                                           
      yy1 = [yy1; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % demeaned regressors and dep vars
      
      yy2 = [yy2; a(2:end,6), a(1:end-1, 29+h), a(2:end, [1,4])];
    end
      b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
      res = [res; b(1, [1, 3])];      
      bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
      res_sign = [res_sign; bsign(1, [1, 3])];
end

% Back out fixed effect alpha, alpha_sign %
alpha_all = []; alpha_sign_all = [];
for h = 1:12
    alpha_h = []; alpha_sign_h = [];
    for n = Ns:Ne
        y1 = z0(z0(:,4)==n,18);
        ybar = mean(y1(2:end));
        x_n = z0(z0(:,4)==n,17+h);                   % undemeaned x to derive alpha
        xbar_n = mean(x_n(1:end-1));
        alpha_n = ybar - res(h,1)*xbar_n;
        alpha_h = [alpha_h; alpha_n];
        
        xsign_n = z0(z0(:,4)==n,41+h);           % undemeaned xsign to derive alphasign
        xbar_sign_s = mean(xsign_n(1:end-1));
        alpha_sign_s = ybar - res_sign(h,1)*xbar_sign_s;
        alpha_sign_h = [alpha_sign_h; alpha_sign_s];
    end 
    alpha_all = [alpha_all, alpha_h];
    alpha_sign_all = [alpha_sign_all, alpha_sign_h];
end
% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),3); Exhat_sign = nan(length(z0),3);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = z0(i,5+h)*res(h,1);      % demeaned x 
        Exhat_sign(i+1,1) = z0(i,29+h)*res_sign(h,1);
        Exhat(i+1,2) = z0(i,17+h)*res(h,1);    % undemeaned x for beta_FE*x 
        Exhat_sign(i+1,2) = z0(i,41+h)*res_sign(h,1);
    end
    for i = 2:length(z0)
        Exhat(i,3) = z0(i,6) - Exhat(i,1);           % residual epsilon_hat                     
        Exhat_sign(i,3) = z0(i,6) - Exhat_sign(i,1);
    end
% Exhat = [beta_FE*dmx, beta_FE*undmx, residuals]
Reshat = [Reshat, Exhat(:,3)];                      
Reshat_sign = [Reshat_sign, Exhat_sign(:,3)];
Pedhat = [Pedhat, Exhat(:,2)];
Pedhat_sign = [Pedhat_sign, Exhat_sign(:,2)];
end
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_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(z0(:,4)==s);
    z0(line3(1),:) = [];
    line4 = find(Pedhat(:,13)==s);
    Pedhat(line4(1),:) = [];
    line5 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line5(1),:) = [];  
end

% Fixed-design 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);
        tY = z0(z0(:,4)==s,1);                             % z0 = [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol]    
        T = length(tY);
        
        Reshat_s = Reshat(Reshat(:,13)==s,:);
        Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);
        Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
        Pedhat_sign_s = Pedhat_sign(Pedhat_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,:) = alpha_all(s,:) + Pedhat_s(t, 1:12) + Reshat_s(t,1:12).*vt;                    
            y_boot_all_sign(t,:) = alpha_sign_all(s,:) + Pedhat_sign_s(t, 1:12)+ Reshat_sign_s(t,1:12).*vt;    
        end
        % Generate demeaned pseudo samples %
        y_dm1 = y_boot_all - mean(y_boot_all);
        y_dm2 = y_boot_all_sign - mean(y_boot_all_sign);
        y_boot1 = [y_boot1; y_dm1];
        y_boot2 = [y_boot2; y_dm2];
    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     dm-boot_y]  
            
          ta = find(z0(:,4)==n);                           %     1          2          3         4         5         6-17                                               
          d = z0(ta,:);                                            %  [ym       xret      rf        ID     Vol  dm-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,29+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
        end
        b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
        res_i = [res_i; b(1, [1, 3])];
  
        bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
        res_sign_i = [res_sign_i; bsign(1, [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_FixCom.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,:);
Return(:,6)=Return(:,6)-24;
Ns = 1; Ne = 9;

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

% Generate demeaned vol-scaled cumulative returns %
Rcumvol_dm = []; Rcumvol_sign_dm = [];
for m = Ns:Ne
    Rcumvol_m = Rcumvol(Rcumvol(:,13)==m,1:12);
    Rcumvol_n = Rcumvol_m-mean(Rcumvol_m);
    Rcumvol_dm = [Rcumvol_dm; Rcumvol_n];
    
    Rcumvol_sign = sign(Rcumvol_m) - mean(sign(Rcumvol_m));
    Rcumvol_sign_dm = [Rcumvol_sign_dm; Rcumvol_sign];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                                                                                                          %            1     2       3   4    5         6-17                  18-29           30-41                    42-53
z0 = [z0, Rcumvol_dm, Rcumvol(:,1:12), Rcumvol_sign_dm, sign(Rcumvol(:,1:12))];      % z0 =  [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol, sign-dm-xcum/vol, sign-xcum/vol]

% Run predictive pool regressions with demeaned values %
res = []; res_sign = [];
for h = 1:12
      yy1 = []; yy2 = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % a = [ym, xret, rf, id, vol, dm-xcum/vol]                                                                           
      yy1 = [yy1; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % demeaned regressors and dep vars
      
      yy2 = [yy2; a(2:end,6), a(1:end-1, 29+h), a(2:end, [1,4])];
    end
      b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
      res = [res; b(1, [1, 3])];      
      bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
      res_sign = [res_sign; bsign(1, [1, 3])];
end

% Back out fixed effect alpha, alpha_sign %
alpha_all = []; alpha_sign_all = [];
for h = 1:12
    alpha_h = []; alpha_sign_h = [];
    for n = Ns:Ne
        y1 = z0(z0(:,4)==n,18);
        ybar = mean(y1(2:end));
        x_n = z0(z0(:,4)==n,17+h);                   % undemeaned x to derive alpha
        xbar_n = mean(x_n(1:end-1));
        alpha_n = ybar - res(h,1)*xbar_n;
        alpha_h = [alpha_h; alpha_n];
        
        xsign_n = z0(z0(:,4)==n,41+h);           % undemeaned xsign to derive alphasign
        xbar_sign_s = mean(xsign_n(1:end-1));
        alpha_sign_s = ybar - res_sign(h,1)*xbar_sign_s;
        alpha_sign_h = [alpha_sign_h; alpha_sign_s];
    end 
    alpha_all = [alpha_all, alpha_h];
    alpha_sign_all = [alpha_sign_all, alpha_sign_h];
end
% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),3); Exhat_sign = nan(length(z0),3);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = z0(i,5+h)*res(h,1);      % demeaned x 
        Exhat_sign(i+1,1) = z0(i,29+h)*res_sign(h,1);
        Exhat(i+1,2) = z0(i,17+h)*res(h,1);    % undemeaned x for beta_FE*x 
        Exhat_sign(i+1,2) = z0(i,41+h)*res_sign(h,1);
    end
    for i = 2:length(z0)
        Exhat(i,3) = z0(i,6) - Exhat(i,1);           % residual epsilon_hat                     
        Exhat_sign(i,3) = z0(i,6) - Exhat_sign(i,1);
    end
% Exhat = [beta_FE*dmx, beta_FE*undmx, residuals]
Reshat = [Reshat, Exhat(:,3)];                      
Reshat_sign = [Reshat_sign, Exhat_sign(:,3)];
Pedhat = [Pedhat, Exhat(:,2)];
Pedhat_sign = [Pedhat_sign, Exhat_sign(:,2)];
end
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_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(z0(:,4)==s);
    z0(line3(1),:) = [];
    line4 = find(Pedhat(:,13)==s);
    Pedhat(line4(1),:) = [];
    line5 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line5(1),:) = [];  
end

% Fixed-design 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);
        tY = z0(z0(:,4)==s,1);                             % z0 = [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol]    
        T = length(tY);
        
        Reshat_s = Reshat(Reshat(:,13)==s,:);
        Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);
        Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
        Pedhat_sign_s = Pedhat_sign(Pedhat_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,:) = alpha_all(s,:) + Pedhat_s(t, 1:12) + Reshat_s(t,1:12).*vt;                    
            y_boot_all_sign(t,:) = alpha_sign_all(s,:) + Pedhat_sign_s(t, 1:12)+ Reshat_sign_s(t,1:12).*vt;    
        end
        % Generate demeaned pseudo samples %
        y_dm1 = y_boot_all - mean(y_boot_all);
        y_dm2 = y_boot_all_sign - mean(y_boot_all_sign);
        y_boot1 = [y_boot1; y_dm1];
        y_boot2 = [y_boot2; y_dm2];
    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     dm-boot_y]  
            
          ta = find(z0(:,4)==n);                           %     1          2          3         4         5         6-17                                               
          d = z0(ta,:);                                            %  [ym       xret      rf        ID     Vol  dm-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,29+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
        end
        b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
        res_i = [res_i; b(1, [1, 3])];
  
        bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
        res_sign_i = [res_sign_i; bsign(1, [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_FixEqt.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,:);
Return(:,6)=Return(:,6)-33;
Ns = 1; Ne = 13;

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

% Generate demeaned vol-scaled cumulative returns %
Rcumvol_dm = []; Rcumvol_sign_dm = [];
for m = Ns:Ne
    Rcumvol_m = Rcumvol(Rcumvol(:,13)==m,1:12);
    Rcumvol_n = Rcumvol_m-mean(Rcumvol_m);
    Rcumvol_dm = [Rcumvol_dm; Rcumvol_n];
    
    Rcumvol_sign = sign(Rcumvol_m) - mean(sign(Rcumvol_m));
    Rcumvol_sign_dm = [Rcumvol_sign_dm; Rcumvol_sign];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                                                                                                          %            1     2       3   4    5         6-17                  18-29           30-41                    42-53
z0 = [z0, Rcumvol_dm, Rcumvol(:,1:12), Rcumvol_sign_dm, sign(Rcumvol(:,1:12))];      % z0 =  [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol, sign-dm-xcum/vol, sign-xcum/vol]

% Run predictive pool regressions with demeaned values %
res = []; res_sign = [];
for h = 1:12
      yy1 = []; yy2 = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % a = [ym, xret, rf, id, vol, dm-xcum/vol]                                                                           
      yy1 = [yy1; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % demeaned regressors and dep vars
      
      yy2 = [yy2; a(2:end,6), a(1:end-1, 29+h), a(2:end, [1,4])];
    end
      b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
      res = [res; b(1, [1, 3])];      
      bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
      res_sign = [res_sign; bsign(1, [1, 3])];
end

% Back out fixed effect alpha, alpha_sign %
alpha_all = []; alpha_sign_all = [];
for h = 1:12
    alpha_h = []; alpha_sign_h = [];
    for n = Ns:Ne
        y1 = z0(z0(:,4)==n,18);
        ybar = mean(y1(2:end));
        x_n = z0(z0(:,4)==n,17+h);                   % undemeaned x to derive alpha
        xbar_n = mean(x_n(1:end-1));
        alpha_n = ybar - res(h,1)*xbar_n;
        alpha_h = [alpha_h; alpha_n];
        
        xsign_n = z0(z0(:,4)==n,41+h);           % undemeaned xsign to derive alphasign
        xbar_sign_s = mean(xsign_n(1:end-1));
        alpha_sign_s = ybar - res_sign(h,1)*xbar_sign_s;
        alpha_sign_h = [alpha_sign_h; alpha_sign_s];
    end 
    alpha_all = [alpha_all, alpha_h];
    alpha_sign_all = [alpha_sign_all, alpha_sign_h];
end
% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),3); Exhat_sign = nan(length(z0),3);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = z0(i,5+h)*res(h,1);      % demeaned x 
        Exhat_sign(i+1,1) = z0(i,29+h)*res_sign(h,1);
        Exhat(i+1,2) = z0(i,17+h)*res(h,1);    % undemeaned x for beta_FE*x 
        Exhat_sign(i+1,2) = z0(i,41+h)*res_sign(h,1);
    end
    for i = 2:length(z0)
        Exhat(i,3) = z0(i,6) - Exhat(i,1);           % residual epsilon_hat                     
        Exhat_sign(i,3) = z0(i,6) - Exhat_sign(i,1);
    end
% Exhat = [beta_FE*dmx, beta_FE*undmx, residuals]
Reshat = [Reshat, Exhat(:,3)];                      
Reshat_sign = [Reshat_sign, Exhat_sign(:,3)];
Pedhat = [Pedhat, Exhat(:,2)];
Pedhat_sign = [Pedhat_sign, Exhat_sign(:,2)];
end
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_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(z0(:,4)==s);
    z0(line3(1),:) = [];
    line4 = find(Pedhat(:,13)==s);
    Pedhat(line4(1),:) = [];
    line5 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line5(1),:) = [];  
end

% Fixed-design 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);
        tY = z0(z0(:,4)==s,1);                             % z0 = [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol]    
        T = length(tY);
        
        Reshat_s = Reshat(Reshat(:,13)==s,:);
        Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);
        Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
        Pedhat_sign_s = Pedhat_sign(Pedhat_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,:) = alpha_all(s,:) + Pedhat_s(t, 1:12) + Reshat_s(t,1:12).*vt;                    
            y_boot_all_sign(t,:) = alpha_sign_all(s,:) + Pedhat_sign_s(t, 1:12)+ Reshat_sign_s(t,1:12).*vt;    
        end
        % Generate demeaned pseudo samples %
        y_dm1 = y_boot_all - mean(y_boot_all);
        y_dm2 = y_boot_all_sign - mean(y_boot_all_sign);
        y_boot1 = [y_boot1; y_dm1];
        y_boot2 = [y_boot2; y_dm2];
    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     dm-boot_y]  
            
          ta = find(z0(:,4)==n);                           %     1          2          3         4         5         6-17                                               
          d = z0(ta,:);                                            %  [ym       xret      rf        ID     Vol  dm-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,29+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
        end
        b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
        res_i = [res_i; b(1, [1, 3])];
  
        bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
        res_sign_i = [res_sign_i; bsign(1, [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_FixBnd.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,:);
Return(:,6)=Return(:,6)-46;
Ns = 1; Ne = 9;

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

% Generate demeaned vol-scaled cumulative returns %
Rcumvol_dm = []; Rcumvol_sign_dm = [];
for m = Ns:Ne
    Rcumvol_m = Rcumvol(Rcumvol(:,13)==m,1:12);
    Rcumvol_n = Rcumvol_m-mean(Rcumvol_m);
    Rcumvol_dm = [Rcumvol_dm; Rcumvol_n];
    
    Rcumvol_sign = sign(Rcumvol_m) - mean(sign(Rcumvol_m));
    Rcumvol_sign_dm = [Rcumvol_sign_dm; Rcumvol_sign];
end

z0 = Return(:, [1, 2, 4, 6, 3]);                                                                                                          %            1     2       3   4    5         6-17                  18-29           30-41                    42-53
z0 = [z0, Rcumvol_dm, Rcumvol(:,1:12), Rcumvol_sign_dm, sign(Rcumvol(:,1:12))];      % z0 =  [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol, sign-dm-xcum/vol, sign-xcum/vol]

% Run predictive pool regressions with demeaned values %
res = []; res_sign = [];
for h = 1:12
      yy1 = []; yy2 = [];
    for n = Ns:Ne 
      tt = find(z0(:,4) == n);
      a = z0(tt,:);                                                 % a = [ym, xret, rf, id, vol, dm-xcum/vol]                                                                           
      yy1 = [yy1; a(2:end,6), a(1:end-1,5+h), a(2:end,[1, 4])];    % demeaned regressors and dep vars
      
      yy2 = [yy2; a(2:end,6), a(1:end-1, 29+h), a(2:end, [1,4])];
    end
      b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
      res = [res; b(1, [1, 3])];      
      bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
      res_sign = [res_sign; bsign(1, [1, 3])];
end

% Back out fixed effect alpha, alpha_sign %
alpha_all = []; alpha_sign_all = [];
for h = 1:12
    alpha_h = []; alpha_sign_h = [];
    for n = Ns:Ne
        y1 = z0(z0(:,4)==n,18);
        ybar = mean(y1(2:end));
        x_n = z0(z0(:,4)==n,17+h);                   % undemeaned x to derive alpha
        xbar_n = mean(x_n(1:end-1));
        alpha_n = ybar - res(h,1)*xbar_n;
        alpha_h = [alpha_h; alpha_n];
        
        xsign_n = z0(z0(:,4)==n,41+h);           % undemeaned xsign to derive alphasign
        xbar_sign_s = mean(xsign_n(1:end-1));
        alpha_sign_s = ybar - res_sign(h,1)*xbar_sign_s;
        alpha_sign_h = [alpha_sign_h; alpha_sign_s];
    end 
    alpha_all = [alpha_all, alpha_h];
    alpha_sign_all = [alpha_sign_all, alpha_sign_h];
end
% Generate yhat and epsilonhat %
Reshat = []; Reshat_sign = []; Pedhat = []; Pedhat_sign = [];
for h = 1:12
    Exhat = nan(length(z0),3); Exhat_sign = nan(length(z0),3);
    for i = 1:length(z0)-1
        Exhat(i+1,1) = z0(i,5+h)*res(h,1);      % demeaned x 
        Exhat_sign(i+1,1) = z0(i,29+h)*res_sign(h,1);
        Exhat(i+1,2) = z0(i,17+h)*res(h,1);    % undemeaned x for beta_FE*x 
        Exhat_sign(i+1,2) = z0(i,41+h)*res_sign(h,1);
    end
    for i = 2:length(z0)
        Exhat(i,3) = z0(i,6) - Exhat(i,1);           % residual epsilon_hat                     
        Exhat_sign(i,3) = z0(i,6) - Exhat_sign(i,1);
    end
% Exhat = [beta_FE*dmx, beta_FE*undmx, residuals]
Reshat = [Reshat, Exhat(:,3)];                      
Reshat_sign = [Reshat_sign, Exhat_sign(:,3)];
Pedhat = [Pedhat, Exhat(:,2)];
Pedhat_sign = [Pedhat_sign, Exhat_sign(:,2)];
end
Reshat = [Reshat, z0(:,4)];
Reshat_sign = [Reshat_sign, z0(:,4)];
Pedhat = [Pedhat, z0(:,4)];
Pedhat_sign = [Pedhat_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(z0(:,4)==s);
    z0(line3(1),:) = [];
    line4 = find(Pedhat(:,13)==s);
    Pedhat(line4(1),:) = [];
    line5 = find(Pedhat_sign(:,13)==s);
    Pedhat_sign(line5(1),:) = [];  
end

% Fixed-design 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);
        tY = z0(z0(:,4)==s,1);                             % z0 = [ym, xret, rf, id, vol, dm-xcum/vol, xcum/vol]    
        T = length(tY);
        
        Reshat_s = Reshat(Reshat(:,13)==s,:);
        Reshat_sign_s = Reshat_sign(Reshat_sign(:,13)==s,:);
        Pedhat_s = Pedhat(Pedhat(:,13)==s,:);
        Pedhat_sign_s = Pedhat_sign(Pedhat_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,:) = alpha_all(s,:) + Pedhat_s(t, 1:12) + Reshat_s(t,1:12).*vt;                    
            y_boot_all_sign(t,:) = alpha_sign_all(s,:) + Pedhat_sign_s(t, 1:12)+ Reshat_sign_s(t,1:12).*vt;    
        end
        % Generate demeaned pseudo samples %
        y_dm1 = y_boot_all - mean(y_boot_all);
        y_dm2 = y_boot_all_sign - mean(y_boot_all_sign);
        y_boot1 = [y_boot1; y_dm1];
        y_boot2 = [y_boot2; y_dm2];
    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     dm-boot_y]  
            
          ta = find(z0(:,4)==n);                           %     1          2          3         4         5         6-17                                               
          d = z0(ta,:);                                            %  [ym       xret      rf        ID     Vol  dm-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,29+h), c(2:end,[1, 2])];   % yy = [y, x,  ym, ID]; 
        end
        b = clusterreg(yy1(:,1), yy1(:,2), yy1(:,3));
        res_i = [res_i; b(1, [1, 3])];
  
        bsign = clusterreg(yy2(:,1), yy2(:,2), yy2(:,3));
        res_sign_i = [res_sign_i; bsign(1, [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_FixCur.mat','beta_all', 'beta_sign_all', 'tstat_all', 'tstat_sign_all')
