% This script takes the results of a Bayesian meta-analysis and estimates
% household and total WTP for lost wetland acres at the state level. 
%
% Two sets of estimates are produced.  Main results apply household WTP to
% in-state residents only.  A sensitivity analysis applies household WTP to
% counties adjacent to the affected stated and includes those in the
% totals.


clear;

rand('state',37); % set arbitrary seed for uniform draws
randn('state',37); % set arbitrary seed for normal draws

load M3thin betamatthin gammatthin sig2matthin;
load stage2_wetland_metatransfer_data;

betamat=betamatthin;
gammat=gammatthin;
sig2mat=sig2matthin;

% Preliminary stuff
%%%%%%%%%%%%%%%%%%%%%%%%
k1 = 11;
k2 = 2;
r2=size(betamat,2);
rp=10; %draws of stochastic terms per draw of parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% settings for BT scenarios
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1   const           1
% 2   lnyear          log(2018-1988+1) = 3.4011
% 3   lninc           state-specific median income (lninc_adj includes adjacent counties weighted by households)
% 4   sagulf          state-specific
% 5   nema            state-specific
% 6   nmw             state-specific
% 7   local           0 (non-local)
% 8   prov            varies by scenario 
% 9   reg             varies by scenario
% 10  cult            varies by scenario
% 11  forest          state-specific-proportion of wetlands that are forested

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set ecosystem services values: proportions equal to study counts in
% meta-data
local   = 0;
prov    = 0;
reg     = 0;
cult    = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare methodological combos - only need to do this once
Jmat=[0 0; 1 0];    % all permissible combos 
kj=size(Jmat,1);
wj=1/kj;            % neutral weight for each combo

