clear all
clc
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  First Load and Get Data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
format long
load regm.raw

state=regm(:,7);
chst=regm(:,8);
[nn kk]=size(regm);

coll=regm(:,1);
merit=regm(:,2);
x=regm(:,3:5); 
year=regm(:,6);
year=year-1988;
state=regm(:,7);
stateo=state;
chst=regm(:,8);
mxst=max(state);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Set up state and year effects
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stid=zeros(nn,1);
nst=0;
nchst=0;
stchid=[];
stnchid=[];
for  ist=min(state):max(state)
 ttt=state(state==ist);
 if length(ttt)>0
  nst=nst+1;
  stid(state==ist)=nst;
  if mean(chst(state==ist))==1
   nchst=nchst+1;
   stchid=[stchid ist];
  else
   stnchid=[stnchid ist];
  end;
 end;
end;
nnchst=nst-nchst;


minyr=min(year);
maxyr=max(year);
nyrs=maxyr-minyr+1;
stdum=zeros(nn,nst);
yrdum=zeros(nn,nyrs);
for ii=1:nn
 stdum(ii,stid(ii))=1;
 yrdum(ii,year(ii)-minyr+1)=1;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Run main difference in differences regression
%%%   and form residuals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

XX=[merit x yrdum stdum(:,2:nst) ];
b=inv(XX'*XX)*XX'*coll;
alpha=b(1)

bt=b;
bt(1)=0;
eta=coll-XX*bt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  form djt which keeps track of years for which merit 
%%%   program was in place
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

djt=zeros(nchst,nyrs);
wt=zeros(nchst,nyrs);
wtdtil=zeros(nchst,nyrs);
for ist=1:nchst
 for iyr=1:nyrs
  djt(ist,iyr)=mean(merit((year==minyr+iyr-1)&(state==stchid(ist))));
  wt(ist,iyr)=sum(((year==minyr+iyr-1)&(state==stchid(ist))));
 end;
 yrt=year(state==stchid(ist));
 dtil=djt(ist,yrt-minyr+1);
 djt(ist,:)=djt(ist,:)-mean(dtil);
 for iyr=1:nyrs,
  wtdtil(ist,iyr)=wt(ist,iyr)*djt(ist,iyr);
 end;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  form etamn which keeps track of stateXyear mean
%%%   in residual
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
etamn=zeros(nnchst+nchst,nyrs);
for ist=1:nnchst
    for it=1:nyrs
        etamn(ist,it)=mean(eta(state==stnchid(ist)&year==minyr+it-1));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Finding critical value - upper bound for 95% CI
% set initial values and dimensions
countmax    = 10^6;     % max number of iterations
tol         = 10^-3;    % tolerance level
gamma       = 0.3;      % convergence parameter

numbsim=3000;
i025=floor(0.025*(numbsim));
i975=ceil(0.975*(numbsim));
i05=floor(0.050*(numbsim));
i95=ceil(0.950*(numbsim));

% randomly draw 10 states out of 51 (3000 times)
stsim=zeros(nchst,numbsim);
for isim=1:numbsim
    [gar idraw]=sort(rand(nnchst+nchst,1));
    stsim(:,isim)=idraw(1:nchst);
end

% start iterations
% hypothesize alpha=0 in the first step
alpha0 = 0;
count  = 0;
dif    = 1;
while dif > tol & count<= countmax
    count = count + 1;
    eta1=eta-alpha0.*merit;
    for ist=1:nchst
        for it=1:nyrs
            etamn(nnchst+ist,it)=mean(eta1(state==stchid(ist)&year==minyr+it-1));
        end
    end
    asim=zeros(numbsim,1);
    for isim=1:numbsim
        snum=0;
        sden=0;
        for ist=1:nchst
            snum=snum+wtdtil(ist,:)*etamn(stsim(ist,isim),:)';
            sden=sden+wtdtil(ist,:)*djt(ist,:)';
        end;
        asim(isim)=snum/sden;
    end;
    asim=sort(asim);
    alpha975=alpha-asim(i025);
    alpha1=alpha975;
    dif = abs((alpha0-alpha1)/(alpha0+0.01));
    alpha0 = gamma*alpha0+(1-gamma)*alpha1;
end
if count >= countmax
    display('maximum number of iterations reached')
    display('results displayed are from last iteration')
end

% Finding critical value - lower bound for 95% CI
% hypothesize alpha=0 in the first step
alpha0=0;
count   = 0;
dif      = 1;
while dif > tol & count<= countmax
    count = count + 1;
    eta1=eta-alpha0.*merit;
    for ist=1:nchst
        for it=1:nyrs
            etamn(nnchst+ist,it)=mean(eta1(state==stchid(ist)&year==minyr+it-1));
        end
    end
    asim=zeros(numbsim,1);
    for isim=1:numbsim
        snum=0;
        sden=0;
        for ist=1:nchst
            snum=snum+wtdtil(ist,:)*etamn(stsim(ist,isim),:)';
            sden=sden+wtdtil(ist,:)*djt(ist,:)';
        end;
        asim(isim)=snum/sden;
    end;
    asim=sort(asim);
    alpha025=alpha-asim(i975);
    alpha1=alpha025;
    dif = abs((alpha0-alpha1)/(alpha0+0.01));
    alpha0 = gamma*alpha0+(1-gamma)*alpha1;
end
if count >= countmax
    display('maximum number of iterations reached')
    display('results displayed are from last iteration')
end

% Finding critical value - upper bound for 90% CI
% hypothesize alpha=0 in the first step
alpha0 = 0;
count  = 0;
dif    = 1;
while dif > tol & count<= countmax
    count = count + 1;
    eta1=eta-alpha0.*merit;
    for ist=1:nchst
        for it=1:nyrs
            etamn(nnchst+ist,it)=mean(eta1(state==stchid(ist)&year==minyr+it-1));
        end
    end
    asim=zeros(numbsim,1);
    for isim=1:numbsim
        snum=0;
        sden=0;
        for ist=1:nchst
            snum=snum+wtdtil(ist,:)*etamn(stsim(ist,isim),:)';
            sden=sden+wtdtil(ist,:)*djt(ist,:)';
        end;
        asim(isim)=snum/sden;
    end;
    asim=sort(asim);
    alpha95=alpha-asim(i05);
    alpha1=alpha95;
    dif = abs((alpha0-alpha1)/(alpha0+0.01));
    alpha0 = gamma*alpha0+(1-gamma)*alpha1;
end
if count >= countmax
    display('maximum number of iterations reached')
    display('results displayed are from last iteration')
end

% Finding critical value - lower bound for 90% CI
% hypothesize alpha=0 in the first step
alpha0=0;
count   = 0;
dif      = 1;
while dif > tol & count<= countmax
    count = count + 1;
    eta1=eta-alpha0.*merit;
    for ist=1:nchst
        for it=1:nyrs
            etamn(nnchst+ist,it)=mean(eta1(state==stchid(ist)&year==minyr+it-1));
        end
    end
    asim=zeros(numbsim,1);
    for isim=1:numbsim
        snum=0;
        sden=0;
        for ist=1:nchst
            snum=snum+wtdtil(ist,:)*etamn(stsim(ist,isim),:)';
            sden=sden+wtdtil(ist,:)*djt(ist,:)';
        end;
        asim(isim)=snum/sden;
    end;
    asim=sort(asim);
    alpha05=alpha-asim(i95);
    alpha1=alpha05;
    dif = abs((alpha0-alpha1)/(alpha0+0.01));
    alpha0 = gamma*alpha0+(1-gamma)*alpha1;
end
if count >= countmax
    display('maximum number of iterations reached')
    display('results displayed are from last iteration')
end

disp('95% confidence interval')
[alpha025 alpha975]

disp('90% confidence interval')
[alpha05 alpha95]

toc