clear all
clc
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  First Load and Get Data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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);
clear regm;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  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;

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Now run the regression the second way which puts
%%%   equal weight on all states
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  First form state year dummies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
styr=zeros(nn,nyrs*nst);
for ii=1:nn,
 ist=stid(ii);
 iyr=year(ii)-minyr+1;
 styr(ii,(ist-1)*nyrs+iyr)=1;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  In first state estimate the state year dummies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XX=[x styr];
clear styr;
b1=inv(XX'*XX)*XX'*coll;
b1(1:3)=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  In second stage estimate the effect of hope
%%%    with state and year dummies
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yrdum=kron(ones(nst,1),eye(nyrs));
stdum=kron(eye(nst),ones(nyrs,1));
state=kron((1:nst)',ones(nyrs,1));
merit=zeros(size(b1));
for ist=1:nchst;
 trid=stid(stateo==stchid(ist));
 trid=trid(1);
 merit((trid-1)*nyrs+1:trid*nyrs)=djt(ist,:)';
end;
XX=[merit yrdum stdum(:,2:nst)];
b2=inv(XX'*XX)*XX'*b1;

alpha=b2(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Again form residuals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bt=b2;
bt(1)=0;
eta=b1-XX*bt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  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;
    
    asim=zeros(numbsim,1);
    for isim=1:numbsim
        snum=0;
        sden=0;
        for ist=1:nchst
            trid=stsim(ist,isim);
            etat=eta1((trid-1)*nyrs+1:trid*nyrs);
            dtil=djt(ist,:);
            dtil=dtil-mean(dtil);
            snum=snum+dtil*etat;
            sden=sden+(dtil*dtil');
        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;
    
    asim=zeros(numbsim,1);
    for isim=1:numbsim
        snum=0;
        sden=0;
        for ist=1:nchst
            trid=stsim(ist,isim);
            etat=eta1((trid-1)*nyrs+1:trid*nyrs);
            dtil=djt(ist,:);
            dtil=dtil-mean(dtil);
            snum=snum+dtil*etat;
            sden=sden+(dtil*dtil');
        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

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

toc