% Electric Vehicle, Wave89 Data
% Estimating the Ordinal Probit Model using Random-Walk MH  Algorithm

clear all
clc
% [data]  = xlsread('/Users/arshadrahman/Library/CloudStorage/Dropbox/1-WorkFiles/2 Research/Electric Vehicle/W89_Apr21/ATP_W89.xlsx','CleanData','A2:N11053');
[data]  = xlsread('ATP_W89.xlsx','CleanData','A2:N11053');


%% Reading the data

% Q: The next time you purchase a vehicle, how likely are you to seriously consider purchasing an electric vehicle? 
EVcar2              = data(:,1);        % 1=very likely, 2=somewhat likely, 3=not too likely, 4=not at all likely

% Q: How much, if anything, have you read or heard about EVs?
EVheard             = data(:,2);        % 1=a lot, 2= a little, 3=nothing at all

% Q: Comparing EVs to gas-powered vehicles, in general, would you say EVs are
EnvImpact             = data(:,3);        % 1=better for the env than gas-powered vehicles, 2=worse ..., 3=about the same

% Q: Do you currently have an electric or hybrid vehicle, or not?
EVcar1              = data(:,4);        % 1=yes, 2=no

% Q: How much do you thing the federal govt. is doing to reduce the effects
% of global climate change
envir8e             = data(:,5);        % 1=too much, 2=too little, 3=about the right amount
%---------
metro               = data(:,6);        % 1=metropolitan, 2=non-metropolitan
region              = data(:,7);        % 1=northeast, 2=midwest, 3=south, west=4
ageCat              = data(:,8);        % 1=age 18-29, 2=age 30-49, 3=age 50-64, 4=age 65+ 
gender              = data(:,9);        % 1=male, 2=female
eduCat              = data(:,10);       % 1=graduate & above, 2=some college, 3=HS and less
YearsinUS           = data(:,11);       % 1=born in US, 2=0-10 years, 3=11-20 years, 4=21+ years 
race                = data(:,12);       % 1=White non-Hispanic, 2=Black non-Hispanic, 3=Hispanic, 4=Other, 5=Asian non-Hispanic
party               = data(:,13);       % 1=Republican, 2=Democrat, 3=Independent, 4=Something else
incCat              = data(:,14);       % 1=Less than $30,000, 2=$30,000 to less than $40,000, 3=$40,000 to less than $50,000
                                        % 4=$50,000 to less than $60,000, 5=$60,000 to less than $70,000, 6=$70,000 to less than $80,000
                                        % 7=$80,000 to less than $90,000,
                                        % 8=$90,000 to less than $100,000,
                                        % 9=$100,000 or more
%% Recoding the variables

y = 5-EVcar2;               % Changing the order
n = size(y,1);              % number of observations

% EVheard
%---------------------
disp('tabulate(EVheard)')
tabulate(EVheard);

EVInfo_LN      = zeros(n,1);      % Little and Nothing merged, latter has 402 obs (3.64 percent)
EVInfo_Lot     = zeros(n,1);

EVInfo_LN(EVheard~=1)  = 1;    % base category
EVInfo_Lot(EVheard==1) = 1;


% EVautod
%---------------------
disp('tabulate(EVautod)')
tabulate(EnvImpact);

EnvBetter    = zeros(n,1);
EnvSameWorse = zeros(n,1);    % Worse=2 or same=3, base category
EnvBetter(EnvImpact==1) = 1;

for i=1:n
    if EnvImpact(i,1)==2||EnvImpact(i,1)==3
        EnvSameWorse(i,1)=1;
    end
end


% EVcar1
%---------------------
disp('tabulate(EVcar1)')
tabulate(EVcar1);

OwnEV    = zeros(n,1);
OwnEV(EVcar1==1) = 1;

NotOwnEV = ones(n,1)-OwnEV;

% Federal Government
%---------------------
disp('tabulate(envir8e)')
tabulate(envir8e);

GovtCC_Little  = zeros(n,1);
GovtCC_Right   = zeros(n,1);
GovtCC_TooMuch = zeros(n,1);

GovtCC_Little(envir8e==2)  = 1;     % base category
GovtCC_Right(envir8e==3)   = 1;
GovtCC_TooMuch(envir8e==1) = 1;

% Metropolitan
%---------------------
disp('tabulate(metro)')
tabulate(metro);

Metropolitan = zeros(n,1);
Metropolitan(metro==1) = 1;

%---------------------
disp('tabulate(region)')
tabulate(region)

NorthEast   = zeros(n,1);
MidWest     = zeros(n,1);
South       = zeros(n,1);
West        = zeros(n,1);

NorthEast(region==1) = 1;   % base category            
MidWest(region==2)   = 1;
South(region==3)     = 1;
West(region==4)      = 1;

% age Category
%---------------------
disp("tabulate(ageCat)")
tabulate(ageCat);

Agebelow50   = zeros(n,1);

for i=1:n
    if ageCat(i,1)==1||ageCat(i,1)==2
        Agebelow50(i,1)=1;
    end
end
% Age above 50 is the base category

Ageabove50 = ones(n,1)-Agebelow50;

% Gender
%---------------------
Female = gender - 1;
% male is base category
Male  = ones(n,1) - Female;

% Education
%---------------------
disp("tabulate(eduCat)")
tabulate(eduCat);

GradnAbove  = zeros(n,1);
SomeCollege = zeros(n,1);
HSnBelow    = zeros(n,1);

GradnAbove(eduCat==1)  = 1;
SomeCollege(eduCat==2) = 1;
HSnBelow(eduCat==3)    = 1;

% HS and below is the base category

% YearsinUS
%---------------------
disp("tabulate(YearsinUS)")
tabulate(YearsinUS);

USborn   = zeros(n,1);
USborn(YearsinUS==1) = 1;

% Not born in the US is the base category


% Race
%---------------------
disp("tabulate(race)")
tabulate(race);

White      = zeros(n,1);      % base category
Black      = zeros(n,1);

White(race==1) = 1;
Black(race==2) = 1;
OtherRaces = ones(n,1) - (White + Black);

% Party
%---------------------
disp("tabulate(party)")
tabulate(party);

Republican  = zeros(n,1);     % base category
Democrat    = zeros(n,1);

Republican(party==1) = 1;
Democrat(party==2)   = 1;
IndnElse   = ones(n,1) - (Republican + Democrat);

% Income
%---------------------
Income = zeros(n,1);

for i=1:n
    if incCat(i,1)  == 1
        Income(i,1) = 20000;
    elseif incCat(i,1) == 2
        Income(i,1) = 35000;
    elseif incCat(i,1) == 3
        Income(i,1) = 45000;
    elseif incCat(i,1) == 4
        Income(i,1) = 55000;
    elseif incCat(i,1) == 5
        Income(i,1) = 65000;
    elseif incCat(i,1) == 6
        Income(i,1) = 75000;
    elseif incCat(i,1) == 7
        Income(i,1) = 85000;
    elseif incCat(i,1) == 8
        Income(i,1) = 95000;
    else 
        Income(i,1) = 110000;
    end
end

%% Data Summary, Table 1

meanIncome = mean(Income/10000);
stdIncome  = std(Income/10000);

countEVInfo_Lot = sum(EVInfo_Lot);
perEVInfo_Lot   = (sum(EVInfo_Lot)/n)*100;

countEVInfo_LN  = sum(EVInfo_LN);
perEVInfo_LN    = (sum(EVInfo_LN)/n)*100;

countEnvBetter   = sum(EnvBetter);
perEnvBetter     = (sum(EnvBetter)/n)*100;

countEnvSameWorse = sum(EnvSameWorse);
perEnvSameWorse   = (sum(EnvSameWorse)/n)*100;

countOwnEV        = sum(OwnEV);
perOwnEV          = (sum(OwnEV)/n)*100;

countGovtCC_TooMuch = sum(GovtCC_TooMuch);
perGovtCC_TooMuch   = (sum(GovtCC_TooMuch)/n)*100;

countGovtCC_Right = sum(GovtCC_Right);
perGovtCC_Right   = (sum(GovtCC_Right)/n)*100;

countGovtCC_Little = sum(GovtCC_Little);
perGovtCC_Little   = (sum(GovtCC_Little)/n)*100;

countAgebelow50  = sum(Agebelow50);
perAgebelow50    = (sum(Agebelow50)/n)*100;

countFemale      = sum(Female);
perFemale        = (sum(Female)/n)*100;

countGradnAbove  = sum(GradnAbove);
perGradnAbove    = (sum(GradnAbove)/n)*100;

countSomeCollege = sum(SomeCollege);
perSomeCollege   = (sum(SomeCollege)/n)*100;

countHSnBelow    = sum(HSnBelow);
perHSnBelow      = (sum(HSnBelow)/n)*100;

countMetropolitan = sum(Metropolitan);
perMetropolitan   = (sum(Metropolitan)/n)*100;

countNorthEast   = sum(NorthEast);
perNorthEast     = (sum(NorthEast)/n)*100;

countMidWest     = sum(MidWest);
perMidWest       = (sum(MidWest)/n)*100;

countSouth       = sum(South);
perSouth         = (sum(South)/n)*100;

countWest        = sum(West);
perWest          = (sum(West)/n)*100;

countUSborn      = sum(USborn);
perUSborn        = (sum(USborn)/n)*100;

countWhite       = sum(White);
perWhite         = (sum(White)/n)*100;

countBlack       = sum(Black);
perBlack         = (sum(Black)/n)*100;

countOtherRaces  = sum(OtherRaces);
perOtherRaces    = (sum(OtherRaces)/n)*100;

countDemocrat    = sum(Democrat);
perDemocrat      = (sum(Democrat)/n)*100;

countRepublican  = sum(Republican);
perRepublican    = (sum(Republican)/n)*100;

countIndnElse    = sum(IndnElse);
perIndnElse      = (sum(IndnElse)/n)*100;

% Data Summary Table
%----------------------------

MeanCounts = [meanIncome countEVInfo_Lot countEVInfo_LN countEnvBetter ...
     countEnvSameWorse countOwnEV countGovtCC_TooMuch countGovtCC_Right countGovtCC_Little ...
     countAgebelow50 countFemale countGradnAbove countSomeCollege countHSnBelow ...
     countMetropolitan countNorthEast countMidWest countSouth countWest countUSborn ...
     countWhite countBlack countOtherRaces countDemocrat countRepublican countIndnElse];

StdPer = [stdIncome perEVInfo_Lot perEVInfo_LN perEnvBetter ...
     perEnvSameWorse perOwnEV perGovtCC_TooMuch perGovtCC_Right perGovtCC_Little ...
     perAgebelow50 perFemale perGradnAbove perSomeCollege perHSnBelow ...
     perMetropolitan perNorthEast perMidWest perSouth perWest perUSborn ...
     perWhite perBlack perOtherRaces perDemocrat perRepublican perIndnElse];

rowlabelsDS = char('Income','EVInfo_Lot','EVInfo_LN','EnvBetter','EnvSameWorse',...
    'OwnEV','GovtCC_TooMuch', 'GovtCC_Right','GovtCC_Little','Agebelow50','Female',...
    'GradnAbove','SomeCollege','HSnBelow', 'Metropolitan','NorthEast','MidWest', ...
    'South','West','USborn','White','Black','OtherRaces','Democrat','Republican',...
     'IndnElse');

    fprintf('____________________________________________ \n')

    fprintf('           Data Summary, April 2021               \n')
    fprintf('____________________________________________ \n')
    fprintf('                 Mean     Std       \n')
    fprintf('                 _____   ______     \n')
      
    for i=1:length(MeanCounts)
        fprintf('%-10s % 7.0f   %4.2f  \n', rowlabelsDS(i,:),...
            MeanCounts(i),StdPer(i))
    end
    fprintf('___________________________________________________ \n')

