%% Fixed-design bootstrap with volatility scaling (Table 6, 1985 - 2009) %%

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

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

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

% Generate demeaned vol-scaled returns %
Retvol_dm = []; Retvol_sign_dm = [];
for nn = Ns:Ne
    Retvol_n = Return(Return(:,6)==nn,7);
    Retvol_m = Retvol_n - mean(Retvol_n);
    Retvol_dm = [Retvol_dm; Retvol_m];
    
    Retvol_sign = sign(Retvol_n) - mean(sign(Retvol_n));
    Retvol_sign_dm = [Retvol_sign_dm; Retvol_sign];
end
z0 = Return(:, [1, 2, 4, 6, 3]);                                                                                      %            1      2       3   4    5         6                 7              8                           9
z0 = [z0, Retvol_dm, Return(:,7), Retvol_sign_dm, sign(Return(:,7))];            % z0 =  [ym, xret, rf, id, vol, dm-xr/vol, xr/vol, dm-sign-xr/vol, sign-xr/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-xr/vol, xr/vol, dm-sign-xr/vol, sign-xr/vol]                                                                          
      yy1 = [yy1; a(h+1:end,6), a(1:end-h,6), a(h+1:end,[1, 4])];    
      
      yy2 = [yy2; a(h+1:end,6), a(1:end-h,8), a(h+1: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,7);                            % xr/vol
        ybar = mean(y1(h+1:end));
        x_n = z0(z0(:,4)==n,7);                          % undemeaned x to derive alpha
        xbar_n = mean(x_n(1:end-h));
        alpha_n = ybar - res(h,1)*xbar_n;
        alpha_h = [alpha_h; alpha_n];
        
        xsign_n = z0(z0(:,4)==n,9);                  % undemeaned xsign to derive alpha_sign
        xbar_sign_s = mean(xsign_n(1:end-h));
        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                                                      % z0 =  [ym, xret, rf, id, vol, dm-xr/vol, xr/vol, dm-sign-xr/vol, sign-xr/vol]
    Exhat = nan(length(z0),3); Exhat_sign = nan(length(z0),3);
    for i = 1:length(z0)-h
        Exhat(i+h,1) = z0(i,6)*res(h,1);           % demeaned x 
        Exhat_sign(i+h,1) = z0(i,8)*res_sign(h,1);
        Exhat(i+h,2) = z0(i,7)*res(h,1);           % undemeaned x for beta_FE*x 
        Exhat_sign(i+h,2) = z0(i,9)*res_sign(h,1);
    end
    for i = 1:length(z0)-h
        Exhat(i+h,3) = z0(i+h,6) - Exhat(i+h,1);      % residuals epsilon_hat                     
        Exhat_sign(i+h,3) = z0(i+h,6) - Exhat_sign(i+h,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-xr/vol, xr/vol, dm-sign-xr/vol, sign-xr/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);
        y_dm1 = nan(T,12);
        y_dm2 = nan(T,12);
        % Construct pseudo sample path %
        for h1 = 1:12
            for t = h1:T
            vt = ((rand_b(t+100,:)>=0.5)*(1) + (rand_b(t+100,:)<0.5)*(-1));
            y_boot_all(t,h1) = alpha_all(s,h1) + Pedhat_s(t, h1) + Reshat_s(t,h1)*vt(:,h1);                    
            y_boot_all_sign(t,h1) = alpha_sign_all(s,h1) + Pedhat_sign_s(t, h1)+ Reshat_sign_s(t,h1)*vt(:,h1);      
            end
        % Generate demeaned pseudo samples %
        y_dm1(:,h1) = y_boot_all(:,h1) - mean(y_boot_all(h1:end, h1));
        y_dm2(:,h1) = y_boot_all_sign(:,h1) - mean(y_boot_all_sign(h1:end, h1));
        end
        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               7                 8                     9
          d = z0(ta,:);                                            % z0 = [ym, xret, rf, id,  vol, dm-xr/vol, xr/vol, dm-sign-xr/vol, sign-xr/vol]
          yy1 = [yy1; a(1+h:end,2+h), d(1:end-h,6), a(h+1:end,[1, 2])];      %  use original independent variables and bootstrapped  dependent variables
            
          td = find(z2(:,2)==n);
          c = z2(td,:);
          yy2 = [yy2; c(1+h:end,2+h), d(1:end-h,8), c(h+1: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\Tab6_FixVol2009.mat','beta_all', 'beta_sign_all', 'tstat_all', 'tstat_sign_all')

