
%% Covariate Effects, Income
%---------------------------------
IncomeX0       = [ones(n,1) EVInfo_Lot EnvBetter OwnEV GovtCC_TooMuch GovtCC_Right   ...
           Agebelow50 Female GradnAbove SomeCollege Metropolitan MidWest South West ...
           USborn Black OtherRaces Democrat IndnElse Income/10000];

IncomeX1       = [ones(n,1) EVInfo_Lot EnvBetter OwnEV GovtCC_TooMuch GovtCC_Right  ...
           Agebelow50 Female GradnAbove SomeCollege Metropolitan MidWest South West ...
           USborn Black OtherRaces Democrat IndnElse (Income+10000)/10000];

% Outcome 1
%------------
oldProb1Income = zeros(n,mc);
newProb1Income = zeros(n,mc);

tic
wbar = waitbar(0,'Covariate Effect, Income, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Income(i,j)    = fast_normcdf(gammaCp(1,j) - IncomeX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Income(i,j)    = fast_normcdf(gammaCp(1,j) - IncomeX1(i,:)*betaStore(:,burn+j));
    end
waitbar(j/mc);    
end
close(wbar);
toc

diffIncomeProb1 = newProb1Income - oldProb1Income;
ceIncomeProb1   = mean(mean(diffIncomeProb1));
sprintf('%.04f',ceIncomeProb1)

% Outcome 2
%------------
oldProb2Income    = zeros(n,mc);        
newProb2Income    = zeros(n,mc);

oldComp2Income    = zeros(n,mc);
newComp2Income    = zeros(n,mc);

tic
wbar = waitbar(0,'Covariate Effect, Income, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Income(i,j)    = fast_normcdf(gammaCp(2,j)- IncomeX0(i,:)*betaStore(:,burn+j));
    oldProb2Income(i,j)    = oldComp2Income(i,j) - oldProb1Income(i,j);
    
    newComp2Income(i,j)    = fast_normcdf(gammaCp(2,j)- IncomeX1(i,:)*betaStore(:,burn+j));
    newProb2Income(i,j)    = newComp2Income(i,j) - newProb1Income(i,j);
    end
waitbar(j/mc);    
end
close(wbar);
toc

diffIncomeProb2 = newProb2Income - oldProb2Income;
ceIncomeProb2   = mean(mean(diffIncomeProb2));
sprintf('%.04f',ceIncomeProb2)

% Outcome 3
%------------
oldProb3Income    = zeros(n,mc);        
newProb3Income    = zeros(n,mc);

oldComp3Income    = zeros(n,mc);
newComp3Income    = zeros(n,mc);

tic
wbar = waitbar(0,'Covariate Effect, Income, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Income(i,j)    = fast_normcdf(gammaCp(3,j)- IncomeX0(i,:)*betaStore(:,burn+j));
    oldProb3Income(i,j)    = oldComp3Income(i,j) - oldComp2Income(i,j);

    newComp3Income(i,j)    = fast_normcdf(gammaCp(3,j)- IncomeX1(i,:)*betaStore(:,burn+j));
    newProb3Income(i,j)    = newComp3Income(i,j) - newComp2Income(i,j);
    end
waitbar(j/mc);    
end
close(wbar);
toc

diffIncomeProb3 = newProb3Income - oldProb3Income;
ceIncomeProb3   = mean(mean(diffIncomeProb3));
sprintf('%.04f',ceIncomeProb3)

% Outcome 4
%------------
oldProb4Income    = zeros(n,mc);        
newProb4Income    = zeros(n,mc);

tic
wbar = waitbar(0,'Covariate Effect, Income, Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4Income(i,j)    = 1 - oldComp3Income(i,j);
    
    newProb4Income(i,j)    = 1 - newComp3Income(i,j);
    end
waitbar(j/mc);    
end
close(wbar);
toc

diffIncomeProb4 = newProb4Income - oldProb4Income;
ceIncomeProb4   = mean(mean(diffIncomeProb4));
sprintf('%.04f',ceIncomeProb4)

disp('CE, Income')
disp([ceIncomeProb1 ceIncomeProb2 ceIncomeProb3 ceIncomeProb4])


%% Table for Covariate Effects

rowlabelsCE = char('EVInfo_Lot','EnvBetter','OwnEV','GovtCC_TooMuch', ...
                 'GovtCC_Right','Agebelow50','Female','GradnAbove','SomeCollege', ... 
                 'Metropolitan','MidWest', 'South','West','USborn','Black','OtherRaces', ...
                 'Democrat','IndnElse','Inc=10k');

ceProb1  = [ceEVInfoProb1 ceEnvImpactProb1 ceOwnEVProb1 ceGACCTMProb1 ceGACCRAProb1 ...
            ceAgeb50Prob1 ceFemaleProb1 ceGradnAboveProb1 ceSomeCollegeProb1 ...
            ceMetropolitanProb1 ceMidwestProb1 ceSouthProb1 ceWestProb1 ceUSbornProb1 ...
            ceBlackProb1 ceOtherRacesProb1 ceDemocratProb1 ceIndnElseProb1 ceIncomeProb1];

ceProb2  = [ceEVInfoProb2 ceEnvImpactProb2 ceOwnEVProb2 ceGACCTMProb2 ceGACCRAProb2 ...
            ceAgeb50Prob2 ceFemaleProb2 ceGradnAboveProb2 ceSomeCollegeProb2 ...
            ceMetropolitanProb2 ceMidwestProb2 ceSouthProb2 ceWestProb2 ceUSbornProb2 ...
            ceBlackProb2 ceOtherRacesProb2 ceDemocratProb2 ceIndnElseProb2 ceIncomeProb2];

ceProb3  = [ceEVInfoProb3 ceEnvImpactProb3 ceOwnEVProb3 ceGACCTMProb3 ceGACCRAProb3 ...
            ceAgeb50Prob3 ceFemaleProb3 ceGradnAboveProb3 ceSomeCollegeProb3 ...
            ceMetropolitanProb3 ceMidwestProb3 ceSouthProb3 ceWestProb3 ceUSbornProb3 ...
            ceBlackProb3 ceOtherRacesProb3 ceDemocratProb3 ceIndnElseProb3 ceIncomeProb3];

ceProb4  = [ceEVInfoProb4 ceEnvImpactProb4 ceOwnEVProb4 ceGACCTMProb4 ceGACCRAProb4 ...
            ceAgeb50Prob4 ceFemaleProb4 ceGradnAboveProb4 ceSomeCollegeProb4 ...
            ceMetropolitanProb4 ceMidwestProb4 ceSouthProb4 ceWestProb4 ceUSbornProb4 ...
            ceBlackProb4 ceOtherRacesProb4 ceDemocratProb4 ceIndnElseProb4 ceIncomeProb4];


    fprintf('    Covariate Effects, Ordinal Probit, April 2021          \n')
    fprintf('___________________________________________________ \n')
    fprintf('                 P(y=1)   P(y=2)   P(y=3)   P(y=4)     \n')
    fprintf('                 _____   ______   ______   ______ \n')
      
    for i=1:length(ceProb1)
        fprintf('%-10s % 7.4f   %7.4f  %7.4f  %7.4f  \n', rowlabelsCE(i,:),...
            ceProb1(i), ceProb2(i), ceProb3(i), ceProb4(i) )
    end
    fprintf('___________________________________________________ \n')

