% Table 6 Conference Draft
% 
%
% 
clear all;
small = 1.0e-10;
big = 1.0e+6;
this_date = datestr(now,'yyyymmdd');
rng(964530);

% -----  File Directories -- overhead, etc.
 outdir = 'out/';
 figdir = 'fig/';
 matdir = 'mat/';
 rtsdir = '../rts/';
 p_initial = path;
 vardef_dir = 'variable_definitions/';            % Directory for necessary matlab_procs and functions                         
 p1=path(vardef_dir,p_initial);
 utility_dir = 'm_utilities';
 path(utility_dir,p1);

% -- Load Data
load_data=0;         % Set = 1 first time to read data from Excel files (then they are saved as MATLAB file)
fhsw_bpea_data;      % fhsw_bpea_data reads in the full dataset .. all variables, etc. saved in datain.xx

% Get Supply Variables
supply_variables;

% Cyclcical Adjustment Parameters
% Run This Program Twice .. once for 1956-2016 and once for 1981 2016
% --------- RUN 1
ismpl = smpl(calvec,[1956 3],[2016 2],4);  % Sample period for estimating Cyclical adjustment coefficients
outfile_name = [outdir ['Table_6_Panel_A_' this_date '.asc'] ];
fileID = fopen(outfile_name,'w');
% --------- RUN 2 
%ismpl = smpl(calvec,[1981 3],[2016 2],4);  % Sample period for estimating Cyclical adjustment coefficients
%outfile_name = [outdir ['Table_6_Panel_B_' this_date '.asc'] ];
%fileID = fopen(outfile_name,'w');


du_nlead = 2;                   % number of leads
du_nlag = 2;                    % number of lags  
% -- Set up regressors for cylical adjustment
d_urate = dif(cps_urate,1);     % first difference of unemployment rate
xca = lag(d_urate,-du_nlead);   % xca is cyclical adjustment variable
for i = (-du_nlead+1):du_nlag;
 xca = [xca lag(d_urate,i)];
end;

dy = d_bus_tfp;
bw_bw = 60;
[dy_cycle,dy_ca,dy_mu_60,dy_irregular,var_dy_mu_60,b_okun_sum,se_b_okun_sum ] = cyclical_adjustment(dy,xca,ismpl,bw_bw);
x = dy_ca(ismpl==1);
cal = calvec(ismpl==1);

T = size(x,1);
mx = mean(x);
sx = std(x);
xm = x - mx;
acov = zeros(0,5);
for h = 0:4;
 acov(h+1) = xm(1:T-h)'*xm(1+h:T)/T;
end;
acor = acov/acov(1);
% Compute NW Standard Error
nw_cov = acov(1);
for h = 1:4;
    nw_cov = nw_cov + 2*((5-h)/5)*acov(h+1);
end;

% Compute Nyblom Statistic
% Several random draws to get p-value
ndraw = 5000;
Y = randn(T,ndraw);
Ym = Y - repmat(mean(Y),T,1);
Ym_ps = cumsum(Ym);
Ym_ps_sq = Ym_ps.^2;
Qdraw = (sum(Ym_ps_sq)/(T^2))';

% Compute Q for these data
Q = sum(cumsum(xm).^2)/(T^2);
Q = Q/nw_cov;
Q_pvalue = mean(Qdraw>Q);

% Compute LFST using Cosine Transforms
q = 12;
[X_mean,X_cos] = xp_compute(x,q);
[lfst,lfst_pvalue] = lfst_pv(X_cos,10);

% QLR Tests
qlr_1break = NaN;
pv_1break = NaN;
bd_1break = NaN;
qlr_2break = NaN;
pv_2break = NaN;
bd_2break = NaN(1,2);
qlr_3break = NaN;
pv_3break = NaN;
bd_3break = NaN(1,3);