disp('Response Variable')
tabulate(y)

%% Bar Graphs (Input)

% These codes are commented. The output from these codes are used in
% BarGraphEV to create the bar graphs.

% % EVInfo
% IndicEVInfo1 = find(EVInfo_LN);
% IndicEVInfo2 = find(EVInfo_Lot);
% 
% ResponseEVInfo1  = y(IndicEVInfo1);
% ResponseEVInfo2  = y(IndicEVInfo2);
% 
% % Env Better
% IndicEnv1  = find(EnvSameWorse);
% IndicEnv2  = find(EnvBetter);
% 
% ResponseEnv1 = y(IndicEnv1);
% ResponseEnv2 = y(IndicEnv2);
% 
% % EVOwner and NotOwnEV
% IndicNotOwnEV = find(NotOwnEV);
% IndicOwnEV    = find(OwnEV);
% 
% ResponseNotOwnEV = y(IndicNotOwnEV);
% ResponseOwnEV    = y(IndicOwnEV);
% 
% % Agebelow 50
% IndicAgeb50  = find(Agebelow50);
% 
% ResponseAgeb50 = y(IndicAgeb50);
% 
% % Female
% IndicFemale  = find(Female);
% %IndicMale    = find(Male);
% 
% ResponseFemale = y(IndicFemale);
% %ResponseMale   = y(IndicMale);
% 
% %Tabulation
% 
% tabulate(ResponseEVInfo1)
% tabulate(ResponseEVInfo2)
% 
% tabulate(ResponseEnv1)
% tabulate(ResponseEnv2)
% 
% tabulate(ResponseNotOwnEV)
% tabulate(ResponseOwnEV)
% 
% tabulate(ResponseAgeb50)
% 
% tabulate(ResponseFemale)
% % tabulate(ResponseMale)


%% Preparing for Regression

% y = 5-EVcar2;               % Defined earlier

x       = [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];

[~, k]  = size(x);
J       = size(unique(y),1);       % number of outcomes


% Getting an initial estimate of delta
yprob   = zeros(n,J);
for i=1:n
yprob(i,y(i))=1;
end
yprob   = sum(yprob)/n;
gam     = norminv(cumsum(yprob(1:J-1)));
deltaIn = log(gam(2:J-1) - gam(1:J-2))';   % Rough starting value.
deltaIn; 
tune    = 2;

% Prior specifications
%--------------------------------------------------------------------------
% beta
b0      = zeros(k,1);    B0 = eye(k);     
invB0b0 = B0\b0;
invB0   = inv(B0);
Btilde  = (invB0 + x'*x)\eye(k); % same as inv(invB0+x'*x)   

% delta
d0      = [0 0]';       
D0      = 0.25*eye(J-2);         % prior specification for delta

% MCMC and burnins
%------------------------------
mc      = 10000;
burn    = 0.25*mc;
nsim       = burn + mc;

% Storage matrices
%------------------------------
% z         = zeros(N,n);         % only if i want store the latent values
betaStore   = zeros(k,nsim);
deltaStore  = zeros(J-2,nsim);
acceptStore = zeros(nsim,1); 

% Initial values
%------------------------------
ytemp            = y - 1.5;
betaStore(:,1)   = (x'*x)\(x'*ytemp);
deltaStore(:,1)  = deltaIn;


% Maximum Likelihood Estimates
%--------------------------------------------------------------------------
cri0     = 1;
cri1     = 0.001;
stepsize = 1;
maxiter  = 10;
h        = 0.002;
dh       = 0.0002;
sw       = 20;
            % Optimization using the Fundamental Theorem of Calculus
             [deltamin,negsumlogL,logL,G,H, checkOutput] = OPminFundTheorem(deltaIn,y,x,betaStore(:,1),cri0,cri1,stepsize,maxiter,h,dh,sw);
            disp('Negative log-likelihood minimized');
            
            deltamin;
            Dhat = -inv(H);

%% MCMC Simulation

% For posterior ordinate: pi(betaStar|y)
%----------------------------------------
postOrdBetaStore   = zeros(nsim-burn,1);
btildeStore        = zeros(k,nsim);
%----------------------------------------

tic
wbar = waitbar(0,'MCMC in progress');

for i=2:nsim
   
    % Sampling delta
    %----------------------------------------------------------------------
     delta0   = deltaStore(:,i-1);
    
     [deltareturn,accept] = OPdrawdelta(y,x,betaStore(:,i-1),delta0,d0,D0,tune,Dhat);
     deltaStore(:,i)      = deltareturn;
     acceptStore(i,1)     = accept;


     % Sampling latent variable
     %----------------------------------------------------------------------
     [z] = OPdrawlatent(y,x,betaStore(:,i-1),deltaStore(:,i));
     
     % Sample beta
     %----------------------------------------------------------------------
     [betaStore(:,i), btildeStore(:,i)] = OPdrawbeta(z,x,invB0b0,Btilde);
        
waitbar(i/nsim);
end
close(wbar);
toc



%% Posterior Mean and Standard Deviation Computation

postMeanBeta   = mean(betaStore(:, burn+1:nsim),2);
postMeanDelta  = mean(deltaStore(:,burn+1:nsim),2);

postStdBeta    = std(betaStore(:, burn+1:nsim),0,2);
postStdDelta   = std(deltaStore(:,burn+1:nsim),0,2);

