% OLSTable5.m
% OLS project, sept. 2007
%
% Runs a Monte Carlo Experiment for adj. OLS x ML 
% (results reported in the paper in Table 5)
% First Price Sealed Bid Auction with values i.i.d. Uniform, mean mu std.
% sigma
% number of bidders: i.i.d. uniform {2,...,maxn}

mu = 3;  %true value of mu
sig = 1; %true value of sigma

S = 1000; %number of samples
maxn = 6; %maximum number of bidders in each auction

opts = optimset('MaxFunEvals',600);
randn('state',1234)
rand('twister',1234)

load as
an = as(1:maxn-1,2);   % this assumes the 2nd column in as is a(n) for uniform 

for T = [50 100 200];   %number of auctions in each sample

    betasML = zeros(S,2);
    betasLeo = zeros(S,2);
    DD=ones(T,1)*(2:maxn);
    a = mu-sqrt(3)*sig; %lower bound of supp(V)

    for s = 1:S,
       n = 2+floor((maxn-1)*rand(T,1));
       Eps = sqrt(12)*(rand(T,maxn)-.5); % Eps ~U[-sqrt(3),sqrt(3)]
       V = mu+sig*Eps;

       p = zeros(T,1);
       for i = 1:T,
          v = V(i,1:n(i))';
          bids = zeros(size(v));
          for j = 1:n(i),
              bids(j) = v(j)-(v(j)-a)/n(i);
          end          
          b = sort(bids,1,'descend');
          p(i) = b(1);
       end

       % estimation by Rezende
       An = (DD==n*ones(1,maxn-1))*an;
       br=regress(p,[ones(T,1),An]);
       betasLeo(s,:) = br';

       % estimation by ML
       bml = fminsearch(@(m) OLSSealUniflike(m,p,n),[3,1],opts);
       betasML(s,:) = bml;
    end
    % betasLeo: sample of adjusted OLS estimates
    % betasML: sample of ML estimates
    % first column: estimates for mu
    % second column: estimates for sigma

    % calculation of statistics
    means = mean([betasLeo,betasML]);
    vars = var([betasLeo,betasML]);
    mse = (means-[3 1 3 1]).^2+vars;
    skew = skewness([betasLeo,betasML]);
    kurt = kurtosis([betasLeo,betasML]);
    JB = S/6*(skew.^2+(kurt-3).^2/4);
    pJB = cdf('chi2',JB,2);

    sb = sort([betasLeo,betasML]);
    median = (sb(S/2,:)+sb(S/2+1,:))/2;
    LQ = (sb(S/4,:)+sb(S/4+1,:))/2;
    UQ = (sb(3*S/4,:)+sb(3*S/4+1,:))/2;

    T
    stats = [means;vars;mse;JB;pJB]
end