Contents
- probit.m
- Create an in-line function for the standard normal c.d.f.
- Eliminate all sates with near-zero probability.
- Set up the binary likelihood function and its derivatives.
- Maximize the binary likelihood functions.
- Set up the full multinomial likelihood function.
- Maximize the ordered probit likelihood function.
- Report results if this is a test run.
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
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.
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