% Inefficiency factors using batch-means method
format short
[~,~,~,~, BMinfactor1] = ConvDiagnostics(betaStore(1,burn+1:nsim)');
[~,~,~,~, BMinfactor2] = ConvDiagnostics(betaStore(2,burn+1:nsim)');
[~,~,~,~, BMinfactor3] = ConvDiagnostics(betaStore(3,burn+1:nsim)');
[~,~,~,~, BMinfactor4] = ConvDiagnostics(betaStore(4,burn+1:nsim)');
[~,~,~,~, BMinfactor5] = ConvDiagnostics(betaStore(5,burn+1:nsim)');

[~,~,~,~, BMinfactor6] = ConvDiagnostics(betaStore(6,burn+1:nsim)');
[~,~,~,~, BMinfactor7] = ConvDiagnostics(betaStore(7,burn+1:nsim)');
[~,~,~,~, BMinfactor8] = ConvDiagnostics(betaStore(8,burn+1:nsim)');
[~,~,~,~, BMinfactor9] = ConvDiagnostics(betaStore(9,burn+1:nsim)');
[~,~,~,~, BMinfactor10] = ConvDiagnostics(betaStore(10,burn+1:nsim)');

[~,~,~,~, BMinfactor11] = ConvDiagnostics(betaStore(11,burn+1:nsim)');
[~,~,~,~, BMinfactor12] = ConvDiagnostics(betaStore(12,burn+1:nsim)');
[~,~,~,~, BMinfactor13] = ConvDiagnostics(betaStore(13,burn+1:nsim)');
[~,~,~,~, BMinfactor14] = ConvDiagnostics(betaStore(14,burn+1:nsim)');
[~,~,~,~, BMinfactor15] = ConvDiagnostics(betaStore(15,burn+1:nsim)');

[~,~,~,~, BMinfactor16] = ConvDiagnostics(betaStore(16,burn+1:nsim)');
[~,~,~,~, BMinfactor17] = ConvDiagnostics(betaStore(17,burn+1:nsim)');
[~,~,~,~, BMinfactor18] = ConvDiagnostics(betaStore(18,burn+1:nsim)');
[~,~,~,~, BMinfactor19] = ConvDiagnostics(betaStore(19,burn+1:nsim)');
[~,~,~,~, BMinfactor20] = ConvDiagnostics(betaStore(20,burn+1:nsim)');

[~,~,~,~, BMinfactor21] = ConvDiagnostics(deltaStore(1,burn+1:nsim)');
[~,~,~,~, BMinfactor22] = ConvDiagnostics(deltaStore(2,burn+1:nsim)');


postMeans  = [postMeanBeta; postMeanDelta]; 
postStds   = [postStdBeta; postStdDelta ];
BMInFactors  = [BMinfactor1 BMinfactor2 BMinfactor3 BMinfactor4 BMinfactor5 BMinfactor6 ...
              BMinfactor7 BMinfactor8 BMinfactor9 BMinfactor10 BMinfactor11 BMinfactor12...
              BMinfactor13 BMinfactor14 BMinfactor15 BMinfactor16 BMinfactor17 ...
              BMinfactor18 BMinfactor19 BMinfactor20 BMinfactor21 BMinfactor22]';


rowlabels = char('Intercept','EVInfo_Lot','EnvBetter','OwnEV','GovtCC_TooMuch', ...
                 'GovtCC_Right','Agebelow50','Female','GradnAbove','SomeCollege', ... 
                 'Metropolitan','MidWest', 'South','West','USborn','Black','OtherRaces', ...
                 'Democrat','IndnElse','Income/10000', 'delta1', 'delta2');

    fprintf('___________________________________________________ \n')
    fprintf('    Results from Ordinal Probit Model, April 2021          \n')
    fprintf('___________________________________________________ \n')
    fprintf('                 Mean     Std      BMIF     \n')
    fprintf('                 _____   ______   ______   \n')
      
    for i=1:length(postMeans)
        fprintf('%-10s % 7.3f   %4.3f  %7.3f \n', rowlabels(i,:),...
            postMeans(i),postStds(i),BMInFactors(i))
    end
    fprintf('___________________________________________________ \n')


%----------- Acceptance rates

disp(' Acceptance rate of (delta)')
disp( mean(acceptStore(burn+1:nsim)) )
            
% Cutpoints
%--------------------------
gammaCp       = zeros(3,1);
expdelta      = exp(postMeanDelta);
for j = 2:J-1
    gammaCp(j)   = sum(expdelta(1:j-1,1));
end

disp('    Estimated cutpoints')
disp('    gamma2    gamma3')
disp([ gammaCp(2)  gammaCp(3) ])

%% Computing the logarithm of Marginal Likelihood

% Computing posterior ordinate: pi(betaStar|y)
%---------------------------------------------
for i=1:(nsim-burn)
    postOrdBetaStore(i,1) = mvNormalPdf(postMeanBeta,btildeStore(:,burn+i),Btilde);
end

postOrdBeta = mean(postOrdBetaStore);


% Working on posterior ordinate of delta: Pi(deltastar|y,betastar)
%-----------------------------------------------------------------

mc2                = mc;
burn2              = burn;
nsim2              = burn2+mc2;        % number of iterations in the reduced run

deltaStoreProp     = zeros(J-2,nsim2);
deltaStoreRedRun   = zeros(J-2,nsim2);
%deltaStoreRedRun(:,1) = postMeanDelta;

postOrdDeltanum    = zeros(nsim2-burn2,1);
postOrdDeltaden    = zeros(nsim2-burn2,1);


tic
wbar = waitbar(0,'Reduced Run in progress');
for i = 2:nsim2  
    
 % Sampling delta
   %----------------------------------------------------------------------
    deltaProp = mvnrnd(postMeanDelta,(tune^2)*Dhat)';
    deltaStoreProp(:,i) = deltaProp;

    delta0   = deltaProp;
    
    [deltareturn,~] = OPdrawdelta(y,x,postMeanBeta,delta0,d0,D0,tune,Dhat);
    deltaStoreRedRun(:,i)      = deltareturn;

    % Sampling latent variable with beta=postMeanBeta
    %----------------------------------------------------------------------
    %[z] = OPdrawlatent(y,x,postMeanBeta,deltaStore(:,i));
    % z is not used in sampling delta, so no need to sample.
     
    % beta fixed at posterior mean
    %----------------------------------------------------------------------   
    waitbar(i/nsim2);
end
close(wbar);
toc

% Estimating the posterior ordinate: Pi(deltastar|y,betastar)
%------------------------------------------------------------
% E1alphaMH_logNum  = zeros(nsim2-burn2,1);
% E1alphaMH_logDen  = zeros(nsim2-burn2,1);

tic
wbar = waitbar(0,'Estimating Pi(deltastar|y,betastar)');

for i=1:nsim2-burn2

    % Working on the numerator
    %-------------------------
    E1alphaMH_loglikeNum = OPnegSumloglike(postMeanDelta,y,x,postMeanBeta);
    E1alphaMH_loglikeDen = OPnegSumloglike(deltaStoreRedRun(:,burn2+i),y,x,postMeanBeta);

    E1alphaMH_logNum = - E1alphaMH_loglikeNum + log(mvNormalPdf(postMeanDelta,d0,D0));
    E1alphaMH_logDen = - E1alphaMH_loglikeDen + log(mvNormalPdf(deltaStoreRedRun(:,burn2+i),d0,D0));

    E1alphaMH = min(1, exp((E1alphaMH_logNum - E1alphaMH_logDen)));
    qpdf      = mvNormalPdf(postMeanDelta, deltaStoreRedRun(:,burn2+i),(tune^2)*Dhat);

    postOrdDeltanum(i,1) = E1alphaMH*qpdf;

    % Working on the denominator
    %---------------------------
    E2alphaMH_loglikeNum = OPnegSumloglike(deltaStoreProp(:,burn2+i),y,x,postMeanBeta);
    E2alphaMH_loglikeDen = OPnegSumloglike(postMeanDelta,y,x,postMeanBeta);

    E2alphaMH_logNum = - E2alphaMH_loglikeNum + log(mvNormalPdf(deltaStoreProp(:,burn2+i),d0,D0));
    E2alphaMH_logDen = - E2alphaMH_loglikeDen + log(mvNormalPdf(postMeanDelta,d0,D0));

    postOrdDeltaden(i,1) = min(1,exp(E2alphaMH_logNum - E2alphaMH_logDen));

 
    waitbar(i/(nsim2-burn2));
end
close(wbar)

toc


postOrdDelta = mean(postOrdDeltanum)/mean(postOrdDeltaden);

priorContBeta  = mvNormalPdf(postMeanBeta,b0,B0);
priorContDelta = mvNormalPdf(postMeanDelta,d0,D0);

logLikeCont    = - OPnegSumloglike(postMeanDelta,y,x,postMeanBeta);
logPriorCont   = log(priorContBeta*priorContDelta);
logPostCont    = log(postOrdBeta*postOrdDelta);

logMarglike    = logLikeCont + logPriorCont - logPostCont

%% Computing the Covariate Effects

% Converting the deltas to gamma's    
% mc = nsim - burn;   % post burn-in draws
expdeltaBurnt  = exp(deltaStore(:,burn+1:nsim));
gammaCp        = zeros(3,mc);
% J              = 4;  

for i=1:mc
    for j = 2:J-1
        gammaCp(j,i)   = sum(expdeltaBurnt(1:j-1,i));
    end
end

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

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

% Outcome 1
%------------
oldProb1EVInfo = zeros(n,mc);
newProb1EVInfo = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Info, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1EVInfo(i,j)    = fast_normcdf(gammaCp(1,j) - EVInfoX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1EVInfo(i,j)    = fast_normcdf(gammaCp(1,j)- EVInfoX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEVInfoProb1 = newProb1EVInfo - oldProb1EVInfo;
ceEVInfoProb1   = mean(mean(diffEVInfoProb1));
sprintf('%.04f',ceEVInfoProb1)

% Outcome 2
%------------
oldProb2EVInfo    = zeros(n,mc);        
newProb2EVInfo    = zeros(n,mc);

oldComp2EVInfo    = zeros(n,mc);
newComp2EVInfo    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Info, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2EVInfo(i,j)    = fast_normcdf(gammaCp(2,j)- EVInfoX0(i,:)*betaStore(:,burn+j));
    oldProb2EVInfo(i,j)    = oldComp2EVInfo(i,j) - oldProb1EVInfo(i,j);
    % fast_normcdf(gammaCp(2,j)- EVInfox0(i,:)*betaStore(:,burn+j)) - fast_normcdf(gammaCp(1,j) - EVInfox0(i,:)*betaStore(:,burn+j));  ; 
    
    newComp2EVInfo(i,j)    = fast_normcdf(gammaCp(2,j)- EVInfoX1(i,:)*betaStore(:,burn+j));
    newProb2EVInfo(i,j)    = newComp2EVInfo(i,j) - newProb1EVInfo(i,j);
    %fast_normcdf(gammaCp(2,j)- EVInfox1(i,:)*betaStore(:,burn+j)) - fast_normcdf(gammaCp(1,j)- EVInfox1(i,:)*betaStore(:,burn+j));;
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEVInfoProb2      = newProb2EVInfo - oldProb2EVInfo;
ceEVInfoProb2   = mean(mean(diffEVInfoProb2));
sprintf('%.04f',ceEVInfoProb2)

% Outcome 3
%------------
oldProb3EVInfo    = zeros(n,mc);        
newProb3EVInfo    = zeros(n,mc);

oldComp3EVInfo    = zeros(n,mc);
newComp3EVInfo    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Info, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3EVInfo(i,j)    = fast_normcdf(gammaCp(3,j)- EVInfoX0(i,:)*betaStore(:,burn+j));
    oldProb3EVInfo(i,j)    = oldComp3EVInfo(i,j) - oldComp2EVInfo(i,j);
    % fast_normcdf(gammaCp(3,j)- EVInfoX0(i,:)*betaStore(:,burn+j)) - fast_normcdf(gammaCp(2,j)- EVInfoX0(i,:)*betaStore(:,burn+j));

    newComp3EVInfo(i,j)    = fast_normcdf(gammaCp(3,j)- EVInfoX1(i,:)*betaStore(:,burn+j));
    newProb3EVInfo(i,j)    = newComp3EVInfo(i,j) - newComp2EVInfo(i,j);
    %fast_normcdf(gammaCp(3,j)- EVInfoX1(i,:)*betaStore(:,burn+j)) - fast_normcdf(gammaCp(2,j)- EVInfoX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEVInfoProb3 = newProb3EVInfo - oldProb3EVInfo;
ceEVInfoProb3   = mean(mean(diffEVInfoProb3));
sprintf('%.04f',ceEVInfoProb3)

% Outcome 4
%------------
oldProb4EVInfo    = zeros(n,mc);        
newProb4EVInfo    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Info, Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4EVInfo(i,j)    = 1 - oldComp3EVInfo(i,j);
    %1 - fast_normcdf(gammaCp(3,j)- EVInfoX0(i,:)*betaStore(:,burn+j));
    
    newProb4EVInfo(i,j)    = 1 - newComp3EVInfo(i,j);
    %1 - normcdf(gammaCp(3,j)- modx(i,:)*betaBurnt(:,j), mu, sigma);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEVInfoProb4 = newProb4EVInfo - oldProb4EVInfo;
ceEVInfoProb4   = mean(mean(diffEVInfoProb4));
sprintf('%.04f',ceEVInfoProb4)

disp('CE, EVInfo')
disp([ceEVInfoProb1 ceEVInfoProb2 ceEVInfoProb3 ceEVInfoProb4])

%--------------------------------------------------------------------------
clearvars oldProb1EVInfo oldComp2EVInfo oldProb2EVInfo oldComp3EVInfo oldProb3EVInfo oldProb4EVInfo
clearvars newProb1EVInfo newComp2EVInfo newProb2EVInfo newComp3EVInfo newProb3EVInfo newProb4EVInfo
clearvars diffEVInfoProb1 diffEVInfoProb2 diffEVInfoProb3 diffEVInfoProb4


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

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

% Outcome 1
%------------
oldProb1EnvImpact = zeros(n,mc);
newProb1EnvImpact = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Env Impact, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1EnvImpact(i,j)    = fast_normcdf(gammaCp(1,j) - EnvImpactX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1EnvImpact(i,j)    = fast_normcdf(gammaCp(1,j)- EnvImpactX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEnvImpactProb1 = newProb1EnvImpact - oldProb1EnvImpact;
ceEnvImpactProb1   = mean(mean(diffEnvImpactProb1));
sprintf('%.04f',ceEnvImpactProb1)

% Outcome 2
%------------
oldProb2EnvImpact    = zeros(n,mc);        
newProb2EnvImpact    = zeros(n,mc);

oldComp2EnvImpact    = zeros(n,mc);
newComp2EnvImpact    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Env Impact, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2EnvImpact(i,j)    = fast_normcdf(gammaCp(2,j)- EnvImpactX0(i,:)*betaStore(:,burn+j));
    oldProb2EnvImpact(i,j)    = oldComp2EnvImpact(i,j) - oldProb1EnvImpact(i,j);
    
    newComp2EnvImpact(i,j)    = fast_normcdf(gammaCp(2,j)- EnvImpactX1(i,:)*betaStore(:,burn+j));
    newProb2EnvImpact(i,j)    = newComp2EnvImpact(i,j) - newProb1EnvImpact(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEnvImpactProb2      = newProb2EnvImpact - oldProb2EnvImpact;
ceEnvImpactProb2   = mean(mean(diffEnvImpactProb2));
sprintf('%.04f',ceEnvImpactProb2)

% Outcome 3
%------------
oldProb3EnvImpact    = zeros(n,mc);        
newProb3EnvImpact    = zeros(n,mc);

oldComp3EnvImpact    = zeros(n,mc);
newComp3EnvImpact    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Env Impact, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3EnvImpact(i,j)    = fast_normcdf(gammaCp(3,j)- EnvImpactX0(i,:)*betaStore(:,burn+j));
    oldProb3EnvImpact(i,j)    = oldComp3EnvImpact(i,j) - oldComp2EnvImpact(i,j);

    newComp3EnvImpact(i,j)    = fast_normcdf(gammaCp(3,j)- EnvImpactX1(i,:)*betaStore(:,burn+j));
    newProb3EnvImpact(i,j)    = newComp3EnvImpact(i,j) - newComp2EnvImpact(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEnvImpactProb3 = newProb3EnvImpact - oldProb3EnvImpact;
ceEnvImpactProb3   = mean(mean(diffEnvImpactProb3));
sprintf('%.04f',ceEnvImpactProb3)

% Outcome 4
%------------
oldProb4EnvImpact    = zeros(n,mc);        
newProb4EnvImpact    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Env Impact, Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4EnvImpact(i,j)    = 1 - oldComp3EnvImpact(i,j);
    
    newProb4EnvImpact(i,j)    = 1 - newComp3EnvImpact(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffEnvImpactProb4 = newProb4EnvImpact - oldProb4EnvImpact;
ceEnvImpactProb4   = mean(mean(diffEnvImpactProb4));
sprintf('%.04f',ceEnvImpactProb4)

disp('CE, EnvImpact')
disp([ceEnvImpactProb1 ceEnvImpactProb2 ceEnvImpactProb3 ceEnvImpactProb4])

%--------------------------------------------------------------------------
clearvars oldProb1EnvImpact oldComp2EnvImpact oldProb2EnvImpact oldComp3EnvImpact oldProb3EnvImpact oldProb4EnvImpact
clearvars newProb1EnvImpact newComp2EnvImpact newProb2EnvImpact newComp3EnvImpact newProb3EnvImpact newProb4EnvImpact
clearvars diffEnvImpactProb1 diffEnvImpactProb2 diffEnvImpactProb3 diffEnvImpactProb4

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

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

% Outcome 1
%------------
oldProb1OwnEV = zeros(n,mc);
newProb1OwnEV = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Owner, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1OwnEV(i,j)    = fast_normcdf(gammaCp(1,j) - OwnEVX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1OwnEV(i,j)    = fast_normcdf(gammaCp(1,j)- OwnEVX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOwnEVProb1 = newProb1OwnEV - oldProb1OwnEV;
ceOwnEVProb1   = mean(mean(diffOwnEVProb1));
sprintf('%.04f',ceOwnEVProb1)

% Outcome 2
%------------
oldProb2OwnEV    = zeros(n,mc);        
newProb2OwnEV    = zeros(n,mc);

oldComp2OwnEV    = zeros(n,mc);
newComp2OwnEV    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Owner, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2OwnEV(i,j)    = fast_normcdf(gammaCp(2,j)- OwnEVX0(i,:)*betaStore(:,burn+j));
    oldProb2OwnEV(i,j)    = oldComp2OwnEV(i,j) - oldProb1OwnEV(i,j);
    
    newComp2OwnEV(i,j)    = fast_normcdf(gammaCp(2,j)- OwnEVX1(i,:)*betaStore(:,burn+j));
    newProb2OwnEV(i,j)    = newComp2OwnEV(i,j) - newProb1OwnEV(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOwnEVProb2      = newProb2OwnEV - oldProb2OwnEV;
ceOwnEVProb2   = mean(mean(diffOwnEVProb2));
sprintf('%.04f',ceOwnEVProb2)

% Outcome 3
%------------
oldProb3OwnEV    = zeros(n,mc);        
newProb3OwnEV    = zeros(n,mc);

oldComp3OwnEV    = zeros(n,mc);
newComp3OwnEV    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Owner, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3OwnEV(i,j)    = fast_normcdf(gammaCp(3,j)- OwnEVX0(i,:)*betaStore(:,burn+j));
    oldProb3OwnEV(i,j)    = oldComp3OwnEV(i,j) - oldComp2OwnEV(i,j);

    newComp3OwnEV(i,j)    = fast_normcdf(gammaCp(3,j)- OwnEVX1(i,:)*betaStore(:,burn+j));
    newProb3OwnEV(i,j)    = newComp3OwnEV(i,j) - newComp2OwnEV(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOwnEVProb3 = newProb3OwnEV - oldProb3OwnEV;
ceOwnEVProb3   = mean(mean(diffOwnEVProb3));
sprintf('%.04f',ceOwnEVProb3)

% Outcome 4
%------------
oldProb4OwnEV    = zeros(n,mc);        
newProb4OwnEV    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, EV Owner, Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4OwnEV(i,j)    = 1 - oldComp3OwnEV(i,j);
    
    newProb4OwnEV(i,j)    = 1 - newComp3OwnEV(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOwnEVProb4 = newProb4OwnEV - oldProb4OwnEV;
ceOwnEVProb4   = mean(mean(diffOwnEVProb4));
sprintf('%.04f',ceOwnEVProb4)

disp('CE, OwnEV')
disp([ceOwnEVProb1 ceOwnEVProb2 ceOwnEVProb3 ceOwnEVProb4])

%--------------------------------------------------------------------------
clearvars oldProb1OwnEV oldComp2OwnEV oldProb2OwnEV oldComp3OwnEV oldProb3OwnEV oldProb4OwnEV
clearvars newProb1OwnEV newComp2OwnEV newProb2OwnEV newComp3OwnEV newProb3OwnEV newProb4OwnEV
clearvars diffOwnEVProb1 diffOwnEVProb2 diffOwnEVProb3 diffOwnEVProb4

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

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

% Outcome 1
%------------
oldProb1GACCTM = zeros(n,mc);
newProb1GACCTM = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC(TM), Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1GACCTM(i,j)    = fast_normcdf(gammaCp(1,j) - GACCTMX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1GACCTM(i,j)    = fast_normcdf(gammaCp(1,j)- GACCTMX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCTMProb1 = newProb1GACCTM - oldProb1GACCTM;
ceGACCTMProb1   = mean(mean(diffGACCTMProb1));
sprintf('%.04f',ceGACCTMProb1)

% Outcome 2
%------------
oldProb2GACCTM    = zeros(n,mc);        
newProb2GACCTM    = zeros(n,mc);

oldComp2GACCTM    = zeros(n,mc);
newComp2GACCTM    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC (TM), Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2GACCTM(i,j)    = fast_normcdf(gammaCp(2,j)- GACCTMX0(i,:)*betaStore(:,burn+j));
    oldProb2GACCTM(i,j)    = oldComp2GACCTM(i,j) - oldProb1GACCTM(i,j);
    
    newComp2GACCTM(i,j)    = fast_normcdf(gammaCp(2,j)- GACCTMX1(i,:)*betaStore(:,burn+j));
    newProb2GACCTM(i,j)    = newComp2GACCTM(i,j) - newProb1GACCTM(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCTMProb2 = newProb2GACCTM - oldProb2GACCTM;
ceGACCTMProb2   = mean(mean(diffGACCTMProb2));
sprintf('%.04f',ceGACCTMProb2)

% Outcome 3
%------------
oldProb3GACCTM    = zeros(n,mc);        
newProb3GACCTM    = zeros(n,mc);

oldComp3GACCTM    = zeros(n,mc);
newComp3GACCTM    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC (TM), Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3GACCTM(i,j)    = fast_normcdf(gammaCp(3,j)- GACCTMX0(i,:)*betaStore(:,burn+j));
    oldProb3GACCTM(i,j)    = oldComp3GACCTM(i,j) - oldComp2GACCTM(i,j);

    newComp3GACCTM(i,j)    = fast_normcdf(gammaCp(3,j)- GACCTMX1(i,:)*betaStore(:,burn+j));
    newProb3GACCTM(i,j)    = newComp3GACCTM(i,j) - newComp2GACCTM(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCTMProb3 = newProb3GACCTM - oldProb3GACCTM;
ceGACCTMProb3   = mean(mean(diffGACCTMProb3));
sprintf('%.04f',ceGACCTMProb3)

% Outcome 4
%------------
oldProb4GACCTM    = zeros(n,mc);        
newProb4GACCTM    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC (TM), Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4GACCTM(i,j)    = 1 - oldComp3GACCTM(i,j);
    
    newProb4GACCTM(i,j)    = 1 - newComp3GACCTM(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCTMProb4 = newProb4GACCTM - oldProb4GACCTM;
ceGACCTMProb4   = mean(mean(diffGACCTMProb4));
sprintf('%.04f',ceGACCTMProb4)

disp('CE, GACCTM')
disp([ceGACCTMProb1 ceGACCTMProb2 ceGACCTMProb3 ceGACCTMProb4])

%--------------------------------------------------------------------------
clearvars oldProb1GACCTM oldComp2GACCTM oldProb2GACCTM oldComp3GACCTM oldProb3GACCTM oldProb4GACCTM
clearvars newProb1GACCTM newComp2GACCTM newProb2GACCTM newComp3GACCTM newProb3GACCTM newProb4GACCTM
clearvars diffGACCTMProb1 diffGACCTMProb2 diffGACCTMProb3 diffGACCTMProb4

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

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

% Outcome 1
%------------
oldProb1GACCRA = zeros(n,mc);
newProb1GACCRA = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC(RA), Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1GACCRA(i,j)    = fast_normcdf(gammaCp(1,j) - GACCRAX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1GACCRA(i,j)    = fast_normcdf(gammaCp(1,j)- GACCRAX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCRAProb1 = newProb1GACCRA - oldProb1GACCRA;
ceGACCRAProb1   = mean(mean(diffGACCRAProb1));
sprintf('%.04f',ceGACCRAProb1)

% Outcome 2
%------------
oldProb2GACCRA    = zeros(n,mc);        
newProb2GACCRA    = zeros(n,mc);

oldComp2GACCRA    = zeros(n,mc);
newComp2GACCRA    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC (RA), Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2GACCRA(i,j)    = fast_normcdf(gammaCp(2,j)- GACCRAX0(i,:)*betaStore(:,burn+j));
    oldProb2GACCRA(i,j)    = oldComp2GACCRA(i,j) - oldProb1GACCRA(i,j);
    
    newComp2GACCRA(i,j)    = fast_normcdf(gammaCp(2,j)- GACCRAX1(i,:)*betaStore(:,burn+j));
    newProb2GACCRA(i,j)    = newComp2GACCRA(i,j) - newProb1GACCRA(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCRAProb2      = newProb2GACCRA - oldProb2GACCRA;
ceGACCRAProb2   = mean(mean(diffGACCRAProb2));
sprintf('%.04f',ceGACCRAProb2)

% Outcome 3
%------------
oldProb3GACCRA    = zeros(n,mc);        
newProb3GACCRA    = zeros(n,mc);

oldComp3GACCRA    = zeros(n,mc);
newComp3GACCRA    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC (RA), Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3GACCRA(i,j)    = fast_normcdf(gammaCp(3,j)- GACCRAX0(i,:)*betaStore(:,burn+j));
    oldProb3GACCRA(i,j)    = oldComp3GACCRA(i,j) - oldComp2GACCRA(i,j);

    newComp3GACCRA(i,j)    = fast_normcdf(gammaCp(3,j)- GACCRAX1(i,:)*betaStore(:,burn+j));
    newProb3GACCRA(i,j)    = newComp3GACCRA(i,j) - newComp2GACCRA(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCRAProb3 = newProb3GACCRA - oldProb3GACCRA;
ceGACCRAProb3   = mean(mean(diffGACCRAProb3));
sprintf('%.04f',ceGACCRAProb3)

% Outcome 4
%------------
oldProb4GACCRA    = zeros(n,mc);        
newProb4GACCRA    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GACC (TM), Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4GACCRA(i,j)    = 1 - oldComp3GACCRA(i,j);
    
    newProb4GACCRA(i,j)    = 1 - newComp3GACCRA(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGACCRAProb4 = newProb4GACCRA - oldProb4GACCRA;
ceGACCRAProb4   = mean(mean(diffGACCRAProb4));
sprintf('%.04f',ceGACCRAProb4)

disp('CE, GACCRA')
disp([ceGACCRAProb1 ceGACCRAProb2 ceGACCRAProb3 ceGACCRAProb4])

%--------------------------------------------------------------------------
clearvars oldProb1GACCRA oldComp2GACCRA oldProb2GACCRA oldComp3GACCRA oldProb3GACCRA oldProb4GACCRA
clearvars newProb1GACCRA newComp2GACCRA newProb2GACCRA newComp3GACCRA newProb3GACCRA newProb4GACCRA
clearvars diffGACCRAProb1 diffGACCRAProb2 diffGACCRAProb3 diffGACCRAProb4

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

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

% Outcome 1
%------------
oldProb1Ageb50 = zeros(n,mc);
newProb1Ageb50 = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Age<50, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Ageb50(i,j)    = fast_normcdf(gammaCp(1,j) - Ageb50X0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Ageb50(i,j)    = fast_normcdf(gammaCp(1,j)- Ageb50X1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffAgeb50Prob1 = newProb1Ageb50 - oldProb1Ageb50;
ceAgeb50Prob1   = mean(mean(diffAgeb50Prob1));
sprintf('%.04f',ceAgeb50Prob1)

% Outcome 2
%------------
oldProb2Ageb50    = zeros(n,mc);        
newProb2Ageb50    = zeros(n,mc);

oldComp2Ageb50    = zeros(n,mc);
newComp2Ageb50    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Age<50, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Ageb50(i,j)    = fast_normcdf(gammaCp(2,j)- Ageb50X0(i,:)*betaStore(:,burn+j));
    oldProb2Ageb50(i,j)    = oldComp2Ageb50(i,j) - oldProb1Ageb50(i,j);
    
    newComp2Ageb50(i,j)    = fast_normcdf(gammaCp(2,j)- Ageb50X1(i,:)*betaStore(:,burn+j));
    newProb2Ageb50(i,j)    = newComp2Ageb50(i,j) - newProb1Ageb50(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffAgeb50Prob2 = newProb2Ageb50 - oldProb2Ageb50;
ceAgeb50Prob2   = mean(mean(diffAgeb50Prob2));
sprintf('%.04f',ceAgeb50Prob2)

% Outcome 3
%------------
oldProb3Ageb50    = zeros(n,mc);        
newProb3Ageb50    = zeros(n,mc);

oldComp3Ageb50    = zeros(n,mc);
newComp3Ageb50    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Age<50, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Ageb50(i,j)    = fast_normcdf(gammaCp(3,j)- Ageb50X0(i,:)*betaStore(:,burn+j));
    oldProb3Ageb50(i,j)    = oldComp3Ageb50(i,j) - oldComp2Ageb50(i,j);

    newComp3Ageb50(i,j)    = fast_normcdf(gammaCp(3,j)- Ageb50X1(i,:)*betaStore(:,burn+j));
    newProb3Ageb50(i,j)    = newComp3Ageb50(i,j) - newComp2Ageb50(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffAgeb50Prob3 = newProb3Ageb50 - oldProb3Ageb50;
ceAgeb50Prob3   = mean(mean(diffAgeb50Prob3));
sprintf('%.04f',ceAgeb50Prob3)

% Outcome 4
%------------
oldProb4Ageb50    = zeros(n,mc);        
newProb4Ageb50    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Age<50, Outcome 4');

for j =1:mc
    for i=1:n
    oldProb4Ageb50(i,j)    = 1 - oldComp3Ageb50(i,j);
    
    newProb4Ageb50(i,j)    = 1 - newComp3Ageb50(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffAgeb50Prob4 = newProb4Ageb50 - oldProb4Ageb50;
ceAgeb50Prob4   = mean(mean(diffAgeb50Prob4));
sprintf('%.04f',ceAgeb50Prob4)

disp('CE, Age<50')
disp([ceAgeb50Prob1 ceAgeb50Prob2 ceAgeb50Prob3 ceAgeb50Prob4])

%--------------------------------------------------------------------------
clearvars oldProb1Ageb50 oldComp2Ageb50 oldProb2Ageb50 oldComp3Ageb50 oldProb3Ageb50 oldProb4Ageb50
clearvars newProb1Ageb50 newComp2Ageb50 newProb2Ageb50 newComp3Ageb50 newProb3Ageb50 newProb4Ageb50
clearvars diffAgeb50Prob1 diffAgeb50Prob2 diffAgeb50Prob3 diffAgeb50Prob4

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

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

% Outcome 1
%------------
oldProb1Female = zeros(n,mc);
newProb1Female = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Female, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Female(i,j)    = fast_normcdf(gammaCp(1,j) - FemaleX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Female(i,j)    = fast_normcdf(gammaCp(1,j)- FemaleX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffFemaleProb1 = newProb1Female - oldProb1Female;
ceFemaleProb1   = mean(mean(diffFemaleProb1));
sprintf('%.04f',ceFemaleProb1)

% Outcome 2
%------------
oldProb2Female    = zeros(n,mc);        
newProb2Female    = zeros(n,mc);

oldComp2Female    = zeros(n,mc);
newComp2Female    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Female, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Female(i,j)    = fast_normcdf(gammaCp(2,j)- FemaleX0(i,:)*betaStore(:,burn+j));
    oldProb2Female(i,j)    = oldComp2Female(i,j) - oldProb1Female(i,j);
    
    newComp2Female(i,j)    = fast_normcdf(gammaCp(2,j)- FemaleX1(i,:)*betaStore(:,burn+j));
    newProb2Female(i,j)    = newComp2Female(i,j) - newProb1Female(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffFemaleProb2 = newProb2Female - oldProb2Female;
ceFemaleProb2   = mean(mean(diffFemaleProb2));
sprintf('%.04f',ceFemaleProb2)

% Outcome 3
%------------
oldProb3Female    = zeros(n,mc);        
newProb3Female    = zeros(n,mc);

oldComp3Female    = zeros(n,mc);
newComp3Female    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Female, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Female(i,j)    = fast_normcdf(gammaCp(3,j)- FemaleX0(i,:)*betaStore(:,burn+j));
    oldProb3Female(i,j)    = oldComp3Female(i,j) - oldComp2Female(i,j);

    newComp3Female(i,j)    = fast_normcdf(gammaCp(3,j)- FemaleX1(i,:)*betaStore(:,burn+j));
    newProb3Female(i,j)    = newComp3Female(i,j) - newComp2Female(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffFemaleProb3 = newProb3Female - oldProb3Female;
ceFemaleProb3   = mean(mean(diffFemaleProb3));
sprintf('%.04f',ceFemaleProb3)

% Outcome 4
%------------
oldProb4Female    = zeros(n,mc);        
newProb4Female    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4Female(i,j)    = 1 - oldComp3Female(i,j);
    
    newProb4Female(i,j)    = 1 - newComp3Female(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffFemaleProb4 = newProb4Female - oldProb4Female;
ceFemaleProb4   = mean(mean(diffFemaleProb4));
sprintf('%.04f',ceFemaleProb4)

disp('CE, Female')
disp([ceFemaleProb1 ceFemaleProb2 ceFemaleProb3 ceFemaleProb4])

%--------------------------------------------------------------------------
clearvars oldProb1Female oldComp2Female oldProb2Female oldComp3Female oldProb3Female oldProb4Female
clearvars newProb1Female newComp2Female newProb2Female newComp3Female newProb3Female newProb4Female
clearvars diffFemaleProb1 diffFemaleProb2 diffFemaleProb3 diffFemaleProb4


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

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

% Outcome 1
%------------
oldProb1GradnAbove = zeros(n,mc);
newProb1GradnAbove = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GradnAbove, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1GradnAbove(i,j)    = fast_normcdf(gammaCp(1,j) - GradnAboveX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1GradnAbove(i,j)    = fast_normcdf(gammaCp(1,j)- GradnAboveX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGradnAboveProb1 = newProb1GradnAbove - oldProb1GradnAbove;
ceGradnAboveProb1   = mean(mean(diffGradnAboveProb1));
sprintf('%.04f',ceGradnAboveProb1)

% Outcome 2
%------------
oldProb2GradnAbove    = zeros(n,mc);        
newProb2GradnAbove    = zeros(n,mc);

oldComp2GradnAbove    = zeros(n,mc);
newComp2GradnAbove    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GradnAbove, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2GradnAbove(i,j)    = fast_normcdf(gammaCp(2,j)- GradnAboveX0(i,:)*betaStore(:,burn+j));
    oldProb2GradnAbove(i,j)    = oldComp2GradnAbove(i,j) - oldProb1GradnAbove(i,j);
    
    newComp2GradnAbove(i,j)    = fast_normcdf(gammaCp(2,j)- GradnAboveX1(i,:)*betaStore(:,burn+j));
    newProb2GradnAbove(i,j)    = newComp2GradnAbove(i,j) - newProb1GradnAbove(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGradnAboveProb2 = newProb2GradnAbove - oldProb2GradnAbove;
ceGradnAboveProb2   = mean(mean(diffGradnAboveProb2));
sprintf('%.04f',ceGradnAboveProb2)

% Outcome 3
%------------
oldProb3GradnAbove    = zeros(n,mc);        
newProb3GradnAbove    = zeros(n,mc);

oldComp3GradnAbove    = zeros(n,mc);
newComp3GradnAbove    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, GradnAbove, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3GradnAbove(i,j)    = fast_normcdf(gammaCp(3,j)- GradnAboveX0(i,:)*betaStore(:,burn+j));
    oldProb3GradnAbove(i,j)    = oldComp3GradnAbove(i,j) - oldComp2GradnAbove(i,j);

    newComp3GradnAbove(i,j)    = fast_normcdf(gammaCp(3,j)- GradnAboveX1(i,:)*betaStore(:,burn+j));
    newProb3GradnAbove(i,j)    = newComp3GradnAbove(i,j) - newComp2GradnAbove(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGradnAboveProb3 = newProb3GradnAbove - oldProb3GradnAbove;
ceGradnAboveProb3   = mean(mean(diffGradnAboveProb3));
sprintf('%.04f',ceGradnAboveProb3)

% Outcome 4
%------------
oldProb4GradnAbove    = zeros(n,mc);        
newProb4GradnAbove    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4GradnAbove(i,j)    = 1 - oldComp3GradnAbove(i,j);
    
    newProb4GradnAbove(i,j)    = 1 - newComp3GradnAbove(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffGradnAboveProb4 = newProb4GradnAbove - oldProb4GradnAbove;
ceGradnAboveProb4   = mean(mean(diffGradnAboveProb4));
sprintf('%.04f',ceGradnAboveProb4)

disp('CE, GradnAbove')
disp([ceGradnAboveProb1 ceGradnAboveProb2 ceGradnAboveProb3 ceGradnAboveProb4])

%--------------------------------------------------------------------------
clearvars oldProb1GradnAbove oldComp2GradnAbove oldProb2GradnAbove oldComp3GradnAbove oldProb3GradnAbove oldProb4GradnAbove
clearvars newProb1GradnAbove newComp2GradnAbove newProb2GradnAbove newComp3GradnAbove newProb3GradnAbove newProb4GradnAbove
clearvars diffGradnAboveProb1 diffGradnAboveProb2 diffGradnAboveProb3 diffGradnAboveProb4

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

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

% Outcome 1
%------------
oldProb1SomeCollege = zeros(n,mc);
newProb1SomeCollege = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, SomeCollege, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1SomeCollege(i,j)    = fast_normcdf(gammaCp(1,j) - SomeCollegeX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1SomeCollege(i,j)    = fast_normcdf(gammaCp(1,j) - SomeCollegeX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSomeCollegeProb1 = newProb1SomeCollege - oldProb1SomeCollege;
ceSomeCollegeProb1   = mean(mean(diffSomeCollegeProb1));
sprintf('%.04f',ceSomeCollegeProb1)

% Outcome 2
%------------
oldProb2SomeCollege    = zeros(n,mc);        
newProb2SomeCollege    = zeros(n,mc);

oldComp2SomeCollege    = zeros(n,mc);
newComp2SomeCollege    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, SomeCollege, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2SomeCollege(i,j)    = fast_normcdf(gammaCp(2,j)- SomeCollegeX0(i,:)*betaStore(:,burn+j));
    oldProb2SomeCollege(i,j)    = oldComp2SomeCollege(i,j) - oldProb1SomeCollege(i,j);
    
    newComp2SomeCollege(i,j)    = fast_normcdf(gammaCp(2,j)- SomeCollegeX1(i,:)*betaStore(:,burn+j));
    newProb2SomeCollege(i,j)    = newComp2SomeCollege(i,j) - newProb1SomeCollege(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSomeCollegeProb2 = newProb2SomeCollege - oldProb2SomeCollege;
ceSomeCollegeProb2   = mean(mean(diffSomeCollegeProb2));
sprintf('%.04f',ceSomeCollegeProb2)

% Outcome 3
%------------
oldProb3SomeCollege    = zeros(n,mc);        
newProb3SomeCollege    = zeros(n,mc);

oldComp3SomeCollege    = zeros(n,mc);
newComp3SomeCollege    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, SomeCollege, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3SomeCollege(i,j)    = fast_normcdf(gammaCp(3,j)- SomeCollegeX0(i,:)*betaStore(:,burn+j));
    oldProb3SomeCollege(i,j)    = oldComp3SomeCollege(i,j) - oldComp2SomeCollege(i,j);

    newComp3SomeCollege(i,j)    = fast_normcdf(gammaCp(3,j)- SomeCollegeX1(i,:)*betaStore(:,burn+j));
    newProb3SomeCollege(i,j)    = newComp3SomeCollege(i,j) - newComp2SomeCollege(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSomeCollegeProb3 = newProb3SomeCollege - oldProb3SomeCollege;
ceSomeCollegeProb3   = mean(mean(diffSomeCollegeProb3));
sprintf('%.04f',ceSomeCollegeProb3)

% Outcome 4
%------------
oldProb4SomeCollege    = zeros(n,mc);        
newProb4SomeCollege    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4SomeCollege(i,j)    = 1 - oldComp3SomeCollege(i,j);
    
    newProb4SomeCollege(i,j)    = 1 - newComp3SomeCollege(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSomeCollegeProb4 = newProb4SomeCollege - oldProb4SomeCollege;
ceSomeCollegeProb4   = mean(mean(diffSomeCollegeProb4));
sprintf('%.04f',ceSomeCollegeProb4)

disp('CE, SomeCollege')
disp([ceSomeCollegeProb1 ceSomeCollegeProb2 ceSomeCollegeProb3 ceSomeCollegeProb4])

%--------------------------------------------------------------------------
clearvars oldProb1SomeCollege oldComp2SomeCollege oldProb2SomeCollege oldComp3SomeCollege oldProb3SomeCollege oldProb4SomeCollege
clearvars newProb1SomeCollege newComp2SomeCollege newProb2SomeCollege newComp3SomeCollege newProb3SomeCollege newProb4SomeCollege
clearvars diffSomeCollegeProb1 diffSomeCollegeProb2 diffSomeCollegeProb3 diffSomeCollegeProb4

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

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

% Outcome 1
%------------
oldProb1Metropolitan = zeros(n,mc);
newProb1Metropolitan = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Metropolitan, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Metropolitan(i,j)    = fast_normcdf(gammaCp(1,j) - MetropolitanX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Metropolitan(i,j)    = fast_normcdf(gammaCp(1,j) - MetropolitanX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMetropolitanProb1 = newProb1Metropolitan - oldProb1Metropolitan;
ceMetropolitanProb1   = mean(mean(diffMetropolitanProb1));
sprintf('%.04f',ceMetropolitanProb1)

% Outcome 2
%------------
oldProb2Metropolitan    = zeros(n,mc);        
newProb2Metropolitan    = zeros(n,mc);

oldComp2Metropolitan    = zeros(n,mc);
newComp2Metropolitan    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Metropolitan, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Metropolitan(i,j)    = fast_normcdf(gammaCp(2,j)- MetropolitanX0(i,:)*betaStore(:,burn+j));
    oldProb2Metropolitan(i,j)    = oldComp2Metropolitan(i,j) - oldProb1Metropolitan(i,j);
    
    newComp2Metropolitan(i,j)    = fast_normcdf(gammaCp(2,j)- MetropolitanX1(i,:)*betaStore(:,burn+j));
    newProb2Metropolitan(i,j)    = newComp2Metropolitan(i,j) - newProb1Metropolitan(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMetropolitanProb2 = newProb2Metropolitan - oldProb2Metropolitan;
ceMetropolitanProb2   = mean(mean(diffMetropolitanProb2));
sprintf('%.04f',ceMetropolitanProb2)

% Outcome 3
%------------
oldProb3Metropolitan    = zeros(n,mc);        
newProb3Metropolitan    = zeros(n,mc);

oldComp3Metropolitan    = zeros(n,mc);
newComp3Metropolitan    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Metropolitan, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Metropolitan(i,j)    = fast_normcdf(gammaCp(3,j)- MetropolitanX0(i,:)*betaStore(:,burn+j));
    oldProb3Metropolitan(i,j)    = oldComp3Metropolitan(i,j) - oldComp2Metropolitan(i,j);

    newComp3Metropolitan(i,j)    = fast_normcdf(gammaCp(3,j)- MetropolitanX1(i,:)*betaStore(:,burn+j));
    newProb3Metropolitan(i,j)    = newComp3Metropolitan(i,j) - newComp2Metropolitan(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMetropolitanProb3 = newProb3Metropolitan - oldProb3Metropolitan;
ceMetropolitanProb3   = mean(mean(diffMetropolitanProb3));
sprintf('%.04f',ceMetropolitanProb3)

% Outcome 4
%------------
oldProb4Metropolitan    = zeros(n,mc);        
newProb4Metropolitan    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4Metropolitan(i,j)    = 1 - oldComp3Metropolitan(i,j);
    
    newProb4Metropolitan(i,j)    = 1 - newComp3Metropolitan(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMetropolitanProb4 = newProb4Metropolitan - oldProb4Metropolitan;
ceMetropolitanProb4   = mean(mean(diffMetropolitanProb4));
sprintf('%.04f',ceMetropolitanProb4)

disp('CE, Metropolitan')
disp([ceMetropolitanProb1 ceMetropolitanProb2 ceMetropolitanProb3 ceMetropolitanProb4])

%--------------------------------------------------------------------------
clearvars oldProb1Metropolitan oldComp2Metropolitan oldProb2Metropolitan oldComp3Metropolitan oldProb3Metropolitan oldProb4Metropolitan
clearvars newProb1Metropolitan newComp2Metropolitan newProb2Metropolitan newComp3Metropolitan newProb3Metropolitan newProb4Metropolitan
clearvars diffMetropolitanProb1 diffMetropolitanProb2 diffMetropolitanProb3 diffMetropolitanProb4

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

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

% Outcome 1
%------------
oldProb1Midwest = zeros(n,mc);
newProb1Midwest = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Midwest, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Midwest(i,j)    = fast_normcdf(gammaCp(1,j) - MidWestX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Midwest(i,j)    = fast_normcdf(gammaCp(1,j) - MidWestX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMidwestProb1 = newProb1Midwest - oldProb1Midwest;
ceMidwestProb1   = mean(mean(diffMidwestProb1));
sprintf('%.04f',ceMidwestProb1)

% Outcome 2
%------------
oldProb2Midwest    = zeros(n,mc);        
newProb2Midwest    = zeros(n,mc);

oldComp2Midwest    = zeros(n,mc);
newComp2Midwest    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Midwest, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Midwest(i,j)    = fast_normcdf(gammaCp(2,j)- MidWestX0(i,:)*betaStore(:,burn+j));
    oldProb2Midwest(i,j)    = oldComp2Midwest(i,j) - oldProb1Midwest(i,j);
    
    newComp2Midwest(i,j)    = fast_normcdf(gammaCp(2,j)- MidWestX1(i,:)*betaStore(:,burn+j));
    newProb2Midwest(i,j)    = newComp2Midwest(i,j) - newProb1Midwest(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMidwestProb2 = newProb2Midwest - oldProb2Midwest;
ceMidwestProb2   = mean(mean(diffMidwestProb2));
sprintf('%.04f',ceMidwestProb2)

% Outcome 3
%------------
oldProb3Midwest    = zeros(n,mc);        
newProb3Midwest    = zeros(n,mc);

oldComp3Midwest    = zeros(n,mc);
newComp3Midwest    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Midwest, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Midwest(i,j)    = fast_normcdf(gammaCp(3,j)- MidWestX0(i,:)*betaStore(:,burn+j));
    oldProb3Midwest(i,j)    = oldComp3Midwest(i,j) - oldComp2Midwest(i,j);

    newComp3Midwest(i,j)    = fast_normcdf(gammaCp(3,j)- MidWestX1(i,:)*betaStore(:,burn+j));
    newProb3Midwest(i,j)    = newComp3Midwest(i,j) - newComp2Midwest(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMidwestProb3 = newProb3Midwest - oldProb3Midwest;
ceMidwestProb3   = mean(mean(diffMidwestProb3));
sprintf('%.04f',ceMidwestProb3)

% Outcome 4
%------------
oldProb4Midwest    = zeros(n,mc);        
newProb4Midwest    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4Midwest(i,j)    = 1 - oldComp3Midwest(i,j);
    
    newProb4Midwest(i,j)    = 1 - newComp3Midwest(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffMidwestProb4 = newProb4Midwest - oldProb4Midwest;
ceMidwestProb4   = mean(mean(diffMidwestProb4));
sprintf('%.04f',ceMidwestProb4)

disp('CE, Midwest')
disp([ceMidwestProb1 ceMidwestProb2 ceMidwestProb3 ceMidwestProb4])

%--------------------------------------------------------------------------
clearvars oldProb1Midwest oldComp2Midwest oldProb2Midwest oldComp3Midwest oldProb3Midwest oldProb4Midwest
clearvars newProb1Midwest newComp2Midwest newProb2Midwest newComp3Midwest newProb3Midwest newProb4Midwest
clearvars diffMidwestProb1 diffMidwestProb2 diffMidwestProb3 diffMidwestProb4

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

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

% Outcome 1
%------------
oldProb1South = zeros(n,mc);
newProb1South = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, South, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1South(i,j)    = fast_normcdf(gammaCp(1,j) - SouthX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1South(i,j)    = fast_normcdf(gammaCp(1,j) - SouthX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSouthProb1 = newProb1South - oldProb1South;
ceSouthProb1   = mean(mean(diffSouthProb1));
sprintf('%.04f',ceSouthProb1)

% Outcome 2
%------------
oldProb2South    = zeros(n,mc);        
newProb2South    = zeros(n,mc);

oldComp2South    = zeros(n,mc);
newComp2South    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, South, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2South(i,j)    = fast_normcdf(gammaCp(2,j)- SouthX0(i,:)*betaStore(:,burn+j));
    oldProb2South(i,j)    = oldComp2South(i,j) - oldProb1South(i,j);
    
    newComp2South(i,j)    = fast_normcdf(gammaCp(2,j)- SouthX1(i,:)*betaStore(:,burn+j));
    newProb2South(i,j)    = newComp2South(i,j) - newProb1South(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSouthProb2 = newProb2South - oldProb2South;
ceSouthProb2   = mean(mean(diffSouthProb2));
sprintf('%.04f',ceSouthProb2)

% Outcome 3
%------------
oldProb3South    = zeros(n,mc);        
newProb3South    = zeros(n,mc);

oldComp3South    = zeros(n,mc);
newComp3South    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, South, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3South(i,j)    = fast_normcdf(gammaCp(3,j)- SouthX0(i,:)*betaStore(:,burn+j));
    oldProb3South(i,j)    = oldComp3South(i,j) - oldComp2South(i,j);

    newComp3South(i,j)    = fast_normcdf(gammaCp(3,j)- SouthX1(i,:)*betaStore(:,burn+j));
    newProb3South(i,j)    = newComp3South(i,j) - newComp2South(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSouthProb3 = newProb3South - oldProb3South;
ceSouthProb3   = mean(mean(diffSouthProb3));
sprintf('%.04f',ceSouthProb3)

% Outcome 4
%------------
oldProb4South    = zeros(n,mc);        
newProb4South    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4South(i,j)    = 1 - oldComp3South(i,j);
    
    newProb4South(i,j)    = 1 - newComp3South(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffSouthProb4 = newProb4South - oldProb4South;
ceSouthProb4   = mean(mean(diffSouthProb4));
sprintf('%.04f',ceSouthProb4)

disp('CE, South')
disp([ceSouthProb1 ceSouthProb2 ceSouthProb3 ceSouthProb4])

%--------------------------------------------------------------------------
clearvars oldProb1South oldComp2South oldProb2South oldComp3South oldProb3South oldProb4South
clearvars newProb1South newComp2South newProb2South newComp3South newProb3South newProb4South
clearvars diffSouthProb1 diffSouthProb2 diffSouthProb3 diffSouthProb4

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

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

% Outcome 1
%------------
oldProb1West = zeros(n,mc);
newProb1West = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, West, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1West(i,j)    = fast_normcdf(gammaCp(1,j) - WestX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1West(i,j)    = fast_normcdf(gammaCp(1,j) - WestX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffWestProb1 = newProb1West - oldProb1West;
ceWestProb1   = mean(mean(diffWestProb1));
sprintf('%.04f',ceWestProb1)

% Outcome 2
%------------
oldProb2West    = zeros(n,mc);        
newProb2West    = zeros(n,mc);

oldComp2West    = zeros(n,mc);
newComp2West    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, West, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2West(i,j)    = fast_normcdf(gammaCp(2,j)- WestX0(i,:)*betaStore(:,burn+j));
    oldProb2West(i,j)    = oldComp2West(i,j) - oldProb1West(i,j);
    
    newComp2West(i,j)    = fast_normcdf(gammaCp(2,j)- WestX1(i,:)*betaStore(:,burn+j));
    newProb2West(i,j)    = newComp2West(i,j) - newProb1West(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffWestProb2 = newProb2West - oldProb2West;
ceWestProb2   = mean(mean(diffWestProb2));
sprintf('%.04f',ceWestProb2)

% Outcome 3
%------------
oldProb3West    = zeros(n,mc);        
newProb3West    = zeros(n,mc);

oldComp3West    = zeros(n,mc);
newComp3West    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, West, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3West(i,j)    = fast_normcdf(gammaCp(3,j)- WestX0(i,:)*betaStore(:,burn+j));
    oldProb3West(i,j)    = oldComp3West(i,j) - oldComp2West(i,j);

    newComp3West(i,j)    = fast_normcdf(gammaCp(3,j)- WestX1(i,:)*betaStore(:,burn+j));
    newProb3West(i,j)    = newComp3West(i,j) - newComp2West(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffWestProb3 = newProb3West - oldProb3West;
ceWestProb3   = mean(mean(diffWestProb3));
sprintf('%.04f',ceWestProb3)

% Outcome 4
%------------
oldProb4West    = zeros(n,mc);        
newProb4West    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4West(i,j)    = 1 - oldComp3West(i,j);
    
    newProb4West(i,j)    = 1 - newComp3West(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffWestProb4 = newProb4West - oldProb4West;
ceWestProb4   = mean(mean(diffWestProb4));
sprintf('%.04f',ceWestProb4)

disp('CE, West')
disp([ceWestProb1 ceWestProb2 ceWestProb3 ceWestProb4])

%--------------------------------------------------------------------------
clearvars oldProb1West oldComp2West oldProb2West oldComp3West oldProb3West oldProb4West
clearvars newProb1West newComp2West newProb2West newComp3West newProb3West newProb4West
clearvars diffWestProb1 diffWestProb2 diffWestProb3 diffWestProb4

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

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

% Outcome 1
%------------
oldProb1USborn = zeros(n,mc);
newProb1USborn = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, USborn, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1USborn(i,j)    = fast_normcdf(gammaCp(1,j) - USbornX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1USborn(i,j)    = fast_normcdf(gammaCp(1,j) - USbornX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffUSbornProb1 = newProb1USborn - oldProb1USborn;
ceUSbornProb1   = mean(mean(diffUSbornProb1));
sprintf('%.04f',ceUSbornProb1)

% Outcome 2
%------------
oldProb2USborn    = zeros(n,mc);        
newProb2USborn    = zeros(n,mc);

oldComp2USborn    = zeros(n,mc);
newComp2USborn    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, USborn, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2USborn(i,j)    = fast_normcdf(gammaCp(2,j)- USbornX0(i,:)*betaStore(:,burn+j));
    oldProb2USborn(i,j)    = oldComp2USborn(i,j) - oldProb1USborn(i,j);
    
    newComp2USborn(i,j)    = fast_normcdf(gammaCp(2,j)- USbornX1(i,:)*betaStore(:,burn+j));
    newProb2USborn(i,j)    = newComp2USborn(i,j) - newProb1USborn(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffUSbornProb2 = newProb2USborn - oldProb2USborn;
ceUSbornProb2   = mean(mean(diffUSbornProb2));
sprintf('%.04f',ceUSbornProb2)

% Outcome 3
%------------
oldProb3USborn    = zeros(n,mc);        
newProb3USborn    = zeros(n,mc);

oldComp3USborn    = zeros(n,mc);
newComp3USborn    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, USborn, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3USborn(i,j)    = fast_normcdf(gammaCp(3,j)- USbornX0(i,:)*betaStore(:,burn+j));
    oldProb3USborn(i,j)    = oldComp3USborn(i,j) - oldComp2USborn(i,j);

    newComp3USborn(i,j)    = fast_normcdf(gammaCp(3,j)- USbornX1(i,:)*betaStore(:,burn+j));
    newProb3USborn(i,j)    = newComp3USborn(i,j) - newComp2USborn(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffUSbornProb3 = newProb3USborn - oldProb3USborn;
ceUSbornProb3   = mean(mean(diffUSbornProb3));
sprintf('%.04f',ceUSbornProb3)

% Outcome 4
%------------
oldProb4USborn    = zeros(n,mc);        
newProb4USborn    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4USborn(i,j)    = 1 - oldComp3USborn(i,j);
    
    newProb4USborn(i,j)    = 1 - newComp3USborn(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffUSbornProb4 = newProb4USborn - oldProb4USborn;
ceUSbornProb4   = mean(mean(diffUSbornProb4));
sprintf('%.04f',ceUSbornProb4)

disp('CE, USborn')
disp([ceUSbornProb1 ceUSbornProb2 ceUSbornProb3 ceUSbornProb4])

%--------------------------------------------------------------------------
clearvars oldProb1USborn oldComp2USborn oldProb2USborn oldComp3USborn oldProb3USborn oldProb4USborn
clearvars newProb1USborn newComp2USborn newProb2USborn newComp3USborn newProb3USborn newProb4USborn
clearvars diffUSbornProb1 diffUSbornProb2 diffUSbornProb3 diffUSbornProb4

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

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

% Outcome 1
%------------
oldProb1Black = zeros(n,mc);
newProb1Black = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Black, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Black(i,j)    = fast_normcdf(gammaCp(1,j) - BlackX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Black(i,j)    = fast_normcdf(gammaCp(1,j) - BlackX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffBlackProb1 = newProb1Black - oldProb1Black;
ceBlackProb1   = mean(mean(diffBlackProb1));
sprintf('%.04f',ceBlackProb1)

% Outcome 2
%------------
oldProb2Black    = zeros(n,mc);        
newProb2Black    = zeros(n,mc);

oldComp2Black    = zeros(n,mc);
newComp2Black    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Black, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Black(i,j)    = fast_normcdf(gammaCp(2,j)- BlackX0(i,:)*betaStore(:,burn+j));
    oldProb2Black(i,j)    = oldComp2Black(i,j) - oldProb1Black(i,j);
    
    newComp2Black(i,j)    = fast_normcdf(gammaCp(2,j)- BlackX1(i,:)*betaStore(:,burn+j));
    newProb2Black(i,j)    = newComp2Black(i,j) - newProb1Black(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffBlackProb2 = newProb2Black - oldProb2Black;
ceBlackProb2   = mean(mean(diffBlackProb2));
sprintf('%.04f',ceBlackProb2)

% Outcome 3
%------------
oldProb3Black    = zeros(n,mc);        
newProb3Black    = zeros(n,mc);

oldComp3Black    = zeros(n,mc);
newComp3Black    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Black, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Black(i,j)    = fast_normcdf(gammaCp(3,j)- BlackX0(i,:)*betaStore(:,burn+j));
    oldProb3Black(i,j)    = oldComp3Black(i,j) - oldComp2Black(i,j);

    newComp3Black(i,j)    = fast_normcdf(gammaCp(3,j)- BlackX1(i,:)*betaStore(:,burn+j));
    newProb3Black(i,j)    = newComp3Black(i,j) - newComp2Black(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffBlackProb3 = newProb3Black - oldProb3Black;
ceBlackProb3   = mean(mean(diffBlackProb3));
sprintf('%.04f',ceBlackProb3)

% Outcome 4
%------------
oldProb4Black    = zeros(n,mc);        
newProb4Black    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4Black(i,j)    = 1 - oldComp3Black(i,j);
    
    newProb4Black(i,j)    = 1 - newComp3Black(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffBlackProb4 = newProb4Black - oldProb4Black;
ceBlackProb4   = mean(mean(diffBlackProb4));
sprintf('%.04f',ceBlackProb4)

disp('CE, Black')
disp([ceBlackProb1 ceBlackProb2 ceBlackProb3 ceBlackProb4])

%--------------------------------------------------------------------------
clearvars oldProb1Black oldComp2Black oldProb2Black oldComp3Black oldProb3Black oldProb4Black
clearvars newProb1Black newComp2Black newProb2Black newComp3Black newProb3Black newProb4Black
clearvars diffBlackProb1 diffBlackProb2 diffBlackProb3 diffBlackProb4


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

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

% Outcome 1
%------------
oldProb1OtherRaces = zeros(n,mc);
newProb1OtherRaces = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, OtherRaces, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1OtherRaces(i,j)    = fast_normcdf(gammaCp(1,j) - OtherRacesX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1OtherRaces(i,j)    = fast_normcdf(gammaCp(1,j) - OtherRacesX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOtherRacesProb1 = newProb1OtherRaces - oldProb1OtherRaces;
ceOtherRacesProb1   = mean(mean(diffOtherRacesProb1));
sprintf('%.04f',ceOtherRacesProb1)

% Outcome 2
%------------
oldProb2OtherRaces    = zeros(n,mc);        
newProb2OtherRaces    = zeros(n,mc);

oldComp2OtherRaces    = zeros(n,mc);
newComp2OtherRaces    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, OtherRaces, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2OtherRaces(i,j)    = fast_normcdf(gammaCp(2,j)- OtherRacesX0(i,:)*betaStore(:,burn+j));
    oldProb2OtherRaces(i,j)    = oldComp2OtherRaces(i,j) - oldProb1OtherRaces(i,j);
    
    newComp2OtherRaces(i,j)    = fast_normcdf(gammaCp(2,j)- OtherRacesX1(i,:)*betaStore(:,burn+j));
    newProb2OtherRaces(i,j)    = newComp2OtherRaces(i,j) - newProb1OtherRaces(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOtherRacesProb2 = newProb2OtherRaces - oldProb2OtherRaces;
ceOtherRacesProb2   = mean(mean(diffOtherRacesProb2));
sprintf('%.04f',ceOtherRacesProb2)

% Outcome 3
%------------
oldProb3OtherRaces    = zeros(n,mc);        
newProb3OtherRaces    = zeros(n,mc);

oldComp3OtherRaces    = zeros(n,mc);
newComp3OtherRaces    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, OtherRaces, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3OtherRaces(i,j)    = fast_normcdf(gammaCp(3,j)- OtherRacesX0(i,:)*betaStore(:,burn+j));
    oldProb3OtherRaces(i,j)    = oldComp3OtherRaces(i,j) - oldComp2OtherRaces(i,j);

    newComp3OtherRaces(i,j)    = fast_normcdf(gammaCp(3,j)- OtherRacesX1(i,:)*betaStore(:,burn+j));
    newProb3OtherRaces(i,j)    = newComp3OtherRaces(i,j) - newComp2OtherRaces(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOtherRacesProb3 = newProb3OtherRaces - oldProb3OtherRaces;
ceOtherRacesProb3   = mean(mean(diffOtherRacesProb3));
sprintf('%.04f',ceOtherRacesProb3)

% Outcome 4
%------------
oldProb4OtherRaces    = zeros(n,mc);        
newProb4OtherRaces    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4OtherRaces(i,j)    = 1 - oldComp3OtherRaces(i,j);
    
    newProb4OtherRaces(i,j)    = 1 - newComp3OtherRaces(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffOtherRacesProb4 = newProb4OtherRaces - oldProb4OtherRaces;
ceOtherRacesProb4   = mean(mean(diffOtherRacesProb4));
sprintf('%.04f',ceOtherRacesProb4)

disp('CE, OtherRaces')
disp([ceOtherRacesProb1 ceOtherRacesProb2 ceOtherRacesProb3 ceOtherRacesProb4])

%--------------------------------------------------------------------------
clearvars oldProb1OtherRaces oldComp2OtherRaces oldProb2OtherRaces oldComp3OtherRaces oldProb3OtherRaces oldProb4OtherRaces
clearvars newProb1OtherRaces newComp2OtherRaces newProb2OtherRaces newComp3OtherRaces newProb3OtherRaces newProb4OtherRaces
clearvars diffOtherRacesProb1 diffOtherRacesProb2 diffOtherRacesProb3 diffOtherRacesProb4


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

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

% Outcome 1
%------------
oldProb1Democrat = zeros(n,mc);
newProb1Democrat = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Democrat, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1Democrat(i,j)    = fast_normcdf(gammaCp(1,j) - DemocratX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1Democrat(i,j)    = fast_normcdf(gammaCp(1,j) - DemocratX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffDemocratProb1 = newProb1Democrat - oldProb1Democrat;
ceDemocratProb1   = mean(mean(diffDemocratProb1));
sprintf('%.04f',ceDemocratProb1)

% Outcome 2
%------------
oldProb2Democrat    = zeros(n,mc);        
newProb2Democrat    = zeros(n,mc);

oldComp2Democrat    = zeros(n,mc);
newComp2Democrat    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Democrat, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2Democrat(i,j)    = fast_normcdf(gammaCp(2,j)- DemocratX0(i,:)*betaStore(:,burn+j));
    oldProb2Democrat(i,j)    = oldComp2Democrat(i,j) - oldProb1Democrat(i,j);
    
    newComp2Democrat(i,j)    = fast_normcdf(gammaCp(2,j)- DemocratX1(i,:)*betaStore(:,burn+j));
    newProb2Democrat(i,j)    = newComp2Democrat(i,j) - newProb1Democrat(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffDemocratProb2 = newProb2Democrat - oldProb2Democrat;
ceDemocratProb2   = mean(mean(diffDemocratProb2));
sprintf('%.04f',ceDemocratProb2)

% Outcome 3
%------------
oldProb3Democrat    = zeros(n,mc);        
newProb3Democrat    = zeros(n,mc);

oldComp3Democrat    = zeros(n,mc);
newComp3Democrat    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, Democrat, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3Democrat(i,j)    = fast_normcdf(gammaCp(3,j)- DemocratX0(i,:)*betaStore(:,burn+j));
    oldProb3Democrat(i,j)    = oldComp3Democrat(i,j) - oldComp2Democrat(i,j);

    newComp3Democrat(i,j)    = fast_normcdf(gammaCp(3,j)- DemocratX1(i,:)*betaStore(:,burn+j));
    newProb3Democrat(i,j)    = newComp3Democrat(i,j) - newComp2Democrat(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffDemocratProb3 = newProb3Democrat - oldProb3Democrat;
ceDemocratProb3   = mean(mean(diffDemocratProb3));
sprintf('%.04f',ceDemocratProb3)

% Outcome 4
%------------
oldProb4Democrat    = zeros(n,mc);        
newProb4Democrat    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4Democrat(i,j)    = 1 - oldComp3Democrat(i,j);
    
    newProb4Democrat(i,j)    = 1 - newComp3Democrat(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffDemocratProb4 = newProb4Democrat - oldProb4Democrat;
ceDemocratProb4   = mean(mean(diffDemocratProb4));
sprintf('%.04f',ceDemocratProb4)

disp('CE, Democrat')
disp([ceDemocratProb1 ceDemocratProb2 ceDemocratProb3 ceDemocratProb4])

%--------------------------------------------------------------------------
clearvars oldProb1Democrat oldComp2Democrat oldProb2Democrat oldComp3Democrat oldProb3Democrat oldProb4Democrat
clearvars newProb1Democrat newComp2Democrat newProb2Democrat newComp3Democrat newProb3Democrat newProb4Democrat
clearvars diffDemocratProb1 diffDemocratProb2 diffDemocratProb3 diffDemocratProb4

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

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

% Outcome 1
%------------
oldProb1IndnElse = zeros(n,mc);
newProb1IndnElse = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, IndnElse, Outcome 1');
for j =1:mc
    for i=1:n
    oldProb1IndnElse(i,j)    = fast_normcdf(gammaCp(1,j) - IndnElseX0(i,:)*betaStore(:,burn+j));       % NOTE: mu=0, sigma=1           
    newProb1IndnElse(i,j)    = fast_normcdf(gammaCp(1,j) - IndnElseX1(i,:)*betaStore(:,burn+j));
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffIndnElseProb1 = newProb1IndnElse - oldProb1IndnElse;
ceIndnElseProb1   = mean(mean(diffIndnElseProb1));
sprintf('%.04f',ceIndnElseProb1)

% Outcome 2
%------------
oldProb2IndnElse    = zeros(n,mc);        
newProb2IndnElse    = zeros(n,mc);

oldComp2IndnElse    = zeros(n,mc);
newComp2IndnElse    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, IndnElse, Outcome 2');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp2IndnElse(i,j)    = fast_normcdf(gammaCp(2,j)- IndnElseX0(i,:)*betaStore(:,burn+j));
    oldProb2IndnElse(i,j)    = oldComp2IndnElse(i,j) - oldProb1IndnElse(i,j);
    
    newComp2IndnElse(i,j)    = fast_normcdf(gammaCp(2,j)- IndnElseX1(i,:)*betaStore(:,burn+j));
    newProb2IndnElse(i,j)    = newComp2IndnElse(i,j) - newProb1IndnElse(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffIndnElseProb2 = newProb2IndnElse - oldProb2IndnElse;
ceIndnElseProb2   = mean(mean(diffIndnElseProb2));
sprintf('%.04f',ceIndnElseProb2)

% Outcome 3
%------------
oldProb3IndnElse    = zeros(n,mc);        
newProb3IndnElse    = zeros(n,mc);

oldComp3IndnElse    = zeros(n,mc);
newComp3IndnElse    = zeros(n,mc);

tic
% wbar = waitbar(0,'Covariate Effect, IndnElse, Outcome 3');
for j =1:mc                  % mc is number of iterations after burn-in
    for i=1:n
    oldComp3IndnElse(i,j)    = fast_normcdf(gammaCp(3,j)- IndnElseX0(i,:)*betaStore(:,burn+j));
    oldProb3IndnElse(i,j)    = oldComp3IndnElse(i,j) - oldComp2IndnElse(i,j);

    newComp3IndnElse(i,j)    = fast_normcdf(gammaCp(3,j)- IndnElseX1(i,:)*betaStore(:,burn+j));
    newProb3IndnElse(i,j)    = newComp3IndnElse(i,j) - newComp2IndnElse(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffIndnElseProb3 = newProb3IndnElse - oldProb3IndnElse;
ceIndnElseProb3   = mean(mean(diffIndnElseProb3));
sprintf('%.04f',ceIndnElseProb3)

% Outcome 4
%------------
oldProb4IndnElse    = zeros(n,mc);        
newProb4IndnElse    = zeros(n,mc);

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

for j =1:mc
    for i=1:n
    oldProb4IndnElse(i,j)    = 1 - oldComp3IndnElse(i,j);
    
    newProb4IndnElse(i,j)    = 1 - newComp3IndnElse(i,j);
    end
% waitbar(j/mc);    
end
% close(wbar);
toc

diffIndnElseProb4 = newProb4IndnElse - oldProb4IndnElse;
ceIndnElseProb4   = mean(mean(diffIndnElseProb4));
sprintf('%.04f',ceIndnElseProb4)

disp('CE, IndnElse')
disp([ceIndnElseProb1 ceIndnElseProb2 ceIndnElseProb3 ceIndnElseProb4])

%--------------------------------------------------------------------------
clearvars oldProb1IndnElse oldComp2IndnElse oldProb2IndnElse oldComp3IndnElse oldProb3IndnElse oldProb4IndnElse
clearvars newProb1IndnElse newComp2IndnElse newProb2IndnElse newComp3IndnElse newProb3IndnElse newProb4IndnElse
clearvars diffIndnElseProb1 diffIndnElseProb2 diffIndnElseProb3 diffIndnElseProb4

%% 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+50000)/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])

%--------------------------------------------------------------------------
clearvars oldProb1Income oldComp2Income oldProb2Income oldComp3Income oldProb3Income oldProb4Income
clearvars newProb1Income newComp2Income newProb2Income newComp3Income newProb3Income newProb4Income
clearvars diffIncomeProb1 diffIncomeProb2 diffIncomeProb3 diffIncomeProb4


%% 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('___________________________________________________ \n')
    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')

