% Program to simulate the power of Common Regression tests
% in the linear regression model. 
% Reproduces Reiss (2016) Figures 1-3 and Models 1-8
% plus calculates auxiliary estimates
clear all; close all; 
state = 16483;
rand('state',state);
randn('state',state);
%
% Set isims to one if simulations are to be compared to formulaic results...
% sims = # simulated regressions
% sims2 = # of simulated deltas used in Models 7 and 8 to do crude freq.
% simulation of eqn (26)
% sims3 = larger # of simulated deltas used in Models 7 and 8 to do crude freq.
% simulation of eqn (26)
isims=1;   sims=100000;    sims2=500;     sims3=5000;
%
%  Testing whether time trend coef equals zero
%  Reg Model    Y = C + Reg_slope * X + e
%  Trend Model  X = Time_slope * t + n
%
R=[0 1]; r =[0];    % Tests slope on x equal to zero
reg_slope=[0.1 0.1 0.1 0.1 0.0 0.1 1.0 1.0]'; % True coefficient on x
time_slope=[1.0 1.0 1.0 1.0 0.1 .1 .1 .1]';   % Trend coefficient
Tvec=[10 20 50 100 20 20 20 20 ]';   % Sample size       
sigsq_n=[0 0 0 0 0 0 1 4];           % Trend error var
sigsq_e=[1 1 1 1 1 1 1 1 ];          % Reg model error var
%
monte_m=[]; monte_dist_m=[]; imhof_m=[]; anal_m=[]; sim_anal=[];
%
% Initial loop to do Monte Carlo - sims times
%   kk indexes model number
%
disp('kk  reg_slope(kk)  time_slope(kk)  sigsq_e(kk)   sigsq_n(kk)   Delta1   Delta2   sig_ee   theta   rho_sq')
for kk=1:size(Tvec,1)
  T=Tvec(kk);   t=(1:1:T)';   %Set sample size and create time trend
  Q_1=[]; Q_2=[]; b_hat=[];   
  beta=[1 reg_slope(kk)]';          % Regression Coefficients
  if isims==1                 % Equals 1 for monte carlos (not 1 means no monte carlos)
  % sig_n = 0 is deterministic trend; sig_n > 0 is stochastic
  for k=1:sims;
    % Generate Simulated Regression Data
    if sigsq_n(kk)>0   % Stochastic regressor case
        x1=time_slope(kk)*t + normrnd(zeros(T,1),sqrt(sigsq_n(kk)));
    else               % Nonstochastic regressor case
        x1=t;
    end
    X=[ones(T,1) x1];  
    y=X*beta+normrnd(zeros(T,1),sqrt(sigsq_e(kk)));
    % Calculate Noncentrality Parameters
    W=X*inv(X'*X)*(R');  mu=zeros(T,1);
    Delt1(k,kk)=(R*beta-r+W'*mu)'*inv(W'*W)*(R*beta-r+W'*mu)/sigsq_e(kk);
    Delt2(k,kk)=mu'*mu - mu'*X*inv(X'*X)*X'*mu/sigsq_e(kk);
    % Regression coefficients plus m stat = r sqr
    b_hat=[b_hat; (inv(X'*X)*X'*y)'];
    z=(y-W*inv(W'*W)*r)/sqrt(sigsq_e(kk));
    Q_1=[Q_1; z'*W*inv(W'*W)*W'*z];
    Q_2=[Q_2; z'*z-z'*X*inv(X'*X)*X'*z];
    monte_m(k,kk)=Q_1(k)/(Q_1(k)+Q_2(k));
  end
  end % End Monte Carlo Simulation Block

% Loop over models continues
% Here we do formula (26) for the stochastic case or formula (27)
% for the nonstochastic case.

K=1; % Remember constant is not part of K calculation
a=K/2; b=(T-K-1)/2; c=(T-1)/2; 
% Set up grid to evaluate cdf and pdf
r_vec=(0:.005:.995)'; theta=[]; rho_sq=[];
% If stochastic regressor case 
if sigsq_n(kk) > 0 
  theta=((time_slope(kk)/sqrt(sigsq_n(kk)))^2)*((T-1)*T*(T+1)/12);  
  rho_sq=((reg_slope(kk)*sqrt(sigsq_n(kk)))^2)/(sigsq_e(kk)+(reg_slope(kk)*sqrt(sigsq_n(kk)))^2);
  % Leading factor in eqn (26)
  fac=exp(-rho_sq*theta*(1-r_vec)./(2*(1-rho_sq*r_vec))).*gamma(c).*(r_vec.^(a)).*(1-r_vec).^(b-1).*...
    ((1-rho_sq)./(1-rho_sq*r_vec)).^c;
  % fac2 tracks triple sum in eqn (26)
  fac2=0;
  z=(theta/2)*(1-rho_sq)./(1-rho_sq*r_vec);
  % Triple sum
  for l=0:b-1;
    for j=0:b-l-1;
        for i=0:j
            fac2=fac2+(((r_vec)./(1-r_vec)).^l).*(rho_sq*(r_vec)./(1-rho_sq*(r_vec))).^j.*...
                (gamma(c+j)/(gamma(b-l-j)*gamma(a+l+j+1))).*(z.^i)/...
                (gamma(j+1-i)*gamma(c+i)*gamma(i+1));
        end;
    end;
  end;
else    % Nonstochastic Regressor calculation -- equation (27)
  delt=(reg_slope(kk))^2*((T-1)*T*(T+1)/12)/sigsq_e(kk);
  fac=exp(-delt*(1-r_vec)/2).*gamma(c).*(r_vec.^(a)).*(1-r_vec).^(b-1);
  % fac2 tracks double sum in eqn (27)
  fac2=0;
  % Double sum
  for l=0:b-1;
    for j=0:b-l-1;
        fac2=fac2+((r_vec./(1-r_vec)).^l).*(delt*r_vec./2).^j./...
                (gamma(b-l-j)*gamma(a+l+j+1)*gamma(j+1));
    end;
  end;
end
% Putting the sums together with the leading terms in (26) or (27) ...
anal_m(:,kk)=(fac.*fac2);
% If stochastic regressor case, average over Deltas in (27)
if sigsq_n(kk) > 0 
  for k=1:sims3
  delt=Delt1(k,kk);
  fac3=exp(-delt*(1-r_vec)/2).*gamma(c).*(r_vec.^(a)).*(1-r_vec).^(b-1);
  fac4=0;
  for l=0:b-1;
    for j=0:b-l-1;
        fac4=fac4+((r_vec./(1-r_vec)).^l).*(delt*r_vec./2).^j./...
                (gamma(b-l-j)*gamma(a+l+j+1)*gamma(j+1));
    end;
  end;
  if k==1 
      sim_anal(:,kk)=fac3.*fac4/sims2;
      sim_anal(:,8+kk)=fac3.*fac4/sims3;
  else if k <= sims2
      sim_anal(:,kk)=sim_anal(:,kk)+fac3.*fac4/sims2;
  end
      sim_anal(:,8+kk)=sim_anal(:,8+kk)+fac3.*fac4/sims3;
  end
  end
end
%
% Imhof rountine calculation for comparison
%
sig_ee=sqrt(reg_slope(kk)^2*sigsq_n(kk)+sigsq_e(kk));
if sigsq_e(kk) > 0 
  t=(1:1:T)';
  X=[ones(T,1) t];   % Regressors
  beta=[1 reg_slope(kk)*time_slope(kk)]';     % Time Trend Model Parameters: y = b_0 + reg_slope x_t + e_t
  sig_ee=sqrt(reg_slope(kk)^2*sigsq_n(kk)+sigsq_e(kk));
  W=X*inv(X'*X)*(R');
  mu=zeros(T,1);   % Conditional mean of e_t (no misspecification here)
  % Noncentrality parameters
  Delta1=(R*beta-r+W'*mu)'*inv(W'*W)*(R*beta-r+W'*mu)/sig_ee^2;
  Delta2=mu'*mu - mu'*X*inv(X'*X)*X'*mu/sig_ee^2;
  delt=[Delta1; Delta2];              % Noncentrality parameters
  multip=[rank(R); T-size(X,2)];    % Multiplicity of roots
  arg=0;                            % q in Imhof's notation  
  u_bnd=100;                        % Integration area width each time integral is evaluated
  amax=size(r_vec(:),1);            % Number of equal-spaced points to evaluate CDF on [0,1]
  M_chi=[];                         % Equal-spaced points.
  lambd_chi=[];                     % Lambda's in Imhof's notation
  Integs=50;
 for j=1:amax;
    M_chi(j)=r_vec(j);         
    lambd_chi(:,j)=[1-M_chi(j); -M_chi(j)];
  % Can't integrate to infinity. Instead, break up the integral into Integs pieces of width 
  %   u_bnd
    temp_res=0.5;
    for jj=1:Integs
      temp_res=temp_res-(1/pi)*quad(@imhof_integrand,1e-15+(jj-1)*u_bnd,...
          u_bnd*jj,1e-10,0,arg,lambd_chi(:,j),multip,delt);
    end;
    imhof_m(j,kk)=temp_res;
 end;
end  % Imhof Calcs
%
%  Compute EDF for simulations
%
   amax=size(r_vec(:),1);
   for j=1:amax;
      M_chi(j)=r_vec(j);
      lambd_chi(:,j)=[1-M_chi(j); -M_chi(j)];
      monte_dist_m(j,kk)=sum(monte_m(:,kk)<M_chi(j))/sims;
   end;
%
disp([kk reg_slope(kk) time_slope(kk) sigsq_e(kk) sigsq_n(kk) Delta1 Delta2 sig_ee theta rho_sq])
end;

figure(1)
     plot(r_vec,anal_m(:,1),'LineWidth',2,'Color','g')
       hold on
     plot(r_vec,anal_m(:,2),'LineWidth',3,'Color','b')
     plot(r_vec,anal_m(:,3),'LineStyle',':','LineWidth',3,'Color','k')
     plot(r_vec,anal_m(:,4),'LineStyle','-.','LineWidth',3,'Color','r')
   
     set(gca,'Ylim',[0,1],'FontWeight','bold','FontSize',12)
     box off
     legend(' Model 1, N = 10',' Model 2, N = 20',' Model 3, N = 50',' Model 4, N = 100','location','Northwest')
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('CDF(R^2)','FontSize',15,'FontWeight','bold')
     %title('CDF of R^2 for Time Trend Model (Models 5-8)','FontSize',15,'FontWeight','bold')
     set(gcf,'PaperOrientation','portrait')
     set(gcf,'PaperPositionMode','manual')
     set(gcf, 'PaperUnits', 'inches');
     set(gcf, 'PaperSize', [8.25 6]);
     set(gcf, 'renderer', 'painters');
     print ('-djpeg100', 'C:\Users\preiss\Documents\Dropbox\Power\powercalc_fig1');
 figure(11)
     plot(r_vec,anal_m(:,5),'LineWidth',2,'Color','g')
       hold on
     plot(r_vec,anal_m(:,6),'LineWidth',3,'Color','b')
     plot(r_vec,anal_m(:,7),'LineStyle',':','LineWidth',3,'Color','k')
     plot(r_vec,anal_m(:,8),'LineStyle','-.','LineWidth',3,'Color','r')
   
     set(gca,'Ylim',[0,1],'FontWeight','bold','FontSize',12)
     box off
     legend(' Model 5 ',' Model 6 ',' Model 7 ',' Model 8 ','location','North')
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('CDF(R^2)','FontSize',15,'FontWeight','bold')
     %title('CDF of R^2 for Time Trend Model (Models 5-8)','FontSize',15,'FontWeight','bold')
          set(gcf,'PaperOrientation','portrait')
     set(gcf,'PaperPositionMode','manual')
     set(gcf, 'PaperUnits', 'inches');
     set(gcf, 'PaperSize', [8.25 6]);
     set(gcf, 'renderer', 'painters');
     print ('-djpeg100', 'C:\Users\preiss\Documents\Dropbox\Power\powercalc_fig2');

figure(2)   % Numerical derivative approximations to pdf 
     M_chi=M_chi(:);
     plot((M_chi(2:end,1)+M_chi(1:end-1,:))*.5,(anal_m(2:end,1)-anal_m(1:end-1,1))/...
          (M_chi(2,1)-M_chi(1,1)),'LineWidth',2,'Color','g')
     hold on
     plot((M_chi(2:end,1)+M_chi(1:end-1,1))*.5,(anal_m(2:end,2)-anal_m(1:end-1,2))/...
          (M_chi(2,1)-M_chi(1,1)),'LineStyle',':','LineWidth',3,'Color','b')
     plot((M_chi(2:end,1)+M_chi(1:end-1,1))*.5,(anal_m(2:end,3)-anal_m(1:end-1,3))/...
          (M_chi(2,1)-M_chi(1,1)),'LineStyle','-.','LineWidth',3,'Color','k')
     plot((M_chi(2:end,1)+M_chi(1:end-1,1))*.5,(anal_m(2:end,4)-anal_m(1:end-1,4))/...
          (M_chi(2,1)-M_chi(1,1)),'LineStyle','--','LineWidth',3,'Color','r')
     box off
    legend(' Model 1, N = 10',' Model 2, N = 20',' Model 3, N = 50',' Model 4, N = 100','location','Northwest')
     set(gca,'Ylim',[0 25],'FontWeight','bold','FontSize',12)
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('PDF(R^2)','FontSize',15,'FontWeight','bold')
     %title('PDF of R^2 for Time Trend Model (Models 5-8)','FontSize',15,'FontWeight','bold')
          set(gcf,'PaperOrientation','portrait')
     set(gcf,'PaperPositionMode','manual')
     set(gcf, 'PaperUnits', 'inches');
     set(gcf, 'PaperSize', [8.25 6]);
     set(gcf, 'renderer', 'painters');
     print ('-djpeg100', 'C:\Users\preiss\Documents\Dropbox\Power\powercalc_fig3');

 figure(12) % Numerical derivative approximations to pdf 
     M_chi=M_chi(:);
     plot((M_chi(2:end,1)+M_chi(1:end-1,:))*.5,(anal_m(2:end,5)-anal_m(1:end-1,5))/...
          (M_chi(2,1)-M_chi(1,1)),'LineWidth',2,'Color','g')
     hold on
     plot((M_chi(2:end,1)+M_chi(1:end-1,1))*.5,(anal_m(2:end,6)-anal_m(1:end-1,6))/...
          (M_chi(2,1)-M_chi(1,1)),'LineStyle',':','LineWidth',3,'Color','b')
     plot((M_chi(2:end,1)+M_chi(1:end-1,1))*.5,(anal_m(2:end,7)-anal_m(1:end-1,7))/...
          (M_chi(2,1)-M_chi(1,1)),'LineStyle','-.','LineWidth',3,'Color','k')
     plot((M_chi(2:end,1)+M_chi(1:end-1,1))*.5,(anal_m(2:end,8)-anal_m(1:end-1,8))/...
          (M_chi(2,1)-M_chi(1,1)),'LineStyle','--','LineWidth',3,'Color','r')
     box off
    legend(' Model 5 ',' Model 6 ',' Model 7 ',' Model 8 ','location','North')
     set(gca,'Ylim',[0 11],'FontWeight','bold','FontSize',12)
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('PDF(R^2)','FontSize',15,'FontWeight','bold')
     %title('PDF of R^2 for Time Trend Model (Models 5-8)','FontSize',15,'FontWeight','bold')
     set(gcf,'PaperOrientation','portrait')
     set(gcf,'PaperPositionMode','manual')
     set(gcf, 'PaperUnits', 'inches');
     set(gcf, 'PaperSize', [8.25 6]);
     set(gcf, 'renderer', 'painters');
     print ('-djpeg100', 'C:\Users\preiss\Documents\Dropbox\Power\powercalc_fig4');

 if isims==1
     figure(3)
     plot(M_chi,anal_m(:,1)-monte_dist_m(:,1),'LineWidth',2,'Color','g')
     hold on
     plot(M_chi,anal_m(:,2)-monte_dist_m(:,2),'LineWidth',2,'Color','b')
     plot(M_chi,anal_m(:,3)-monte_dist_m(:,3),'LineStyle',':','LineWidth',2,'Color','k')
     plot(M_chi,anal_m(:,4)-monte_dist_m(:,4),'LineStyle','-.','LineWidth',2,'Color','r')
     set(gca,'Ylim',[-.005 .005],'FontWeight','bold','FontSize',12)
     box off
     legend(' Model 1 ',' Model 2 ',' Model 3 ',' Model 4 ')
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('Approx Error','FontSize',15,'FontWeight','bold')
     title('CDF Approximation Error: Series Formula - Monte Carlo','FontSize',15,'FontWeight','bold')
     figure(13)
     plot(M_chi,anal_m(:,5)-monte_dist_m(:,5),'LineWidth',2,'Color','g')
     hold on
     plot(M_chi,anal_m(:,6)-monte_dist_m(:,6),'LineWidth',2,'Color','b')
     plot(M_chi,anal_m(:,7)-monte_dist_m(:,7),'LineStyle',':','LineWidth',2,'Color','k')
     plot(M_chi,anal_m(:,8)-monte_dist_m(:,8),'LineStyle','-.','LineWidth',2,'Color','r')
     set(gca,'Ylim',[-.005 .005],'FontWeight','bold','FontSize',12)
     box off
     legend(' Model 5 ',' Model 6 ',' Model 7 ',' Model 8 ')
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('Approx Error','FontSize',15,'FontWeight','bold')
     title('CDF Approximation Error: Series Formula - Monte Carlo','FontSize',15,'FontWeight','bold')
 end
 
figure(4)
     plot(M_chi,anal_m(:,1)-imhof_m(:,1),'LineWidth',2,'Color','g')
     hold on
     plot(M_chi,anal_m(:,2)-imhof_m(:,2),'LineWidth',2,'Color','b')
     set(gca,'Ylim',[-.0005 .0005],'FontWeight','bold','FontSize',12)
     box off
     legend(['Model 5:   \beta = ' num2str(reg_slope(1)) ', \alpha = ' num2str(time_slope(1)) ...
             ', \sigma_\epsilon^2 = ' num2str(sigsq_e(kk)) ', \sigma_\eta^2 = ' num2str(sigsq_n(1))],...
            ['Model 6:   \beta = ' num2str(reg_slope(2)) ', \alpha = ' num2str(time_slope(2)) ...
             ', \sigma_\epsilon^2 = ' num2str(sigsq_e(kk)) ', \sigma_\eta^2 = ' num2str(sigsq_n(2))])
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('Approx Error','FontSize',15,'FontWeight','bold')
     title('CDF Approximation Error: Series Formula - Monte Carlo','FontSize',15,'FontWeight','bold')
     
figure(5)
     plot(M_chi,anal_m(:,7)-sim_anal(:,7),'LineWidth',2,'Color','g')
     hold on
     plot(M_chi,anal_m(:,7)-sim_anal(:,15),'LineWidth',2,'Color','r')
     plot(M_chi,anal_m(:,8)-sim_anal(:,8),'LineWidth',2,'Color','b')
     plot(M_chi,anal_m(:,8)-sim_anal(:,16),'LineWidth',2,'Color','k')
     set(gca,'Ylim',[-.025 .025],'FontWeight','bold','FontSize',12)
     box off
     legend(['Model 7:   Average over 500 Deltas'],...
            ['Model 7:   Average over 5000 Deltas'], ...
            ['Model 8:   Average over 500 Deltas'], ...
            ['Model 8:   Average over 5000 Deltas'])
     xlabel('R^2','FontSize',15,'FontWeight','bold')
     ylabel('Approx Error','FontSize',15,'FontWeight','bold')
     title('CDF Approximation Error: Series Formula - Monte Carlo','FontSize',15,'FontWeight','bold')
     
figure(6)
     ksdensity(monte_m(:,1),'support',[0 1])
     hold on
     ksdensity(monte_m(:,2),'support',[0 1])
     ksdensity(monte_m(:,3),'support',[0 1])
     ksdensity(monte_m(:,4),'support',[0 1])
     legend('1','2','3','4')
     set(gca,'Ylim',[0 10])