% Single Break
wstat = zeros(T,2);
Tmin = ceil(0.15*T);
for t = 1:T;
    t1 = t;
    t2 = T-t;
    if min([t1;t2]) >= Tmin;
      m1 = mean(x(1:t));
      m2 = mean(x(t+1:T));
      t1 = t;
      t2 = T-t;
      var_hac = nw_cov*(1/t1 + 1/t2);
      wstat(t,1) = ((m1-m2)^2)/var_hac;
      wstat(t,2) = cal(t);
    end;
end;

[qlr imax]=max(wstat(:,1));

% Compute pvalue using simulations
ndraw = 10000;
xdraw = randn(T,ndraw);
wdraw = zeros(T,ndraw);
for t = 1:T;
    t1 = t;
    t2 = T-t;
    if min([t1;t2]) >= Tmin;
      m1 = mean(xdraw(1:t,:));
      m2 = mean(xdraw(t+1:T,:));
      t1 = t;
      t2 = T-t;
      var = (1/t1 + 1/t2);
      wdraw(t,:) = ((m1-m2).^2)/var;
    end;
end;
qlrdraw = max(wdraw)';
pv = mean(qlrdraw > qlr);
qlr_1break = qlr;
pv_1break = pv;
bd_1break = cal(imax);

% 2 breaks 
Tmin = ceil(0.15*T);
R = zeros(2,3);
R(1,1) = 1;
R(1,2) = -1;
R(2,1) = 1;
R(2,3) = -1;
ik = 0;
for t1 = 1:T;
    for t2 = t1+1:T;
      n1 = t1;
      n2 = t2-t1;
      n3 = T-t2;
      if min([n1;n2;n3]) >= Tmin;
       m1 = mean(x(1:t1));
       m2 = mean(x(t1+1:t2));
       m3 = mean(x(t2+1:T));
       m = [m1;m2;m3];
       v = zeros(3,3);
       v(1,1) = 1/n1;
       v(2,2) = 1/n2;
       v(3,3) = 1/n3;
       v = nw_cov*v;
       RvR_inv = inv(R*v*R');
       ch = chol(RvR_inv);
       Rm = R*m;
       CRm = ch*Rm;
       wstat = sum(CRm.^2);
       tmp = [wstat t1 t2];
       ik = ik+1;
       if ik == 1;
             wstat_save = tmp;
       else
             wstat_save = [wstat_save;tmp];
       end
      end;
    end;
end;

[qlr imax]=max(wstat_save(:,1));


ndraw = 1000;
qlrdraw = NaN(ndraw,1);
for idraw = 1:ndraw;
z = randn(T,1);
ik = 0;
for t1 = 1:T;
    for t2 = t1+1:T;
      n1 = t1;
      n2 = t2-t1;
      n3 = T-t2;
      if min([n1;n2;n3]) >= Tmin;
       m1 = mean(z(1:t1));
       m2 = mean(z(t1+1:t2));
       m3 = mean(z(t2+1:T));
       m = [m1;m2;m3];
       v = zeros(3,3);
       v(1,1) = 1/n1;
       v(2,2) = 1/n2;
       v(3,3) = 1/n3;
       RvR_inv = inv(R*v*R');
       ch = chol(RvR_inv);
       Rm = R*m;
       CRm = ch*Rm;
       wstat = sum(CRm.^2);
       ik = ik+1;
       if ik == 1;
             wdraw = wstat;
       else
             wdraw = [wdraw;wstat];
       end
      end;
    end;
end;
qlrdraw(idraw) = max(wdraw);
end;
qlr_2break = qlr;
pv = mean(qlrdraw > qlr);
pv_2break = pv;
bd_2break = [cal(wstat_save(imax,2)) cal(wstat_save(imax,3))];
% 
% % 3 breaks 
Tmin = ceil(0.15*T);
R = zeros(3,4);
R(1,1) = 1;
R(1,2) = -1;
R(2,1) = 1;
R(2,3) = -1;
R(3,1) = 1;
R(3,4) = -1;
ik = 0;
for t1 = 1:T;
    for t2 = t1+1:T;
        for t3 = t2+1:T;
          n1 = t1;
          n2 = t2-t1;
          n3 = t3-t2;
          n4 = T-t3;
          if min([n1;n2;n3;n4]) >= Tmin;
            m1 = mean(x(1:t1));
            m2 = mean(x(t1+1:t2));
            m3 = mean(x(t2+1:t3));
            m4 = mean(x(t3+1:T));
            m = [m1;m2;m3;m4];
            v = zeros(4,4);
            v(1,1) = 1/n1;
            v(2,2) = 1/n2;
            v(3,3) = 1/n3;
            v(4,4) = 1/n4;
            v = nw_cov*v;
            RvR_inv = inv(R*v*R');
            ch = chol(RvR_inv);
            Rm = R*m;
            CRm = ch*Rm;
            wstat = sum(CRm.^2);
            tmp = [wstat t1 t2 t3];
            ik = ik+1;
            if ik == 1;
              wstat_save = tmp;
            else
             wstat_save = [wstat_save;tmp];
            end
          end;
        end;
    end;
end;

[qlr imax]=max(wstat_save(:,1));

ndraw = 1000;
qlrdraw = NaN(ndraw,1);
for idraw = 1:ndraw;
idraw
tic;
z = randn(T,1);
ik = 0;
tic;
for t1 = 1:T;
    for t2 = t1+1:T;
        for t3 = t2+1:T;
          n1 = t1;
          n2 = t2-t1;
          n3 = t3-t2;
          n4 = T-t3;
          if min([n1;n2;n3;n4]) >= Tmin;
            m1 = mean(z(1:t1));
            m2 = mean(z(t1+1:t2));
            m3 = mean(z(t2+1:t3));
            m4 = mean(z(t3+1:T));
            m = [m1;m2;m3;m4];
            v = zeros(4,4);
            v(1,1) = 1/n1;
            v(2,2) = 1/n2;
            v(3,3) = 1/n3;
            v(4,4) = 1/n4;
            RvR_inv = inv(R*v*R');
            ch = chol(RvR_inv);
            Rm = R*m;
            CRm = ch*Rm;
            wstat = sum(CRm.^2);
            ik = ik+1;
            if ik == 1;
              wdraw = wstat;
            else
              wdraw = [wdraw;wstat];
            end
          end;
        end;
    end;
end;
qlrdraw(idraw) = max(wdraw);
toc;
end;

pv = mean(qlrdraw > qlr);
qlr_3break = qlr;
pv_3break = pv;
bd_3break = [cal(wstat_save(imax,2)) cal(wstat_save(imax,3)) cal(wstat_save(imax,4))];

 % Output file
  fprintf(fileID,'Sample Period: %5.2f - %5.2f \n\n',[cal(1) cal(end)]);
  fprintf(fileID,'Newey-West Long-run Standard Deviation: %5.4f\n',sqrt(nw_cov));
  fprintf(fileID,'Nyblom Statistic and Pvalue: %5.4f  %5.4f \n',[Q Q_pvalue]);
  fprintf(fileID,'LFST Statistic and Pvalue: %5.4f  %5.4f \n',[lfst lfst_pvalue]);
  fprintf(fileID,'QLR 1 Break \n');
  fprintf(fileID,'    Stat and Pvalue %5.4f  %5.4f \n',[qlr_1break pv_1break]);
  fprintf(fileID,'    Break date: %5.4f \n',bd_1break);
  fprintf(fileID,'QLR 2 Breaks \n');
  fprintf(fileID,'    Stat and Pvalue %5.4f  %5.4f \n',[qlr_2break pv_2break]);
  fprintf(fileID,'    Break date: %5.4f  %5.4f \n',bd_2break);
  fprintf(fileID,'QLR 3 Breaks \n');
  fprintf(fileID,'    Stat and Pvalue %5.4f  %5.4f \n',[qlr_3break pv_3break]);
  fprintf(fileID,'    Break date: %5.4f  %5.4f  %5.4f\n',bd_3break);
 
path(p_initial);  % Reset path