Convolutions Method POE

[Federica Di Marcantonio]

Dears

I am trying to create a code in Stata that allows me to perform the estimation of the difference of two distributions starting from the mean and se of estimated willingness to pay (WTP).
In theory I have to make one random draw from the m ∗ n, Xˆ i − Yˆ j combinations. In this paper the authors suggest that the significance of the difference be calculated a large number (L) of times and averaged across the L sample repetitions, using standard bootstrapping techniques of repeated sampling with replacement from the two distributions. Letting BL (Xˆ ) indicate the Lth random sampling with replacement from distribution Xˆ , and BL (Yˆ ) be the corresponding sample from Yˆ , then estimating the significance of the difference could be accomplished by computing the number of negative difference values as a proportion of all BL i (Xˆ ) − BL i (Yˆ ) paired differences, repeating this process for L bootstrap-replicated samples, and calculating the average proportion of negative differences across the L samples.

I copy here below the code used in Matlab which I was not able to re-tranlate in stata. WOuld any of you be so kind to help me (and probably many others in my situation) to construct a code for stata?

% Test differences treatment A vs B
clear
clc
close all
 
H = 1000;    %produce H distribution, from each of them I draw one time
HH = 1000;    %number of draws from the mean 
p_valuesLess = 999*ones(18,2);
p_valuesMore = 999*ones(18,2);
p_values2Side = 999*ones(18,2);
r = 1;
t=1;
 
% ------------------------------------------------------------------------
% load data without and with made for information
Data =[
 
];
 
 
Data_I =[
 
 
];
 
for i=1:18
    mC1=Data(r,1)
    mC2=Data(r+1,1)
    stC1=Data(r,2)
    stC2=Data(r+1,2)
    mC1_I=Data_I(r,1)
    mC2_I=Data_I(r+1,1)
    stC1_I=Data_I(r,2);
    stC2_I=Data_I(r+1,2);
 
% ------------------------------------------------------------------------
% Calculations
R_1  = mvnrnd(mC1,stC1.^2,HH)';
R_2  = mvnrnd(mC2,stC2.^2,HH)';
R_1_I  = mvnrnd(mC1_I,stC1_I.^2,HH)';
R_2_I  = mvnrnd(mC2_I,stC2_I.^2,HH)';
 
 
 
% G&P test for mean values
vectdiff_b1 = zeros(H,1);
i=1;
while i <= H
    for j=1:H
         vectdiff_b1(j,:)= R_1_I(1,i)' - R_1(1,j)' ;
    end
    if i <2
        vectdiff_old = vectdiff_b1;
    else
        vectdiff_1 = vertcat(vectdiff_b1,vectdiff_old);
        vectdiff_old = vectdiff_1; 
    end
    i= i+1;
end
 
vectdiff_b2 = zeros(H,1);
i=1;
while i <= H
    for j=1:H
         vectdiff_b2(j,:)= R_2_I(1,i)' - R_2(1,j)';
    end
    if i <2
        vectdiff_old = vectdiff_b2;
    else
        vectdiff_2 = vertcat(vectdiff_b2,vectdiff_old);
        vectdiff_old = vectdiff_2;
    end
    i= i+1;
end
 
[value_1, pvalue_1]= min(abs(sort(vectdiff_1)));
[value_2, pvalue_2]= min(abs(sort(vectdiff_2)));
 
%one side tests
pvalue_Less = 999*ones(1,2);
pvalue_Less(1,1) =  1-(pvalue_1/(H*H));
pvalue_Less(1,2) =  1-(pvalue_2/(H*H));
 
pvalue_More = 999*ones(1,2);
pvalue_More(1,1) =  (pvalue_1/(H*H));
pvalue_More(1,2) =  (pvalue_2/(H*H));
 
% 2-tail test
pvalue_2Side = 999*ones(1,2);
if pvalue_1 > .5*(H*H)
   pvalue_2Side(1,1) = 2*(1- (pvalue_1(1)/(H*H)));
else
  pvalue_2Side(1,1) =  2*(pvalue_1(1)/(H*H));
end
 
if pvalue_2 > .5*(H*H)
   pvalue_2Side(1,2) = 2*(1- (pvalue_2(1)/(H*H)));
else
  pvalue_2Side(1,2) =  2*(pvalue_2(1)/(H*H));
end
 
%home country proud
p_valuesLess(t,:) = pvalue_Less
%inferiority complex
p_valuesMore(t,:) = pvalue_More
p_values2Side(t,:) = pvalue_2Side
 
r=r+2
t=t+1
end

Thanks

Federica

[Federica Di Marcantonio]

Thanks to @Clyde here you can find a solution.
F