Contents

function [mProbitC,mProbitSigma] = probit(maxN,minN,omega,nCstates,overlineC,underlineC,v,varargin)

probit.m

This script calculates population values of Bresnahan and Reiss's estimates of

$$ \pi(N)/\pi(1) $$

based on their ordered probit model of entry.

CN=[kron(ones(maxN-minN+1,1),omega') kron((minN:1:maxN)',ones(nCstates,1))];

CData = exp(CN(:,1));
NData = CN(:,2);
Input argument "maxN" is undefined.

Error in ==> probit at 11
CN=[kron(ones(maxN-minN+1,1),omega') kron((minN:1:maxN)',ones(nCstates,1))];

Create an in-line function for the standard normal c.d.f.

Phi = @(x) 0.5*(x>=0).*(erf(abs(x)/sqrt(2))+1)+0.5*(x<0).*erfc(abs(x)/sqrt(2));

Eliminate all sates with near-zero probability.

This avoids numerical errors arising from multiplying very, very small probabilities by ``infinite'' log likelihood evaluations.

eps=1e-10;
CData=CData(v>eps); NData=NData(v>eps); v=v(v>eps); v=v/sum(v);

Set up the binary likelihood function and its derivatives.

We first estimate the entry thresholds with standard binary probits. These estimates then serve as starting values for the ordered probit estimation. So that it can be used for all of the thresholds, the rank under consideration is read as one of the likelihood function's inputs.

$$ L = I(N\geq R)\ln(\Phi(\ln(C/\overline{C}_R)/\sigma))+I(N<R)\ln(\Phi(\ln(\overline{C}_R/C)/\sigma))$$

lf = @(x) v'*((NData>=x(3)).*log(Phi(log(CData/x(1))/x(2))) + (NData<x(3)).*log(Phi(log(x(1)./CData)/x(2))));

% Calculate the minimum and maximum values of N consistent with the ergodic distribution.
Nlow=min(NData);
Nhigh=max(NData);

Maximize the binary likelihood functions.

This uses the optimization toolbox for its minimization.

bProbitC=zeros(Nhigh-Nlow-1,1);
bProbitSigma=bProbitC;
for R=Nlow+1:Nhigh;
    f = @(x) -lf([x R]);
    x=fminsearch(f,[0.5*(exp(overlineC(R))+exp(underlineC(R))) 0.2]);
    bProbitC(R)=x(1);
    bProbitSigma(R)=x(2);
end

Set up the full multinomial likelihood function.

As written, x must be a row vector with [sigma epsilon Cbar(Nlow) Cbar(Nlow+1) ... Cbar(Nhigh) infinity].

logprobs = @(x) log(Phi(log(CData./x(NData-Nlow+2)')/x(1))...
    -Phi(log(CData./x(NData-Nlow+3)')/x(1)));

%Likelihood function. First argument is standard deviation, and remaining arguments are thresholds. (``zero'' and
%``infinity'' coming first and last).
lf = @(x) v'*logprobs([ x(1) 1e-10 x(2:end) 1e10]);

Maximize the ordered probit likelihood function.

f= @(x) -lf(x);
xm=fminsearch(f,[mean(bProbitSigma) bProbitC']);
mProbitC=xm(2:end);
mProbitSigma=xm(1);

Report results if this is a test run.

if nargin > 7
   disp(mProbitC);
   disp(mProbitSigma);
   disp(f(xm));
end