Meb=zeros(1,r2);
for j=1:r2
    int=repmat(betamat(k1+1:k1+k2,j)',kj,1); %kj by k2
    int=int.*Jmat;  % will zero out betas not needed for given methodology combos
    int=wj*sum(sum(int),2); % can do this since wj unchanged over combos; 
    Meb(j)=int;
end
Meb=reshape(repmat(Meb,rp,1),1,rp*r2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare stochastic terms - only needed to do this once
% only used for mean estimates.  zero stochastic terms out for medians
int=normrnd(0,repmat(sqrt(sig2mat),rp,1)); %rp by r2
epmat=reshape(int,1,r2*rp); %r2*rp as needed


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate hosuehold WTP for in-state households
for i = 1:length(state_name)
    xp = [1 3.4011 lninc(i,:) sagulf(i,:) nema(i,:) nmw(i,:) 0 prov reg cult forest(i,:)]';
    f_s=log((1./gammat).*(exp(gammat*q1(i,:))-exp(gammat*q0(i,:))));
    
    Xpmat   = repmat(xp,1,r2);          % k1 by r2
    bmat    = betamat(1:k1,:);          % k1 by r2
    %Xb      = sum(Xpmat.*bmat) + f_s;   % 1  by r2 - add nonlinear part
    Xb      = sum(Xpmat.*bmat) + f_s + betamat(13,:);  %1 by r2 - add nonlinear part AND DRAWS OF LUMP SUM COEFFICIENT
    Xb      = reshape(repmat(Xb,rp,1),1,rp*r2);
    WTP(:,i) = exp(Xb+Meb+epmat); %%%%%%%%%% For Mean %%%%%%%%%% 
    %WTP(:,i) = exp(Xb+Meb);     %%%%%%%%%% For Median %%%%%%%%%%
    clc
    fprintf('In-State WTP: Result %2.0f out of 51\n',i)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate hosuehold WTP for in-state households AND ADJACENT COUNTIES
for i = 1:length(state_name)
    xp = [1 3.4011 lninc_adj(i,:) sagulf(i,:) nema(i,:) nmw(i,:) 0 prov reg cult forest(i,:)]';
    f_s=log((1./gammat).*(exp(gammat*q1(i,:))-exp(gammat*q0(i,:))));
    
    Xpmat   = repmat(xp,1,r2);          % k1 by r2
    bmat    = betamat(1:k1,:);          % k1 by r2
    %Xb      = sum(Xpmat.*bmat) + f_s;   % 1  by r2 - add nonlinear part
    Xb      = sum(Xpmat.*bmat) + f_s + betamat(13,:);  %1 by r2 - add nonlinear part AND DRAWS OF LUMP SUM COEFFICIENT
    Xb      = reshape(repmat(Xb,rp,1),1,rp*r2);
    WTP_adj(:,i) = exp(Xb+Meb+epmat); %%%%%%%%%% For Mean %%%%%%%%%% 
    %WTP_adj(:,i) = exp(Xb+Meb);     %%%%%%%%%% For Median %%%%%%%%%%
    clc
    fprintf('Including Adjacent Counties: Result %2.0f out of 51\n',i)
end

% Replace states with no adjacent counties (AK & HI) with in-state WTP results 
[row,col] = find(adj_hh==0);
WTP_adj(:,row) = WTP(:,row);

tp = 99; % Truncation point 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prepare results for reporting 

for i = 1:length(state_name)
    if sum(WTP(:,i)>0)==0       % required for truncation because states with zero affected acres will have no truncation point
        mean_wtp(:,i) = 0;       mean_wtp_adj(:,i) = 0;
        std_wtp(:,i) = 0;        std_wtp_adj(:,i) = 0;
        lower_5(:,i) = 0;        lower_5_adj(:,i) = 0;
        upper_95(:,i) = 0;       upper_95_adj(:,i) = 0;
        wtp_per_acre(:,i) = 0;   wtp_per_acre_adj(:,i) = 0;
        state_wtp_dist(:,i) = 0; state_wtp_dist_adj(:,i) =0;
    else
        mxm = prctile(WTP(:,i),tp);                         mxm_adj = prctile(WTP_adj(:,i),tp);
        index = WTP(:,i)<=mxm;                              index_adj = WTP_adj(:,i)<=mxm_adj;
        mean_wtp(:,i) = mean(WTP(index,i));                 mean_wtp_adj(:,i) = mean(WTP_adj(index_adj,i));
        std_wtp(:,i) = std(WTP(index,i));                   std_wtp_adj(:,i) = std(WTP_adj(index_adj,i));
        lower_5(:,i) = prctile(WTP(index,i),5);             lower_5_adj(:,i) = prctile(WTP_adj(index_adj,i),5);
        upper_95(:,i) = prctile(WTP(index,i),95);           upper_95_adj(:,i) = prctile(WTP_adj(index_adj,i),95);
        wtp_per_acre(:,i) = WTP(index,i)/(delta(i)*1000);   wtp_per_acre_adj(:,i) = WTP_adj(index_adj,i)/(delta(i)*1000);
        state_wtp_dist(:,i) = WTP(index,i)*n_HH(i);         state_wtp_dist_adj(:,i) = WTP_adj(index_adj,i)*adj_hh(i);
    end
    clc
    fprintf('Preparing Results: %2.0f out of 51\n',i)
end 
clc
mean_wtp_per_acre = mean(wtp_per_acre); mean_wtp_per_acre(isnan(mean_wtp_per_acre)) = 0;
mean_state_wtp = mean(state_wtp_dist);

mean_wtp_per_acre_adj = mean(wtp_per_acre_adj); mean_wtp_per_acre_adj(isnan(mean_wtp_per_acre_adj)) = 0;
mean_state_wtp_adj = mean(state_wtp_dist_adj);

lower5_wtp_per_acre  = prctile(wtp_per_acre,5,1);    lower5_wtp_per_acre_adj  = prctile(wtp_per_acre_adj,5,1);
upper95_wtp_per_acre = prctile(wtp_per_acre,95,1);   upper95_wtp_per_acre_adj = prctile(wtp_per_acre_adj,95,1);
lower5_state_wtp     = prctile(state_wtp_dist,5,1);  lower5_state_wtp_adj     = prctile(state_wtp_dist_adj,5,1);
upper95_state_wtp    = prctile(state_wtp_dist,95,1); upper95_state_wtp_adj    = prctile(state_wtp_dist_adj,95,1);

save('wtp_dist.mat','WTP','WTP_adj','state_name','delta','n_HH')

% %%%%%%%%%%%%% Display in-state results %%%%%%%%%%%%%%%%
% disp('In-State Wetlands Benefits')
%     disp(   'State  n_HH    Acres   HH WTP/Acre  Total WTP   HH WTP/Acre 5th   Total WTP 5th   HH WTP/Acre 95th   Total WTP 95th')
% for i = 1:length(state_name)
%     fprintf('%5s   %6.0f   %4.1f       %4.4f       %4.4f         %4.4f           %4.4f             %4.4f              %4.4f\n',...
%     state_name(i,:), n_HH(i), delta(i)*1000, mean_wtp_per_acre(i), mean_state_wtp(i), lower5_wtp_per_acre(i), ...
%     lower5_state_wtp(i), upper95_wtp_per_acre(i), upper95_state_wtp(i))
% end 
% results = [n_HH, delta*1000, mean_wtp_per_acre', mean_state_wtp', lower5_wtp_per_acre', lower5_state_wtp', upper95_wtp_per_acre', upper95_state_wtp'];
% xlswrite('final WTP results_stage2.xlsx',results,'state-level in-state wtp only','B2:I53')
% 
% %%%%%%%%%%%%% Display state and adjacent county results %%%%%%%%%%%%%%%%
% disp(' ')
% disp('State and Adjacent Counties Wetlands Benefits')
%     disp(   'State  n_HH    Acres   HH WTP/Acre  Total WTP   HH WTP/Acre 5th   Total WTP 5th   HH WTP/Acre 95th   Total WTP 95th')
% for i = 1:length(state_name)
%     fprintf('%5s   %6.0f   %4.1f       %4.4f       %4.4f         %4.4f           %4.4f             %4.4f              %4.4f\n',...
%     state_name(i,:), adj_hh(i), delta(i)*1000, mean_wtp_per_acre_adj(i), mean_state_wtp_adj(i), lower5_wtp_per_acre_adj(i), ...
%     lower5_state_wtp_adj(i), upper95_wtp_per_acre_adj(i), upper95_state_wtp_adj(i))
% end 
% results_adj = [adj_hh, delta*1000, mean_wtp_per_acre_adj', mean_state_wtp_adj', lower5_wtp_per_acre_adj', lower5_state_wtp_adj', upper95_wtp_per_acre_adj', upper95_state_wtp_adj'];
%xlswrite('final WTP results_stage2.xlsx',results_adj,'state-level with adj counties','B2:I53')

% Summary statistics for totals
total_wtp_dist = sum(state_wtp_dist,2);
total_wtp_dist_adj = sum(state_wtp_dist_adj,2);
summary_stats = [mean(total_wtp_dist )*10^-6, prctile(total_wtp_dist ,5)*10^-6, prctile(total_wtp_dist ,95)*10^-6;
                 mean(total_wtp_dist_adj )*10^-6, prctile(total_wtp_dist_adj ,5)*10^-6, prctile(total_wtp_dist_adj ,95)*10^-6];
             
xlswrite('final WTP results_stage2.xlsx',summary_stats,'summary stats','B3:D4')
sim_settings = [prov;reg;cult;tp];
xlswrite('final WTP results_stage2.xlsx',sim_settings,'summary stats','B7:B10')

fprintf('\n')
disp('                         Mean total benefits    Lower 5th    Upper 95th')
fprintf('In-State Only:               %10.3f          %10.3f   %10.3f\n',mean(total_wtp_dist )*10^-6, prctile(total_wtp_dist ,5)*10^-6, prctile(total_wtp_dist ,95)*10^-6)
fprintf('State and Adjacent Counties: %10.3f          %10.3f   %10.3f\n',mean(total_wtp_dist_adj )*10^-6, prctile(total_wtp_dist_adj ,5)*10^-6, prctile(total_wtp_dist_adj ,95)*10^-6)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Federalism scenarios
% Summary statistics for in-state totals

% Federalism scenarios
wtp_scenario1 = state_wtp_dist;
for i = 1:length(state_name)
    if fed_cat(i,:) > 2
        wtp_scenario1(:,i)=0;
    end
end 
total_wtp_scenario1 = sum(wtp_scenario1,2);

wtp_scenario3 = state_wtp_dist;
for i = 1:length(state_name)
    if fed_cat(i,:) ~= 1
        wtp_scenario3(:,i)=0;
    end
end
total_wtp_scenario3 = sum(wtp_scenario3,2);

fprintf('\n')
disp('Federalism Scenarios - in-state Only')
disp('Scenario            Mean Total WTP   Standard Dev   Lower 5th   Upper 95th')
fprintf('Full Deregulation     %10.2f   %10.2f   %10.2f   %10.2f\n',mean(total_wtp_dist ), std(total_wtp_dist), prctile(total_wtp_dist ,5), prctile(total_wtp_dist ,95))
fprintf('Secenario 1 & 2       %10.2f   %10.2f   %10.2f   %10.2f\n',mean(total_wtp_scenario1 ), std(total_wtp_scenario1), prctile(total_wtp_scenario1,5), prctile(total_wtp_scenario1,95))
fprintf('Secenario 3           %10.2f   %10.2f   %10.2f   %10.2f\n',mean(total_wtp_scenario3 ), std(total_wtp_scenario3), prctile(total_wtp_scenario3,5), prctile(total_wtp_scenario3,95))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Summary statistics for state+adjacent county totals

% Federalism scenarios
wtp_scenario1_adj = state_wtp_dist_adj;
for i = 1:length(state_name)
    if fed_cat(i,:) > 2
        wtp_scenario1_adj(:,i)=0;
    end
end 
total_wtp_scenario1_adj = sum(wtp_scenario1_adj,2);

wtp_scenario3_adj = state_wtp_dist_adj;
for i = 1:length(state_name)
    if fed_cat(i,:) ~= 1
        wtp_scenario3_adj(:,i)=0;
    end
end
total_wtp_scenario3_adj = sum(wtp_scenario3_adj,2);

fprintf('\n')
disp('Federalism Scenarios - Including Adjacent Counties')
disp('Scenario            Mean Total WTP   Standard Dev   Lower 5th   Upper 95th')
fprintf('Full Deregulation     %10.2f   %10.2f   %10.2f   %10.2f\n',mean(total_wtp_dist_adj ), std(total_wtp_dist_adj), prctile(total_wtp_dist_adj ,5), prctile(total_wtp_dist_adj ,95))
fprintf('Secenario 1 & 2       %10.2f   %10.2f   %10.2f   %10.2f\n',mean(total_wtp_scenario1_adj ), std(total_wtp_scenario1_adj), prctile(total_wtp_scenario1_adj,5), prctile(total_wtp_scenario1_adj,95))
fprintf('Secenario 3           %10.2f   %10.2f   %10.2f   %10.2f\n',mean(total_wtp_scenario3_adj ), std(total_wtp_scenario3_adj), prctile(total_wtp_scenario3_adj,5), prctile(total_wtp_scenario3_adj,95))

