clc

omega=0;
ResFileName='Results.xlsx';
TabNames={'Table 1-1','Table 1-2','Table 2-1','Table 3-1','Table 4-1','Table 4-2','Table 4-3','Table 5-1','Table 5-2','Table 5-3'};
Boot=true;
CaseCalc=true;
Weighted=true;
Nrep=200;
Itmax=20;
WeightVar=39;
RightVars=[22,23,21,17,18,19,20,12,13,14,28,29,31,32,33,34,35,36,37,38];
Ncases=size(Cases,1);
pCases=zeros(Ncases,Nrep);
N=numel(RightVars)+1;
Extremes=[];
xExtremes=[];
xlswrite(ResFileName,0,TabNames{1},'a99');  %  check to see if we can write on the output file before investing in long calculations
for LeftVar=1:10
    
    %  Data for this equation
    
    TabName=TabNames{LeftVar}
    y=Data(:,LeftVar+1);
    Nobs=size(Data,1);
    weight=Data(:,WeightVar);
    Nob=size(Data,1);
    betas=[];
    pCases=[];
    X=zeros(Nob,N);
    X(:,1)=1;
    for var=2:N;
        X(:,var)=Data(:,RightVars(var-1));
    end
    
    %  Remove bad observations
    
    good=0;
    for obs=1:Nob
        check=[y(obs),X(obs,2:N)];
        bad=false;
        for var=1:N
            if isnan(check(var))
                bad=true;
                break
            end
        end
        % bad=bad && (isnan(weight(obs))||weight(obs)<0);
        if ~bad
            good=good+1;
            y(good)=y(obs);
            X(good,:)=X(obs,:);
            weight(good)=weight(obs)^.5;
        end
    end
    yS=y(1:good);
    XS=X(1:good,:);
    weightS=weight(1:good);
    Nob=good;
    X=zeros(Nob,N);
    y=zeros(Nob,1);
    beta=zeros(N,1);
    beta(1)=0.5;
    yt=zeros(Nob,1);
    Xt=zeros(Nob,N);
    w=zeros(Nob,1);
    p=zeros(Nob,1);
    
    if Boot
        for rep=1:Nrep
            
            %  Draw bootstrap sample
            
            for obs=1:Nob
                use=ceil(rand*Nob);
                y(obs)=yS(use);
                X(obs,:)=XS(use,:);
                weight(obs)=weightS(use);
            end
            
            %  Estimate
            
            for it=1:Itmax
                for obs=1:Nob
                    p(obs)=1/(1+exp(-X(obs,:)*beta));
                    w(obs)=1-omega +omega*(p(obs)*(1-p(obs)))^.5;
                    if w(obs)<10^-10
                        forget=true;
                        break
                    else
                        forget=false;
                    end
                    if Weighted
                        wgt=weight(obs);
                    else
                        wgt=1;
                    end
                    yt(obs)=wgt*(y(obs)-p(obs))/w(obs);
                    Xt(obs,:)=wgt*X(obs,:)*w(obs);
                end
                if forget
                    break
                end
                XX=Xt'*Xt;
                if rcond(XX)<10^-12
                    disp 'Ill-conditioned RHVs in bootstrap'
                    return
                end
                delta=XX\(Xt'*yt);
                if isnan(delta(1))
                    disp 'bad delta'
                    return
                end
                beta=beta+delta;
                for var=1:N
                    if beta(var)>5  %  runaway coefficient corresponding to p=1
                        beta(var)=5;
                    end
                    if beta(var)<-5  % runaway coefficient corresponding to p=0
                        beta(var)=-5;
                    end
                end
                Norm=(delta'*delta)^.5;
                if Norm < 10^(-10)
                    break
                end
            end
            if forget
                break
            end
            betas=[betas,beta];
            
            %  Probabilities for cases
            
            if CaseCalc
                for icase=1:Ncases
                    pCase(icase,1)=1/(1+exp(-Cases(icase,:)*beta));
                end
                pCases=[pCases,pCase];
            end
            
        end  %  rep
        
        betaSDs=std(betas')';
        xlswrite(ResFileName,betaSDs,TabName,'c2');
        NrepGood=size(pCases,2)
        pCasesD=zeros(Ncases,NrepGood);
        if CaseCalc
            pCasesSDs=std(pCases')';
            for icase=2:Ncases
                pCasesD(icase,:)=pCases(icase,:)-pCases(1,:);
            end
            xlswrite(ResFileName,pCasesSDs,TabName,'b30');
            pCasesDSDs=std(pCasesD')';
            xlswrite(ResFileName,pCasesDSDs,TabName,'e30');
        end
    end
    
    %  Standard estimation
    
    beta=zeros(N,1);
    for it=1:Itmax
        Iteration=it
        for obs=1:Nob
            p(obs)=1/(1+exp(-XS(obs,:)*beta));
            w(obs)=(p(obs)*(1-p(obs)))^.5;
            if Weighted
                wgt=weightS(obs);
            else
                wgt=1;
            end
            yt(obs)=wgt*(yS(obs)-p(obs))/w(obs);
            Xt(obs,:)=wgt*XS(obs,:)*w(obs);
        end
        XX=Xt'*Xt;
        if rcond(XX)<10^-12
            disp 'Ill-conditioned RHVs in standard estimation'
            return
        end
        delta=XX\(Xt'*yt);
        if isnan(delta(1))
            disp 'bad delta in standard'
            return
        end
        beta=beta+delta;
        for var=1:N
            if beta(var)>5  %  runaway coefficient corresponding to p=1
                beta(var)=5;
            end
            if beta(var)<-5  % runaway coefficient corresponding to p=0
                beta(var)=-5;
            end
        end
        Norm=(delta'*delta)^.5;
        if Norm < 10^(-10)
            break
        end
    end
    disp 'Write to Excel'
    xlswrite(ResFileName,RightVars',TabName,'a3');
    xlswrite(ResFileName,beta,TabName,'b2');
    xlswrite(ResFileName,mean(XS)',TabName,'e2');
    
    %  Probabilities for cases
    
    if CaseCalc
        pCases=zeros(Ncases,1);
        for icase=1:Ncases
            pCase=1/(1+exp(-Cases(icase,:)*beta));
            pCases(icase)=pCase;
        end
        xlswrite(ResFileName,pCases,TabName,'a30');
        pCasesD=zeros(Ncases,1);
        for icase=2:Ncases
            pCasesD(icase,:)=pCases(icase,:)-pCases(1,:);
        end
        
        %  Put in table format
        
        %  Levels
        
        loc=1;
        for icase=1:Ncases
            Table(loc,1)=pCases(icase);
            loc=loc+1;
            Table(loc,1)=-pCasesSDs(icase);
            loc=loc+1;
        end
        
        %  Differences
        
        loc=3;
        for icase=2:Ncases
            Table(loc,2)=pCasesD(icase);
            loc=loc+1;
            Table(loc,2)=-pCasesDSDs(icase);
            loc=loc+1;
        end
        disp 'Write table data'
        xlswrite(ResFileName,Table,TabName,'I2');
    end
    
    %  Extreme cases of fitted probability
    
    pMax=-1;
    pMin=2;
    for obs=1:Nob
        if p(obs)>pMax
            pMax=p(obs);
            obsMax=obs;
        end
        if p(obs)<pMin
            pMin=p(obs);
            obsMin=obs;
        end
    end
    Extremes=[Extremes,[p(obsMax);XS(obsMax,:)'],[p(obsMin);XS(obsMin,:)']];
    
    % Extreme values of RHVs
    
    xMax(1:N,1)=-10^10;
    xMin(1:N,1)=10^10;
    for obs=1:Nob
        for var=1:N
            if XS(obs,var)>xMax(var)
                xMax(var)=XS(obs,var);
            end
            if XS(obs,var)<xMin(var)
                xMin(var)=XS(obs,var);
            end
        end
    end
    xExtremes=[xExtremes,xMax,xMin];
end  %  LeftVar

xlswrite(ResFileName,Extremes,'Extremes','b2');
xlswrite(ResFileName,xExtremes,'xExtremes','b2');
